BACKGROUND
rna secondary structure prediction is still an important problem in computational biology.
with the advent of next generation sequencing and rna-seq technologies, many rna structural changes are being found to play important roles in regulating gene expression  <cit> .
gene regulation studies can now be done on a genome-wide scale.
in some cases rna secondary structures can be experimentally determined on a genome-wide level  <cit> , but these methods require rna isolation and many not preserve in vivo structures.
rna secondary structure prediction programs are still often used to predict structures across the genome  <cit> .
the predicted secondary structures, and predicted structural changes, are being used to find relationships and suggest mechanisms in gene regulatory networks.
some methods for rna secondary structure prediction only consider a single sequence as input.
however, prediction quality can be improved by using multiple sequences, assuming that rna secondary structure is conserved through evolution.
even without a complex evolutionary model, these additional structural constraints provide valuable information on folding.
comparative methods for rna secondary structure prediction are based on this observation, and use evolutionary information from multiple alignments to improve prediction quality.
methods for rna secondary structure prediction generally fall into two categories.
thermodynamic models make use of free-energy functions, which take experimentally determined energy parameters for individual structural elements.
dynamic programming is then used to find the secondary structure with the minimum free energy, which is reported as the predicted structure.
this has been successfully implemented in programs such as rnafold  <cit>  and unafold  <cit> .
thermodynamic methods typically deal with the single-sequence prediction problem, but extensions such as rnaalifold  <cit>  and petfold  <cit>  allow for comparative prediction.
stochastic context-free grammars , on the other hand, define a probability distribution over the space of rna secondary structures.
posterior decoding techniques are typically used to determine, for example, the maximum expected accuracy structure  <cit> .
scfgs have been used for the single-sequence prediction problem  <cit> , but their advantage comes through coupling with a molecular evolution model.
the first comparative scfg-based approach was developed in pfold  <cit> , where alignment column probabilities were determined through single and paired column evolution models, calculated via the felsenstein pruning algorithm  <cit> .
for a more complete review on rna secondary structure prediction, see  <cit> .
in genome-wide predictions of rna secondary structure, the accuracy of the secondary structure pre-diction program is rarely factored into analysis.
typically only the mean accuracy is reported in rna secondary structure prediction benchmarks  <cit> , with the variance in the accuracy given little thought.
variance in accuracy is particularly problematic in the case of single-sequence prediction.
figure  <dig> shows the cumulative density function of predictive accuracy for two single-sequence applications of rna secondary structure prediction, rnafold and ppfold , on  <dig> sequences taken from the rnastrand database  <cit> .
additionally, a uniform  cumulative density function is shown for comparison.
the figure illustrates that for sequences in this data set, the predictive accuracy of rnafold and ppfold is not very much different from a random number between  <dig> and  <dig>  when genome-wide rna secondary structure prediction is done on a large number of single sequences, many predictions will be poor ones.
consequently, it is important to understand more about the variability of rna secondary structure prediction programs.
in comparative prediction, there are many sources of variability: alignment quality, the number of sequences chosen and which alignment samples have been taken, the evolutionary relationships between sequences, as well as the ill-conditioned nature of the folding model itself.
understanding and quantifying these variances is key for biological applications that rely on these folding programs.
additionally, other bioinformatics software that utilize these folding programs– for example inverse rna folding algorithms  <cit> – may often experience a fundamental limitation in performance due to variance in structure prediction quality.
sequence alignment is a fundamental step in most comparative sequence analysis pipelines.
the typical approach is to create a single, trusted multiple alignment of the sequences using methods based on an artificial scoring scheme and heuristics to find a highly scoring alignment  <cit> .
although this methodology is successful when the alignment is well resolved, it has been shown in the context of downstream analyses that the end result can be highly sensitive to the choice of alignment  <cit> .
rna secondary structure prediction methods take a variety of approaches with respect to possible errors in rna alignments.
some methods  invoke a fold-and-align approach directly, where alignment is done simultaneously with structure prediction.
pfold, instead, takes a fixed alignment as input, but allocates a nite probability to a nucleotide being any other nucleotide; this makes the model more robust to poor alignment quality.
most modern methods  still consider prediction from a single, fixed, alignment.
recently, alignment-free methods have also been proposed  <cit> .
however, even after considering poor alignment quality, there are many additional variances associated with poor comparative rna secondary structure prediction.
sequence selection is another important variable.
alignment methods may produce poor alignments due to poor individual sequences, which will in turn produce poor structure predictions.
there have been methods developed to select homologous sequences, particularly  <cit> , which is based on evolutionary models and structural constraints.
this is implemented in  <cit> , which shows strong results in rna secondary structure prediction can be found by selecting useful sequences.
in this paper we consider the problem of variances in comparative rna secondary structure prediction.
we present statistical analyses of different variances including the relationship between structure prediction quality and chosen alignment distance from the reference alignment, and a predictive algorithm for accuracy is provided.
factors like the number of sequences in the alignment and the evolutionary distance of the sequences are considered.
finally, a novel method is presented which extends information entropy for stochastic context-free grammars  <cit>  to consider variation over alignments.
statistical alignment
statistical alignment  <cit>  provides a solution to many of the issues encountered with the traditional approach of sequence alignment.
it models sequence evolution as a stochastic process involving sequence insertions, deletions and character substitutions, which defines a probability distribution over the alignments of the sequences.
using techniques such as expectation maximisation or markov chain monte carlo , it is possible to either maximize the likelihood of the alignment, or generate a representative set of putative alignments by sampling from the alignment distribution.
several statistical alignment implementations have emerged in past years  <cit> , some of which allow co-sampling other entities such as the evolutionary tree or the locations of cis-regulatory motifs.
such methods can highlight homoplasy and alignment uncertainty with high accuracy or can be used to decrease alignment uncertainty effects in downstream analyses, for instance protein secondary structure prediction  <cit> .
statalign is a statistical alignment software package  <cit>  that allows joint bayesian analysis of multiple alignments, phylogenetic trees and evolutionary parameters.
the background model for insertions and deletions is a modified version of the tkf <dig> model  <cit>  as described in  <cit> .
the indel model can be coupled with an arbitrary substitution model .
the bayesian analysis is based on mcmc, the transition kernels are improved versions of those in  <cit> .
statalign generates random samples from the joint posterior distribution of sequence alignments, evolutionary trees and model parameters.
this high-dimensional joint distribution can be analysed in several ways, ranging from the simple statistics of marginalised single dimensions  to the application of other tools to the alignment samples.
alignment and rna secondary structure accuracy metrics
to analyse variances of rna secondary structure as alignment quality varies, we calculate a similarity score that measures how close a sample alignment is to the reference alignment.
we use an alignment metric, taken from  <cit> , which is generalised to an alignment method in  <cit> .
let as <dig> s <dig> be an alignment of a sequence s <dig> of length n to a sequence s <dig> of length m. each column of as <dig> s <dig> can be expressed as pairs of the form s1i,s2j, s1i,−, and −,s2i.
we define
has <dig> s2=i,j|s1i,s2j∈as <dig> s <dig> ias <dig> s2=j|−,s2j∈as <dig> s <dig> anddas <dig> s2=i|s1i,−∈as <dig> s <dig>
sets which represent ‘homology’, ‘insertion’ and ‘deletion’ respectively.
given these sets, we define the distance between two alignments as <dig> s2k and as <dig> s2ℓ of two sequences s <dig> and s <dig> to be
das <dig> s2k,as <dig> s2ℓ=n+m−2has <dig> s2k∩has <dig> s2ℓ−das <dig> s2k∩das <dig> s2ℓ−ias <dig> s2k∩ias <dig> s2ℓ.
for example, consider the case as <dig> s2k=as <dig> s2ℓ.
then we have n=m,has <dig> s2k∩has <dig> s2ℓ=n,das <dig> s2k∩das <dig> s2ℓ= <dig> as there are no deletions, and ias <dig> s2k∩ias <dig> s2ℓ= <dig> as there are no insertions.
this gives the distance between the alignments as zero, as would be expected.
equation  <dig> can be generalised to sequence alignments with more than two sequences.
assuming now that ak and aℓ are alignments of m sequences s <dig>  …, sm of lengths n <dig>  …, nm, we have
dak,aℓ=∑p=1m−1∑q=p+1mdasp,sqk,asp,sqℓ,
that is, summing all the pairwise alignment distances from equation  <dig>  this alignment metric if is then normalized and subtracted from  <dig> to produce a similarity score
ssak,aℓ=1−dak,aℓm−1∑t=1mnt
the denominator of the fraction, m−1∑t=1mnt, is the normalizing constant, the maximum that the alignment distance dak,aℓ can be.
the similarity score is bounded by  <dig> and  <dig>  with  <dig> indicating that the sample alignment is identical to the reference alignment.
rna secondary structure metrics
there are a wealth of available metrics on rna secondary structure  <cit> .
here we use sensitivity, positive predictive value , and f-score .
defining true positives  as the number of base-pairs correctly predicted, false positives  as the number of true base-pairs not predicted, and false negatives  as the number of base-pairs predicted which are incorrect, we have
sensitivity=tptp+fnppv=tptp+fpf−score=2×tp2×tp+fn+fp
the strength of these rna secondary structure accuracy metrics is that they are easy to interpret, and make it straightforward to compare methods across different datasets.
an f-score of  <dig> would represent a structure prediction that was completely correct and an f-score of  <dig> a structure prediction that only predicted incorrect base-pairs.
information entropy
as we later develop calculations for information entropy for a set of alignments, here we outline the computation of information entropy for a single alignment.
the information entropy h of a probability distribution p with a set of events x is defined as:
hp=−∑x∈xpxlog2px.
information entropy is a measure for the “spread” of the probability distribution, and has well-defined lower and upper bounds.
the minimum entropy of  <dig> occurs when there is only one outcome with probability  <dig>  and the maximum entropy of log <dig>  occurs when there are n possible outcomes, each with probability 1/n, that is the uniform distribution.
when the base of the logarithm is  <dig>  the entropy is measured in bits.
for a probability distribution, an entropy of k bits indicates that the expected value of the information content of observing a single outcome is k bits.
in the context of secondary structure prediction, a low entropy therefore indicates that few secondary structures dominate the probability space, whereas a high entropy indicates a more even probability distribution over possible secondary structures.
thus, information entropy is a useful single quantity to characterize the underlying probability distribution of secondary structures.
the information entropy of the probability distribution over the possible secondary structures generated by a phylo-scfg can be obtained using expected rule frequencies, which can be computed using the inside-outside algorithm  <cit> .
this is outlined here.
let the set of all derivations for the input alignment be Φ. since the probability of a derivation d can be written as the product of the scfg production rule probabilities and the phylogenetic probabilities, we can write the total probability t of the grammar producing the input string as
t=∑d∈Φpd=∑d∈Φpgdptd,
where pg denotes the prior probabilities obtained from the scfg part of the model, and pt  are the likelihood factors obtained from the phylogenetic model.
conditioning on producing the input string, the normalized probability of a derivation d is pΦd=1tpd=1tpgdptd.
consequently, we have that the information entropy of the input alignment under the phylo-scfg model is
hΦg=−∑d∈ΦpΦdlog2pΦd,
which can be written using equation  <dig> as
hΦg=log2t−1t∑d∈Φpdlog2pgd−1t∑d∈Φpdlog2ptd,
that is, separating out the scfg contribution and the phylogenetic contribution.
to calculate the entropy in practice, firstly we use a simplified form of the scfg contribution.
for a scfg with set of production rules r, we can write the scfg contribution in terms of the expected production rule frequency,
∑d∈Φpdlog2pgd=∑r∈rlog2pgreusesofr.
secondly, we can simplify the phylogenetic contribution.
let ra ∈ r be a scfg rule which produces base pairs, and rb ∈ r a scfg rule which produces unpaired bases.
define 1d , the indicator function for whether the column pair  is emitted from a rule ra , and 1s , the indicator function for whether column i is emitted from a rule rb .
finally, define fd  as the frequency that rule ra is used in derivation d. then
∑d∈Φpdlog2ptd=∑d∈Φpd∑ralog2ptc|cunpairedfdra+∑rblog2ptc,c′|c,c′fdrb=∑ilog2pti|iunpaired∑d∈Φ1sipd+∑i,jlog2pti,j|i,jpaired∑d∈Φ1di,jpd
as ∑d∈Φ1di,jpd is the total probability under the model that positions i and j are emitted from a rule ra, and ∑d∈Φ1sipd is just the total probability under the model that position i is emitted from a rule rb, the quantity ∑d∈Φpdlog2ptd can be computed using the expected rule frequencies obtained through the inside-outside algorithm  <cit> .
methods
data
statalign dataset
to test factors relating to alignment quality and secondary structure prediction quality, a large number of alignment samples from trusted reference alignments with known secondary structures are needed.
we have created a curated rna dataset based on the rfam database  <cit>  for the purposes of evaluating the framework.
alignments of homologous rna sequences with known consensus secondary structure were extracted from rfam seed alignments.
from these,  <dig> rna families with at least  <dig> sequences were randomly selected  in the section “statalign dataset”.
from each family, in a pre-filtering step we removed divergent sequences with long indels, as follows.
we defined insertion as consecutive non-gap characters of a sequence in the reference alignment which appear in columns where over 80% of the sequences have gaps.
deletions were defined analogously.
columns with fraction of gaps between 20% and 80% were regarded ambiguous and ignored.
to over-penalize long indels, we applied the super linear score function l × log <dig>  for indels of length l, indels being defined as above.
then, a sequence was removed from a family if its total indel score was beyond  <dig> and the difference between its indel score and the mean indel score in the family was beyond  <dig>  times the standard deviation of the indel scores in family, i.e.
if the sequence had significantly more and/or longer indels than what is representative of the rest of the family.
then,  <dig> sequences were selected at random, and further random selection was done to get pairs, triplets etc.
up to  <dig> sequences in alignments.
from these samples of known reference alignment, we could produce many different alignment samples using statalign  <cit> .
for each rna alignment,  <dig> alignment samples were taken, and the reference alignments were also kept to for comparison.
we refer to this dataset throughout as the statalign dataset.
random alignment data
we also wanted to measure the effect of alignment accuracy on secondary structure independently of alignment method.
therefore we created a dataset based on the rna families of the statalign dataset, where alignments were sampled uniformly at various fixed distances from the reference alignment.
using the alignment distance measure in eq.
<dig>  we created a metropolis-coupled mcmc framework that runs several parallel mcmc chains to take alignment samples from the target distribution
πa=exp−2da,ar−dt
where ar is the reference alignment, d is the target distance to get samples from and t is the temperature of the chain.
to improve the mixing properties of the chains we allowed each chain to explore alignments that do not exactly match the specified target distance  but then rejected non-exact matches when taking samples from the cold chain .
the state space of the markov chains is the set of all possible multiple alignments of the input sequences.
alignments that represent the same set of homology statements, and only differ by the order of the alignment columns, are treated as different .
the following basic alignment rearrangement moves are iterated:
<dig>  breaking an alignment column into two columns by moving one of its characters into a new column, placing gaps in every other row:
cccaaaggg⇌ <dig> .ccc−aa−a−ggg
<dig>  the exact reverse of the previous, i.e.
joining two compatible adjacent columns – one having gaps in all but one row, the other having a gap in that row – to form a single column
<dig>  relocating one character of a column to a gap position of an adjoining column:
cccagaga−ttt⇔ <dig> cccagag−attt
as the space of alignments is vast and the moves are very local, to get a good approximation of uniform sampling, the number of steps that have to be done is on the order of millions, even for small inputs when a single chain is run.
to this end, we created a special alignment representation where the above rearrangements can be carried out very efficiently .
for moderate input sizes  this representation allowed us to take  <dig>  million rearrangement steps in a second .
a single chain is sufficient for small alignments, but very slow mixing becomes an issue for practical alignment sizes.
we have found that standard parallel tempering techniques effectively speed up mixing if the chain temperatures are chosen correctly.
because the optimal temperatures vary significantly depending on reference alignment characteristics and target distribution, a simple acceptance optimization routine was added that tunes chain temperatures to achieve chain swap acceptance ratios between neighboring chains around a pre-set value.
we found ratios  <dig> – <dig>  to be the most effective.
with this framework, running 8– <dig> parallel chains was best to maximize the speed of convergence to the uniform distribution as compared to a single chain with sampling times 8–10-fold.
the above described framework was utilized to create a dataset consisting of the rna families of statalign dataset, where for each family and each selection of  <dig> representative sequences from the family,  <dig> samples were taken at a distance corresponding to a similarity of  <dig>  to the reference alignment , then  <dig> samples at a similarity of  <dig>  etc., down to a similarity ratio of  <dig> .
we refer to this dataset as the random alignment dataset.
extending information entropy to alignment space
the information entropy defined for a single alignment contains the length of the alignment as a parameter.
attempting to extend the measure to the probability mass over rna secondary structure space, variable alignment length is a concern.
for example, if we have two alignments and corresponding secondary structures
))))ccccaaaa−ggggccccaaaaggggccccaaa−aggggccccaaaagggg
we would not want to suggest that these alignments give two different secondary structures.
consequently, we use a projection method to give alignment column pairing probability matrices the same dimension, so that the matrices can be averaged.
for a given set of input sequences, the sequence containing the greatest number of non-gap characters was chosen as the reference sequence.
each pairing probability matrix is projected by deleting columns and rows of the matrix corresponding to gap positions in the reference sequence, thus ensuring each matrix corresponding to an alignment sample has the same dimensions .
for example, we might start with an   matrix, delete row i and column i due to a gap in position i in the reference sequence, then end up with an n × n matrix as required.
to calculate information entropy over alignments, we need to be able to calculate the probability of each alignment.
however, we cannot calculate the information entropy explicitly, since the probability of a given secondary structure ss is
pss=∑a∈alignmentspss|apa,
and there is no known efficient way to recurse over all possible alignments.
instead, we create an information entropy measure based on samples from the alignment space, and show that, in the sample-size limit, the alignment-sample information entropy tends to the true information entropy.
consider alignment samples a <dig>  …, an from the space of all alignments of m sequences, sampled according to their probability.
if we are using statistical alignment, as in statalign, we will be sampling alignments in this fashion.
then we have for a column c, once alignments have been projected to the same length, the probability of being unpaired
pc|cunpaired=∑i=1npc|cunpairedinalignmentaipai,
with an analogous result holding for paired columns.
we now define a sample phylogenetic probability ps as the average of the sample phylogenetic probabilities:
psc|cunpaired=∑i=1n1npc|cunpairedinalignmentai,
to show this sample probability converges to the true probability as sample size tends to infinity, we first note that, rearranging equation 12:
psc|cunpaired=∑i=1n1n∑a∈alignmentspc|cunpairedinalignmenta1a=ai,
with  <dig> being the indicator function.
hence
psc|cunpaired=∑a∈alignmentspc|cunpairedinalignmenta∑i=1n1n1a=ai.
taking the limit n → ∞, by the weak law of large numbers
psc|cunpaired→∑a∈alignmentspc|cunpairedinalignmentapaasn→∞
=pc|cunpaired
as required.
now, we have from above that the entropy hΦ  of grammar derivations Φ of a grammar g is
hΦp=log2t−1t∑d∈Φpgdptdlog2pgd−1t∑d∈Φpgdptdlog2ptd
since tree probabilities are the product of unpaired and paired column probabilities in the derivation, the tree probabilities can be recalculated from the sample of alignments.
these can then be substituted into the above equation to get an approximation to the information entropy over the space of sampled alignments as well.
we refer to this entropy of more than one alignment as the alignment consensus entropy.
RESULTS
alignment quality and predictive accuracy
a common question in comparative rna secondary structure prediction is how many sequences are required to get a good structure prediction.
this is briefly addressed in  <cit> , but the sample size is quite small, and only total accuracy is considered.
with  <dig> sequences in the alignment, we assume that no more evolutionary information can be gained by adding further sequences, but when fewer sequences are present, lack of information might yield a poorer structure prediction.
instead of considering total accuracy, we wanted to quantify relatively how much accuracy is lost when fewer sequences are present.
for example, if an alignment with  <dig> sequences is predicted with an average f-score of  <dig> , and an alignment of the same family with  <dig> sequences is predicted with an average f-score of  <dig> , then 80% of the accuracy will have been retained, that is the alignment with  <dig> sequences has a relative f-score of  <dig> .
to investigate how many sequences are needed for a good structure prediction, we took the statalign dataset and considered the relative f-score for each family.
almost 100% of the possible f-score was achieved when the alignment contained  <dig> sequences, for both ppfold and rnaalifold.
interestingly, the accuracy of rnaalifold decreased slightly as the number of sequences was increased.
this is due to the increased number of non-canonical base-pairs in the alignments, which the thermodynamic method could not predict.
ppfold, on the other hand, has a small probability for non-canonical base-pairs, so is not affected by these in the same way.
overall these results suggest that  <dig> sequences are sufficient for approaching maximal predictive accuracy.
to consider the effect alignment quality has on rna secondary structure prediction, we took the statalign and random alignment datasets and measured their similarity to the reference alignment using the similarity score above.
again, percentage of accuracy retained was calculated by normalizing against the accuracy achieved on the reference alignment.
log-scale heatmaps showing the accuracy and percentage of accuracy retained for the statalign dataset  and random alignment dataset  for ppfold and rnaalifold can be seen in figure  <dig>  as expected, decreasing alignment quality decreases the accuracy of structure predictions.
however, other patterns also emerge from these graphs.
firstly, in the statalign dataset , we observe a weaker correlation between alignment distance and accuracy than for the random alignment dataset .
this suggests that predictions are better for the statalign dataset on alignments far from the reference alignment.
statalign looks to produce alignments with a high likelihood under an evolutionary model, which the random alignments do not consider, and so statalign's alignments could be considered more realistic.
this con rms the expectation that statalign's alignments are more useful for rna secondary structure prediction than random alignments.
secondly, for the statalign dataset, we see a much higher correlation with the f-score than with the relative f-score.
some families of rna consistently produce the same alignment, which skews the graphs.
for example, an alignment which consistently has similarity to the reference alignment of  <dig> , and f-score  <dig>  would give a relative f-score of  <dig> every time, and would support the correlations seen on each graph.
because we can control the spread of distances in the random alignment data set, we don't see this behaviour.
as expected, variation comes more with families that produce more variable alignments.
in the statalign dataset, this is obscured by those families which produce consistent alignments.
lastly, for the random alignment dataset , we see many more zero quality predictions in the case of rnaalifold than in the case of ppfold.
this is most easily seen by the larger intensity of red in the main body of the heatmap for ppfold.
this suggests that ppfold functions better than rnaalifold when given a low-quality alignment, likely due to its more complete model for molecular evolution.
evolutionary distance
we also consider the effects of evolutionary distance on rna secondary structure prediction quality.
one might expect that there is a “sweet spot” for evolutionary distance– sequences too close to each other do not display enough co-variation to benefit the evolutionary model, but the evolutionary signal might be lost if the distance is too large.
to investigate this, we measured evolutionary distance in the phylogenetic trees predicted by ppfold using four different measures:
measure 1– mean of all the evolutionary distances,
measure 2– standard deviation of all the evolutionary distances,
measure 3– maximum evolutionary distance,
measure 4– maximum difference between evolutionary distances.
all four measures would be expected to be correlated with sequence length, which is well known to correlate with predictive accuracy.
to account for this, we considered relative evolutionary distance, similar to the relative predictive accuracy above.
a measure was normalized by the average for that family and number of sequences, so that it could be seen whether an alignment had greater or less evolutionary distance than might be expected.
we then looked at the correlation between relative evolutionary distance and predictive accuracy.
all methods correlated extremely poorly with predictive accuracy and relative predictive accuracy, measures  <dig> to  <dig> having r <dig> correlations with relative predictive accuracy of  <dig> ,  <dig> ,  <dig> , and  <dig>  respectively .
this suggests that evolutionary distance is not an underlying factor for variation in rna secondary structure prediction, and those other factors, such as those seen in  <cit> , play a more important role in determining predictive accuracy.
alignment distances and maximum posterior decoding
given the results concerning accuracy lost as alignment quality decreases, it would be desirable to be able to predict alignment quality, with the hope of predicting structure prediction quality.
this has previously been attempted in  <cit> .
first, the sequences were aligned with clustalw  <cit> .
the sequences were then re-aligned using  <dig> other programs  and the similarity between the alignment generated using clustalw to each of the  <dig> other alignments was measured.
the maximum of the  <dig> similarities, max , was chosen as a predictor of alignment quality.
the authors of  <cit>  detected a strong correlation between the true similarity  and max .
we implemented a modified version of this method.
for a given set of input sequences we aligned with both amap  <cit>  and with statalign, obtaining the maximum posterior decoding alignment  from statalign.
the similarity between the amap alignment and the mpd alignment was used as our predicted similarity measure.
this produced a strong correlation between our predicted similarity and true similarity, with an adjusted r-squared value of  <dig> .
we also implemented another method, which calculates an estimate of the expected similarity score using posterior probabilities from the mpd alignment.
for each column, we might expect that a posterior probability close to  <dig> would contribute a score of close to  <dig> to the similarity measure.
so our predicted alignment distance is just the average of the column posterior probabilities.
figure  <dig> shows an example of the correlation between predicted similarity and true similarity, here giving an adjusted r-squared value of  <dig> .
our new predicted similarity can be calculated efficiently, and is a strong predictor of true similarity to the reference alignment.
extending information entropy to alignment space
to test the information entropy extension developed above, we calculated the alignment consensus entropy for samples of alignments from the statalign dataset.
figure  <dig> gives information entropy for  <dig> different representative rna families from the statalign dataset.
on each graph, the information entropy for each of  <dig> statistical alignment samples is given, as well as the alignment consensus entropy.
the leftmost figure is one where the alignment samples were very similar, the rightmost figure where the alignment samples were very different, and the middle figure closer to the median value of alignment sample similarity.
for family  <dig> , we see that the alignment consensus entropy is comparative to the mean entropy of the individual samples.
this is expected, as it indicates there is little uncertainty in the alignment.
on the other hand, the high-diversity family has much higher alignment consensus entropy than for each individual sample.
this is again expected, as the difference in entropies indicates a high uncertainty in alignment.
in this way, we can incorporate alignment uncertainty into our understanding of comparative rna secondary structure prediction.
predicting secondary structure accuracy
given strong correlations between the alignment quality and the structure prediction quality, we might expect that we could find a predictor of structure prediction accuracy.
by “integrating out” alignment uncertainty, we may find a reliability score which is more reliable than the one currently reported by the ppfold.
to test this, we predicted accuracy for one of the five-sequence alignments of each family and then tested the predicted accuracy against the true accuracy.
the ppfold reliability score produced an adjusted r <dig> score of  <dig>  when considering correlation with the true f-scores.
for our new reliability score, we adjusted the ppfold reliability score to consider only base-pairs, as the f-score considers only base-pairs .
we then performed linear regression with the average of the mpd column probabilities, the information entropy of the alignment space, and this pairs-only reliability score against the known f-score measure.
this multiple regression improved the reliability score significantly.
figure  <dig> shows the predicted f-score against the true f-score, for a randomly chosen five-sequence alignments from each family of the statalign dataset.
the adjusted r <dig> value with the new reliability measure improved to  <dig> .
these results seem to indicate that while alignment quality does affect structure prediction quality, the actual structure prediction model still plays a great role in the overall prediction accuracy.
consequently, improving these models, possibly by incorporating other kinds of information , is an area where new research e orts are still needed in rna secondary structure prediction.
CONCLUSIONS
in this paper we have explored a number of factors which can contribute to poor rna secondary structure prediction quality.
we established a relationship between alignment quality and expected loss of accuracy.
furthermore, we provided a method to predict alignment quality based only on statistical alignment samples.
while our predictor of accuracy improves on the ppfold reliability score, there is still a large amount of variability unaccounted for, which we therefore suggest comes from the predictive model itself.
to consider this further, we extended the information entropy measure for scfgs to consider uncertainty in alignments.
the fact that our accuracy predictor did not account for all the variances associated with rna secondary structure prediction, despite good predictors being found for alignment quality and a strong correlation between alignment quality and predictive accuracy, suggests that whilst alignment quality is an important factor, the predictive model itself determines plays a large part in the quality of prediction.
given what is shown in figure  <dig> for single sequence predictions, that the accuracy of ppfold and rnafold is very variable, it is unsurprising that variances remain.
clearly then, further efforts should be put into creating stronger single-sequence models, and then the advantages of evolutionary modelling and additional structural constraints will benefit further.
the use of experimental data from new probing experiments as well as more biologically realistic constraints, such as kinetic or co-transcriptional folding, may improve the results of rna secondary structure prediction.
competing interests
the authors declare they have no competing interests.
authors’ contributions
jwja, an, and zs formulated the ideas presented here, carried out some of the analyses, and wrote the manuscript.
an generated the statalign and random alignment datasets, and mg, pa, and ie generated some of the results for the reliability scores, alignment distances, and predictive algorithm.
jwja and zs developed the theoretical framework for extending information entropy.
all authors read and approved the final manuscript.
supplementary material
additional file 1
information on the  <dig> rna families randomly selected from rfam.
click here for file
acknowledgements
in part, this work was carried out as part of the oxford summer school in computational biology,  <dig>  in conjunction with the department of plant sciences and the department of zoology.
funding was provided by the eu cogangs grant.
Án would like to thank bbsrc for continued funding through ocisb.
jwja would like to thank the epsrc for funding.

 rna secondary structure prediction is still an important problem in computational biology. with the advent of next generation sequencing and rna-seq technologies, many rna structural changes are being found to play important roles in regulating gene expression. gene regulation studies can now be done on a genome-wide scale. in some cases rna secondary structures can be experimentally determined on a genome-wide level, but these methods require rna isolation and many not preserve in vivo structures. rna secondary structure prediction programs are still often used to predict structures across the genome. the predicted secondary structures, and predicted structural changes, are being used to find relationships and suggest mechanisms in gene regulatory networks. some methods for rna secondary structure prediction only consider a single sequence as input. however, prediction quality can be improved by using multiple sequences, assuming that rna secondary structure is conserved through evolution. even without a complex evolutionary model, these additional structural constraints provide valuable information on folding. comparative methods for rna secondary structure prediction are based on this observation, and use evolutionary information from multiple alignments to improve prediction quality. methods for rna secondary structure prediction generally fall into two categories. thermodynamic models make use of free-energy functions, which take experimentally determined energy parameters for individual structural elements. dynamic programming is then used to find the secondary structure with the minimum free energy, which is reported as the predicted structure. this has been successfully implemented in programs such as rnafold and unafold. thermodynamic methods typically deal with the single-sequence prediction problem, but extensions such as rnaalifold and petfold allow for comparative prediction. stochastic context-free grammars, on the other hand, define a probability distribution over the space of rna secondary structures. posterior decoding techniques are typically used to determine, for example, the maximum expected accuracy structure. scfgs have been used for the single-sequence prediction problem, but their advantage comes through coupling with a molecular evolution model. the first comparative scfg-based approach was developed in pfold, where alignment column probabilities were determined through single and paired column evolution models, calculated via the felsenstein pruning algorithm. for a more complete review on rna secondary structure prediction, see. in genome-wide predictions of rna secondary structure, the accuracy of the secondary structure pre-diction program is rarely factored into analysis. typically only the mean accuracy is reported in rna secondary structure prediction benchmarks, with the variance in the accuracy given little thought. variance in accuracy is particularly problematic in the case of single-sequence prediction. figure  shows the cumulative density function of predictive accuracy for two single-sequence applications of rna secondary structure prediction, rnafold and ppfold, on sequences taken from the rnastrand database. additionally, a uniform cumulative density function is shown for comparison. the figure illustrates that for sequences in this data set, the predictive accuracy of rnafold and ppfold is not very much different from a random number between and when genome-wide rna secondary structure prediction is done on a large number of single sequences, many predictions will be poor ones. consequently, it is important to understand more about the variability of rna secondary structure prediction programs. in comparative prediction, there are many sources of variability: alignment quality, the number of sequences chosen and which alignment samples have been taken, the evolutionary relationships between sequences, as well as the ill-conditioned nature of the folding model itself. understanding and quantifying these variances is key for biological applications that rely on these folding programs. additionally, other bioinformatics software that utilize these folding programs– for example inverse rna folding algorithms – may often experience a fundamental limitation in performance due to variance in structure prediction quality. sequence alignment is a fundamental step in most comparative sequence analysis pipelines. the typical approach is to create a single, trusted multiple alignment of the sequences using methods based on an artificial scoring scheme and heuristics to find a highly scoring alignment. although this methodology is successful when the alignment is well resolved, it has been shown in the context of downstream analyses that the end result can be highly sensitive to the choice of alignment. rna secondary structure prediction methods take a variety of approaches with respect to possible errors in rna alignments. some methods invoke a fold-and-align approach directly, where alignment is done simultaneously with structure prediction. pfold, instead, takes a fixed alignment as input, but allocates a nite probability to a nucleotide being any other nucleotide; this makes the model more robust to poor alignment quality. most modern methods still consider prediction from a single, fixed, alignment. recently, alignment-free methods have also been proposed. however, even after considering poor alignment quality, there are many additional variances associated with poor comparative rna secondary structure prediction. sequence selection is another important variable. alignment methods may produce poor alignments due to poor individual sequences, which will in turn produce poor structure predictions. there have been methods developed to select homologous sequences, particularly, which is based on evolutionary models and structural constraints. this is implemented in, which shows strong results in rna secondary structure prediction can be found by selecting useful sequences. in this paper we consider the problem of variances in comparative rna secondary structure prediction. we present statistical analyses of different variances including the relationship between structure prediction quality and chosen alignment distance from the reference alignment, and a predictive algorithm for accuracy is provided. factors like the number of sequences in the alignment and the evolutionary distance of the sequences are considered. finally, a novel method is presented which extends information entropy for stochastic context-free grammars to consider variation over alignments. statistical alignment statistical alignment provides a solution to many of the issues encountered with the traditional approach of sequence alignment. it models sequence evolution as a stochastic process involving sequence insertions, deletions and character substitutions, which defines a probability distribution over the alignments of the sequences. using techniques such as expectation maximisation or markov chain monte carlo, it is possible to either maximize the likelihood of the alignment, or generate a representative set of putative alignments by sampling from the alignment distribution. several statistical alignment implementations have emerged in past years, some of which allow co-sampling other entities such as the evolutionary tree or the locations of cis-regulatory motifs. such methods can highlight homoplasy and alignment uncertainty with high accuracy or can be used to decrease alignment uncertainty effects in downstream analyses, for instance protein secondary structure prediction. statalign is a statistical alignment software package that allows joint bayesian analysis of multiple alignments, phylogenetic trees and evolutionary parameters. the background model for insertions and deletions is a modified version of the tkf model as described in. the indel model can be coupled with an arbitrary substitution model. the bayesian analysis is based on mcmc, the transition kernels are improved versions of those in. statalign generates random samples from the joint posterior distribution of sequence alignments, evolutionary trees and model parameters. this high-dimensional joint distribution can be analysed in several ways, ranging from the simple statistics of marginalised single dimensions to the application of other tools to the alignment samples. alignment and rna secondary structure accuracy metrics to analyse variances of rna secondary structure as alignment quality varies, we calculate a similarity score that measures how close a sample alignment is to the reference alignment. we use an alignment metric, taken from, which is generalised to an alignment method in. let as s be an alignment of a sequence s of length n to a sequence s of length m. each column of as s can be expressed as pairs of the form s1i,s2j, s1i,−, and −,s2i. we define has s2=i,j|s1i,s2j∈as s ias s2=j|−,s2j∈as s anddas s2=i|s1i,−∈as s sets which represent ‘homology’, ‘insertion’ and ‘deletion’ respectively. given these sets, we define the distance between two alignments as s2k and as s2ℓ of two sequences s and s to be das s2k,as s2ℓ=n+m−2has s2k∩has s2ℓ−das s2k∩das s2ℓ−ias s2k∩ias s2ℓ. for example, consider the case as s2k=as s2ℓ. then we have n=m,has s2k∩has s2ℓ=n,das s2k∩das s2ℓ= as there are no deletions, and ias s2k∩ias s2ℓ= as there are no insertions. this gives the distance between the alignments as zero, as would be expected. equation  can be generalised to sequence alignments with more than two sequences. assuming now that ak and aℓ are alignments of m sequences s  …, sm of lengths n  …, nm, we have dak,aℓ=∑p=1m−1∑q=p+1mdasp,sqk,asp,sqℓ, that is, summing all the pairwise alignment distances from equation  this alignment metric if is then normalized and subtracted from to produce a similarity score ssak,aℓ=1−dak,aℓm−1∑t=1mnt the denominator of the fraction, m−1∑t=1mnt, is the normalizing constant, the maximum that the alignment distance dak,aℓ can be. the similarity score is bounded by and with indicating that the sample alignment is identical to the reference alignment. rna secondary structure metrics there are a wealth of available metrics on rna secondary structure. here we use sensitivity, positive predictive value, and f-score. defining true positives as the number of base-pairs correctly predicted, false positives as the number of true base-pairs not predicted, and false negatives as the number of base-pairs predicted which are incorrect, we have sensitivity=tptp+fnppv=tptp+fpf−score=2×tp2×tp+fn+fp the strength of these rna secondary structure accuracy metrics is that they are easy to interpret, and make it straightforward to compare methods across different datasets. an f-score of would represent a structure prediction that was completely correct and an f-score of a structure prediction that only predicted incorrect base-pairs. information entropy as we later develop calculations for information entropy for a set of alignments, here we outline the computation of information entropy for a single alignment. the information entropy h of a probability distribution p with a set of events x is defined as: hp=−∑x∈xpxlog2px. information entropy is a measure for the “spread” of the probability distribution, and has well-defined lower and upper bounds. the minimum entropy of occurs when there is only one outcome with probability and the maximum entropy of log occurs when there are n possible outcomes, each with probability 1/n, that is the uniform distribution. when the base of the logarithm is the entropy is measured in bits. for a probability distribution, an entropy of k bits indicates that the expected value of the information content of observing a single outcome is k bits. in the context of secondary structure prediction, a low entropy therefore indicates that few secondary structures dominate the probability space, whereas a high entropy indicates a more even probability distribution over possible secondary structures. thus, information entropy is a useful single quantity to characterize the underlying probability distribution of secondary structures. the information entropy of the probability distribution over the possible secondary structures generated by a phylo-scfg can be obtained using expected rule frequencies, which can be computed using the inside-outside algorithm. this is outlined here. let the set of all derivations for the input alignment be Φ. since the probability of a derivation d can be written as the product of the scfg production rule probabilities and the phylogenetic probabilities, we can write the total probability t of the grammar producing the input string as t=∑d∈Φpd=∑d∈Φpgdptd, where pg denotes the prior probabilities obtained from the scfg part of the model, and pt are the likelihood factors obtained from the phylogenetic model. conditioning on producing the input string, the normalized probability of a derivation d is pΦd=1tpd=1tpgdptd. consequently, we have that the information entropy of the input alignment under the phylo-scfg model is hΦg=−∑d∈ΦpΦdlog2pΦd, which can be written using equation  as hΦg=log2t−1t∑d∈Φpdlog2pgd−1t∑d∈Φpdlog2ptd, that is, separating out the scfg contribution and the phylogenetic contribution. to calculate the entropy in practice, firstly we use a simplified form of the scfg contribution. for a scfg with set of production rules r, we can write the scfg contribution in terms of the expected production rule frequency, ∑d∈Φpdlog2pgd=∑r∈rlog2pgreusesofr. secondly, we can simplify the phylogenetic contribution. let ra ∈ r be a scfg rule which produces base pairs, and rb ∈ r a scfg rule which produces unpaired bases. define 1d, the indicator function for whether the column pair is emitted from a rule ra, and 1s, the indicator function for whether column i is emitted from a rule rb. finally, define fd as the frequency that rule ra is used in derivation d. then ∑d∈Φpdlog2ptd=∑d∈Φpd∑ralog2ptc|cunpairedfdra+∑rblog2ptc,c′|c,c′fdrb=∑ilog2pti|iunpaired∑d∈Φ1sipd+∑i,jlog2pti,j|i,jpaired∑d∈Φ1di,jpd as ∑d∈Φ1di,jpd is the total probability under the model that positions i and j are emitted from a rule ra, and ∑d∈Φ1sipd is just the total probability under the model that position i is emitted from a rule rb, the quantity ∑d∈Φpdlog2ptd can be computed using the expected rule frequencies obtained through the inside-outside algorithm. methods data statalign dataset to test factors relating to alignment quality and secondary structure prediction quality, a large number of alignment samples from trusted reference alignments with known secondary structures are needed. we have created a curated rna dataset based on the rfam database for the purposes of evaluating the framework. alignments of homologous rna sequences with known consensus secondary structure were extracted from rfam seed alignments. from these, rna families with at least sequences were randomly selected in the section “statalign dataset”. from each family, in a pre-filtering step we removed divergent sequences with long indels, as follows. we defined insertion as consecutive non-gap characters of a sequence in the reference alignment which appear in columns where over 80% of the sequences have gaps. deletions were defined analogously. columns with fraction of gaps between 20% and 80% were regarded ambiguous and ignored. to over-penalize long indels, we applied the super linear score function l × log for indels of length l, indels being defined as above. then, a sequence was removed from a family if its total indel score was beyond and the difference between its indel score and the mean indel score in the family was beyond times the standard deviation of the indel scores in family, i.e. if the sequence had significantly more and/or longer indels than what is representative of the rest of the family. then, sequences were selected at random, and further random selection was done to get pairs, triplets etc. up to sequences in alignments. from these samples of known reference alignment, we could produce many different alignment samples using statalign. for each rna alignment, alignment samples were taken, and the reference alignments were also kept to for comparison. we refer to this dataset throughout as the statalign dataset. random alignment data we also wanted to measure the effect of alignment accuracy on secondary structure independently of alignment method. therefore we created a dataset based on the rna families of the statalign dataset, where alignments were sampled uniformly at various fixed distances from the reference alignment. using the alignment distance measure in eq.  we created a metropolis-coupled mcmc framework that runs several parallel mcmc chains to take alignment samples from the target distribution πa=exp−2da,ar−dt where ar is the reference alignment, d is the target distance to get samples from and t is the temperature of the chain. to improve the mixing properties of the chains we allowed each chain to explore alignments that do not exactly match the specified target distance but then rejected non-exact matches when taking samples from the cold chain. the state space of the markov chains is the set of all possible multiple alignments of the input sequences. alignments that represent the same set of homology statements, and only differ by the order of the alignment columns, are treated as different. the following basic alignment rearrangement moves are iterated: breaking an alignment column into two columns by moving one of its characters into a new column, placing gaps in every other row: cccaaaggg⇌.ccc−aa−a−ggg the exact reverse of the previous, i.e. joining two compatible adjacent columns – one having gaps in all but one row, the other having a gap in that row – to form a single column relocating one character of a column to a gap position of an adjoining column: cccagaga−ttt⇔ cccagag−attt as the space of alignments is vast and the moves are very local, to get a good approximation of uniform sampling, the number of steps that have to be done is on the order of millions, even for small inputs when a single chain is run. to this end, we created a special alignment representation where the above rearrangements can be carried out very efficiently. for moderate input sizes this representation allowed us to take million rearrangement steps in a second. a single chain is sufficient for small alignments, but very slow mixing becomes an issue for practical alignment sizes. we have found that standard parallel tempering techniques effectively speed up mixing if the chain temperatures are chosen correctly. because the optimal temperatures vary significantly depending on reference alignment characteristics and target distribution, a simple acceptance optimization routine was added that tunes chain temperatures to achieve chain swap acceptance ratios between neighboring chains around a pre-set value. we found ratios – to be the most effective. with this framework, running 8– parallel chains was best to maximize the speed of convergence to the uniform distribution as compared to a single chain with sampling times 8–10-fold. the above described framework was utilized to create a dataset consisting of the rna families of statalign dataset, where for each family and each selection of representative sequences from the family, samples were taken at a distance corresponding to a similarity of to the reference alignment, then samples at a similarity of etc., down to a similarity ratio of. we refer to this dataset as the random alignment dataset. extending information entropy to alignment space the information entropy defined for a single alignment contains the length of the alignment as a parameter. attempting to extend the measure to the probability mass over rna secondary structure space, variable alignment length is a concern. for example, if we have two alignments and corresponding secondary structures ))))ccccaaaa−ggggccccaaaaggggccccaaa−aggggccccaaaagggg we would not want to suggest that these alignments give two different secondary structures. consequently, we use a projection method to give alignment column pairing probability matrices the same dimension, so that the matrices can be averaged. for a given set of input sequences, the sequence containing the greatest number of non-gap characters was chosen as the reference sequence. each pairing probability matrix is projected by deleting columns and rows of the matrix corresponding to gap positions in the reference sequence, thus ensuring each matrix corresponding to an alignment sample has the same dimensions. for example, we might start with an matrix, delete row i and column i due to a gap in position i in the reference sequence, then end up with an n × n matrix as required. to calculate information entropy over alignments, we need to be able to calculate the probability of each alignment. however, we cannot calculate the information entropy explicitly, since the probability of a given secondary structure ss is pss=∑a∈alignmentspss|apa, and there is no known efficient way to recurse over all possible alignments. instead, we create an information entropy measure based on samples from the alignment space, and show that, in the sample-size limit, the alignment-sample information entropy tends to the true information entropy. consider alignment samples a …, an from the space of all alignments of m sequences, sampled according to their probability. if we are using statistical alignment, as in statalign, we will be sampling alignments in this fashion. then we have for a column c, once alignments have been projected to the same length, the probability of being unpaired pc|cunpaired=∑i=1npc|cunpairedinalignmentaipai, with an analogous result holding for paired columns. we now define a sample phylogenetic probability ps as the average of the sample phylogenetic probabilities: psc|cunpaired=∑i=1n1npc|cunpairedinalignmentai, to show this sample probability converges to the true probability as sample size tends to infinity, we first note that, rearranging equation 12: psc|cunpaired=∑i=1n1n∑a∈alignmentspc|cunpairedinalignmenta1a=ai, with being the indicator function. hence psc|cunpaired=∑a∈alignmentspc|cunpairedinalignmenta∑i=1n1n1a=ai. taking the limit n → ∞, by the weak law of large numbers psc|cunpaired→∑a∈alignmentspc|cunpairedinalignmentapaasn→∞ =pc|cunpaired as required. now, we have from above that the entropy hΦ of grammar derivations Φ of a grammar g is hΦp=log2t−1t∑d∈Φpgdptdlog2pgd−1t∑d∈Φpgdptdlog2ptd since tree probabilities are the product of unpaired and paired column probabilities in the derivation, the tree probabilities can be recalculated from the sample of alignments. these can then be substituted into the above equation to get an approximation to the information entropy over the space of sampled alignments as well. we refer to this entropy of more than one alignment as the alignment consensus entropy.  alignment quality and predictive accuracy a common question in comparative rna secondary structure prediction is how many sequences are required to get a good structure prediction. this is briefly addressed in, but the sample size is quite small, and only total accuracy is considered. with sequences in the alignment, we assume that no more evolutionary information can be gained by adding further sequences, but when fewer sequences are present, lack of information might yield a poorer structure prediction. instead of considering total accuracy, we wanted to quantify relatively how much accuracy is lost when fewer sequences are present. for example, if an alignment with sequences is predicted with an average f-score of, and an alignment of the same family with sequences is predicted with an average f-score of, then 80% of the accuracy will have been retained, that is the alignment with sequences has a relative f-score of. to investigate how many sequences are needed for a good structure prediction, we took the statalign dataset and considered the relative f-score for each family. almost 100% of the possible f-score was achieved when the alignment contained sequences, for both ppfold and rnaalifold. interestingly, the accuracy of rnaalifold decreased slightly as the number of sequences was increased. this is due to the increased number of non-canonical base-pairs in the alignments, which the thermodynamic method could not predict. ppfold, on the other hand, has a small probability for non-canonical base-pairs, so is not affected by these in the same way. overall these results suggest that sequences are sufficient for approaching maximal predictive accuracy. to consider the effect alignment quality has on rna secondary structure prediction, we took the statalign and random alignment datasets and measured their similarity to the reference alignment using the similarity score above. again, percentage of accuracy retained was calculated by normalizing against the accuracy achieved on the reference alignment. log-scale heatmaps showing the accuracy and percentage of accuracy retained for the statalign dataset and random alignment dataset for ppfold and rnaalifold can be seen in figure  as expected, decreasing alignment quality decreases the accuracy of structure predictions. however, other patterns also emerge from these graphs. firstly, in the statalign dataset, we observe a weaker correlation between alignment distance and accuracy than for the random alignment dataset. this suggests that predictions are better for the statalign dataset on alignments far from the reference alignment. statalign looks to produce alignments with a high likelihood under an evolutionary model, which the random alignments do not consider, and so statalign's alignments could be considered more realistic. this con rms the expectation that statalign's alignments are more useful for rna secondary structure prediction than random alignments. secondly, for the statalign dataset, we see a much higher correlation with the f-score than with the relative f-score. some families of rna consistently produce the same alignment, which skews the graphs. for example, an alignment which consistently has similarity to the reference alignment of, and f-score would give a relative f-score of every time, and would support the correlations seen on each graph. because we can control the spread of distances in the random alignment data set, we don't see this behaviour. as expected, variation comes more with families that produce more variable alignments. in the statalign dataset, this is obscured by those families which produce consistent alignments. lastly, for the random alignment dataset, we see many more zero quality predictions in the case of rnaalifold than in the case of ppfold. this is most easily seen by the larger intensity of red in the main body of the heatmap for ppfold. this suggests that ppfold functions better than rnaalifold when given a low-quality alignment, likely due to its more complete model for molecular evolution. evolutionary distance we also consider the effects of evolutionary distance on rna secondary structure prediction quality. one might expect that there is a “sweet spot” for evolutionary distance– sequences too close to each other do not display enough co-variation to benefit the evolutionary model, but the evolutionary signal might be lost if the distance is too large. to investigate this, we measured evolutionary distance in the phylogenetic trees predicted by ppfold using four different measures: measure 1– mean of all the evolutionary distances, measure 2– standard deviation of all the evolutionary distances, measure 3– maximum evolutionary distance, measure 4– maximum difference between evolutionary distances. all four measures would be expected to be correlated with sequence length, which is well known to correlate with predictive accuracy. to account for this, we considered relative evolutionary distance, similar to the relative predictive accuracy above. a measure was normalized by the average for that family and number of sequences, so that it could be seen whether an alignment had greater or less evolutionary distance than might be expected. we then looked at the correlation between relative evolutionary distance and predictive accuracy. all methods correlated extremely poorly with predictive accuracy and relative predictive accuracy, measures to having r correlations with relative predictive accuracy of,,, and respectively. this suggests that evolutionary distance is not an underlying factor for variation in rna secondary structure prediction, and those other factors, such as those seen in, play a more important role in determining predictive accuracy. alignment distances and maximum posterior decoding given the results concerning accuracy lost as alignment quality decreases, it would be desirable to be able to predict alignment quality, with the hope of predicting structure prediction quality. this has previously been attempted in. first, the sequences were aligned with clustalw. the sequences were then re-aligned using other programs and the similarity between the alignment generated using clustalw to each of the other alignments was measured. the maximum of the similarities, max, was chosen as a predictor of alignment quality. the authors of detected a strong correlation between the true similarity and max. we implemented a modified version of this method. for a given set of input sequences we aligned with both amap and with statalign, obtaining the maximum posterior decoding alignment from statalign. the similarity between the amap alignment and the mpd alignment was used as our predicted similarity measure. this produced a strong correlation between our predicted similarity and true similarity, with an adjusted r-squared value of. we also implemented another method, which calculates an estimate of the expected similarity score using posterior probabilities from the mpd alignment. for each column, we might expect that a posterior probability close to would contribute a score of close to to the similarity measure. so our predicted alignment distance is just the average of the column posterior probabilities. figure  shows an example of the correlation between predicted similarity and true similarity, here giving an adjusted r-squared value of. our new predicted similarity can be calculated efficiently, and is a strong predictor of true similarity to the reference alignment. extending information entropy to alignment space to test the information entropy extension developed above, we calculated the alignment consensus entropy for samples of alignments from the statalign dataset. figure  gives information entropy for different representative rna families from the statalign dataset. on each graph, the information entropy for each of statistical alignment samples is given, as well as the alignment consensus entropy. the leftmost figure is one where the alignment samples were very similar, the rightmost figure where the alignment samples were very different, and the middle figure closer to the median value of alignment sample similarity. for family, we see that the alignment consensus entropy is comparative to the mean entropy of the individual samples. this is expected, as it indicates there is little uncertainty in the alignment. on the other hand, the high-diversity family has much higher alignment consensus entropy than for each individual sample. this is again expected, as the difference in entropies indicates a high uncertainty in alignment. in this way, we can incorporate alignment uncertainty into our understanding of comparative rna secondary structure prediction. predicting secondary structure accuracy given strong correlations between the alignment quality and the structure prediction quality, we might expect that we could find a predictor of structure prediction accuracy. by “integrating out” alignment uncertainty, we may find a reliability score which is more reliable than the one currently reported by the ppfold. to test this, we predicted accuracy for one of the five-sequence alignments of each family and then tested the predicted accuracy against the true accuracy. the ppfold reliability score produced an adjusted r score of when considering correlation with the true f-scores. for our new reliability score, we adjusted the ppfold reliability score to consider only base-pairs, as the f-score considers only base-pairs. we then performed linear regression with the average of the mpd column probabilities, the information entropy of the alignment space, and this pairs-only reliability score against the known f-score measure. this multiple regression improved the reliability score significantly. figure  shows the predicted f-score against the true f-score, for a randomly chosen five-sequence alignments from each family of the statalign dataset. the adjusted r value with the new reliability measure improved to. these results seem to indicate that while alignment quality does affect structure prediction quality, the actual structure prediction model still plays a great role in the overall prediction accuracy. consequently, improving these models, possibly by incorporating other kinds of information, is an area where new research e orts are still needed in rna secondary structure prediction.  in this paper we have explored a number of factors which can contribute to poor rna secondary structure prediction quality. we established a relationship between alignment quality and expected loss of accuracy. furthermore, we provided a method to predict alignment quality based only on statistical alignment samples. while our predictor of accuracy improves on the ppfold reliability score, there is still a large amount of variability unaccounted for, which we therefore suggest comes from the predictive model itself. to consider this further, we extended the information entropy measure for scfgs to consider uncertainty in alignments. the fact that our accuracy predictor did not account for all the variances associated with rna secondary structure prediction, despite good predictors being found for alignment quality and a strong correlation between alignment quality and predictive accuracy, suggests that whilst alignment quality is an important factor, the predictive model itself determines plays a large part in the quality of prediction. given what is shown in figure  for single sequence predictions, that the accuracy of ppfold and rnafold is very variable, it is unsurprising that variances remain. clearly then, further efforts should be put into creating stronger single-sequence models, and then the advantages of evolutionary modelling and additional structural constraints will benefit further. the use of experimental data from new probing experiments as well as more biologically realistic constraints, such as kinetic or co-transcriptional folding, may improve the results of rna secondary structure prediction. competing interests the authors declare they have no competing interests. authors’ contributions jwja, an, and zs formulated the ideas presented here, carried out some of the analyses, and wrote the manuscript. an generated the statalign and random alignment datasets, and mg, pa, and ie generated some of the results for the reliability scores, alignment distances, and predictive algorithm. jwja and zs developed the theoretical framework for extending information entropy. all authors read and approved the final manuscript. supplementary material additional file 1 information on the rna families randomly selected from rfam. click here for file acknowledgements in part, this work was carried out as part of the oxford summer school in computational biology, in conjunction with the department of plant sciences and the department of zoology. funding was provided by the eu cogangs grant. Án would like to thank bbsrc for continued funding through ocisb. jwja would like to thank the epsrc for funding.
in this paper we explore a number of factors which can contribute to poor rna secondary structure prediction quality. however, current approaches to rna secondary structure modelling produce predictions with a high variance in predictive accuracy, and we have little quantifiable knowledge about the reasons for these variances. our predictive accuracy improves on the ppfold reliability score. these results often rely on rna secondary structure predictions. however, there is still variability unaccounted for, which we therefore suggest comes from the rna secondary structure predictive model itself. furthermore, we define two new measures to quantify uncertainty in alignment-based structure predictions. the other measure considers the information entropy for scfgs over a space of input alignments. we establish a quantified relationship between alignment quality and loss of accuracy.
 with the advancement of next-generation sequencing and transcriptomics technologies, regulatory effects involving rna, in particular rna structural changes are being detected. these results often rely on rna secondary structure predictions. however, current approaches to rna secondary structure modelling produce predictions with a high variance in predictive accuracy, and we have little quantifiable knowledge about the reasons for these variances.  in this paper we explore a number of factors which can contribute to poor rna secondary structure prediction quality. we establish a quantified relationship between alignment quality and loss of accuracy. furthermore, we define two new measures to quantify uncertainty in alignment-based structure predictions. one of the measures improves on the “reliability score” reported by ppfold, and considers alignment uncertainty as well as base-pair probabilities. the other measure considers the information entropy for scfgs over a space of input alignments.  our predictive accuracy improves on the ppfold reliability score. we can successfully characterize many of the underlying reasons for and variances in poor prediction. however, there is still variability unaccounted for, which we therefore suggest comes from the rna secondary structure predictive model itself. 