BACKGROUND
the advent of novel high-throughput technologies such as next-generation sequencing  has empowered the gathering of transcriptional data of unprecedented quantity and quality.
it has lead to the discovery that ~75% of the human genome is transcribed to rna at some point in certain cell types  <cit> , with other estimates being as high as 90%  <cit> .
only ~1-2% of these transcripts are translated into protein form  <cit> , while the remainder are classified as non-coding rnas   <cit> .
the fact that ncrnas are relatively abundant, expressed in a developmentally regulated fashion  <cit> , and often exhibit precise sub-cellular localization  <cit>  supports the notion that they play important biological roles.
in particular, their ability to base-pair with other transcripts and regions suggests they may be responsible for a variety of regulatory functions  <cit> .
attempts to untangle the complex landscape of ncrnas have led to crude classification of ncrnas based on their length   <cit> , function , transfer rnas ), regulatory potential , long non-coding rnas )  <cit> , and subcellular localization , small nucleolar rnas , cytoplasm-located piwi-interacting rnas , and short interfering rnas ).
other unusual ncrna species such as trans-spliced transcripts, macrornas that encompass enormous genomic distances, and multi-gene transcripts that encompass several genes or even the whole chromosome further confound efforts for systematic classification  <cit> .
in reality, however, clear categorization of ncrna classes has been quite difficult, as many ncrna transcripts often share the properties of multiple categories.
ncrnas can also be broadly categorized by genomic properties such as 1) sense or antisense transcripts, when the ncrna transcript co-localizes with exons of another transcript on the same or the opposite strand, respectively, 2) bidirectional, when the expression of the ncrna and of a neighboring coding transcript in the opposite strand are initiated in close proximity, 3) intronic, when the ncrna arises from an intron of another transcript, and 4) intergenic, when the ncrna is localized between two coding transcripts  <cit> .
other ncrna classes, such as promoter-associated rnas  and enhancer rnas  are also being characterized  <cit> .
these considerations alone suggest the role of genomic organization in the biogenesis of ncrnas is more complex than previously thought  <cit> .
recent years have seen a rapid growth of publicly available data on genome organization and functional annotation.
the encyclopedia of dna elements  project has been actively cataloging functional elements in the human genome, such as cell type-specific histone modification profiles, chromatin states, and transcription factor binding sites  <cit> .
in this study we will refer to these functional and regulatory regions as epigenomic elements, i.e., genomic data other than nucleotide sequence that describes functions, properties, or experimental values associated with genome regions  <cit> .
although the precise definition of the word "epigenomics" is hotly debated  <cit>  and is often narrower in scope, we feel our broader definition of epigenomic elements is suitable to convey the concepts of our research to the readers.
the aim of the current study was to use the encode data, such as cell type-specific histone modification profiles, chromatin states, and experimentally validated transcription factor binding sites to examine how epigenomic elements statistically significantly associated with different classes of ncrnas may reveal functional roles and relationships among ncrna classes.
we investigated whether the genomic coordinates of  <dig> ncrna classes significantly overlapped with or were depleted in a total of  <dig> annotated epigenomic elements.
we found that pirna and trna ncrnas were strongly under- and over-represented with the majority of epigenomic elements examined, respectively.
other ncrna classes showed less extreme enrichments, but we were able to confirm known and identify new statistically significant ncrna class-specific epigenomic associations.
clustering ncrna classes by the significance of their epigenomic associations captures their known hierarchy, providing a means to utilize epigenomic background for classification purposes.
in summary, our study demonstrates a means to use genome annotation data to identify regulatory commonalities and differences among genomic and epigenomic elements.
methods
ncrna data
grch37/hg <dig> genome assembly coordinates for different classes of non-coding rnas were extracted from mirbase , ucsc genome database , and rfam  using custom scripts available at git://wrenlab.org/ncrna-loci.
pirnas were extracted from the ucsc knowngenes table by string matching on "*pirna*", and rfam categories were used for rfam-derived elements.
we ignored ncrna classes with less than  <dig> members, as their analysis did not show any statistically significant associations .
the relationship among ncrna classes is shown on figure  <dig>  a total of  <dig>  ncrnas were analyzed .
functional/regulatory elements data
as part of our interest in automating the search for biologically meaningful correlations among high-throughput and high-information data  <cit> , we developed genomerunner , a software program that searches for statistically significant co-localization between a set of genomic regions of interest  and sets of annotated genomic features   <cit> .
we tested each different class of ncrnas to see if they were statistically significantly associated with or depleted in three groups of epigenomic elements: chromatin state segmentation by hmm from encode/broad, histone modifications by chip-seq from encode/broad institute, and experimentally validated transcription factor chip-seq from encode .
a total of  <dig> epigenomic elements were obtained from the ucsc genome database  <cit> , and stored in a local mysql database accessible for the genomerunner community.
enrichment analysis
briefly, sets of ncrnas were tested for genomic overlap of at least  <dig> nucleotide with  <dig> epigenomic elements.
we first calculated the total number of members from an ncrna class that overlapped with an epigenomic element.
if a given ncrna overlaps with > <dig> epigenomic element of the same type, it is only counted once, to reflect the fact of overlap.
in other words, overlap counting is ncrna-centric, e.g., if a trna overlaps with  <dig> nfkb binding sites, only one is added to the trna-nfkb overlap counter.
we then performed random sampling from the total pool of all  <dig> classes of ncrnas, selecting the same number of random ncrnas as in the class being analyzed.
the use of all ncrnas as a background enables us to search, specifically, for associations that might distinguish one ncrna class from the rest.
this background selection also restricts random sampling to avoid low complexity genomic regions, such as transposable elements, duplications and inversions, repeats, comprising up to 47% of the human genome  <cit> .
we selected a random number of ncrna elements according to the size of each ncrna class and performed  <dig> such random samplings, estimating the average number and variance of random overlapping.
a chi-square test was used to determine whether there was a significant difference between the total number of co-localization for the ncrna class as compared with what could be expected by random chance.
transformation of p-values
to emphasize epigenomic elements most significantly associated with an ncrna class, p-values were adjusted for multiple testing using the benjamini-hochberg procedure  <cit> .
for easier visualization and comparison, we converted p-values into decimal scale by -log10-transformation.
a "-"sign was added if a p-value represents an underrepresented association.
this allows representation of significant associations in an intuitive format - larger numbers equal more statistically significant overrepresented associations of a ncrna class with an epigenomic element, while smaller negative numbers represent more statistically significant under-represented associations.
clustering and visualization
the transformed p-values can be directly visualized with a blue-yellow gradient representing under- and over-represented associations, respectively.
a n × m matrix of transformed p-values, where n is the number of ncrna classes and m is the number of epigenomic elements, was assembled .
epigenomic elements showing no statistically significant associations  with at least one ncrna class were removed.
further, epigenomic elements showing consistent associations across the ncrna classes  were also filtered.
these filtering steps simplify visualization and allow a reader to focus on the most significant epigenomic associations differentially enriched among ncrna classes.
hierarchical clustering was performed using "maximum" distance to measure dissimilarity between rows and columns, and the "ward" agglomeration method  <cit> .
clustering and visualization were performed within r computing environment  <cit> .
pair-wise correlation analysis
we compared epigenomic similarities and differences among ncrna classes.
that is, pair-wise correlations between all ncrna class-specific transformed p-values were measured using pearson's correlation coefficient.
we expect ncrna classes showing similar enrichment patterns with epigenomic elements to correlate with positive pearson's correlation coefficient close to  <dig>  while ncrna classes differentially associated with epigenomic elements would show negative pearson's correlation coefficient.
combining these pearson's correlation coefficients into n × n matrix, where n is the number of ncrna classes, allows clustering and visualization of epigenomic similarities among ncrna classes as describe above, using blue/yellow gradient to highlight negative/positive correlations, respectively.
RESULTS
the number of members per ncrna class varies considerably, ranging from  <dig> scarnas to  <dig>  small nuclear rnas .
this variability in class size was accompanied by variability in their statistically significant enrichments versus all  <dig> epigenomic features.
the trna and pirna classes showed the largest variability of the enrichment/depletion p-values, after correction for multiple testing .
there was a positive correlation between the size of ncrna classes and the standard deviation of the p-value , although it was not significant , likely due to the extreme p-value variability in trna and pirna classes.
removing those classes increased pearson's correlation coefficient to  <dig>  .
notable, that both classes did not show extreme size or length abnormalities, being comparable in size with snorna and mirna classes , suggesting strong epigenomic associations driving p-value variability.
pair-wise correlation analysis of epigenomic associations reveal known and novel relationships among ncrna classes
we investigated how similar ncrna classes were across all their associations with epigenomic elements.
that is, p-values for each ncrna-epigenomic element were calculated and their similarities across all  <dig> elements analyzed were quantified by pearson's correlation coefficient .
out of  <dig> elements analyzed,  <dig> did not show statistically significant associations with any of the ncrna classes.
therefore,  <dig> were used for the global correlation analysis .
expectedly, snorna, c/d-box and h/aca-box ncrna classes shared similar epigenomic associations, being subgroups of a larger snrna class .
yet scarnas, a subgroup of snornas, did not share similar significant associations.
this may be, in part, due to size differences in the groups .
along with trna class, scarnas had epigenomic associations distinct from other ncrna classes.
spliceosomal, rrna, 7sk, pirna, srp-rna and y-rna classes formed the largest group with similar epigenomic associations.
while ncrna classes that cluster together by their epigenomic associations provide a good means to assess potential functional/regulatory similarities, each class can also have anti-correlated epigenomic associations.
pairs of ncrna classes best  correlating with each other are listed in table  <dig>
pirnas are depleted, while trnas are enriched in histone modification marks
combinations of histone modifications mark the boundaries of transcriptionally active regions, and are associated with distinct regulatory sites  <cit> .
we expected to observe both enriched and depleted associations of ncrna classes with histone marks.
surprisingly, ncrna classes showed either preferentially enriched or depleted statistically significant associations with histone modification marks.
pirnas, spliceosomal, srp-rnas, rrnas, y-rnas were predominantly depleted in histone modification marks from the  <dig> cell lines examined.
trnas, h/aca-box, pri-mirnas, mirnas, lncrnas, c/d-box and snornas were enriched in histone modification marks.
notably, pirnas and trnas showed the largest variability in enrichment .
transfer rnas and piwi-interacting rnas showed the strongest anti-correlated associations with epigenomic elements .
the trnas were strongly enriched in the actively transcribed genomic regions, permissive and activation-related h3k <dig> di- and trimethylation, h3k <dig> and h3k <dig> acetylation marks, and h2az histone variant.
contrary to the overall picture of active transcription regions associated with trnas, they were also associated with the transcription insulator ctcf binding sites  that blocks communication between promoters and downstream genes  <cit> .
although we identified that trnas tend to be enriched in regions with ctcf insulators, whether or not ctcf binds in these regions remains to be determined.
conversely, pirnas were associated with heterochromatin regions in the genome and underrepresented in active promoter regions .
pirnas were depleted in the same types of histone marks the trnas were enriched in.
that is, pirnas were depleted in activation-related h3k <dig> mono-, di- and trimethylation marks, h3k9ac, h3k27ac and h3k36me <dig> marks, and the h2az histone variant.
yet, the association of pirnas with transcriptionally silent regions of the genome is blurred by their depletion in h4k20me <dig>  a marker of silent chromatin, and in h3k27me <dig> marker of gene repression.
heatmap visualization helps summarize over- and underrepresentation trends for other ncrna classes .
the h3k <dig> di- and trimethylation marks were the most frequently differentially enriched, followed by h3k9ac.
we did not observe bias in the cell type-specific origin of these histone marks, although k <dig>  nhek, h1hesc, hela-s <dig> and hepg <dig> specific marks dominated the heatmap.
the driving force behind the clustering appeared to be the histone marks themselves.
this can be seen from groups of h3k79me <dig>  h3k4me <dig> and h3k4me <dig> showing similar enrichment patterns, thus clustered together, on the heatmap.
interestingly, even though c/d box and h/aca box rnas are subclasses of snornas, the h/aca box subclass was not enriched for most of the associations that either snornas or c/d box rnas were, in general.
this is in contrast to the general similarity of associated features shown in figure 5b and indicates that one of the factors that may separate the two snorna subclasses is their regulatory locations.
the only mark h/aca box and c/d boxes seemed to both be enriched in was h3k79me <dig> which has been associated with the cell cycle.
this suggests that h/aca and c/d box rnas may be orchestrated similarly by histone modification marks during the cell cycle, but not beyond.
pirnas are enriched in the repressed chromatin states, trnas are associated with transcriptionally active regions
the chromatin state segmentation by hmm from encode/broad covers a spectrum of chromatin states, ranging from active promoters and strong enhancers to heterochromatin and repetitive/copy number variation regions.
expectedly, we found nearly symmetrical distribution of enrichment/depletion of ncrna associations in opposing chromatin states, reflecting their presence in either active or repressed regions and, conversely, absence in the region with opposite regulatory role.
again, pirna and trna classes showed greater range of p-values .
pirnas were strongly enriched in the repressed chromatin , heterochromatin  and insulator  regions, as well as in the repetitive/cnv regions in the h1hesc cell line.
conversely, they were depleted in the active, weak and poised promoters , transcription elongation and transition regions , strong and weak enhancers  .
although these findings suggest pirnas are generally located at the inactive regions of the genome, cell type-specific chromatin states also showed an exception from these observations.
for example, pirnas were depleted in the heterochromatin regions in only the k <dig> cell line.
trnas showed greater diversity in cell type-specific enriched and depleted associations with chromatin states.
they were enriched in the active, weak, and poised promoters , and depleted in heterochromatin regions in  <dig> cell lines.
although this suggests active transcription processes associated with trnas, they were depleted in transcription elongation, transition, weak transcription , and in weak enhancers in the hmec and huvec cell lines.
trnas were also strongly overrepresented in the repetitive/cnv regions in  <dig> cell lines .
other ncrna classes showed similar enrichment/depletion patterns in the antagonistic chromatin state regions.
similar to the histone marks, chromatin states were the major driving force behind clustering.
for example, pri-mirna, mirna, and lncrna classes were depleted in heterochromatin regions in the majority of the cell lines, while rrna, 7sk and spliceosomal classes were enriched in these regions .
the snorna class, and its members c/d box, h/aca box and scarnas were enriched in the "transcription elongation" and "transcription transition" regions in h1hesc cells, and depleted in the "repressed" regions in nhek and hepg <dig> cells.
these classes were also enriched in the "strong enhancers" and depleted in "weak transcription" chromatin regions in the k <dig> cell line.
the y-rna and srp-rna classes, on the contrary, were enriched in weak transcription regions in k <dig> cells .
pirnas are depleted in, while trnas are strongly enriched in, transcription factor binding sites
trnas showed strong enrichment in their tendencies to co-localize with  <dig> tfbss  .
not unexpectedly, the strongest association was rpc <dig> , a catalytic core and the largest  component of rna polymerase iii , which synthesizes trnas as well as small rnas, such as 5s rrna.
other pol iii components were also enriched, such as bdp <dig> , tfiiic- <dig> , brf <dig> , pol <dig> .
tbp, tata-binding protein, a subunit of the transcription factor tfiid  has also been shown to participate in the initiation complex assembly on trna  <cit> .
a general co-activator protein p <dig>  was also strongly associated with trnas.
pol iii can sense non-self dsdna that serves as template for transcription into dsrna.
such non-self pol iii transcripts induce type i interferon and nf-kb through rig-i pathway  <cit> .
we observed strong enrichment of immune-related transcription factor binding sites associated with trnas.
notably, ets <dig>  a transcription factor involved in hematopoietic cell differentiation and the development of lymphoid tissues, was the most strongly enriched , accompanied by elf <dig> .
other well known immune-related transcription factors include nfkb, irf <dig>  irf <dig>  and stat <dig> .
pirnas showed completely the opposite picture , being statistically significantly depleted in  <dig> tfbss at p-value adjusted for multiple testing < <dig>  .
pirnas were strongly depleted in pol <dig> binding sites assessed by four different antibodies , p.adj =  <dig> e-09; pol <dig>  p.adj =  <dig> e-06).
we also observed the opposite behavior of the pol iii-related proteins, with tbp , rcp <dig> , bdp <dig> , tfiiic- <dig> , pol <dig> , and p <dig>  binding sites.
the aforementioned rna polymerase ii- and iii-related tfbss were also statistically significantly associated with other ncrna classes .
a general trend was that an ncrna class shows opposite associations with rna polymerases ii and iii, being generally associated with pol ii sites and depleted in pol iii-related sites.
mirnas were found enriched in pol ii sites, which corroborates earlier reports  <cit> .
discussion
we comprehensively investigated statistically significant genomic associations between a total of  <dig> ncrna classes versus  <dig> cell type-specific epigenomic elements.
although the encode project provides cell type-specific information about epigenomic elements, we generally did not observe preferential cell type specificity of the enrichments.
epigenomic elements associated with different ncrna classes were generally similarly enriched across multiple cell lines and were the driving force behind clustering patterns confirming that epigenomic relationships hold across multiple cell lines  <cit> .
this can be illustrated by similar patterns of enrichment among spliceosomal, rrna and 7sk ncrnas in the heterochromatin regions in multiple cell lines , association of c/d box and snornas with h3k79me <dig> mark across several cell lines , and by similar observations in other ncrna classes .
there were exceptions from this general trend; for example, pirnas were strongly enriched in the heterochromatin regions in all but the k <dig> cell line, which instead showed statistically significant underrepresentation.
y-rna and srp-rna classes were strongly associated with weak transcription regions exclusively in the k <dig> cell line.
therefore, we kept the cell type-specific results of our analysis  while describing our results from the epigenomic point of view, mentioning cell line specificity only where appropriate.
our study revealed both expected  and novel relationships among ncrna classes and epigenomic elements.
for example, c/d box, h/acabox, and snorna were grouped by their common epigenomic background, being subsets of a larger snrna class.
yet, the relationship between long non-coding rnas and mirnas is interesting, as they are located in different genomic regions but shared similar epigenomic background.
both lncrna and mirna classes were depleted in the heterochromatin regions in multiple cell lines, were located in the poised promoter regions in the h1hesc cell line, and showed similar enrichment in several histone modification marks .
yet they were differentially associated with active promoters, with only lncrnas being associated with these regions in the hsmm and huvec cell lines.
the lncrna class was also enriched in the repressed regions of the genome in gm <dig> and nhek cell lines, which is consistent with an observation that >20% of lncrnas are bound by the polycomb repressive complex  <dig>   <cit> , with mirnas showing no measurable enrichment.
mirnas, on the contrary, were enriched in the transcription elongation regions in hmec, huvec, and nhlf cell lines.
mirnas were depleted in the pol iii transcription factor binding regions, and enriched only in the pol <dig> binding sites, while lncrnas did not show any associations with any of these regions.
instead, lncrnas were enriched in suz <dig> binding sites, consistent with an observation that prc2-transcribed smaller ~50- <dig> nt rnas that interact with suz <dig> to mediate gene repression  <cit> .
we observed strong statistically significant associations of the trna and pirna classes of non-coding rnas with epigenomic elements.
this can't be attributed solely to the size of these classes, since ncrna classes having more members did not show such enrichment diversity.
the mean length of trnas and pirnas were among the smallest among all ncrna classes , thus diminishing the probability of overlapping with epigenomic elements on a per-transcript basis.
yet, the significance of the identified associations suggests strong epigenomic regulation of trna and pirna.
we expect our analyses to produce more focused results as the ncrna classes will become better defined with more members showing clearer enrichment patterns.
this is especially important for small ncrna classes, like scarnas, which did not have sufficient numbers to demonstrate statistically significant enrichments.
advancements of the encode project with more epigenomic data being available will further aid analysis of cell type-specific associations.
our enrichment analysis has limitations.
we define ncrna classes by the genomic coordinates of ncrna elements, and test these locations for the enriched co-localizations with epigenomic elements.
thus, our enrichment analysis of transcription factor binding sites does not consider the promoter regions of the ncrnas or, if they don't have promoters, then their local cis-regulatory genomic environment.
yet, the transcription of ncrnas can be different from that of genes, with ncrnas often originating from exons and introns  <cit> , promoters of the genes  <cit> , or intergenic regions by mechanisms apparently different from canonical transcription factor binding  <cit> .
transcription factors, on the other hand, show a great variety of binding patterns, being located within gene bodies, overlapping first exons, and occurring in the intergenic regions  <cit> .
our analysis revealed tf binding preference  within or partially overlapping ncrna bodies without making assumptions about the mechanism of ncrna biogenesis.
while the full potential of transcription factor co-localization with ncrna classes remains to be tested experimentally, we are confident our findings of differential transcription factor co-localization with ncrnas will reveal novel relationships among non-coding and coding parts of the genome.
in the current study we did not consider the amount of overlap of ncrnas with epigenomic elements.
while this is not critical for the analysis of, for example, single nucleotide polymorphisms with length equal to  <dig> nt, ncrnas may show a variety of co-localization patterns, ranging from being located within an epigenomic element to overlapping it completely or partially.
we restricted our analysis by considering any overlap as co-localization for computational efficiency and ease of interpretation.
however, our future work will include comparison of the percent of overlap observed for an ncrna class vs. that of for a random selection.
this "overlap enrichment" analysis may provide better resolution of spatial relationships of ncrna classes with epigenomic elements.
CONCLUSIONS
we present the use of our method for finding statistically significant relationships between experimental regions of interest  with cell type-specific epigenomic data from the encode project .
we demonstrate the utility of statistically significant epigenomic associations to classify and contrast ncrna classes, and to outline their potential functional roles and mechanisms of biogenesis from epigenomic perspective.
our method, implemented as the open source software genomerunner   <cit>  can provide interpretation of any experimental genome-wide data within the growing amount of epigenomic data from the encode project  <cit> .
availability
genomerunner, the main method used in this work, is available at http://sourceforge.net/projects/genomerunner
competing interests
the authors declare that they have no competing interests.
authors' contributions
mgd conceived the idea of using epigenomic elements for systematic ncrna classification, designed and conducted the experiments, interpreted the results, drafted the manuscript and prepared it for submission.
cbg helped with data analysis and generation of the results.
kak helped with generation of the results.
jdw helped with interpretation of the results and writing the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
the encode data summary used in the current study
click here for file
additional file 2
cell types used to obtain cell type-specific epigenomic elements.
the vocabulary is extracted from the encode consortium data coordination center at ucsc .
m/f/u indicate cell/tissue donor was of male, female, or unknown gender, respectively.
click here for file
additional file 3
brief description of the known roles of histone modification marks
click here for file
additional file 4
enrichment matrix for the different ncrna classes.
rows show cell type specific epigenomic marks, except for the transcription factor binding sites .
-log <dig> values are shown.
red font highlights epigenomic marks enriched at p.adj <  <dig> .
green font and a "-" sign highlight depleted epigenomic marks at p.adj <  <dig> .
rows containing no statistically significant values in at least one ncrna class are removed for clarity.
click here for file
additional file 5
epigenomic elements statistically significantly associated with trna and pirna classes.
regular & red font indicates enriched associations, italics & gree indicates depleted associations.
common transcription factors for both trna and pirna classes are not colored, as they show enriched associations with trna and depleted associations with pirna classes.
click here for file
additional file 6
epigenomic elements statistically significantly associated with ncrna classes other that trna and pirna classes.
regular & red font indicates enriched associations, italics & green indicates depleted associations.
click here for file
declarations
publication of this article is funded in part by the national institutes of health , the national institute of arthritis and musculoskeletal and skin diseases , and an institutional development award  from the national institute of general medical sciences .
the content of this article is solely the responsibility of the authors and does not necessarily represent the official views of the national institutes of health.
this article has been published as part of bmc bioinformatics volume  <dig> supplement  <dig>  2013: proceedings of the tenth annual mcbios conference.
discovery in a sea of data.
the full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/14/s <dig>

 the advent of novel high-throughput technologies such as next-generation sequencing has empowered the gathering of transcriptional data of unprecedented quantity and quality. it has lead to the discovery that ~75% of the human genome is transcribed to rna at some point in certain cell types, with other estimates being as high as 90%. only ~1-2% of these transcripts are translated into protein form, while the remainder are classified as non-coding rnas. the fact that ncrnas are relatively abundant, expressed in a developmentally regulated fashion, and often exhibit precise sub-cellular localization supports the notion that they play important biological roles. in particular, their ability to base-pair with other transcripts and regions suggests they may be responsible for a variety of regulatory functions. attempts to untangle the complex landscape of ncrnas have led to crude classification of ncrnas based on their length, function, transfer rnas ), regulatory potential, long non-coding rnas ), and subcellular localization, small nucleolar rnas, cytoplasm-located piwi-interacting rnas, and short interfering rnas ). other unusual ncrna species such as trans-spliced transcripts, macrornas that encompass enormous genomic distances, and multi-gene transcripts that encompass several genes or even the whole chromosome further confound efforts for systematic classification. in reality, however, clear categorization of ncrna classes has been quite difficult, as many ncrna transcripts often share the properties of multiple categories. ncrnas can also be broadly categorized by genomic properties such as 1) sense or antisense transcripts, when the ncrna transcript co-localizes with exons of another transcript on the same or the opposite strand, respectively, 2) bidirectional, when the expression of the ncrna and of a neighboring coding transcript in the opposite strand are initiated in close proximity, 3) intronic, when the ncrna arises from an intron of another transcript, and 4) intergenic, when the ncrna is localized between two coding transcripts. other ncrna classes, such as promoter-associated rnas and enhancer rnas are also being characterized. these considerations alone suggest the role of genomic organization in the biogenesis of ncrnas is more complex than previously thought. recent years have seen a rapid growth of publicly available data on genome organization and functional annotation. the encyclopedia of dna elements project has been actively cataloging functional elements in the human genome, such as cell type-specific histone modification profiles, chromatin states, and transcription factor binding sites. in this study we will refer to these functional and regulatory regions as epigenomic elements, i.e., genomic data other than nucleotide sequence that describes functions, properties, or experimental values associated with genome regions. although the precise definition of the word "epigenomics" is hotly debated and is often narrower in scope, we feel our broader definition of epigenomic elements is suitable to convey the concepts of our research to the readers. the aim of the current study was to use the encode data, such as cell type-specific histone modification profiles, chromatin states, and experimentally validated transcription factor binding sites to examine how epigenomic elements statistically significantly associated with different classes of ncrnas may reveal functional roles and relationships among ncrna classes. we investigated whether the genomic coordinates of ncrna classes significantly overlapped with or were depleted in a total of annotated epigenomic elements. we found that pirna and trna ncrnas were strongly under- and over-represented with the majority of epigenomic elements examined, respectively. other ncrna classes showed less extreme enrichments, but we were able to confirm known and identify new statistically significant ncrna class-specific epigenomic associations. clustering ncrna classes by the significance of their epigenomic associations captures their known hierarchy, providing a means to utilize epigenomic background for classification purposes. in summary, our study demonstrates a means to use genome annotation data to identify regulatory commonalities and differences among genomic and epigenomic elements. methods ncrna data grch37/hg genome assembly coordinates for different classes of non-coding rnas were extracted from mirbase, ucsc genome database, and rfam using custom scripts available at git://wrenlab.org/ncrna-loci. pirnas were extracted from the ucsc knowngenes table by string matching on "*pirna*", and rfam categories were used for rfam-derived elements. we ignored ncrna classes with less than members, as their analysis did not show any statistically significant associations. the relationship among ncrna classes is shown on figure a total of ncrnas were analyzed. functional/regulatory elements data as part of our interest in automating the search for biologically meaningful correlations among high-throughput and high-information data, we developed genomerunner, a software program that searches for statistically significant co-localization between a set of genomic regions of interest and sets of annotated genomic features. we tested each different class of ncrnas to see if they were statistically significantly associated with or depleted in three groups of epigenomic elements: chromatin state segmentation by hmm from encode/broad, histone modifications by chip-seq from encode/broad institute, and experimentally validated transcription factor chip-seq from encode. a total of epigenomic elements were obtained from the ucsc genome database, and stored in a local mysql database accessible for the genomerunner community. enrichment analysis briefly, sets of ncrnas were tested for genomic overlap of at least nucleotide with epigenomic elements. we first calculated the total number of members from an ncrna class that overlapped with an epigenomic element. if a given ncrna overlaps with > epigenomic element of the same type, it is only counted once, to reflect the fact of overlap. in other words, overlap counting is ncrna-centric, e.g., if a trna overlaps with nfkb binding sites, only one is added to the trna-nfkb overlap counter. we then performed random sampling from the total pool of all classes of ncrnas, selecting the same number of random ncrnas as in the class being analyzed. the use of all ncrnas as a background enables us to search, specifically, for associations that might distinguish one ncrna class from the rest. this background selection also restricts random sampling to avoid low complexity genomic regions, such as transposable elements, duplications and inversions, repeats, comprising up to 47% of the human genome. we selected a random number of ncrna elements according to the size of each ncrna class and performed such random samplings, estimating the average number and variance of random overlapping. a chi-square test was used to determine whether there was a significant difference between the total number of co-localization for the ncrna class as compared with what could be expected by random chance. transformation of p-values to emphasize epigenomic elements most significantly associated with an ncrna class, p-values were adjusted for multiple testing using the benjamini-hochberg procedure. for easier visualization and comparison, we converted p-values into decimal scale by -log10-transformation. a "-"sign was added if a p-value represents an underrepresented association. this allows representation of significant associations in an intuitive format - larger numbers equal more statistically significant overrepresented associations of a ncrna class with an epigenomic element, while smaller negative numbers represent more statistically significant under-represented associations. clustering and visualization the transformed p-values can be directly visualized with a blue-yellow gradient representing under- and over-represented associations, respectively. a n × m matrix of transformed p-values, where n is the number of ncrna classes and m is the number of epigenomic elements, was assembled. epigenomic elements showing no statistically significant associations with at least one ncrna class were removed. further, epigenomic elements showing consistent associations across the ncrna classes were also filtered. these filtering steps simplify visualization and allow a reader to focus on the most significant epigenomic associations differentially enriched among ncrna classes. hierarchical clustering was performed using "maximum" distance to measure dissimilarity between rows and columns, and the "ward" agglomeration method. clustering and visualization were performed within r computing environment. pair-wise correlation analysis we compared epigenomic similarities and differences among ncrna classes. that is, pair-wise correlations between all ncrna class-specific transformed p-values were measured using pearson's correlation coefficient. we expect ncrna classes showing similar enrichment patterns with epigenomic elements to correlate with positive pearson's correlation coefficient close to while ncrna classes differentially associated with epigenomic elements would show negative pearson's correlation coefficient. combining these pearson's correlation coefficients into n × n matrix, where n is the number of ncrna classes, allows clustering and visualization of epigenomic similarities among ncrna classes as describe above, using blue/yellow gradient to highlight negative/positive correlations, respectively.  the number of members per ncrna class varies considerably, ranging from scarnas to small nuclear rnas. this variability in class size was accompanied by variability in their statistically significant enrichments versus all epigenomic features. the trna and pirna classes showed the largest variability of the enrichment/depletion p-values, after correction for multiple testing. there was a positive correlation between the size of ncrna classes and the standard deviation of the p-value, although it was not significant, likely due to the extreme p-value variability in trna and pirna classes. removing those classes increased pearson's correlation coefficient to. notable, that both classes did not show extreme size or length abnormalities, being comparable in size with snorna and mirna classes, suggesting strong epigenomic associations driving p-value variability. pair-wise correlation analysis of epigenomic associations reveal known and novel relationships among ncrna classes we investigated how similar ncrna classes were across all their associations with epigenomic elements. that is, p-values for each ncrna-epigenomic element were calculated and their similarities across all elements analyzed were quantified by pearson's correlation coefficient. out of elements analyzed, did not show statistically significant associations with any of the ncrna classes. therefore, were used for the global correlation analysis. expectedly, snorna, c/d-box and h/aca-box ncrna classes shared similar epigenomic associations, being subgroups of a larger snrna class. yet scarnas, a subgroup of snornas, did not share similar significant associations. this may be, in part, due to size differences in the groups. along with trna class, scarnas had epigenomic associations distinct from other ncrna classes. spliceosomal, rrna, 7sk, pirna, srp-rna and y-rna classes formed the largest group with similar epigenomic associations. while ncrna classes that cluster together by their epigenomic associations provide a good means to assess potential functional/regulatory similarities, each class can also have anti-correlated epigenomic associations. pairs of ncrna classes best correlating with each other are listed in table pirnas are depleted, while trnas are enriched in histone modification marks combinations of histone modifications mark the boundaries of transcriptionally active regions, and are associated with distinct regulatory sites. we expected to observe both enriched and depleted associations of ncrna classes with histone marks. surprisingly, ncrna classes showed either preferentially enriched or depleted statistically significant associations with histone modification marks. pirnas, spliceosomal, srp-rnas, rrnas, y-rnas were predominantly depleted in histone modification marks from the cell lines examined. trnas, h/aca-box, pri-mirnas, mirnas, lncrnas, c/d-box and snornas were enriched in histone modification marks. notably, pirnas and trnas showed the largest variability in enrichment. transfer rnas and piwi-interacting rnas showed the strongest anti-correlated associations with epigenomic elements. the trnas were strongly enriched in the actively transcribed genomic regions, permissive and activation-related h3k di- and trimethylation, h3k and h3k acetylation marks, and h2az histone variant. contrary to the overall picture of active transcription regions associated with trnas, they were also associated with the transcription insulator ctcf binding sites that blocks communication between promoters and downstream genes. although we identified that trnas tend to be enriched in regions with ctcf insulators, whether or not ctcf binds in these regions remains to be determined. conversely, pirnas were associated with heterochromatin regions in the genome and underrepresented in active promoter regions. pirnas were depleted in the same types of histone marks the trnas were enriched in. that is, pirnas were depleted in activation-related h3k mono-, di- and trimethylation marks, h3k9ac, h3k27ac and h3k36me marks, and the h2az histone variant. yet, the association of pirnas with transcriptionally silent regions of the genome is blurred by their depletion in h4k20me a marker of silent chromatin, and in h3k27me marker of gene repression. heatmap visualization helps summarize over- and underrepresentation trends for other ncrna classes. the h3k di- and trimethylation marks were the most frequently differentially enriched, followed by h3k9ac. we did not observe bias in the cell type-specific origin of these histone marks, although k nhek, h1hesc, hela-s and hepg specific marks dominated the heatmap. the driving force behind the clustering appeared to be the histone marks themselves. this can be seen from groups of h3k79me h3k4me and h3k4me showing similar enrichment patterns, thus clustered together, on the heatmap. interestingly, even though c/d box and h/aca box rnas are subclasses of snornas, the h/aca box subclass was not enriched for most of the associations that either snornas or c/d box rnas were, in general. this is in contrast to the general similarity of associated features shown in figure 5b and indicates that one of the factors that may separate the two snorna subclasses is their regulatory locations. the only mark h/aca box and c/d boxes seemed to both be enriched in was h3k79me which has been associated with the cell cycle. this suggests that h/aca and c/d box rnas may be orchestrated similarly by histone modification marks during the cell cycle, but not beyond. pirnas are enriched in the repressed chromatin states, trnas are associated with transcriptionally active regions the chromatin state segmentation by hmm from encode/broad covers a spectrum of chromatin states, ranging from active promoters and strong enhancers to heterochromatin and repetitive/copy number variation regions. expectedly, we found nearly symmetrical distribution of enrichment/depletion of ncrna associations in opposing chromatin states, reflecting their presence in either active or repressed regions and, conversely, absence in the region with opposite regulatory role. again, pirna and trna classes showed greater range of p-values. pirnas were strongly enriched in the repressed chromatin, heterochromatin and insulator regions, as well as in the repetitive/cnv regions in the h1hesc cell line. conversely, they were depleted in the active, weak and poised promoters, transcription elongation and transition regions, strong and weak enhancers. although these findings suggest pirnas are generally located at the inactive regions of the genome, cell type-specific chromatin states also showed an exception from these observations. for example, pirnas were depleted in the heterochromatin regions in only the k cell line. trnas showed greater diversity in cell type-specific enriched and depleted associations with chromatin states. they were enriched in the active, weak, and poised promoters, and depleted in heterochromatin regions in cell lines. although this suggests active transcription processes associated with trnas, they were depleted in transcription elongation, transition, weak transcription, and in weak enhancers in the hmec and huvec cell lines. trnas were also strongly overrepresented in the repetitive/cnv regions in cell lines. other ncrna classes showed similar enrichment/depletion patterns in the antagonistic chromatin state regions. similar to the histone marks, chromatin states were the major driving force behind clustering. for example, pri-mirna, mirna, and lncrna classes were depleted in heterochromatin regions in the majority of the cell lines, while rrna, 7sk and spliceosomal classes were enriched in these regions. the snorna class, and its members c/d box, h/aca box and scarnas were enriched in the "transcription elongation" and "transcription transition" regions in h1hesc cells, and depleted in the "repressed" regions in nhek and hepg cells. these classes were also enriched in the "strong enhancers" and depleted in "weak transcription" chromatin regions in the k cell line. the y-rna and srp-rna classes, on the contrary, were enriched in weak transcription regions in k cells. pirnas are depleted in, while trnas are strongly enriched in, transcription factor binding sites trnas showed strong enrichment in their tendencies to co-localize with tfbss. not unexpectedly, the strongest association was rpc, a catalytic core and the largest component of rna polymerase iii, which synthesizes trnas as well as small rnas, such as 5s rrna. other pol iii components were also enriched, such as bdp, tfiiic-, brf, pol. tbp, tata-binding protein, a subunit of the transcription factor tfiid has also been shown to participate in the initiation complex assembly on trna. a general co-activator protein p was also strongly associated with trnas. pol iii can sense non-self dsdna that serves as template for transcription into dsrna. such non-self pol iii transcripts induce type i interferon and nf-kb through rig-i pathway. we observed strong enrichment of immune-related transcription factor binding sites associated with trnas. notably, ets a transcription factor involved in hematopoietic cell differentiation and the development of lymphoid tissues, was the most strongly enriched, accompanied by elf. other well known immune-related transcription factors include nfkb, irf irf and stat. pirnas showed completely the opposite picture, being statistically significantly depleted in tfbss at p-value adjusted for multiple testing <. pirnas were strongly depleted in pol binding sites assessed by four different antibodies, p.adj = e-09; pol p.adj = e-06). we also observed the opposite behavior of the pol iii-related proteins, with tbp, rcp, bdp, tfiiic-, pol, and p binding sites. the aforementioned rna polymerase ii- and iii-related tfbss were also statistically significantly associated with other ncrna classes. a general trend was that an ncrna class shows opposite associations with rna polymerases ii and iii, being generally associated with pol ii sites and depleted in pol iii-related sites. mirnas were found enriched in pol ii sites, which corroborates earlier reports. discussion we comprehensively investigated statistically significant genomic associations between a total of ncrna classes versus cell type-specific epigenomic elements. although the encode project provides cell type-specific information about epigenomic elements, we generally did not observe preferential cell type specificity of the enrichments. epigenomic elements associated with different ncrna classes were generally similarly enriched across multiple cell lines and were the driving force behind clustering patterns confirming that epigenomic relationships hold across multiple cell lines. this can be illustrated by similar patterns of enrichment among spliceosomal, rrna and 7sk ncrnas in the heterochromatin regions in multiple cell lines, association of c/d box and snornas with h3k79me mark across several cell lines, and by similar observations in other ncrna classes. there were exceptions from this general trend; for example, pirnas were strongly enriched in the heterochromatin regions in all but the k cell line, which instead showed statistically significant underrepresentation. y-rna and srp-rna classes were strongly associated with weak transcription regions exclusively in the k cell line. therefore, we kept the cell type-specific results of our analysis while describing our results from the epigenomic point of view, mentioning cell line specificity only where appropriate. our study revealed both expected and novel relationships among ncrna classes and epigenomic elements. for example, c/d box, h/acabox, and snorna were grouped by their common epigenomic background, being subsets of a larger snrna class. yet, the relationship between long non-coding rnas and mirnas is interesting, as they are located in different genomic regions but shared similar epigenomic background. both lncrna and mirna classes were depleted in the heterochromatin regions in multiple cell lines, were located in the poised promoter regions in the h1hesc cell line, and showed similar enrichment in several histone modification marks. yet they were differentially associated with active promoters, with only lncrnas being associated with these regions in the hsmm and huvec cell lines. the lncrna class was also enriched in the repressed regions of the genome in gm and nhek cell lines, which is consistent with an observation that >20% of lncrnas are bound by the polycomb repressive complex, with mirnas showing no measurable enrichment. mirnas, on the contrary, were enriched in the transcription elongation regions in hmec, huvec, and nhlf cell lines. mirnas were depleted in the pol iii transcription factor binding regions, and enriched only in the pol binding sites, while lncrnas did not show any associations with any of these regions. instead, lncrnas were enriched in suz binding sites, consistent with an observation that prc2-transcribed smaller ~50- nt rnas that interact with suz to mediate gene repression. we observed strong statistically significant associations of the trna and pirna classes of non-coding rnas with epigenomic elements. this can't be attributed solely to the size of these classes, since ncrna classes having more members did not show such enrichment diversity. the mean length of trnas and pirnas were among the smallest among all ncrna classes, thus diminishing the probability of overlapping with epigenomic elements on a per-transcript basis. yet, the significance of the identified associations suggests strong epigenomic regulation of trna and pirna. we expect our analyses to produce more focused results as the ncrna classes will become better defined with more members showing clearer enrichment patterns. this is especially important for small ncrna classes, like scarnas, which did not have sufficient numbers to demonstrate statistically significant enrichments. advancements of the encode project with more epigenomic data being available will further aid analysis of cell type-specific associations. our enrichment analysis has limitations. we define ncrna classes by the genomic coordinates of ncrna elements, and test these locations for the enriched co-localizations with epigenomic elements. thus, our enrichment analysis of transcription factor binding sites does not consider the promoter regions of the ncrnas or, if they don't have promoters, then their local cis-regulatory genomic environment. yet, the transcription of ncrnas can be different from that of genes, with ncrnas often originating from exons and introns, promoters of the genes, or intergenic regions by mechanisms apparently different from canonical transcription factor binding. transcription factors, on the other hand, show a great variety of binding patterns, being located within gene bodies, overlapping first exons, and occurring in the intergenic regions. our analysis revealed tf binding preference within or partially overlapping ncrna bodies without making assumptions about the mechanism of ncrna biogenesis. while the full potential of transcription factor co-localization with ncrna classes remains to be tested experimentally, we are confident our findings of differential transcription factor co-localization with ncrnas will reveal novel relationships among non-coding and coding parts of the genome. in the current study we did not consider the amount of overlap of ncrnas with epigenomic elements. while this is not critical for the analysis of, for example, single nucleotide polymorphisms with length equal to nt, ncrnas may show a variety of co-localization patterns, ranging from being located within an epigenomic element to overlapping it completely or partially. we restricted our analysis by considering any overlap as co-localization for computational efficiency and ease of interpretation. however, our future work will include comparison of the percent of overlap observed for an ncrna class vs. that of for a random selection. this "overlap enrichment" analysis may provide better resolution of spatial relationships of ncrna classes with epigenomic elements.  we present the use of our method for finding statistically significant relationships between experimental regions of interest with cell type-specific epigenomic data from the encode project. we demonstrate the utility of statistically significant epigenomic associations to classify and contrast ncrna classes, and to outline their potential functional roles and mechanisms of biogenesis from epigenomic perspective. our method, implemented as the open source software genomerunner can provide interpretation of any experimental genome-wide data within the growing amount of epigenomic data from the encode project. availability genomerunner, the main method used in this work, is available at http://sourceforge.net/projects/genomerunner competing interests the authors declare that they have no competing interests. authors' contributions mgd conceived the idea of using epigenomic elements for systematic ncrna classification, designed and conducted the experiments, interpreted the results, drafted the manuscript and prepared it for submission. cbg helped with data analysis and generation of the results. kak helped with generation of the results. jdw helped with interpretation of the results and writing the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 the encode data summary used in the current study click here for file additional file 2 cell types used to obtain cell type-specific epigenomic elements. the vocabulary is extracted from the encode consortium data coordination center at ucsc. m/f/u indicate cell/tissue donor was of male, female, or unknown gender, respectively. click here for file additional file 3 brief description of the known roles of histone modification marks click here for file additional file 4 enrichment matrix for the different ncrna classes. rows show cell type specific epigenomic marks, except for the transcription factor binding sites. -log values are shown. red font highlights epigenomic marks enriched at p.adj <. green font and a "-" sign highlight depleted epigenomic marks at p.adj <. rows containing no statistically significant values in at least one ncrna class are removed for clarity. click here for file additional file 5 epigenomic elements statistically significantly associated with trna and pirna classes. regular & red font indicates enriched associations, italics & gree indicates depleted associations. common transcription factors for both trna and pirna classes are not colored, as they show enriched associations with trna and depleted associations with pirna classes. click here for file additional file 6 epigenomic elements statistically significantly associated with ncrna classes other that trna and pirna classes. regular & red font indicates enriched associations, italics & green indicates depleted associations. click here for file declarations publication of this article is funded in part by the national institutes of health, the national institute of arthritis and musculoskeletal and skin diseases, and an institutional development award from the national institute of general medical sciences. the content of this article is solely the responsibility of the authors and does not necessarily represent the official views of the national institutes of health. this article has been published as part of bmc bioinformatics volume supplement 2013: proceedings of the tenth annual mcbios conference. discovery in a sea of data. the full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/14/s 
conversely, transfer rnas were depleted in heterochromatin regions and strongly enriched in regions containing activating h3k di- and trimethylation marks, h2az histone variant, and a variety of tfbss. we examined different ncrna classes for statistically significant genomic co-localizations with cell type-specific chromatin segmentation states, transcription factor binding sites, and histone modification marks using genomerunner. we clustered and visualized the ncrna classes by the strength of their statistical enrichments and depletions. searching for statistically significant associations between ncrna classes and epigenomic elements permits detection of potential functional and/or regulatory relationships among ncrna classes, and suggests cell type-specific biological roles of ncrnas. we found piwi-interacting rnas to be depleted in regions containing activating histone modification marks, such as h3k mono-, di- and trimethylation, h3k acetylation, as well as certain tfbss. pirnas were further depleted in active promoters, weak transcription, and transcription elongation regions, and enriched in repressed and heterochromatic regions. ncrnanon-coding rnaepigeneticsgenomeencodegenomerunner5- april tenth annual mcbios conference. we hypothesized that examining how different classes of ncrnas co-localized with annotated epigenomic elements could help understand the functions, regulatory mechanisms, and relationships among ncrna families. p-values were obtained using a chi-square test and corrected for multiple testing using the benjamini-hochberg procedure. even though only % of the human genome is translated into proteins, recent reports indicate that most of it is transcribed into non-coding rnas, which are becoming the subject of increased scientific interest.
 even though only % of the human genome is translated into proteins, recent reports indicate that most of it is transcribed into non-coding rnas, which are becoming the subject of increased scientific interest. we hypothesized that examining how different classes of ncrnas co-localized with annotated epigenomic elements could help understand the functions, regulatory mechanisms, and relationships among ncrna families.  we examined different ncrna classes for statistically significant genomic co-localizations with cell type-specific chromatin segmentation states, transcription factor binding sites, and histone modification marks using genomerunner. p-values were obtained using a chi-square test and corrected for multiple testing using the benjamini-hochberg procedure. we clustered and visualized the ncrna classes by the strength of their statistical enrichments and depletions. we found piwi-interacting rnas to be depleted in regions containing activating histone modification marks, such as h3k mono-, di- and trimethylation, h3k acetylation, as well as certain tfbss. pirnas were further depleted in active promoters, weak transcription, and transcription elongation regions, and enriched in repressed and heterochromatic regions. conversely, transfer rnas were depleted in heterochromatin regions and strongly enriched in regions containing activating h3k di- and trimethylation marks, h2az histone variant, and a variety of tfbss. interestingly, regions containing ctcf insulator protein binding sites were associated with trnas. trnas were also enriched in the active, weak and poised promoters and, surprisingly, in regions with repetitive/copy number variations.  searching for statistically significant associations between ncrna classes and epigenomic elements permits detection of potential functional and/or regulatory relationships among ncrna classes, and suggests cell type-specific biological roles of ncrnas. ncrnanon-coding rnaepigeneticsgenomeencodegenomerunner5- april tenth annual mcbios conference. discovery in a sea of data columbia, mo, usa 