BACKGROUND
rna molecules are essential to many functions in the cells of all organisms.
for example, these molecules are involved in gene translation and also act as catalysts and as regulators of gene expression  <cit> .
because function is determined by molecular structure, there is significant investment in computational methods for predicting rna secondary structure, which in turn is useful for inferring tertiary structure  <cit> .
our interest is in assessing the merits of some recent advances in secondary structure prediction in statistically robust ways.
we focus on thermodynamically informed approaches for predicting pseudoknot free secondary structures from the base sequence.
a widely used method finds the minimum free energy  structure with respect to the nearest neighbour thermodynamic model of mathews, turner and colleagues  <cit> .
some recent advances in secondary structure prediction are the new maximum expected accuracy  and maximum pseudo-expected accuracy  methods of lu et al.
<cit>  and hamada et al.
<cit> .
these approaches generally maximize  expected base pair accuracy as a function of base pair probabilities calculated using a partition function method and have higher average accuracy than the mfe algorithm on the turner and andronescu et al.
energy parameters.
knudsen et al.
<cit>  and do et al.
<cit>  also presented some rna secondary structure prediction methods based on probabilistic models of structures.
but since their probabilistic approaches are not determined using a thermodynamic model, we don't include their methods in our later comparisons.
another advance is estimation of new energy parameters from both thermodynamic and structural data using state-of-the-art estimation techniques.
andronescu et al.
<cit>  derived two parameter sets by inference from energies that were derived from optical melting experiments as well as from structural data.
the two energy parameter sets are called bl* and cg*, named after the boltzmann likelihood and constraint generation methods used to infer them.
these parameter sets have yielded significant improvements in prediction accuracy of the mfe method, compared with the turner parameters, with the bl* parameters being slightly better than the cg* parameters.
here and throughout, the accuracy of a prediction refers to its f-measure, which is the harmonic mean of sensitivity and positive predictive value .
all of this work assesses algorithm accuracy on specific classes of rnas, such as introns or transfer rnas, as well as overall average accuracy on rnas taken over all such classes.
this recent work motivates the following questions.
are comparisons of the mea-based and mfe approaches on specific rna classes reliable when the size of available datasets is small?
do the mea- or pseudo-mea-based approaches produce significantly more accurate predictions than mfe on the latest energy parameter sets?
what is the best combination of algorithm and thermodynamic model?
to answer these questions, we report on the accuracy of both mea-based and mfe methods with respect to two versions of the turner parameters as well as the recent bl* and cg* parameters of andronescu et al., on datasets for specific rna classes as well as large datasets that combine multiple rna classes.
we present three main findings.
first, we show that f-measure accuracies on our large datasets are likely to be reliable estimates of accuracy of a population represented by such sets, in the sense that high-confidence interval widths for f-measure obtained using the bootstrap percentile method are within a small, 2% range.
average accuracy on smaller classes is less reliable.
for example, confidence intervals for both mea and mfe have an 8% range on a class of  <dig> group i introns that has been used previously in benchmarking algorithms.
second, there is a clear "winner" in terms of overall prediction accuracy, namely the pseudo-mea-based method of hamada et al.
<cit> .
however, the relative accuracy of the mfe and mea-based approaches depends on the underlying energy parameters: using a permutation test we find that, at a statistically significant level, the accuracy of mfe-based prediction on our large datasets is better on two of the four energy parameter sets that we consider, while mea-based prediction is better than mfe-based prediction on a third parameter set.
finally, both mea-based and mfe methods achieve the highest accuracy when using the fourth parameter set we consider, namely the bl* energy parameters of andronescu et al.
<cit> .
methods
in this section we first describe the datasets, thermodynamic models, and algorithms considered in this paper.
we then describe the accuracy measures and statistical methods used in our analyses.
datasets
we use three datasets, as follows
• s-full is a comprehensive set of  <dig>  rna sequences and their secondary structures that has been assembled from numerous reliable databases  <cit> .
sequences in this and our other datasets have length at most  <dig> nucleotides; in some cases these were derived by partitioning larger sequences such as 16s ribosomal rna sequences.
the average length of sequences in s-full is 270nt.
• mt was used by lu et al.
<cit>  in their study of the mea algorithm and contains rnas from the following classes: 16s ribosomal rna, 23s ribosomal rna, 5s ribosomal rna, group i intron, group ii intron, ribonuclease p rna, signal recognition rna and transfer rna.
• ma is a subset of the s-full set, containing exactly those sequence of s-full that are in the rna classes included in mt.
we formed the ma dataset in order to compare our algorithms on the same classes as did lu et al., while using sets of rnas from these classes with as large a size as possible.
an overview of different rna classes in the mt and ma data sets is shown in table  <dig>  this table presents the number of rna sequences, the mean and standard deviation of their lengths and the average of normalized similarity between rna structures for each rna class in these two data sets .
overview of different rna classes in the mt and ma data sets, including the number of rnas, the mean and standard deviation of the length and the average normalized similarity between rna sequences in these two data sets for each rna class.
lu et al.
reported their results on a restricted format of the mt set in which certain bases of some sequences are in lower case, indicating that the base is unpaired in the reference structure.
they have used this structural information in their predictions.
however, we have not employed this information, and therefore our accuracy measures on the mt set are different from those of lu et al.
<cit> .
thermodynamic models
a thermodynamic model consists of features--small structural motifs such as stacked pairs--plus free energy change parameters, one per feature.
the first model we use, the turner model  <cit>  called "turner99", is the most widely used energy model for prediction of rna secondary structures.
the model has over  <dig> features, which are based on turner's nearest neighbor rules and reflect the assumption that the stability of a base pair or loop depends on its sequence and on the adjacent base pair or unpaired bases.
the model is additive in that the overall free energy change of a secondary structure for a given sequence is the sum of the free energy changes for features of the structure.
the parameters of the model were derived from optical melting experiments, the most commonly used experimental approach to determine the free energy change of rna structures.
we also consider variants of the turner model, used by andronescu et al.
<cit> .
the t99-multirnafold  model is derived from the turner <dig> model but includes only  <dig> features.
parameters for features in the turner <dig> model can be obtained by extrapolation from the parameters of the t99-mrf model.
using maximum likelihood and constraint optimization methods, andronescu et al.
<cit>  derived new free energy change parameters for these  <dig> features; the resulting models are called bl*  and cg*  respectively.
we used all three models, namely t99-mrf, bl* and cg*, in this work, in order to assess the dependence of algorithm accuracy on model parameters.
we note that in lu et al.
's work  <cit> , the parameters of a newer version of turner model, called turner <dig>  were used for one of the structural motifs, coaxial stacking.
however, for the rest of structural motifs the parameters of turner <dig> model were engaged.
so, we also call the parameter set used for the lu et al.
's benchmark turner <dig> since most of its parameters are turner <dig> ones and also the coaxial stacking motif is not employed in the other model that we study.
zakov et al.
<cit>  also obtained parameters that improve rna structure prediction.
but we don't consider their parameters in this work, since those are not applicable for the partition function calculation and therefore for the probability calculation required for the mea method.
algorithms
we analyze four rna secondary structure prediction algorithms.
the first predicts secondary structures that have minimum free energy  with respect to a given thermodynamic model.
the second is the maximum expected accuracy  algorithm as introduced by lu et al.
<cit> , which maximizes expected base pair accuracy as a function of base pair probabilities calculated using a partition function method.
we implemented the mea algorithm for use with the multirnafold models.
as a result, we worked not only with the algorithms of lu et al., which we refer to as rsmfe and rsmea, but also with the mfe algorithm of andronescu et al., referred to as ubcmfe, and a new implementation of mea which we developed, called ubcmea.
lu et al.
's benchmark  <cit>  showed that rsmea gives the best prediction accuracy when its γ parameter  is equal to  <dig>  accordingly, we also set γ to be  <dig> in ubcmea.
the third algorithm that we analyze is the generalized centroid estimator method of hamada et al.
<cit> .
this is similar to the mea method, but uses a somewhat different objective function, namely a gamma-centroid estimator, to infer a structure from base pair probabilities.
this method, which we refer to as gc-g <dig>  employs a parameter γ to balance sensitivity and positive predictive value; based on the results of hamada et al., we set γ =  <dig>  the fourth method is another algorithm of hamada et al.
<cit>  that generalizes the centroid estimator by using a pseudo-expected accuracy maximization technique that automatically selects γ to balance sensitivity and positive predictive value for a given input.
we refer to this algorithm as gc-pmfmeas .
the rsmea and rsmfe algorithms always use the parameters of turner <dig> model as their free energy parameter set, while ubcmea and ubcmfe use parameter sets in the multirnafold model format, namely the bl*, cg* and t99-mrf sets.
the gc-g <dig> and gc-pmfmeas algorithms also employ bl* parameters using a turner <dig> format.
accuracy measures
we use three measures for determining the structural prediction accuracy, namely sensitivity , positive predictive value or ppv , and f-measure, which combines the sensitivity and ppv into a single measure.
sensitivity is the ratio of correctly predicted base pairs to the total base pairs in the reference structures.
ppv is the fraction of correctly predicted base pairs, out of all predicted base pairs.
f-measure is the harmonic mean of the sensitivity and ppv.
this value is equal to the arithmetic mean when sensitivity and ppv are equal.
however, f-measure becomes smaller than the arithmetic mean as one of the numbers approaches  <dig> .
the possible values for these three measures are between  <dig> and 1; the closer to  <dig>  the better prediction.
the f-measure is widely used measure in the literature; it is also the common measure in the studies by hamada et al.
and by lu et al., to which refer in our study.
mostly to facilitate comparison of their results to ours, we decided to use the f-measure rather than the matthews correlation coefficient .
sensitivity=numberofcorrectlypredictedbasepairsnumberofbasepairsinthereferencestructureppv=numberofcorrectlypredictedbasepairsnumberofpredictedbasepairsf - measure=2×sensitivity×ppvsensitivity+ppv
throughout this paper, we calculate three types of averages for a given measure "m" , namely unweighted averages, weighted averages and s-weighted averages, defined below.
the weighted average counts each sequence equally, regardless of which class it belongs to.
the unweighted average, on the other hand, counts each class equally and was used by lu et al  <cit> .
a potential problem with the unweighted average is that an rna class with many highly similar sequences can have disproportionate influence on the overall accuracy, relative to its sequence diversity.
therefore, we introduce the s-weighted average, which takes into account the similarities between rna sequences in each rna class and gives a weight to each one according to its average normalized similarity, in such a way that classes with highly similar sequences have lower weight.
the three averages are defined as follows when there are n rna classes, c <dig>  c <dig>  ..., cn, with cardinalities l <dig> l <dig>  ...,ln, respectively.
for the remainder of our study, n =  <dig>  and the classes are those listed in table 1).
unweightedaverageofm=1n∑∀c1i∈c1ml1+∑∀c2i∈c2ml2+⋯+∑∀cni∈cnmln,
weightedaverageofm=∑∀c1i∈c1m+ ∑∀c2i∈c2m+...+ ∑∀cni∈cnml1+l2+...+ln,
and
s - weightedaverageofm=l1-s1⋅∑∀c1i∈c1m+l2-s2⋅∑∀c2i∈c2m+...+ln-sn⋅∑∀cni∈cnml11-s1+l21-s2+...+ln1-sn,
where si is the mean of the normalized similarities measured between the reference structures of any two of rna sequences in the corresponding rna class ci.
our normalized similarities were computed using the simtree procedure by eden et al.
<cit> , which takes into account secondary structure similarities in addition to sequence similarities.
the average normalized similarities are reported in table  <dig> and always fall between zero and one, where an average normalized similarity of one means that all sequences in the set are identical.
the s-weighted average equals the weighted average when si is zero for all rna classes, and it equals the unweighted average when si is one for all classes.
for the remaining cases, the value of the s-weighted average would fall between that of the unweighted and weighted averages.
bootstrap percentile confidence intervals
we calculated the confidence intervals of average f-measure for various rna classes, by using the bootstrap percentile method  <cit>  following the recent work by aghaeepour and hoos  <cit> .
we first took  <dig> resamples with replacement from the original f-measure values obtained by the prediction method under investigation on the given reference dataset, where all resamples had the same size as the original sample.
we then computed the average f-measures of the resamples and employed these into the bootstrap distribution.
finally, we determined the lower and upper limits of the 95% bootstrap percentile confidence interval as the  <dig> th and  <dig> th percentile of the bootstrap distribution, respectively.
these calculations were all performed using the "boot" package of the r statistics software environment  <cit> .
we verified that the bootstrap distributions are close to gaussian, using the anderson-darling test, which indicates that the bootstrap percentile intervals can be expected to be reasonably accurate  <cit> .
permutation test method
to assess the statistical significance of the observed performance differences, we used permutation tests, following aghaeepour and hoos  <cit> .
since lu et al.
<cit>  and hamada et al.
<cit>  reported that the mea-based methods outperform mfe, we used a one-sided permutation test  <cit>  to determine whether mea-based methods have significantly better accuracy than mfe on our parameter sets.
the test that we applied proceeds as follows.
first, we calculated the difference in means between sets of f-measure values obtained by the two given structure prediction procedures, an mea-based method and mfe.
for simplicity, we call these two sets a and b and we denote their sizes as na and nb, respectively.
then we combined the f-measure values of sets a and b. next, we calculated and recorded the difference in sample means for  <dig> randomly chosen ways of dividing these combined values into two sets of size na and nb.
the p-value was then calculated as the proportion of the sample means thus determined whose difference was less than or equal to that of the means of sets a and b. then, if the p-value is less than the 5% significance level, we reject the null hypothesis that mfe and the mea-based method have equal accuracy and thus accept the alternative hypothesis that the mea-based method has significantly better accuracy than mfe.
otherwise, we cannot reject the null hypothesis and therefore we cannot accept the alternative hypothesis.
furthermore, to assess whether the difference in accuracy between the mea-based methods and the mfe method on a given parameter set is significant, we performed a two-sided permutation test.
this test works exactly like the one-sided permutation test, except that its p-value is calculated as the proportion of the sampled permutations where the absolute difference was greater than or equal to that of absolute difference of the means of sets a and b. then, if the p-value of this test is less than the 5% significance level, we reject the null hypothesis that mfe and mea have equal accuracy, otherwise, we cannot reject the null hypothesis.
all of these calculations were performed using the "perm" package of the r statistics software environment.
RESULTS
in this section, we investigate to which degree the prediction accuracy of the energy-based methods studied in this paper is dependent on different datasets and different thermodynamic parameter sets used for prediction.
we start by considering how the size of a dataset will influence the accuracy achieved by rna secondary structure prediction methods.
we then study the dependency of prediction accuracy of the mfe, mea, gc-g <dig>  and gc-pmfmeas algorithms on the thermodynamic parameter sets that they use.
dependency of the energy-based methods on data characteristics
results of the energy-based methods on different data sets
because measures of accuracy on reference datasets are used to assess the quality of prediction achieved by various algorithms, it is important to understand how such measures vary depending on the reference datasets used.
later in this section we also consider how accuracy measures vary depending on the energy model used.
the baseline data in tables  <dig> and  <dig> show that there can be significant differences in accuracy of a given algorithm on rna classes within the ma versus mt sets.
for example, on ribonuclease p rna, when using the bl* parameter set, ubcmea achieves an f-measure of  <dig>  on the mt dataset  and  <dig>  on the ma set , an absolute difference of about 17%.
using data from tables  <dig> and  <dig>  figure  <dig> illustrates the difference between accuracy measures on the mt versus ma datasets on both the ribonuclease p rna and the group i intron classes, for the ubcmea, ubcmfe and gc-pmfmeas algorithms, using the bl* parameters, as well as for rsmea and rsmfe using the turner <dig> parameters.
in the figure, the further the points are from the dotted diagonal line, the larger the difference .
of the ten data points shown, only one represents a difference smaller than  <dig> %.
these data clearly indicate that, based only on prediction accuracy and without further statistical analysis, one cannot draw meaningful conclusions about the average accuracy of a particular algorithm on the overall population of an rna class, based only on the average accuracy on currently available datasets.
bootstrap confidence intervals for the prediction accuracy of the energy-based methods on different rna classes
following the work of aghaeepour and hoos  <cit> , we use bootstrap percentile confidence intervals to assess the accuracy of the average f-measures obtained by a secondary structure prediction procedure on a given rna dataset; these provide a measure of the overall average accuracy on the whole population from which the dataset is drawn .
we chose to calculate confidence intervals using weighted average f-measure.
we note from the data of table  <dig> that on the ma dataset, all three averages  are qualitatively similar in the sense that if the f-measure accuracy of one algorithm is better than another with respect to one average, then the same is true with respect to the other averages.
thus, we would expect the same qualitative conclusions if we had used a different average.
confidence intervals are obtained by the bootstrap percentile method for the ma and s-full sets and different rna classes in ma  when ubcmea and ubcmfe use the bl* parameter set and rsmea and rsmfe use the turner <dig> parameter set and the confidence level  is set to  <dig> .
first, the confidence intervals of all algorithms on the ma and s-full sets have a width of at most  <dig> , indicating that the average accuracy measured on these sets is likely to be a good estimate - within 1% - of the accuracy of a population of rna molecules represented by these sets.
however, the interval widths on individual classes can much higher, e.g.,  <dig>  for ubcmea on the group i intron class, suggesting that average accuracy is not a reliable measure of a method's overall accuracy on such classes.
we note that, without the use of rigorous statistical methods, a  <dig>  difference in accuracy is considered significant in some related work  <cit> .
second, the confidence interval widths of rna classes do not strictly decrease as the size of the class increases .
for example, as shown in table  <dig>  for the ubcmea  algorithm the confidence interval width of the ribonuclease p rna class is  <dig>   less than the interval width for the transfer rna class, even though the transfer rna class contains roughly  <dig>  times as many rnas as are in the ribonuclease p rna class.
thus factors specific to the classes, other than class size, must account for confidence interval widths.
one possibility is that classes with greater similarity among structures would have smaller confidence intervals, since in the limit, if all sequences in a class are identical then the confidence interval has zero width.
however, data from tables  <dig> and  <dig> did not support such a correlation between average normalized similarity and confidence interval width, even for classes not too different in size.
for example, the ribonuclease p and transfer rna classes in the ma set have relatively similar sizes , with the ribonuclease p rna class having lower normalized similarity  than the transfer rna class , as one might expect; yet, for ubcmea, the confidence interval for the ribonuclease p rna class is narrower  than those for the ribonuclease p rna class .
understanding class-specific factors underlying differences in prediction accuracy remains a goal for further study.
third, for any given dataset , the confidence interval width for the ubcmfe algorithm is close to that for ubcmea and gc-pmfmeas, even though the location of the interval  can be quite different.
for example, the average f-measure of ubcmea and ubcmfe on the ribonuclease p rna class differ by  <dig> , but the confidence interval widths are identical.
as another example, the smallest and largest confidence interval widths of the 16s ribosomal rna class among all four algorithms are, respectively,  <dig>  and  <dig> , which indicates a difference of less than  <dig>  in the interval widths for this class .
figure  <dig> shows the width of confidence intervals versus the size of different rna classes for the ubcmea, ubcmfe and gc-pmfmeas methods.
as one can see, almost all of the points for these three methods are coincident, indicating that ubcmea, ubcmfe, and gc-pmfmeas have identical width in the corresponding sizes.
the other important observation from this figure is that, as one could expect, the width of the confidence interval tends to decrease with increasing number of rnas of a given type in the reference set, but  there are notable exceptions to this trend for sets of size below  <dig>
the final interesting observation from figure  <dig> and table  <dig> is that the gc-pmfmeas method using the bl* parameter set outperforms all the other methods on most rna classes and also on our large ma and s-full sets.
hamada et al.
<cit>  also showed that on a small set of  <dig> rnas, gc-pmfmeas achieves better prediction accuracy than mfe and the other mea-based methods when using the bl* parameter set.
our results confirm their finding on our large ma and s-full data sets.
similar graphs for the rsmea and rsmfe algorithms are provided in the additional files.
additional files  <dig> and  <dig>  figures s <dig> and s <dig> support findings similar to those reported for ubcmea, ubcmfe and gc-pmfmeas.
since the overall results on the gc-g <dig> algorithm were roughly the same as the results of the ubcmea method, we report its prediction accuracies only in tables  <dig> and  <dig> and did not study it further.
dependency of the energy-based methods on thermodynamic parameters
the accuracy of algorithmic methods achieved by energy-based secondary structure prediction approaches, such as gc-pmfmeas, mea and mfe, depends on the underlying thermodynamic parameter sets.
lu et al.
showed that the accuracy of mea is better than that of mfe when using the turner <dig> parameter set .
but does mea or gc-pmfmeas outperform mfe on other published parameter sets?
if so, is the difference in accuracy statistically significant?
we address these questions in the following.
the table presents the prediction accuracy of different algorithms with different thermodynamic sets in terms of f-measure.
the 95% percentile confidence intervals of their accuracies are also shown in parentheses.
the parameter set t99-mrf refers to the turner <dig> parameters in multirnafold format.
bl* and cg* are the parameter sets obtained by the bl and cg approaches of andronescu et al.
<cit> , respectively.
also, the turner <dig> parameter set is the parameter set obtained by mathews et al.
<cit> .
"n/a" indicates cases in which a given algorithm is not applicable to a parameter set, since that does not match the energy model underlying the algorithm.
the highest accuracies for mea and mfe are shown in bold.
in contrast, the relative accuracy of mea versus mfe varies, with mfe being slightly better than mea in two cases, mea being slightly better than mfe in one case, and both having approximately the same accuracy in one case, namely bl*.
to assess whether mea has significantly better accuracy than mfe on our parameter sets, we first performed the one-sided permutation test  with the alternative hypothesis that mea yields higher prediction accuracy than mfe.
as recorded in table  <dig>  the resulting p-value falls below the standard significance threshold of  <dig>  only for turner <dig>  indicating that only for this parameter set, mea outperforms mfe, and its performance, only on this parameter set, can be considered better than mfe.
in contrast, the p-values obtained for the other parameter sets, namely cg*, t99-mrf and bl* are above  <dig> , indicating that we cannot conclude that mea outperforms mfe for these parameters sets.
p-value of one-sided permutation test is computed for the mea versus mfe algorithm on different parameter sets; p-values below the standard significance level of  <dig>  are shown in bold face and indicate significantly higher prediction accuracy of mea compared to mfe.
we next performed the two-sided permutation test on all mentioned parameter sets.
the results, as recorded in table  <dig>  indicate that for the bl* parameter set, the performance difference between mfe and mea is not statistically significant , while the performance differences for the other two parameter sets, cg* and t99-mrf, are statistically significant; we note that in both cases, mfe achieved higher prediction accuracy than mea and their confidence intervals don't overlap .
p-value of two-sided permutation test is computed for the mea and mfe algorithms on different parameter sets; p-values below the standard significance level of  <dig>  are shown in bold face and indicate significant differences in prediction accuracy.
furthermore, the results in table  <dig> show that of the parameter sets we considered, bl* gives the highest prediction accuracy, regardless of whether the mea or mfe prediction algorithm is used.
this is consistent with earlier results by andronescu et al.
<cit>  regarding the performance of ubcmfe using the bl* parameter set.
finally, we observe that for a given class of rnas, depending on the energy model, mfe sometimes performs better than mea and vice versa .
for example, when considering our set of 16s ribosomal rnas and using the bl* parameter set, mfe outperforms mea by  <dig> , while for the full turner  <dig> parameters, mea outperforms mfe by  <dig>  .
overall, we conclude that the relative performance of mea versus mfe depends on the thermodynamic parameter set used, and that both ubcmea and ubcmfe with the bl* parameter set have significantly better accuracy than rsmea with the turner <dig> parameter set.
CONCLUSIONS
improvements both in algorithmic methods and in thermodynamic models lead to important advances in secondary structure prediction.
in this work, we showed that gc-pmfmeas with the bl* parameter set significantly outperforms the other mfe and mea-based methods we studied.
however, the relative performance of mea-based and mfe methods vary depending on the thermodynamic parameter set used.
for example, mea-based methods significantly outperform minimum free energy  methods for the turner <dig> model but the opposite is true for other models and in fact the difference in performance between mea-based and mfe methods is negligible on our best thermodynamic model, bl*.
our findings suggest that, as thermodynamic models continue to evolve, a diverse toolbox of algorithmic methods will continue to be important.
we also showed the importance of using large datasets when assessing accuracy of specific algorithms and thermodynamic models.
specifically, we showed that bootstrap percentile confidence intervals for average prediction accuracy on our two largest datasets, ma and s-full, have narrow widths , indicating that the average accuracies measured on these sets are likely to be good estimates of the accuracies of the populations of rna molecules they represent.
in contrast, interval widths for several of the specific rna classes studied in the paper were much larger, with no clear correlation between confidence interval width and either class size or average normalized similarity.
it may be the case that confidence interval widths are correlated with the distribution of evolutionary distances among the class members, rather than on the more simplistic average normalized similarity that we consider in this paper; studying this further would be an interesting direction for future research.
regardless, our analysis shows that larger datasets are needed to obtain reliable accuracy estimates on specific classes of rnas, underscoring the importance of continued expansion of rna secondary structure data sets.
our work illustrates the use of statistical techniques to gauge the reliability and significance of measured accuracies of rna secondary structure prediction algorithms.
we hope that this work will provide a useful basis for ongoing assessment of the merits of rna secondary structure prediction algorithms.
authors' contributions
mh acquired, analyzed and interpreted the data and drafted the manuscript.
ac supervised the research and participated in analysis of data and revising the manuscript.
hh consulted in approaches used and participated in data interpretation and manuscript revision.
all authors read and approved the final manuscript.
supplementary material
additional file 1
95% bootstrap percentile confidence interval graphs for the f-measure average of the rsmea and rsmfe.95% bootstrap percentile confidence intervals are shown for the f-measure average of the rsmea , and rsmfe  algorithms on the ma and s-full sets and also different rna classes in ma.
click here for file
additional file 2
confidence interval width versus rna class size in the ma set for the rsmea and rsmfe methods.
the figure shows the confidence interval width of rna classes in the ma set for the rsmea and rsmfe methods.
click here for file
additional file 3
accuracy comparison of different prediction algorithms with various parameter sets on the s-full-test set.
the table presents the prediction accuracy of different algorithms with different thermodynamic sets in terms of f-measure.
the parameter set t99-mrf refers to the turner <dig> parameters in multirnafold format.
bl* and cg* are the parameter sets obtained by the bl and cg approaches of andronescu et al.
<cit> , respectively.
also, the turner <dig> parameter set is the parameter set obtained by mathews et al.
<cit> .
"n/a" indicates cases in which a given algorithm is not applicable to a parameter set, since that does not match the energy model underlying the algorithm.
the highest accuracies for mea and mfe are shown in bold.
click here for file
acknowledgements
the authors thank and acknowledge dr. mirela andronescu and dr. david h. mathews for providing their data sets and software used in the presented work.
they also thank the anonymous reviewers for their very helpful comments on an earlier draft of the paper.
this research was funded by a grant from the natural sciences and engineering research council of canada .

 rna molecules are essential to many functions in the cells of all organisms. for example, these molecules are involved in gene translation and also act as catalysts and as regulators of gene expression. because function is determined by molecular structure, there is significant investment in computational methods for predicting rna secondary structure, which in turn is useful for inferring tertiary structure. our interest is in assessing the merits of some recent advances in secondary structure prediction in statistically robust ways. we focus on thermodynamically informed approaches for predicting pseudoknot free secondary structures from the base sequence. a widely used method finds the minimum free energy structure with respect to the nearest neighbour thermodynamic model of mathews, turner and colleagues. some recent advances in secondary structure prediction are the new maximum expected accuracy and maximum pseudo-expected accuracy methods of lu et al. and hamada et al.. these approaches generally maximize expected base pair accuracy as a function of base pair probabilities calculated using a partition function method and have higher average accuracy than the mfe algorithm on the turner and andronescu et al. energy parameters. knudsen et al. and do et al. also presented some rna secondary structure prediction methods based on probabilistic models of structures. but since their probabilistic approaches are not determined using a thermodynamic model, we don't include their methods in our later comparisons. another advance is estimation of new energy parameters from both thermodynamic and structural data using state-of-the-art estimation techniques. andronescu et al. derived two parameter sets by inference from energies that were derived from optical melting experiments as well as from structural data. the two energy parameter sets are called bl* and cg*, named after the boltzmann likelihood and constraint generation methods used to infer them. these parameter sets have yielded significant improvements in prediction accuracy of the mfe method, compared with the turner parameters, with the bl* parameters being slightly better than the cg* parameters. here and throughout, the accuracy of a prediction refers to its f-measure, which is the harmonic mean of sensitivity and positive predictive value. all of this work assesses algorithm accuracy on specific classes of rnas, such as introns or transfer rnas, as well as overall average accuracy on rnas taken over all such classes. this recent work motivates the following questions. are comparisons of the mea-based and mfe approaches on specific rna classes reliable when the size of available datasets is small? do the mea- or pseudo-mea-based approaches produce significantly more accurate predictions than mfe on the latest energy parameter sets? what is the best combination of algorithm and thermodynamic model? to answer these questions, we report on the accuracy of both mea-based and mfe methods with respect to two versions of the turner parameters as well as the recent bl* and cg* parameters of andronescu et al., on datasets for specific rna classes as well as large datasets that combine multiple rna classes. we present three main findings. first, we show that f-measure accuracies on our large datasets are likely to be reliable estimates of accuracy of a population represented by such sets, in the sense that high-confidence interval widths for f-measure obtained using the bootstrap percentile method are within a small, 2% range. average accuracy on smaller classes is less reliable. for example, confidence intervals for both mea and mfe have an 8% range on a class of group i introns that has been used previously in benchmarking algorithms. second, there is a clear "winner" in terms of overall prediction accuracy, namely the pseudo-mea-based method of hamada et al.. however, the relative accuracy of the mfe and mea-based approaches depends on the underlying energy parameters: using a permutation test we find that, at a statistically significant level, the accuracy of mfe-based prediction on our large datasets is better on two of the four energy parameter sets that we consider, while mea-based prediction is better than mfe-based prediction on a third parameter set. finally, both mea-based and mfe methods achieve the highest accuracy when using the fourth parameter set we consider, namely the bl* energy parameters of andronescu et al.. methods in this section we first describe the datasets, thermodynamic models, and algorithms considered in this paper. we then describe the accuracy measures and statistical methods used in our analyses. datasets we use three datasets, as follows • s-full is a comprehensive set of rna sequences and their secondary structures that has been assembled from numerous reliable databases. sequences in this and our other datasets have length at most nucleotides; in some cases these were derived by partitioning larger sequences such as 16s ribosomal rna sequences. the average length of sequences in s-full is 270nt. • mt was used by lu et al. in their study of the mea algorithm and contains rnas from the following classes: 16s ribosomal rna, 23s ribosomal rna, 5s ribosomal rna, group i intron, group ii intron, ribonuclease p rna, signal recognition rna and transfer rna. • ma is a subset of the s-full set, containing exactly those sequence of s-full that are in the rna classes included in mt. we formed the ma dataset in order to compare our algorithms on the same classes as did lu et al., while using sets of rnas from these classes with as large a size as possible. an overview of different rna classes in the mt and ma data sets is shown in table this table presents the number of rna sequences, the mean and standard deviation of their lengths and the average of normalized similarity between rna structures for each rna class in these two data sets. overview of different rna classes in the mt and ma data sets, including the number of rnas, the mean and standard deviation of the length and the average normalized similarity between rna sequences in these two data sets for each rna class. lu et al. reported their results on a restricted format of the mt set in which certain bases of some sequences are in lower case, indicating that the base is unpaired in the reference structure. they have used this structural information in their predictions. however, we have not employed this information, and therefore our accuracy measures on the mt set are different from those of lu et al.. thermodynamic models a thermodynamic model consists of features--small structural motifs such as stacked pairs--plus free energy change parameters, one per feature. the first model we use, the turner model called "turner99", is the most widely used energy model for prediction of rna secondary structures. the model has over features, which are based on turner's nearest neighbor rules and reflect the assumption that the stability of a base pair or loop depends on its sequence and on the adjacent base pair or unpaired bases. the model is additive in that the overall free energy change of a secondary structure for a given sequence is the sum of the free energy changes for features of the structure. the parameters of the model were derived from optical melting experiments, the most commonly used experimental approach to determine the free energy change of rna structures. we also consider variants of the turner model, used by andronescu et al.. the t99-multirnafold model is derived from the turner model but includes only features. parameters for features in the turner model can be obtained by extrapolation from the parameters of the t99-mrf model. using maximum likelihood and constraint optimization methods, andronescu et al. derived new free energy change parameters for these features; the resulting models are called bl* and cg* respectively. we used all three models, namely t99-mrf, bl* and cg*, in this work, in order to assess the dependence of algorithm accuracy on model parameters. we note that in lu et al.'s work, the parameters of a newer version of turner model, called turner were used for one of the structural motifs, coaxial stacking. however, for the rest of structural motifs the parameters of turner model were engaged. so, we also call the parameter set used for the lu et al.'s benchmark turner since most of its parameters are turner ones and also the coaxial stacking motif is not employed in the other model that we study. zakov et al. also obtained parameters that improve rna structure prediction. but we don't consider their parameters in this work, since those are not applicable for the partition function calculation and therefore for the probability calculation required for the mea method. algorithms we analyze four rna secondary structure prediction algorithms. the first predicts secondary structures that have minimum free energy with respect to a given thermodynamic model. the second is the maximum expected accuracy algorithm as introduced by lu et al., which maximizes expected base pair accuracy as a function of base pair probabilities calculated using a partition function method. we implemented the mea algorithm for use with the multirnafold models. as a result, we worked not only with the algorithms of lu et al., which we refer to as rsmfe and rsmea, but also with the mfe algorithm of andronescu et al., referred to as ubcmfe, and a new implementation of mea which we developed, called ubcmea. lu et al.'s benchmark showed that rsmea gives the best prediction accuracy when its γ parameter is equal to accordingly, we also set γ to be in ubcmea. the third algorithm that we analyze is the generalized centroid estimator method of hamada et al.. this is similar to the mea method, but uses a somewhat different objective function, namely a gamma-centroid estimator, to infer a structure from base pair probabilities. this method, which we refer to as gc-g employs a parameter γ to balance sensitivity and positive predictive value; based on the results of hamada et al., we set γ = the fourth method is another algorithm of hamada et al. that generalizes the centroid estimator by using a pseudo-expected accuracy maximization technique that automatically selects γ to balance sensitivity and positive predictive value for a given input. we refer to this algorithm as gc-pmfmeas. the rsmea and rsmfe algorithms always use the parameters of turner model as their free energy parameter set, while ubcmea and ubcmfe use parameter sets in the multirnafold model format, namely the bl*, cg* and t99-mrf sets. the gc-g and gc-pmfmeas algorithms also employ bl* parameters using a turner format. accuracy measures we use three measures for determining the structural prediction accuracy, namely sensitivity, positive predictive value or ppv, and f-measure, which combines the sensitivity and ppv into a single measure. sensitivity is the ratio of correctly predicted base pairs to the total base pairs in the reference structures. ppv is the fraction of correctly predicted base pairs, out of all predicted base pairs. f-measure is the harmonic mean of the sensitivity and ppv. this value is equal to the arithmetic mean when sensitivity and ppv are equal. however, f-measure becomes smaller than the arithmetic mean as one of the numbers approaches. the possible values for these three measures are between and 1; the closer to the better prediction. the f-measure is widely used measure in the literature; it is also the common measure in the studies by hamada et al. and by lu et al., to which refer in our study. mostly to facilitate comparison of their results to ours, we decided to use the f-measure rather than the matthews correlation coefficient. sensitivity=numberofcorrectlypredictedbasepairsnumberofbasepairsinthereferencestructureppv=numberofcorrectlypredictedbasepairsnumberofpredictedbasepairsf - measure=2×sensitivity×ppvsensitivity+ppv throughout this paper, we calculate three types of averages for a given measure "m", namely unweighted averages, weighted averages and s-weighted averages, defined below. the weighted average counts each sequence equally, regardless of which class it belongs to. the unweighted average, on the other hand, counts each class equally and was used by lu et al. a potential problem with the unweighted average is that an rna class with many highly similar sequences can have disproportionate influence on the overall accuracy, relative to its sequence diversity. therefore, we introduce the s-weighted average, which takes into account the similarities between rna sequences in each rna class and gives a weight to each one according to its average normalized similarity, in such a way that classes with highly similar sequences have lower weight. the three averages are defined as follows when there are n rna classes, c c..., cn, with cardinalities l l...,ln, respectively. for the remainder of our study, n = and the classes are those listed in table 1). unweightedaverageofm=1n∑∀c1i∈c1ml1+∑∀c2i∈c2ml2+⋯+∑∀cni∈cnmln, weightedaverageofm=∑∀c1i∈c1m+ ∑∀c2i∈c2m+...+ ∑∀cni∈cnml1+l2+...+ln, and s - weightedaverageofm=l1-s1⋅∑∀c1i∈c1m+l2-s2⋅∑∀c2i∈c2m+...+ln-sn⋅∑∀cni∈cnml11-s1+l21-s2+...+ln1-sn, where si is the mean of the normalized similarities measured between the reference structures of any two of rna sequences in the corresponding rna class ci. our normalized similarities were computed using the simtree procedure by eden et al., which takes into account secondary structure similarities in addition to sequence similarities. the average normalized similarities are reported in table and always fall between zero and one, where an average normalized similarity of one means that all sequences in the set are identical. the s-weighted average equals the weighted average when si is zero for all rna classes, and it equals the unweighted average when si is one for all classes. for the remaining cases, the value of the s-weighted average would fall between that of the unweighted and weighted averages. bootstrap percentile confidence intervals we calculated the confidence intervals of average f-measure for various rna classes, by using the bootstrap percentile method following the recent work by aghaeepour and hoos. we first took resamples with replacement from the original f-measure values obtained by the prediction method under investigation on the given reference dataset, where all resamples had the same size as the original sample. we then computed the average f-measures of the resamples and employed these into the bootstrap distribution. finally, we determined the lower and upper limits of the 95% bootstrap percentile confidence interval as the th and th percentile of the bootstrap distribution, respectively. these calculations were all performed using the "boot" package of the r statistics software environment. we verified that the bootstrap distributions are close to gaussian, using the anderson-darling test, which indicates that the bootstrap percentile intervals can be expected to be reasonably accurate. permutation test method to assess the statistical significance of the observed performance differences, we used permutation tests, following aghaeepour and hoos. since lu et al. and hamada et al. reported that the mea-based methods outperform mfe, we used a one-sided permutation test to determine whether mea-based methods have significantly better accuracy than mfe on our parameter sets. the test that we applied proceeds as follows. first, we calculated the difference in means between sets of f-measure values obtained by the two given structure prediction procedures, an mea-based method and mfe. for simplicity, we call these two sets a and b and we denote their sizes as na and nb, respectively. then we combined the f-measure values of sets a and b. next, we calculated and recorded the difference in sample means for randomly chosen ways of dividing these combined values into two sets of size na and nb. the p-value was then calculated as the proportion of the sample means thus determined whose difference was less than or equal to that of the means of sets a and b. then, if the p-value is less than the 5% significance level, we reject the null hypothesis that mfe and the mea-based method have equal accuracy and thus accept the alternative hypothesis that the mea-based method has significantly better accuracy than mfe. otherwise, we cannot reject the null hypothesis and therefore we cannot accept the alternative hypothesis. furthermore, to assess whether the difference in accuracy between the mea-based methods and the mfe method on a given parameter set is significant, we performed a two-sided permutation test. this test works exactly like the one-sided permutation test, except that its p-value is calculated as the proportion of the sampled permutations where the absolute difference was greater than or equal to that of absolute difference of the means of sets a and b. then, if the p-value of this test is less than the 5% significance level, we reject the null hypothesis that mfe and mea have equal accuracy, otherwise, we cannot reject the null hypothesis. all of these calculations were performed using the "perm" package of the r statistics software environment.  in this section, we investigate to which degree the prediction accuracy of the energy-based methods studied in this paper is dependent on different datasets and different thermodynamic parameter sets used for prediction. we start by considering how the size of a dataset will influence the accuracy achieved by rna secondary structure prediction methods. we then study the dependency of prediction accuracy of the mfe, mea, gc-g and gc-pmfmeas algorithms on the thermodynamic parameter sets that they use. dependency of the energy-based methods on data characteristics results of the energy-based methods on different data sets because measures of accuracy on reference datasets are used to assess the quality of prediction achieved by various algorithms, it is important to understand how such measures vary depending on the reference datasets used. later in this section we also consider how accuracy measures vary depending on the energy model used. the baseline data in tables and show that there can be significant differences in accuracy of a given algorithm on rna classes within the ma versus mt sets. for example, on ribonuclease p rna, when using the bl* parameter set, ubcmea achieves an f-measure of on the mt dataset and on the ma set, an absolute difference of about 17%. using data from tables and figure illustrates the difference between accuracy measures on the mt versus ma datasets on both the ribonuclease p rna and the group i intron classes, for the ubcmea, ubcmfe and gc-pmfmeas algorithms, using the bl* parameters, as well as for rsmea and rsmfe using the turner parameters. in the figure, the further the points are from the dotted diagonal line, the larger the difference. of the ten data points shown, only one represents a difference smaller than %. these data clearly indicate that, based only on prediction accuracy and without further statistical analysis, one cannot draw meaningful conclusions about the average accuracy of a particular algorithm on the overall population of an rna class, based only on the average accuracy on currently available datasets. bootstrap confidence intervals for the prediction accuracy of the energy-based methods on different rna classes following the work of aghaeepour and hoos, we use bootstrap percentile confidence intervals to assess the accuracy of the average f-measures obtained by a secondary structure prediction procedure on a given rna dataset; these provide a measure of the overall average accuracy on the whole population from which the dataset is drawn. we chose to calculate confidence intervals using weighted average f-measure. we note from the data of table that on the ma dataset, all three averages are qualitatively similar in the sense that if the f-measure accuracy of one algorithm is better than another with respect to one average, then the same is true with respect to the other averages. thus, we would expect the same qualitative conclusions if we had used a different average. confidence intervals are obtained by the bootstrap percentile method for the ma and s-full sets and different rna classes in ma when ubcmea and ubcmfe use the bl* parameter set and rsmea and rsmfe use the turner parameter set and the confidence level is set to. first, the confidence intervals of all algorithms on the ma and s-full sets have a width of at most, indicating that the average accuracy measured on these sets is likely to be a good estimate - within 1% - of the accuracy of a population of rna molecules represented by these sets. however, the interval widths on individual classes can much higher, e.g., for ubcmea on the group i intron class, suggesting that average accuracy is not a reliable measure of a method's overall accuracy on such classes. we note that, without the use of rigorous statistical methods, a difference in accuracy is considered significant in some related work. second, the confidence interval widths of rna classes do not strictly decrease as the size of the class increases. for example, as shown in table for the ubcmea algorithm the confidence interval width of the ribonuclease p rna class is less than the interval width for the transfer rna class, even though the transfer rna class contains roughly times as many rnas as are in the ribonuclease p rna class. thus factors specific to the classes, other than class size, must account for confidence interval widths. one possibility is that classes with greater similarity among structures would have smaller confidence intervals, since in the limit, if all sequences in a class are identical then the confidence interval has zero width. however, data from tables and did not support such a correlation between average normalized similarity and confidence interval width, even for classes not too different in size. for example, the ribonuclease p and transfer rna classes in the ma set have relatively similar sizes, with the ribonuclease p rna class having lower normalized similarity than the transfer rna class, as one might expect; yet, for ubcmea, the confidence interval for the ribonuclease p rna class is narrower than those for the ribonuclease p rna class. understanding class-specific factors underlying differences in prediction accuracy remains a goal for further study. third, for any given dataset, the confidence interval width for the ubcmfe algorithm is close to that for ubcmea and gc-pmfmeas, even though the location of the interval can be quite different. for example, the average f-measure of ubcmea and ubcmfe on the ribonuclease p rna class differ by, but the confidence interval widths are identical. as another example, the smallest and largest confidence interval widths of the 16s ribosomal rna class among all four algorithms are, respectively, and, which indicates a difference of less than in the interval widths for this class. figure shows the width of confidence intervals versus the size of different rna classes for the ubcmea, ubcmfe and gc-pmfmeas methods. as one can see, almost all of the points for these three methods are coincident, indicating that ubcmea, ubcmfe, and gc-pmfmeas have identical width in the corresponding sizes. the other important observation from this figure is that, as one could expect, the width of the confidence interval tends to decrease with increasing number of rnas of a given type in the reference set, but there are notable exceptions to this trend for sets of size below the final interesting observation from figure and table is that the gc-pmfmeas method using the bl* parameter set outperforms all the other methods on most rna classes and also on our large ma and s-full sets. hamada et al. also showed that on a small set of rnas, gc-pmfmeas achieves better prediction accuracy than mfe and the other mea-based methods when using the bl* parameter set. our results confirm their finding on our large ma and s-full data sets. similar graphs for the rsmea and rsmfe algorithms are provided in the additional files. additional files and figures s and s support findings similar to those reported for ubcmea, ubcmfe and gc-pmfmeas. since the overall results on the gc-g algorithm were roughly the same as the results of the ubcmea method, we report its prediction accuracies only in tables and and did not study it further. dependency of the energy-based methods on thermodynamic parameters the accuracy of algorithmic methods achieved by energy-based secondary structure prediction approaches, such as gc-pmfmeas, mea and mfe, depends on the underlying thermodynamic parameter sets. lu et al. showed that the accuracy of mea is better than that of mfe when using the turner parameter set. but does mea or gc-pmfmeas outperform mfe on other published parameter sets? if so, is the difference in accuracy statistically significant? we address these questions in the following. the table presents the prediction accuracy of different algorithms with different thermodynamic sets in terms of f-measure. the 95% percentile confidence intervals of their accuracies are also shown in parentheses. the parameter set t99-mrf refers to the turner parameters in multirnafold format. bl* and cg* are the parameter sets obtained by the bl and cg approaches of andronescu et al., respectively. also, the turner parameter set is the parameter set obtained by mathews et al.. "n/a" indicates cases in which a given algorithm is not applicable to a parameter set, since that does not match the energy model underlying the algorithm. the highest accuracies for mea and mfe are shown in bold. in contrast, the relative accuracy of mea versus mfe varies, with mfe being slightly better than mea in two cases, mea being slightly better than mfe in one case, and both having approximately the same accuracy in one case, namely bl*. to assess whether mea has significantly better accuracy than mfe on our parameter sets, we first performed the one-sided permutation test with the alternative hypothesis that mea yields higher prediction accuracy than mfe. as recorded in table the resulting p-value falls below the standard significance threshold of only for turner indicating that only for this parameter set, mea outperforms mfe, and its performance, only on this parameter set, can be considered better than mfe. in contrast, the p-values obtained for the other parameter sets, namely cg*, t99-mrf and bl* are above, indicating that we cannot conclude that mea outperforms mfe for these parameters sets. p-value of one-sided permutation test is computed for the mea versus mfe algorithm on different parameter sets; p-values below the standard significance level of are shown in bold face and indicate significantly higher prediction accuracy of mea compared to mfe. we next performed the two-sided permutation test on all mentioned parameter sets. the results, as recorded in table indicate that for the bl* parameter set, the performance difference between mfe and mea is not statistically significant, while the performance differences for the other two parameter sets, cg* and t99-mrf, are statistically significant; we note that in both cases, mfe achieved higher prediction accuracy than mea and their confidence intervals don't overlap. p-value of two-sided permutation test is computed for the mea and mfe algorithms on different parameter sets; p-values below the standard significance level of are shown in bold face and indicate significant differences in prediction accuracy. furthermore, the results in table show that of the parameter sets we considered, bl* gives the highest prediction accuracy, regardless of whether the mea or mfe prediction algorithm is used. this is consistent with earlier results by andronescu et al. regarding the performance of ubcmfe using the bl* parameter set. finally, we observe that for a given class of rnas, depending on the energy model, mfe sometimes performs better than mea and vice versa. for example, when considering our set of 16s ribosomal rnas and using the bl* parameter set, mfe outperforms mea by, while for the full turner parameters, mea outperforms mfe by. overall, we conclude that the relative performance of mea versus mfe depends on the thermodynamic parameter set used, and that both ubcmea and ubcmfe with the bl* parameter set have significantly better accuracy than rsmea with the turner parameter set.  improvements both in algorithmic methods and in thermodynamic models lead to important advances in secondary structure prediction. in this work, we showed that gc-pmfmeas with the bl* parameter set significantly outperforms the other mfe and mea-based methods we studied. however, the relative performance of mea-based and mfe methods vary depending on the thermodynamic parameter set used. for example, mea-based methods significantly outperform minimum free energy methods for the turner model but the opposite is true for other models and in fact the difference in performance between mea-based and mfe methods is negligible on our best thermodynamic model, bl*. our findings suggest that, as thermodynamic models continue to evolve, a diverse toolbox of algorithmic methods will continue to be important. we also showed the importance of using large datasets when assessing accuracy of specific algorithms and thermodynamic models. specifically, we showed that bootstrap percentile confidence intervals for average prediction accuracy on our two largest datasets, ma and s-full, have narrow widths, indicating that the average accuracies measured on these sets are likely to be good estimates of the accuracies of the populations of rna molecules they represent. in contrast, interval widths for several of the specific rna classes studied in the paper were much larger, with no clear correlation between confidence interval width and either class size or average normalized similarity. it may be the case that confidence interval widths are correlated with the distribution of evolutionary distances among the class members, rather than on the more simplistic average normalized similarity that we consider in this paper; studying this further would be an interesting direction for future research. regardless, our analysis shows that larger datasets are needed to obtain reliable accuracy estimates on specific classes of rnas, underscoring the importance of continued expansion of rna secondary structure data sets. our work illustrates the use of statistical techniques to gauge the reliability and significance of measured accuracies of rna secondary structure prediction algorithms. we hope that this work will provide a useful basis for ongoing assessment of the merits of rna secondary structure prediction algorithms. authors' contributions mh acquired, analyzed and interpreted the data and drafted the manuscript. ac supervised the research and participated in analysis of data and revising the manuscript. hh consulted in approaches used and participated in data interpretation and manuscript revision. all authors read and approved the final manuscript. supplementary material additional file 1 95% bootstrap percentile confidence interval graphs for the f-measure average of the rsmea and rsmfe.95% bootstrap percentile confidence intervals are shown for the f-measure average of the rsmea, and rsmfe algorithms on the ma and s-full sets and also different rna classes in ma. click here for file additional file 2 confidence interval width versus rna class size in the ma set for the rsmea and rsmfe methods. the figure shows the confidence interval width of rna classes in the ma set for the rsmea and rsmfe methods. click here for file additional file 3 accuracy comparison of different prediction algorithms with various parameter sets on the s-full-test set. the table presents the prediction accuracy of different algorithms with different thermodynamic sets in terms of f-measure. the parameter set t99-mrf refers to the turner parameters in multirnafold format. bl* and cg* are the parameter sets obtained by the bl and cg approaches of andronescu et al., respectively. also, the turner parameter set is the parameter set obtained by mathews et al.. "n/a" indicates cases in which a given algorithm is not applicable to a parameter set, since that does not match the energy model underlying the algorithm. the highest accuracies for mea and mfe are shown in bold. click here for file acknowledgements the authors thank and acknowledge dr. mirela andronescu and dr. david h. mathews for providing their data sets and software used in the presented work. they also thank the anonymous reviewers for their very helpful comments on an earlier draft of the paper. this research was funded by a grant from the natural sciences and engineering research council of canada.
than when using the parameters of mathews and turner, and there is no significant difference between the accuracy of mea-based methods and mfe when using the bl* parameters. second, on our large datasets, the algorithm with best overall accuracy is a pseudo mea-based algorithm of hamada et al. however, average accuracy on smaller classes of rnas such as a class of group i introns used previously in benchmarking algorithm accuracy is not reliable enough to draw meaningful conclusions about the relative merits of the mfe and mea-based algorithms. that uses a generalized centroid estimator of base pairs. large datasets should be used to obtain reliable measures of the accuracy of rna structure prediction algorithms, and average accuracies on specific classes should be interpreted with caution, considering the relatively small size of currently available datasets for such classes. first, using the bootstrap percentile method, we show that the average f-measure accuracy of the mfe and mea-based algorithms, as measured on our largest datasets with over rnas from diverse families, is a reliable estimate of the accuracy of a population of rna molecules represented by this set. however, between mfe and other mea-based methods, there is no clear winner in the sense that the relative accuracy of the mfe versus mea-based algorithms changes depending on the underlying energy parameters. our work advances such understanding for the mfe and mea-based methods, with respect to the latest datasets and energy parameters. since such benchmarks document progress in improving accuracy of computational prediction methods, it is important to understand how measures of accuracy vary as a function of the reference datasets and whether advances in algorithms or thermodynamic parameters yield statistically significant improvements. advances in prediction methods are typically benchmarked using sensitivity, positive predictive value and their harmonic mean, namely f-measure, on datasets of known reference structures. rna molecules play critical roles in the cells of organisms, including roles in gene regulation, catalysis, and synthesis of proteins. since rna function depends in large part on its folded structures, much effort has been invested in developing accurate methods for prediction of rna secondary structure from the base sequence. the accuracy of the mea-based methods is significantly higher when using the bl* parameter set of andronescu et al. with the bl* parameter set significantly outperforms all other mfe and mea-based algorithms on our large data sets. minimum free energy predictions are widely used, based on nearest neighbor thermodynamic parameters of mathews, turner et al. third, of the four parameter sets we considered, the best accuracy for the mfe-, mea-based, and pseudo-mea-based methods is,, and, respectively and is obtained with a thermodynamic parameter set obtained by andronescu et al.
 rna molecules play critical roles in the cells of organisms, including roles in gene regulation, catalysis, and synthesis of proteins. since rna function depends in large part on its folded structures, much effort has been invested in developing accurate methods for prediction of rna secondary structure from the base sequence. minimum free energy predictions are widely used, based on nearest neighbor thermodynamic parameters of mathews, turner et al. or those of andronescu et al. some recently proposed alternatives that leverage partition function calculations find the structure with maximum expected accuracy or pseudo-expected accuracy methods. advances in prediction methods are typically benchmarked using sensitivity, positive predictive value and their harmonic mean, namely f-measure, on datasets of known reference structures. since such benchmarks document progress in improving accuracy of computational prediction methods, it is important to understand how measures of accuracy vary as a function of the reference datasets and whether advances in algorithms or thermodynamic parameters yield statistically significant improvements. our work advances such understanding for the mfe and mea-based methods, with respect to the latest datasets and energy parameters.  we present three main findings. first, using the bootstrap percentile method, we show that the average f-measure accuracy of the mfe and mea-based algorithms, as measured on our largest datasets with over rnas from diverse families, is a reliable estimate of the accuracy of a population of rna molecules represented by this set. however, average accuracy on smaller classes of rnas such as a class of group i introns used previously in benchmarking algorithm accuracy is not reliable enough to draw meaningful conclusions about the relative merits of the mfe and mea-based algorithms. second, on our large datasets, the algorithm with best overall accuracy is a pseudo mea-based algorithm of hamada et al. that uses a generalized centroid estimator of base pairs. however, between mfe and other mea-based methods, there is no clear winner in the sense that the relative accuracy of the mfe versus mea-based algorithms changes depending on the underlying energy parameters. third, of the four parameter sets we considered, the best accuracy for the mfe-, mea-based, and pseudo-mea-based methods is,, and, respectively and is obtained with a thermodynamic parameter set obtained by andronescu et al. called bl*.  large datasets should be used to obtain reliable measures of the accuracy of rna structure prediction algorithms, and average accuracies on specific classes should be interpreted with caution, considering the relatively small size of currently available datasets for such classes. the accuracy of the mea-based methods is significantly higher when using the bl* parameter set of andronescu et al. than when using the parameters of mathews and turner, and there is no significant difference between the accuracy of mea-based methods and mfe when using the bl* parameters. the pseudo-mea-based method of hamada et al. with the bl* parameter set significantly outperforms all other mfe and mea-based algorithms on our large data sets. 