eukaryotic gene expression is controlled at multiple levels , and splicing of mrna is an important regulatory step in the production of functional proteins . during mrna splicing , portions of the rna , introns ,
are removed by the spliceosome complex , and the remaining protein - coding rnas , exons , are joined together .
alternative splicing , in which different combinations of exons are included in the final protein - coding mrna , is responsible for the diversity of protein - coding mrnas that can be derived from a single open reading frame [ 14 ] .
the location of these splice junctions is generally defined by a g - t dna sequence signal at the splice donor and an a - g dna sequence signal at the splice acceptor [ 1 , 5 ] .
it is likely that more subtle sequence features in addition to the local sequence context of these splice sequences play a role in both constitutive splicing and alternative splicing .
transcription and splicing appear to be mechanistically coupled [ 614 ] , and the precise rules governing recognition and regulation of constitutive and alternative splicing are poorly understood .
several recent papers have focused on the role that the packaging of the template dna into chromatin plays in the splicing process [ 6 , 1523 ] .
the packaging of eukaryotic dna into chromatin is expected to affect all dna templated processes .
the fundamental subunit of chromatin is the nucleosome , 150 base pairs ( bp ) of dna wrapped around a histone octamer .
the position and density of nucleosomes play key regulatory roles and are controlled both by chromatin regulatory complexes and by features intrinsic to the dna sequence [ 2426 ] .
there is limited information on how these two determinants of nucleosome occupancy coordinate to regulate responses , such as transcription in mammalian cells .
nucleosome forming and nucleosome inhibitory properties were derived from first principles more than three decades ago .
more recently , maps of nucleosome distribution have enabled the development of computational models that use dna sequence features to predict the nucleosome forming or inhibitory potential of dna .
position specific scoring matrices   were developed by aligning strong nucleosome forming sequences and identifying the positions of enriched dinucleotide composition .
machine learning tools [ 24 , 28 ] , such as support vector machines ( svms ) , have been used to predict the intrinsic nucleosome occupancy likelihood ( nol ) for any dna sequence .
the support vector machine used in these studies is discriminative , rather than generative , and uses dna sequences representing the strongest and the weakest affinity for forming nucleosomes . in these studies , dna sequences protected from mnase digestion
genomic loci with the highest and the lowest nucleosome occupancy signals were used to train an svm .
the resulting algorithm can be used to classify the nucleosome forming potential of any genomic sequence .
a comparative assessment of available nucleosome occupancy prediction algorithms revealed that the svm trained on human chromatin worked well on related species with relatively large , complex genomes .
these seminal discoveries affirmed an important role for intrinsic dna sequence features influencing chromatin organization and revealed the utility of combined genomics and computational approaches for chromatin research .
these exonic nucleosomes have been proposed to act as  speed bumps  that allow for cotranscriptional splicing [ 32 , 33 ] .
additionally , posttranslational modifications of the nucleosomal histones at these exonic boundaries have been shown to affect splice site usage [ 18 , 19 , 34 ] .
this exonic nucleosomal occupancy is conserved throughout metazoan evolution , indicating a role for genomic sequence in this process .
we have previously described a computational model of nucleosomal occupancy trained on dna sequence content .
a large fraction of nucleosome positions can be accurately predicted based on dna sequence , indicating a significant dna sequence component to nucleosome positioning [ 25 , 26 ] .
nucleosome positioning signals at exon boundaries are conserved throughout evolution , and we reasoned that the analysis of predicted dna - encoded nucleosomal occupancy would improve characterization of regulatory features associated with intron - exon boundaries .
we used predictions of nucleosome occupancy to characterize intron / exon boundaries in five metazoan species and identified a cryptic sequence - based dna property that is specific to a subset of fundamental metabolic genes .
our results suggest that , by defining regulated positions for nucleosomes , dna features other than the consensus splice site sequence play a role in splicing .
the dna sequences for the current builds of all organisms ( human , hg19 ; rat , rn4 ; zebrafish , danrer7 ; fly , dm3 ; worm , ce10 ; yeast , saccer3 ) in this analysis were downloaded from the ucsc genome bioinformatics website ( http://hgdownload.cse.ucsc.edu/downloads.html ) .
annotations for gene structure ( transcription start site , strand , and boundaries between introns and exons ) were acquired from the refgene tables of the ucsc genome browser public mysql database with the exception of yeast for which the sgdgene table was used .
nucleosome occupancy likelihood ( nol ) scores were generated using the support vector machine ( svm ) model derived from the ozsolak a375 dataset ( described in ) .
the svm scores a 50  bp window of dna sequence and uses a 1  bp step sliding window for sequences longer than 50  bp .
the resulting score indicates the likelihood that the associated 50-mer is nucleosome forming ( positive value ) or nucleosome inhibitory ( negative value ) .
gene ontology analysis was completed with the gorilla software   using genes of interest in the target set and the list of all genes in the associated genome as the background set . for calculations of significance across genomic statistics in enriched ontological categories ,
outliers were defined as those values less than the first quartile minus the interquartile range ( iqr ) or greater than the third quartile plus the iqr .
p values were then calculated using a random sample of the same size from the whole genome and all values from the set of interest to perform a two sample t - test with a p value of < 0.05 being considered significant .
we were first interested in discovering if there was any intrinsic nucleosome occupancy information in the regions flanking intron - exon boundaries .
we reasoned that if chromatin structure plays a role in cotranscriptional splicing , then a robust location to store that chromatin structural information would be within the dna sequence itself . to this end , we retrieved all intron - exon boundaries from the refseq annotation of the human genome .
we calculated nucleosome occupancy likelihood ( nol ) scores for sequences spanning 1000  bp centered on intron - exon boundaries .
nol scores were calculated for each 50-mer and the resulting value was assigned to the center nucleotide position .
nol scores were plotted to reveal overarching trends of predicted nucleosome positions at intron - exon boundaries .
consistent with previous observations , we detected characteristic nucleosome positions at these boundaries ( figure 1(a ) ) .
sequence analysis of the regions surrounding these boundaries showed clear consensus splice donor and splice acceptor sequences ( figure 1(a ) ) .
closer inspection of the plots shown in figure 1(a ) revealed a sharp dip in average nol scores for the introns both upstream and downstream of exons .
this sharp dip , identified by the nol , could only result from a set of sequence features that reduce the nucleosome forming potential at that location .
we were next interested in understanding the characteristics of the sequences that contributed to this sharp dip .
we hypothesized that the dip in the nol scores may represent a functional dna - encoded chromatin - regulatory structural element , as this is what the nol plots measure . to investigate this possibility , we aligned the entire dataset to this putative regulatory feature by centering each region on the minimum found in the 100 base pairs centered on the boundary ( figure 1(b ) ) .
the percentage representation of each base at all positions was calculated and graphically represented to identify a previously undetected consensus sequence defining intron - exon boundaries .
the sequence composition of these shifted subsets reflected equal representation of all four dna bases and therefore did not reveal any clear consensus sequence feature ( figure 1(b ) , sequence diagram ) .
this result suggests that a more cryptic combination of sequence features reflecting some physical property of dna may be defining this low scoring element .
different biophysical properties emerge with different dna sequence combinations ( e.g. , dna wedge angles ) .
we were , therefore , interested in determining if the dips upstream and downstream conferred by the dna sequence of the exon might identify a cryptic characteristic of intron boundary architecture .
we calculated the location of the minimum nol score for each boundary window and then represented these data as a histogram of distance of the minima from the intron - exon boundary ( figure 2(a ) ) .
we observed a striking overrepresentation 26 nucleotides ( nt ) from the annotated intron - exon boundary , upstream of the exon ( u  26 ) .
we also found a similar feature + 26  nt from the exon - intron boundary , downstream of the exon ( d + 26 ) .
the discovery of the u  26 and d + 26 features prompted us to investigate how many and what types of genes include these features . in order to understand the numbers and types of genes associated with the u  26 and d + 26 feature , we selected the exons represented in each of these groups , u  26 and d + 26 , for further analysis
. there are 9578 genes represented in the u  26 group and 7360 genes represented in the d + 26 group , representing 24.1% and 18.5% of all open reading frames tested , respectively .
we next wanted to know if the u  26 and d + 26 features were found in the same sets or different sets of genes .
3369 genes , or 19% , overlap between the u  26 and d + 26 groups ( figure 2(b ) ) . as positioned nucleosomes flanking exons are phylogenetically conserved , we were next interested in determining whether the prominent u  26 and the d + 26 signature are conserved in other metazoan species .
conservation of these features would suggest an important role for the u  26 and the d + 26 causing it to be maintained by evolution .
we identified the nucleotide position of the minimum at the boundary between intron and exon for rat , zebrafish , fly , and worm . a conspicuous u  26 and d + 26 signature exists for all of these species ( figures 2(c ) , 2(e ) , 2(g ) , and 2(i ) ) . as with the human example
this result led us to hypothesize that the signature was associated with a particular gene set that is conserved across each of the species tested .
we next wanted to identify the feature that is present in groups of genes with related function . in order to test whether the u  26 or the d + 26 signatures identified groups of genes that share a common function , we searched for ontological enrichment .
the u  26 and the d + 26 signatures both showed enrichment for overlapping groups of gene function ( table 1 ) .
these groups include atpase activity , atp binding , helicase activity , and motor activity . with little exception , these enrichments persist throughout all metazoan species tested ( table 1 ) .
these results led us to hypothesize that these groups of genes contain genomic characteristics that differ from the set of all genes . in order to test whether genes found in the enriched functional categories containing the u  26 and the d + 26 signature had genomic characteristics that varied significantly from the rest of the genome , we compared exon size , intron size , and number of exons for each function category to the same values calculated for the genome as a whole ( figure 3 ) .
the exon sizes of genes included in each of the categories did not differ substantially from the exon sizes encoded by each genome ( supplementary figure  1 , available online at http://dx.doi.org/10.1155/2015/167578 ) .
intron sizes of the genes in the ontological function categories were significantly shorter when compared to the respective whole genomes with few exceptions ( rat : motor activity , zebrafish : atp binding , and worm : atp binding ) . across all organisms , with the exception of atpase activity and helicase activity in rat , all ontological categories across all organisms exhibited significantly higher numbers of exons per gene .
for example , the human genome has a median of seven exons per gene while the subset of human motor activity genes has approximately five times the number of exons with a median of 35 exons per gene .
as our experiments , to this point , had been purely based on in silico experiments and genome sequence data , we wanted to see how our results compared to published in vivo data .
we next were most interested in testing if the u  26 or the d + 26 signatures are observed in vivo .
we used the nucleosome maps generated by the encode consortium for the two cell lines , gm12878 and k562 .
we retrieved all nucleosome measurements for intron - exon boundaries for each of the two cell lines from http://genome.ucsc.edu/   and identified the location of the minimum within the central 100 base pairs surrounding intron - exon boundaries as described previously using the nol scores .
approximately 2.5% of genes contained intron - exon boundaries exhibiting the d  26 and u + 26 characteristic ( table 2 ) .
we wanted to know whether this in vivo signature was associated with the shared common function consistent with those identified from the nol predictions .
three of the four categories identified by the nol score signal , atp binding , atpase activity , and motor activity , were also significantly enriched in the in vivo data .
these findings indicate that sequence - based nucleosome forming signals at the boundaries between introns and exons may play a role in the regulation of these genes .
the upstream intron contains a region between the branchpoint sequence and the 3 splice site that is generally depleted of ag dinucleotides .
this region is generally within 40 nucleotides of the ag splice site and encompasses the u  26 feature .
we were interested in determining if the loci containing the u  26 feature were enriched or depleted for any dinucleotide occurrences relative to the rest of the genome .
we calculated dinucleotide frequency for the ten dinucleotides for the loci containing the u  26 feature and compared that to the dinucleotide frequency for an equal number of other intron - exon boundaries in the genome ( figure 4 ) .
we found that the intronic region off these genes is actually slightly more depleted of ag dinucleotides than other intron - exon boundaries .
we were next interested in determining if the remaining dinucleotides showed differential intron and exon content .
recent work has shown that differential g / c content plays a role in the intron exon definition and splice site selection [ 14 , 38 ] .
cc , cg , and gc dinucleotides were depleted in introns of the loci containing the u  26 feature .
likewise , aa , at , and ta dinucleotides were enriched in the introns of the loci containing the u  26 feature compared to equivalent regions in the rest of the genome ( figure 4 ) . loci containing the u  26 feature have a lower overall intronic g / c content .
the exonic region , however , shows the opposite trend with decreased a / t and increased g / c content .
the cc dinucleotide feature most strongly defines this set of loci and is depleted from introns and enriched in exons .
this would suggest that dinucleotide content plays a larger role in the definition of intron - exon boundary than was previously anticipated .
we have identified a set of conserved genes sharing a common function using a nucleosome positioning signature .
we have further characterized the set of genes as having increased numbers of exons while having average number and length of introns .
thus , we have identified a set of conserved genes with common function and distinguishing features that suggest shared regulation .
our observations and the classification of a particular subset of genes could not have been accomplished through alignment of nucleotide content .
the nol scores point toward a physical property of dna related to the ability of a particular dna sequence to form a nucleosome .
the organization and architecture of dna around the nucleosome may likely play a role in the mechanism of pre - mrna splicing .
 the regulation of metazoan gene expression occurs in part by pre - mrna splicing into mature rnas . 
 signals affecting the efficiency and specificity with which introns are removed have not been completely elucidated . 
 splicing likely occurs cotranscriptionally , with chromatin structure playing a key regulatory role . 
 we calculated dna encoded nucleosome occupancy likelihood ( nol ) scores at the boundaries between introns and exons across five metazoan species . 
 we found that ( i ) nol scores reveal a sequence - based feature at the introns on both sides of the intron - exon boundary ; ( ii ) this feature is not part of any recognizable consensus sequence ; ( iii ) this feature is conserved throughout metazoa ; ( iv ) this feature is enriched in genes sharing similar functions : atpase activity , atp binding , helicase activity , and motor activity ; ( v ) genes with these functions exhibit different genomic characteristics ; 
 ( vi ) in vivo nucleosome positioning data confirm ontological enrichment at this feature ; and ( vii ) genes with this feature exhibit unique dinucleotide distributions at the intron - exon boundary . 
 the nol scores point toward a physical property of dna that may play a role in the mechanism of pre - mrna splicing . 
 these results provide a foundation for identification of a new set of regulatory dna elements involved in splicing regulation . 
1. Background
2. Methods
3. Results and Discussion
4. Conclusion