BACKGROUND
one of the major findings of genome sequencing has been that the primary sequence of roughly 5% of the human and mouse genomes is under purifying selection, indicating functionality  <cit> .
however, less than 2% is accounted for by mrna exons.
the remaining 3% presumably encompasses cis-regulatory sequence, signals for transcriptional initiation, termination, rna processing, chromosomal features such as replication origins, and genes encoding ncrnas such as trna, snorna, mirna, and others.
accurate computational identification of novel ncrna genes and mrna structural elements  in genome sequence has proven to be more difficult than identification of exons, due generally to a limited and highly variable sequence signature  <cit> .
in bacteria, which have compact genomes, searching for transcription initiation signals  <cit> , primary sequence conservation  <cit> , and base composition  <cit>  have been fruitful approaches to de novo ncrna discovery; however, these features alone are unlikely to be sufficiently specific in large eukaryotic genomes.
most, albeit not all, functional ncrna features possess some degree of secondary structure, either as part of the precursor or the functional rna itself.
following the assumption that structural rna sequences should be more thermodynamically stable than random permutations of the same base composition, thermodynamic stability  is an additional feature than can be incorporated into genomic searches for new ncrnas.
major classes of structural rnas have lower Δg than corresponding shuffled sequences.
it has been debated whether Δg is a sufficiently accurate discriminant when only a single  sequence is analyzed  <cit> ; however, Δg has been proposed to be comparable or superior to more sophisticated algorithms  when applied independently to segments of a pairwise alignment  <cit> .
what is clear is that disruption of dinucleotides in the random permutation dramatically affects the perceived precision of predictions  <cit> , presumably because dinucleotide contributions are a key determinant of stability of an rna fold.
covariance  is also now a widely-accepted approach to ncrna discovery.
a variety of recently-described ncrna search algorithms  score for covariance to discriminate structural rna elements .
success of covariance requires that sequences be sufficiently conserved to achieve a correct alignment, yet contain some nucleotide changes in order to assess compensatory mutations.
an advantage of methods that do not utilize covariance is that they can identify structures common to sequences without high sequence similarity  <cit>  and apparently even sequences that fail to align at the primary sequence level  <cit> .
however, even considering both covariance and thermodynamic stability, some classes of ncrnas appear to be more difficult to detect than others  <cit> .
hybrid/oligoarrayaux.php
tools/zmfold.pl
to our knowledge, most ncrna search tools have not been assessed or compared systematically by an independent laboratory.
moreover, previous studies using these tools for genomic scanning have acknowledged that the false-positive rate is high and cannot be determined accurately  <cit> .
here, with the goal of independently evaluating these tools for de novo discovery of ncrna elements in a genomic context, we have compiled an extensive sequence and alignment data set, and developed a new shuffling algorithm for pairwise alignments  that simultaneously preserves key features of the alignment.
we used these sequences, as well as real and shuffled genomic sequence, as input for a panel of published tools for discovery of novel ncrnas.
we also evaluated whether different schemes utilizing only thermodynamic stability could discriminate real from shuffled ncrnas, and examined the output of a subset of tools on a genomic scale.
RESULTS
a test set for comparing ncrna-finding tools
we tested the ncrna search tools  on a test set derived by extracting  <dig> genomic alignments corresponding to known eukaryotic ncrnas and mrna structural elements .
the positive test set is available at  <cit>  and in additional files  <dig> and  <dig>  and consists of  <dig> mirna,  <dig> trna,  <dig> box h/aca snorna,  <dig> box c/d snorna,  <dig> snrna,  <dig> utr regulatory element alignments, in addition to  <dig> other structural ncrna alignments  and  <dig> non-structural rna alignments  for a total of  <dig>  bases.
we generated the alignments by mapping all available ncrna sequences from human, mouse, c. elegans, d. melanogaster, and s. cerevisiae, to available ucsc pairwise genomic alignments .
to obtain negative controls, we created  <dig> sequence shuffles of each sequence alignment using a new shuffling algorithm, shuffle-pair.pl , to simulate the problem of finding bona fide structural elements against a much larger background that resembles true genomic alignments in as many features as possible.
briefly, the shuffle-pair.pl algorithm randomizes the columns of the alignment while preserving conservation, gaps, and dinucleotides in both sequences, all of which may contribute to scoring obtained from at least a subset of algorithms tested  <cit> .
it shuffles by swapping the middle bases between trinucleotides that are identical at positions  <dig> and  <dig> as described  <cit>  with the additional constraints of swapping in aligned trinucleotides and maintaining conservation and gaps of the alignment.
the algorithm finds all swappable positions that preserve the criteria above and then randomly selects two bases to swap.
it keeps track of positions that have been swapped and continues shuffling until it runs out of positions to shuffle.
a detailed description of the algorithm and pseudocode is found in the methods section.
on average, mirna alignments shuffled 57% , trna alignments – 61%, snrna alignments – 62%, snorna alignments – 46%, and regulatory elements – 43%; overall the  <dig> alignments shuffled on average 53%.
we did not observe a correlation between the degree of shuffling and the scores output by tools .
we observed that the degree of shuffling possible depends on sequence length and conservation: increasing length and conservation increased the shufflability of the alignments, presumably due to more shuffling opportunities resulting from more bases and less constraint imposed by preserving local conservation patterns.
conversely, increasing the number of aligned sequences greatly diminishes shuffling potential due the constraint of preserving dinucleotides in all aligned sequences.
consequently, here we consider only findings pertaining to pairwise alignments.
we acknowledge that we are not taking full advantage of msari, ddbrna, evofold, and rnaz, since these tools were designed to handle multiple sequence alignments.
nonetheless, we reasoned that it would be of value to examine the outputs of these tools under a situation that controls for all variables believed to be relevant.
we ran all of the tools with default settings.
all but one of the tools ran without giving error messages.
rnaz produces error messages with sequences shorter than  <dig> and longer than  <dig> nt and for alignments with very low sequence identity, but still produces scores that we retained for the analysis since  only a small proportion of the test set  falls into this category, and  and challenges imposed in scoring these alignments apply to all of the tools.
msari output was excluded from subsequent evaluation due to its extremely poor performance on alignments of  <dig> sequences .
overview of results on test set
for each of the tools, we scored precision ) and recall )  across its range of output scores.
in comparison to area under the roc curve, and specificity ), we view precision and recall as most relevant to real-world genomic screening: precision estimates the false-positive rate one would face if experimentally verifying predicted ncrnas, and recall estimates the proportion of all ncrnas that could be found.
our results are consistent with expectation in several ways.
first, taking an operating point at which each tool yields ~50% overall recall, the recall of different classes is ranked on average mirna > trna > regulatory > box h/aca snorna > snrna > other , roughly as has been described by others .
algorithms utilizing only thermodynamic stability appear to have a relative advantage in detecting mirnas, presumably due to the simple hairpin structure.
at the same time, they have relative difficulty identifying trnas, perhaps due to the more complicated fold.
box c/d snornas and non-structural elements were virtually indistinguishable from shuffles.
we note that some of the tools may have been trained or their parameters optimized using a portion of our test set, and this may slightly bias outcome.
evofold, for example, was trained on a large set of ncrnas including many mirnas  <cit> , and it does perform well on mirnas, even in our pairwise alignments; however, we confirmed that results for evofold are virtually identical if its training set rnas are omitted , presumably because it does not learn any structural or sequence information , and therefore results for the full data set are shown.
nonetheless, on the whole, the tools appear to detect different ncrna classes in a roughly similar manner at the 50% recall operating point .
at lower overall recall  the results become less comparable due to a strong bias towards detection of mirnas by the thermodynamics-only tools.
at higher recall , all tools suffer from low precision .
our results also verify the recent demonstration of uzilov et al.
<cit>  that, on pairwise alignments, z-scores based on Δg scores for both strands yield substantially higher precision at roughly equivalent recall for all ncrna classes tested, taking either an overall average or a median of class averages .
our "zrnafold " implementation is conceptually similar to "dynalign" used by uzilov  <cit>  .
in our hands dynalign was too slow  to be useful for genome-wide scanning  and we thus did not test it.
a less expected but intriguing result was that tools utilizing apparently unrelated algorithms often mutually identified and missed the same ncrnas .
for example, evofold, one of the covariance tools, displayed a very significant overlap in true positives with all of the thermodynamic-only tools ).
importance of dinucleotide frequency
a somewhat discouraging outcome of our test was that all of the tools displayed low apparent precision in our task.
at 50% recall, all of the tools displayed between  <dig> and 20% precision .
in a genomic scan, in which the background sequence most likely exceeds the 100× ratio we employed , this would translate to laboratory confirmations being overwhelmingly negative, in the absence of further filtering.
the results of our analysis must be viewed in light of being based only on pairwise alignments, which substantially handicaps tools such as ddbrna, msari, and evofold that rely solely on covariance and were created solely to scan multiple alignments.
nonetheless, we were surprised that rnaz, which has both covariance and thermodynamic components  <cit> , did not generally outperform zrnafold , which is our relatively simple two-strand z-score implementation of the rnaz thermodynamic input.
we reasoned that retention of dinucleotide counts might account for the high false-positive rate we observed.
to ask whether this is the case, we examined the results obtained using a negative-control test set generated using an alternative shuffling procedure, shuffle-aln.pl  <cit> , which preserves local conservation, gap structure, and base composition, but not dinucleotide content.
indeed, with shuffle-aln.pl, the apparent false-positive rate  decreased  <dig> -fold on average; most drastically for rnaz and qrna, which dropped  <dig>  and  <dig> -fold, respectively, at 50% recall .
this indicates that virtually all tools, including those that do not explicitly have a thermodynamic component, are nonetheless influenced in some way by dinucleotide counts, which it is known should be accounted for in estimation of false-positive rates based on thermodynamic stability  <cit> .
considerations in genomic scanning: dinucleotides and other aspects of sequence content
we next explored whether the dinucleotide issue might influence genomic scans.
we ran five of the eight tools on 120-base tiling windows of human-mouse alignments of chromosome  <dig> , and compared the distributions of scores with and without preservation of dinucleotide content, again by comparing shuffle-pair.pl with shuffle-aln.pl).
, and alifoldz were not run on these sequences).
in this analysis, the estimated false-positive rates drastically increased when using shuffles with conserved dinucleotides: in figure  <dig>  the scores obtained from shuffled sequences with dinucleotide frequency preserved  nearly overlap those obtained from the original sequence , whereas there is a clear separation of both from the scores with mononucleotide shuffling .
we reasoned that insufficient shuffling  might provide a trivial explanation for this observation; however, when we subset the analysis to those shuffles which differed in at least 50% of their nucleotide positions from their original sequences, the trend clearly remained .
even among the  <dig> windows in which, by chance, the degree of shuffling in the dinucleotide-constrained set exceeded that in the mononucleotide-shuffled set, the trend remained largely intact  .
why do dinucleotide counts influence scores from these tools?
for those relying on thermodynamics , the dependence of Δg on dinucleotides provides a likely explanation.
however, qrna does not utilize thermodynamic stability.
indeed, all the tools we examined displayed negligible biases towards specific nucleotides or dinucleotides when run on randomly-generated alignments .
only on ncrnas and genomic sequence , which contain non-random distributions of di- and even mononucleotides , did high-scoring  sequences emerge with enrichment of specific dinucleotides.
for example, in shuffled versions of the test set, the highest-scoring sequences from rnaz and qrna tended to be slightly enriched for cc and gg, and depleted of aa and at .
this trend was even more marked for shuffled versions of chromosome  <dig> alignments .
evofold appeared to be even more sensitive to dinucleotide composition of the shuffled genomic alignments, favoring those with aa, tt, at, and ta .
manual examination of the specific shuffled genomic sequences that scored highly with evofold revealed that many of them contained unusual distributions of homopolymeric or near-homopolymeric runs and simple repeat-like sequence  that are extremely rare in randomly-generated sequence, but very common in the human genome.
because the mono-, di-, and/or tri-nucleotide contents have reduced complexity, these patterns are occasionally preserved in shuffled controls.
in fact, we observed striking correlations between scores for all of the chromosome  <dig> tiling windows before and after shuffling  , suggesting that simple sequence properties can and do influence scores emerging from these tools.
in an effort to gain further insight, we examined the highest-scoring ~ <dig>  real genomic sequences from chromosome  <dig> for each of three tools utilizing covariance .
these are posted in the supplementary data .
manual examination of all three lists revealed the presence of many homopolymeric or near-homopolymeric runs, repetitive short motifs that escape repeatmasking, and nonuniform distribution of nucleotides.
the sequences that rnaz and qrna scored most highly were clearly enriched for g and c  while those scored most highly by evofold appeared enriched for homopolymeric stretches of a and t  .
anecdotally, these non-repeat-masked features do appear to contribute to the predicted rna structures.
while we cannot rule out that a subset of these are bona fide undiscovered functional ncrna features, it is also possible that the non-random distribution of nucleotides and short sequence features in genomes  could lead to a higher degree of apparent secondary structure than would be expected from a completely random evolutionary process.
we currently have no unifying explanation for why these tools appear to favour specific types of sequences, and suggest that this subject might form the basis for further study.
these trends do seem to suggest that caution should be exercised in drawing conclusions about the number of ncrna features in the human genome, on the basis of results of computational genomic scans.
from our analyses, we conclude that an accurate assessment of state-of-the-art ncrna search tools on multiple alignments, which was our initial objective, is not realistic in the absence of a shuffling algorithm that retains at least an approximation of both dinucleotide frequencies and real alignment structure.
moreover, it may be worth considering the fact that genome sequences contain a non-random distribution of mononucleotides and other short sequence features.
genomic scanning for novel conserved rna structural elements reveals enrichment for stable secondary structure in intergenic regions, introns, and utrs, and selection against structure in coding regions
our analysis above indicates that with any of the tools we tested, the vast majority of high-scoring windows in a genome scan using pairwise alignments would be false-positives.
nonetheless, we next asked whether use of shuffle-pair.pl, which we believe provides a realistic estimate of the false-positive rate, could aid in making new observations regarding frequency of structural elements across genomic features.
we used a tiling strategy similar to that employed previously  <cit> .
we used a window size of  <dig> nt  because this window size resulted in the greatest sensitivity at little cost to specificity relative to other tiling window sizes .
this window size is within range of previously computed genomic scans:  <dig>  <cit>  and  <dig>  <cit>  used in bacterial scans, and  <dig>  <cit>  used for the human genome.
on average, 23% of nucleotides in a  <dig> base window of human-mouse pairwise alignments change identity during shuffling; however, 13% of sliding windows have >50% shuffling.
we scanned pairwise genomic alignments  <cit>  of human chromosome  <dig> and genomic regions corresponding to alignable portions of  <dig>  refseq  <cit>  genes with either zmfold or evofold, as practical representatives of the thermodynamic-stability-only regime and the covariance regime.
we scanned real and shuffled counterparts  of each window.
we found that the distributions of the real and  shuffled scores largely overlap, with few or no real scores completely exceeding the shuffled distribution .
however, the real and shuffled score distributions were significantly different .
a histogram of the difference in number of real and shuffled sequences at different scores is shown at the right in figure 7a.
these results obviously bring into question the significance of any individual element, but also suggest a selection for some level of conserved secondary structure across the genome.
to ask whether results obtained from genomic scanning appear to be biologically relevant on the whole, and to gain evidence that we have not mis-calibrated the negative controls , we examined the overall distribution of score differences between real and shuffled sequences obtained from different types of genomic features , which might be expected to display differing propensities towards secondary structure.
using either zmfold or evofold, we found that striking differences exist among results obtained from these features .
in particular, introns, intergenic sequence and 3' utrs tend to have a higher distribution of scores, indicating enrichment for structural features, relative to shuffled controls.
in contrast, we observed marginal enrichment within 5' utrs.
strikingly, we found that coding regions have a lower score distribution than shuffled controls, indicating an overall depletion of rna structures in open reading frames, which is consistent with the possibility that these elements would inhibit translation  <cit> .
these results indicate that our analysis is not, on average, biased towards falsely reporting structures, since some genomic features appear to be enriched for structures while others are depleted or neutral.
these results also provide evidence for the presence of biologically functional structural features, although specificity limitations of the individual tools preclude direct prediction of individual structural elements.
to further investigate whether the tools are in fact identifying a biologically meaningful selection for structural elements we asked whether zmfold and evofold, which represent very different approaches, tend to give high scores to the same utr segments.
we found that in 3'utrs, where an  obvious enrichment for structural elements exists, the top zmfold and evofold scores are correlated .
furthermore, 48% of evofold and zmfold highest scores within each utr correspond to the same region of the utr  which is approximately two-fold higher than expected by chance .
discussion
to our knowledge, the panel of ncrna search tools examined here have not been examined systematically by an independent laboratory.
with the goal of evaluating how useful these tools are for de novo discovery of ncrna elements in a genomic context, in which the anticipated sparseness of real features makes it particularly important to control false discovery, we compiled an extensive sequence and alignment data set comprised of  <dig> fragments of genomic alignments  <cit>  corresponding to ncrnas from a variety of eukaryotes.
we also developed a shuffling algorithm for pairwise alignments  that maintains dinucleotide frequencies, in addition to other features typically maintained in alignment shuffling   <cit> .
due to shuffling constraints, which we show are important, we limited our test to pairwise alignments.
in our test, we found that failure to account for background sequence properties leads to over-estimation of precision  for most tools.
we were surprised to find that tools with a strong or exclusive covariance component were sensitive, in at least some situations, to  dinucleotide counts, which should primarily impact thermodynamic stability.
estimation of background rates from shuffling procedures that do not preserve dinucleotides in our study resulted in an underestimate of the false-positive rate of up to 3-fold in our analysis.
we also confirmed that, for pairwise alignments, z-scores derived from assessing thermodynamic stability  contained comparable discriminatory potential to the best alignment-based tools applied to pairwise alignments, on the basis of precision at comparable recall.
although our shuffling constraints precluded analysis of more than pairwise alignments, it is possible that the false-positive rates of these tools on multiple alignments are higher than has been previously estimated, as their false-positive rates are usually estimated by shuffling mononucleotides  <cit>  or by simulating evolution without conserving dinucleotides  <cit> .
for rnaz, a possible explanation is that the score incorporates thermodynamic stability which is calculated by sampling from precomputed stabilities from sequences with similar length and mononucleotide composition .
the svm sampling could be extended to include a dinucleotide sampling which would likely alleviate any dinucleotide bias.
we are less certain of the explanation for the apparent dinucleotide biases evident for qrna  and evofold .
it seems unlikely to be a result of dinucleotide bias per se although qrna may have learned a gc bias from trnas in its test set.
we propose that these algorithms may suffer from not having been exposed to a large amount of nonselected mammalian genomic sequence as negative controls, as genomic sequence is known to contain substantial local variation in nucleotide content as well as a non-random distribution of other simple sequence features  <cit> .
our results underscore the importance of background modelling, which we believe to be one of the major unresolved issues in de novo ncrna searching in practice.
although our demonstration is based on pairwise alignments, multiple alignments should not be immune from these same trends.
unfortunately, the difficulty of preserving dinucleotides while shuffling increases substantially as the number of sequences in the alignment increases, even in the absence of other shuffling constraints.
we see no clean solution to this dilemma, although general strategies can be envisioned.
first, the shuffling problem itself deserves more detailed investigation.
our algorithm is only a simple heuristic and it might be improved by finding more efficient ways to explore the space of permutations.
another possibility might be to relax the constraint that every dinucleotide count in every alignment must be preserved.
this could be done by repeatedly shuffling the multiple alignment while keeping track of dinucleotide content until a shuffled variant sufficiently close to the native alignment is attained.
pollard et al.
<cit>  generated alignments in silico by simulating evolutionary sequence divergence; this method might be adapted to produce similar overall dinucleotide counts and alignment structures to test sequences.
an additional possibility is to use longer shuffling windows to shuffle multiple alignments.
however, most known structural ncrnas are relatively short , and precision and recall decrease substantially with tiling window sizes exceeding  <dig> nt .
unfortunately, altschul-and-erickson-type shuffling  <cit> , which retains dinucleotide counts in individual sequences but not the alignment structure, is inappropriate even for thermodynamic tools that do not depend on conserved structure for their overall score .
this is because independently shuffling the sequences within the alignment eliminates the dependency of the score contributions from the individual sequences to the overall score, and lowers the range of background scores .
independent shuffling is clearly inappropriate for tools that require conservation of structure.
an additional issue is that covariance-based tools are reliant on correct alignment as it is assumed that the covarying bases are structurally paired.
unlike protein coding regions, where independent structural evidence often exists to validate sequence alignments, aligning non-coding regions presents additional challenges.
pollard et al.
<cit>  reported that local alignment tools such as blastz  had higher specificity but lower sensitivity than global alignment tools  again using simulated evolutionary sequence divergence as a control for background.
this implies that although full genome alignments represent orthologous sequences, the exact nucleotide positions may not always be accurate, and thus may not provide an accurate covariance signature.
selection pressure to maintain or avoid conserved secondary structure depends on genomic location
despite these caveats, comparisons of the aggregate scores for real and shuffled alignments among different types of genomic features, including coding and non-coding parts of the genome, revealed striking differences with consequent biological implications.
in general, coding sequences appear to possess lower propensity for conserved secondary structure than would be expected by chance.
this is not an artefact of the high conservation of coding sequences; highly conserved sequences in fact tend to shuffle more completely  which would aid in discriminating bona fide structures.
on the basis of this result, it seems unlikely that we are dramatically overestimating the degree of selected secondary structure in the genome because at least one type of genomic feature scores as containing less than would be expected by chance.
since the random alignments exactly match the background expectation , we conclude that there is likely to be an evolutionary selection against formation of secondary structures in open reading frames.
to our knowledge, this trend has not yet been demonstrated on a genome-wide scale, although previous experimental data has shown that the introduction of stable secondary structure can indeed impact on translation efficiency  <cit> .
the most obvious rationale for such a selection is that inhibition of translation is, on the whole, detrimental to fitness of the organism.
in contrast, we find evidence that 3'utrs, introns, and intergenic regions on the whole contain slightly higher structural scores than corresponding shuffles.
although the biological significance of all but a small fraction of the highest-scoring individual sequences is questionable, our results does nevertheless support the existence of many additional rna structural elements in utrs and introns, and potentially additional classes of ncrnas in intergenic regions.
CONCLUSIONS
accurate prediction of rna structural elements in aligned genome sequence remains a difficult problem.
further work is required to properly assess background detection rates.
nonetheless, the general trends we observed in genomic scans appear to roughly reflect expectations, supporting the presence of secondary structural elements in utrs, and providing compelling evidence for evolutionary selection against secondary structures in coding regions.

 one of the major findings of genome sequencing has been that the primary sequence of roughly 5% of the human and mouse genomes is under purifying selection, indicating functionality. however, less than 2% is accounted for by mrna exons. the remaining 3% presumably encompasses cis-regulatory sequence, signals for transcriptional initiation, termination, rna processing, chromosomal features such as replication origins, and genes encoding ncrnas such as trna, snorna, mirna, and others. accurate computational identification of novel ncrna genes and mrna structural elements in genome sequence has proven to be more difficult than identification of exons, due generally to a limited and highly variable sequence signature. in bacteria, which have compact genomes, searching for transcription initiation signals, primary sequence conservation, and base composition have been fruitful approaches to de novo ncrna discovery; however, these features alone are unlikely to be sufficiently specific in large eukaryotic genomes. most, albeit not all, functional ncrna features possess some degree of secondary structure, either as part of the precursor or the functional rna itself. following the assumption that structural rna sequences should be more thermodynamically stable than random permutations of the same base composition, thermodynamic stability is an additional feature than can be incorporated into genomic searches for new ncrnas. major classes of structural rnas have lower Δg than corresponding shuffled sequences. it has been debated whether Δg is a sufficiently accurate discriminant when only a single sequence is analyzed; however, Δg has been proposed to be comparable or superior to more sophisticated algorithms when applied independently to segments of a pairwise alignment. what is clear is that disruption of dinucleotides in the random permutation dramatically affects the perceived precision of predictions, presumably because dinucleotide contributions are a key determinant of stability of an rna fold. covariance is also now a widely-accepted approach to ncrna discovery. a variety of recently-described ncrna search algorithms score for covariance to discriminate structural rna elements. success of covariance requires that sequences be sufficiently conserved to achieve a correct alignment, yet contain some nucleotide changes in order to assess compensatory mutations. an advantage of methods that do not utilize covariance is that they can identify structures common to sequences without high sequence similarity and apparently even sequences that fail to align at the primary sequence level. however, even considering both covariance and thermodynamic stability, some classes of ncrnas appear to be more difficult to detect than others. hybrid/oligoarrayaux.php tools/zmfold.pl to our knowledge, most ncrna search tools have not been assessed or compared systematically by an independent laboratory. moreover, previous studies using these tools for genomic scanning have acknowledged that the false-positive rate is high and cannot be determined accurately. here, with the goal of independently evaluating these tools for de novo discovery of ncrna elements in a genomic context, we have compiled an extensive sequence and alignment data set, and developed a new shuffling algorithm for pairwise alignments that simultaneously preserves key features of the alignment. we used these sequences, as well as real and shuffled genomic sequence, as input for a panel of published tools for discovery of novel ncrnas. we also evaluated whether different schemes utilizing only thermodynamic stability could discriminate real from shuffled ncrnas, and examined the output of a subset of tools on a genomic scale.  a test set for comparing ncrna-finding tools we tested the ncrna search tools on a test set derived by extracting genomic alignments corresponding to known eukaryotic ncrnas and mrna structural elements. the positive test set is available at and in additional files and and consists of mirna, trna, box h/aca snorna, box c/d snorna, snrna, utr regulatory element alignments, in addition to other structural ncrna alignments and non-structural rna alignments for a total of bases. we generated the alignments by mapping all available ncrna sequences from human, mouse, c. elegans, d. melanogaster, and s. cerevisiae, to available ucsc pairwise genomic alignments. to obtain negative controls, we created sequence shuffles of each sequence alignment using a new shuffling algorithm, shuffle-pair.pl, to simulate the problem of finding bona fide structural elements against a much larger background that resembles true genomic alignments in as many features as possible. briefly, the shuffle-pair.pl algorithm randomizes the columns of the alignment while preserving conservation, gaps, and dinucleotides in both sequences, all of which may contribute to scoring obtained from at least a subset of algorithms tested. it shuffles by swapping the middle bases between trinucleotides that are identical at positions and as described with the additional constraints of swapping in aligned trinucleotides and maintaining conservation and gaps of the alignment. the algorithm finds all swappable positions that preserve the criteria above and then randomly selects two bases to swap. it keeps track of positions that have been swapped and continues shuffling until it runs out of positions to shuffle. a detailed description of the algorithm and pseudocode is found in the methods section. on average, mirna alignments shuffled 57%, trna alignments – 61%, snrna alignments – 62%, snorna alignments – 46%, and regulatory elements – 43%; overall the alignments shuffled on average 53%. we did not observe a correlation between the degree of shuffling and the scores output by tools. we observed that the degree of shuffling possible depends on sequence length and conservation: increasing length and conservation increased the shufflability of the alignments, presumably due to more shuffling opportunities resulting from more bases and less constraint imposed by preserving local conservation patterns. conversely, increasing the number of aligned sequences greatly diminishes shuffling potential due the constraint of preserving dinucleotides in all aligned sequences. consequently, here we consider only findings pertaining to pairwise alignments. we acknowledge that we are not taking full advantage of msari, ddbrna, evofold, and rnaz, since these tools were designed to handle multiple sequence alignments. nonetheless, we reasoned that it would be of value to examine the outputs of these tools under a situation that controls for all variables believed to be relevant. we ran all of the tools with default settings. all but one of the tools ran without giving error messages. rnaz produces error messages with sequences shorter than and longer than nt and for alignments with very low sequence identity, but still produces scores that we retained for the analysis since only a small proportion of the test set falls into this category, and and challenges imposed in scoring these alignments apply to all of the tools. msari output was excluded from subsequent evaluation due to its extremely poor performance on alignments of sequences. overview of results on test set for each of the tools, we scored precision ) and recall ) across its range of output scores. in comparison to area under the roc curve, and specificity ), we view precision and recall as most relevant to real-world genomic screening: precision estimates the false-positive rate one would face if experimentally verifying predicted ncrnas, and recall estimates the proportion of all ncrnas that could be found. our results are consistent with expectation in several ways. first, taking an operating point at which each tool yields ~50% overall recall, the recall of different classes is ranked on average mirna > trna > regulatory > box h/aca snorna > snrna > other, roughly as has been described by others. algorithms utilizing only thermodynamic stability appear to have a relative advantage in detecting mirnas, presumably due to the simple hairpin structure. at the same time, they have relative difficulty identifying trnas, perhaps due to the more complicated fold. box c/d snornas and non-structural elements were virtually indistinguishable from shuffles. we note that some of the tools may have been trained or their parameters optimized using a portion of our test set, and this may slightly bias outcome. evofold, for example, was trained on a large set of ncrnas including many mirnas, and it does perform well on mirnas, even in our pairwise alignments; however, we confirmed that results for evofold are virtually identical if its training set rnas are omitted, presumably because it does not learn any structural or sequence information, and therefore results for the full data set are shown. nonetheless, on the whole, the tools appear to detect different ncrna classes in a roughly similar manner at the 50% recall operating point. at lower overall recall the results become less comparable due to a strong bias towards detection of mirnas by the thermodynamics-only tools. at higher recall, all tools suffer from low precision. our results also verify the recent demonstration of uzilov et al. that, on pairwise alignments, z-scores based on Δg scores for both strands yield substantially higher precision at roughly equivalent recall for all ncrna classes tested, taking either an overall average or a median of class averages. our "zrnafold " implementation is conceptually similar to "dynalign" used by uzilov. in our hands dynalign was too slow to be useful for genome-wide scanning and we thus did not test it. a less expected but intriguing result was that tools utilizing apparently unrelated algorithms often mutually identified and missed the same ncrnas. for example, evofold, one of the covariance tools, displayed a very significant overlap in true positives with all of the thermodynamic-only tools ). importance of dinucleotide frequency a somewhat discouraging outcome of our test was that all of the tools displayed low apparent precision in our task. at 50% recall, all of the tools displayed between and 20% precision. in a genomic scan, in which the background sequence most likely exceeds the 100× ratio we employed, this would translate to laboratory confirmations being overwhelmingly negative, in the absence of further filtering. the results of our analysis must be viewed in light of being based only on pairwise alignments, which substantially handicaps tools such as ddbrna, msari, and evofold that rely solely on covariance and were created solely to scan multiple alignments. nonetheless, we were surprised that rnaz, which has both covariance and thermodynamic components, did not generally outperform zrnafold, which is our relatively simple two-strand z-score implementation of the rnaz thermodynamic input. we reasoned that retention of dinucleotide counts might account for the high false-positive rate we observed. to ask whether this is the case, we examined the results obtained using a negative-control test set generated using an alternative shuffling procedure, shuffle-aln.pl, which preserves local conservation, gap structure, and base composition, but not dinucleotide content. indeed, with shuffle-aln.pl, the apparent false-positive rate decreased -fold on average; most drastically for rnaz and qrna, which dropped and -fold, respectively, at 50% recall. this indicates that virtually all tools, including those that do not explicitly have a thermodynamic component, are nonetheless influenced in some way by dinucleotide counts, which it is known should be accounted for in estimation of false-positive rates based on thermodynamic stability. considerations in genomic scanning: dinucleotides and other aspects of sequence content we next explored whether the dinucleotide issue might influence genomic scans. we ran five of the eight tools on 120-base tiling windows of human-mouse alignments of chromosome, and compared the distributions of scores with and without preservation of dinucleotide content, again by comparing shuffle-pair.pl with shuffle-aln.pl)., and alifoldz were not run on these sequences). in this analysis, the estimated false-positive rates drastically increased when using shuffles with conserved dinucleotides: in figure the scores obtained from shuffled sequences with dinucleotide frequency preserved nearly overlap those obtained from the original sequence, whereas there is a clear separation of both from the scores with mononucleotide shuffling. we reasoned that insufficient shuffling might provide a trivial explanation for this observation; however, when we subset the analysis to those shuffles which differed in at least 50% of their nucleotide positions from their original sequences, the trend clearly remained. even among the windows in which, by chance, the degree of shuffling in the dinucleotide-constrained set exceeded that in the mononucleotide-shuffled set, the trend remained largely intact. why do dinucleotide counts influence scores from these tools? for those relying on thermodynamics, the dependence of Δg on dinucleotides provides a likely explanation. however, qrna does not utilize thermodynamic stability. indeed, all the tools we examined displayed negligible biases towards specific nucleotides or dinucleotides when run on randomly-generated alignments. only on ncrnas and genomic sequence, which contain non-random distributions of di- and even mononucleotides, did high-scoring sequences emerge with enrichment of specific dinucleotides. for example, in shuffled versions of the test set, the highest-scoring sequences from rnaz and qrna tended to be slightly enriched for cc and gg, and depleted of aa and at. this trend was even more marked for shuffled versions of chromosome alignments. evofold appeared to be even more sensitive to dinucleotide composition of the shuffled genomic alignments, favoring those with aa, tt, at, and ta. manual examination of the specific shuffled genomic sequences that scored highly with evofold revealed that many of them contained unusual distributions of homopolymeric or near-homopolymeric runs and simple repeat-like sequence that are extremely rare in randomly-generated sequence, but very common in the human genome. because the mono-, di-, and/or tri-nucleotide contents have reduced complexity, these patterns are occasionally preserved in shuffled controls. in fact, we observed striking correlations between scores for all of the chromosome tiling windows before and after shuffling, suggesting that simple sequence properties can and do influence scores emerging from these tools. in an effort to gain further insight, we examined the highest-scoring ~ real genomic sequences from chromosome for each of three tools utilizing covariance. these are posted in the supplementary data. manual examination of all three lists revealed the presence of many homopolymeric or near-homopolymeric runs, repetitive short motifs that escape repeatmasking, and nonuniform distribution of nucleotides. the sequences that rnaz and qrna scored most highly were clearly enriched for g and c while those scored most highly by evofold appeared enriched for homopolymeric stretches of a and t. anecdotally, these non-repeat-masked features do appear to contribute to the predicted rna structures. while we cannot rule out that a subset of these are bona fide undiscovered functional ncrna features, it is also possible that the non-random distribution of nucleotides and short sequence features in genomes could lead to a higher degree of apparent secondary structure than would be expected from a completely random evolutionary process. we currently have no unifying explanation for why these tools appear to favour specific types of sequences, and suggest that this subject might form the basis for further study. these trends do seem to suggest that caution should be exercised in drawing conclusions about the number of ncrna features in the human genome, on the basis of results of computational genomic scans. from our analyses, we conclude that an accurate assessment of state-of-the-art ncrna search tools on multiple alignments, which was our initial objective, is not realistic in the absence of a shuffling algorithm that retains at least an approximation of both dinucleotide frequencies and real alignment structure. moreover, it may be worth considering the fact that genome sequences contain a non-random distribution of mononucleotides and other short sequence features. genomic scanning for novel conserved rna structural elements reveals enrichment for stable secondary structure in intergenic regions, introns, and utrs, and selection against structure in coding regions our analysis above indicates that with any of the tools we tested, the vast majority of high-scoring windows in a genome scan using pairwise alignments would be false-positives. nonetheless, we next asked whether use of shuffle-pair.pl, which we believe provides a realistic estimate of the false-positive rate, could aid in making new observations regarding frequency of structural elements across genomic features. we used a tiling strategy similar to that employed previously. we used a window size of nt because this window size resulted in the greatest sensitivity at little cost to specificity relative to other tiling window sizes. this window size is within range of previously computed genomic scans: and used in bacterial scans, and used for the human genome. on average, 23% of nucleotides in a base window of human-mouse pairwise alignments change identity during shuffling; however, 13% of sliding windows have >50% shuffling. we scanned pairwise genomic alignments of human chromosome and genomic regions corresponding to alignable portions of refseq genes with either zmfold or evofold, as practical representatives of the thermodynamic-stability-only regime and the covariance regime. we scanned real and shuffled counterparts of each window. we found that the distributions of the real and shuffled scores largely overlap, with few or no real scores completely exceeding the shuffled distribution. however, the real and shuffled score distributions were significantly different. a histogram of the difference in number of real and shuffled sequences at different scores is shown at the right in figure 7a. these results obviously bring into question the significance of any individual element, but also suggest a selection for some level of conserved secondary structure across the genome. to ask whether results obtained from genomic scanning appear to be biologically relevant on the whole, and to gain evidence that we have not mis-calibrated the negative controls, we examined the overall distribution of score differences between real and shuffled sequences obtained from different types of genomic features, which might be expected to display differing propensities towards secondary structure. using either zmfold or evofold, we found that striking differences exist among results obtained from these features. in particular, introns, intergenic sequence and 3' utrs tend to have a higher distribution of scores, indicating enrichment for structural features, relative to shuffled controls. in contrast, we observed marginal enrichment within 5' utrs. strikingly, we found that coding regions have a lower score distribution than shuffled controls, indicating an overall depletion of rna structures in open reading frames, which is consistent with the possibility that these elements would inhibit translation. these results indicate that our analysis is not, on average, biased towards falsely reporting structures, since some genomic features appear to be enriched for structures while others are depleted or neutral. these results also provide evidence for the presence of biologically functional structural features, although specificity limitations of the individual tools preclude direct prediction of individual structural elements. to further investigate whether the tools are in fact identifying a biologically meaningful selection for structural elements we asked whether zmfold and evofold, which represent very different approaches, tend to give high scores to the same utr segments. we found that in 3'utrs, where an obvious enrichment for structural elements exists, the top zmfold and evofold scores are correlated. furthermore, 48% of evofold and zmfold highest scores within each utr correspond to the same region of the utr which is approximately two-fold higher than expected by chance. discussion to our knowledge, the panel of ncrna search tools examined here have not been examined systematically by an independent laboratory. with the goal of evaluating how useful these tools are for de novo discovery of ncrna elements in a genomic context, in which the anticipated sparseness of real features makes it particularly important to control false discovery, we compiled an extensive sequence and alignment data set comprised of fragments of genomic alignments corresponding to ncrnas from a variety of eukaryotes. we also developed a shuffling algorithm for pairwise alignments that maintains dinucleotide frequencies, in addition to other features typically maintained in alignment shuffling. due to shuffling constraints, which we show are important, we limited our test to pairwise alignments. in our test, we found that failure to account for background sequence properties leads to over-estimation of precision for most tools. we were surprised to find that tools with a strong or exclusive covariance component were sensitive, in at least some situations, to dinucleotide counts, which should primarily impact thermodynamic stability. estimation of background rates from shuffling procedures that do not preserve dinucleotides in our study resulted in an underestimate of the false-positive rate of up to 3-fold in our analysis. we also confirmed that, for pairwise alignments, z-scores derived from assessing thermodynamic stability contained comparable discriminatory potential to the best alignment-based tools applied to pairwise alignments, on the basis of precision at comparable recall. although our shuffling constraints precluded analysis of more than pairwise alignments, it is possible that the false-positive rates of these tools on multiple alignments are higher than has been previously estimated, as their false-positive rates are usually estimated by shuffling mononucleotides or by simulating evolution without conserving dinucleotides. for rnaz, a possible explanation is that the score incorporates thermodynamic stability which is calculated by sampling from precomputed stabilities from sequences with similar length and mononucleotide composition. the svm sampling could be extended to include a dinucleotide sampling which would likely alleviate any dinucleotide bias. we are less certain of the explanation for the apparent dinucleotide biases evident for qrna and evofold. it seems unlikely to be a result of dinucleotide bias per se although qrna may have learned a gc bias from trnas in its test set. we propose that these algorithms may suffer from not having been exposed to a large amount of nonselected mammalian genomic sequence as negative controls, as genomic sequence is known to contain substantial local variation in nucleotide content as well as a non-random distribution of other simple sequence features. our results underscore the importance of background modelling, which we believe to be one of the major unresolved issues in de novo ncrna searching in practice. although our demonstration is based on pairwise alignments, multiple alignments should not be immune from these same trends. unfortunately, the difficulty of preserving dinucleotides while shuffling increases substantially as the number of sequences in the alignment increases, even in the absence of other shuffling constraints. we see no clean solution to this dilemma, although general strategies can be envisioned. first, the shuffling problem itself deserves more detailed investigation. our algorithm is only a simple heuristic and it might be improved by finding more efficient ways to explore the space of permutations. another possibility might be to relax the constraint that every dinucleotide count in every alignment must be preserved. this could be done by repeatedly shuffling the multiple alignment while keeping track of dinucleotide content until a shuffled variant sufficiently close to the native alignment is attained. pollard et al. generated alignments in silico by simulating evolutionary sequence divergence; this method might be adapted to produce similar overall dinucleotide counts and alignment structures to test sequences. an additional possibility is to use longer shuffling windows to shuffle multiple alignments. however, most known structural ncrnas are relatively short, and precision and recall decrease substantially with tiling window sizes exceeding nt. unfortunately, altschul-and-erickson-type shuffling, which retains dinucleotide counts in individual sequences but not the alignment structure, is inappropriate even for thermodynamic tools that do not depend on conserved structure for their overall score. this is because independently shuffling the sequences within the alignment eliminates the dependency of the score contributions from the individual sequences to the overall score, and lowers the range of background scores. independent shuffling is clearly inappropriate for tools that require conservation of structure. an additional issue is that covariance-based tools are reliant on correct alignment as it is assumed that the covarying bases are structurally paired. unlike protein coding regions, where independent structural evidence often exists to validate sequence alignments, aligning non-coding regions presents additional challenges. pollard et al. reported that local alignment tools such as blastz had higher specificity but lower sensitivity than global alignment tools again using simulated evolutionary sequence divergence as a control for background. this implies that although full genome alignments represent orthologous sequences, the exact nucleotide positions may not always be accurate, and thus may not provide an accurate covariance signature. selection pressure to maintain or avoid conserved secondary structure depends on genomic location despite these caveats, comparisons of the aggregate scores for real and shuffled alignments among different types of genomic features, including coding and non-coding parts of the genome, revealed striking differences with consequent biological implications. in general, coding sequences appear to possess lower propensity for conserved secondary structure than would be expected by chance. this is not an artefact of the high conservation of coding sequences; highly conserved sequences in fact tend to shuffle more completely which would aid in discriminating bona fide structures. on the basis of this result, it seems unlikely that we are dramatically overestimating the degree of selected secondary structure in the genome because at least one type of genomic feature scores as containing less than would be expected by chance. since the random alignments exactly match the background expectation, we conclude that there is likely to be an evolutionary selection against formation of secondary structures in open reading frames. to our knowledge, this trend has not yet been demonstrated on a genome-wide scale, although previous experimental data has shown that the introduction of stable secondary structure can indeed impact on translation efficiency. the most obvious rationale for such a selection is that inhibition of translation is, on the whole, detrimental to fitness of the organism. in contrast, we find evidence that 3'utrs, introns, and intergenic regions on the whole contain slightly higher structural scores than corresponding shuffles. although the biological significance of all but a small fraction of the highest-scoring individual sequences is questionable, our results does nevertheless support the existence of many additional rna structural elements in utrs and introns, and potentially additional classes of ncrnas in intergenic regions.  accurate prediction of rna structural elements in aligned genome sequence remains a difficult problem. further work is required to properly assess background detection rates. nonetheless, the general trends we observed in genomic scans appear to roughly reflect expectations, supporting the presence of secondary structural elements in utrs, and providing compelling evidence for evolutionary selection against secondary structures in coding regions. 
current algorithms identify and score potential rna secondary structures on the basis of thermodynamic stability, conservation, and/or covariance in sequence alignments. accurate identification of novel, functional noncoding rna features in genome sequence has proven more difficult than for exons. accurate prediction of novel rna structural elements in genome sequence remains a difficult problem, and development of an appropriate negative-control strategy for multiple alignments is an important practical challenge. we developed a shuffling algorithm, shuffle-pair.pl, that simultaneously preserves dinucleotide frequency, gaps, and local conservation in pairwise sequence alignments. nonetheless, the general trends we observed for the distributions of predicted ncrnas across genomic features are biologically meaningful, supporting the presence of secondary structural elements in many 3' utrs, and providing evidence for evolutionary selection against secondary structures in coding regions. we used shuffle-pair.pl to assess precision and recall of six ncrna search tools was significantly higher for real than shuffled sequence, while the distribution for coding sequences was lower than that of corresponding shuffles. furthermore, due to issues in modelling background signal, it has been difficult to gauge the precision of these algorithms on a genomic scale, in which even a seemingly small false-positive rate can result in a vast excess of false discoveries.
 accurate identification of novel, functional noncoding rna features in genome sequence has proven more difficult than for exons. current algorithms identify and score potential rna secondary structures on the basis of thermodynamic stability, conservation, and/or covariance in sequence alignments. neither the algorithms nor the information gained from the individual inputs have been independently assessed. furthermore, due to issues in modelling background signal, it has been difficult to gauge the precision of these algorithms on a genomic scale, in which even a seemingly small false-positive rate can result in a vast excess of false discoveries.  we developed a shuffling algorithm, shuffle-pair.pl, that simultaneously preserves dinucleotide frequency, gaps, and local conservation in pairwise sequence alignments. we used shuffle-pair.pl to assess precision and recall of six ncrna search tools was significantly higher for real than shuffled sequence, while the distribution for coding sequences was lower than that of corresponding shuffles.  accurate prediction of novel rna structural elements in genome sequence remains a difficult problem, and development of an appropriate negative-control strategy for multiple alignments is an important practical challenge. nonetheless, the general trends we observed for the distributions of predicted ncrnas across genomic features are biologically meaningful, supporting the presence of secondary structural elements in many 3' utrs, and providing evidence for evolutionary selection against secondary structures in coding regions. 