Characterising RNA secondary structure space using information entropy
 Zsuzsanna Sükösd^{1, 2, 3}Email author,
 Bjarne Knudsen^{4},
 James WJ Anderson^{5},
 Ádám Novák^{5, 6},
 Jørgen Kjems^{2, 3} and
 Christian NS Pedersen^{1, 7}
https://doi.org/10.1186/1471210514S2S22
© Sükösd et al.; licensee BioMed Central Ltd. 2013
Published: 21 January 2013
Abstract
Comparative methods for RNA secondary structure prediction use evolutionary information from RNA alignments to increase prediction accuracy. The model is often described in terms of stochastic contextfree grammars (SCFGs), which generate a probability distribution over secondary structures. It is, however, unclear how this probability distribution changes as a function of the input alignment. As prediction programs typically only return a single secondary structure, better characterisation of the underlying probability space of RNA secondary structures is of great interest. In this work, we show how to efficiently compute the information entropy of the probability distribution over RNA secondary structures produced for RNA alignments by a phyloSCFG, and implement it for the PPfold model. We also discuss interpretations and applications of this quantity, including how it can clarify reasons for low prediction reliability scores. PPfold and its source code are available from http://birc.au.dk/software/ppfold/.
Background
The function of RNA molecules is known to depend on their threedimensional structure, which is stabilized by a secondary structure scaffold of basepairing. The secondary structure is defined by hydrogen bonds between nucleotides, which form across the structure for thermodynamic stability and molecular function. Despite its importance, the accurate prediction of RNA secondary structure remains an unsolved challenge in computational biology.
With the advent of nextgeneration sequencing technologies and new methods in transcriptomics, an explosively growing amount of biological RNA data is available in public databases such as Rfam [1] and RNA STRAND [2]. This makes it possible to acquire a large number of RNA alignments to be used in comparative RNA secondary structure predictions. This is especially significant in the case of long RNAs such as RNA viral genomes and long genomic introns, many of which are known to have functional, conserved secondary structures.
Several methods have been established to predict RNA secondary structures from nucleotide sequences. In this paper, we focus entirely on nonpseudoknotted secondary structure prediction. Thermodynamic optimisation based on minimising freeenergy functions has been used to great effect in algorithms such as mfold [3], UNAFold [4] and RNAfold [5]. In a different approach, stochastic contextfree grammars (SCFGs) have also been successfully used to model RNA secondary structure. The Pfold program [6, 7], for example, combines molecular evolution with a lightweight SCFG model (known as a phylogrammar) to predict the consensus secondary structure of RNA alignments, and has in the past shown to be highly accurate for structural alignments [8]. PPfold is a recent multithreaded reimplementation of Pfold [9].
Common to these methods is that they produce a probability distribution over all possible nested secondary structures for the input sequences, but usually only a single, optimal secondary structure is reported to the user. A particularly interesting question is how the underlying distribution changes as a function of input data. Due to the large space of possible secondary structures, however, it is difficult to report useful quantities to describe this. Information entropy is one such measure.
Entropy computations in the context of RNA secondary structure prediction have been considered previously from a thermodynamic perspective, calculating the thermodynamic entropy over both secondary structure space [10, 11] and tertiary structure space [12]. Positional thermodynamic entropy [5] as well as thermodynamic entropy changes in response to basepair mutations [13] have also been computed. Additionally, SCFGbased methods have recently been utilised for calculating both the information entropy of individual basepairs in a singlesequence context [14], and the thermodynamic entropy changes in response to basepair mutations. However, no form of entropy has been computed previously in the case of comparative RNA secondary structure prediction.
The information entropy is a measure for the "spread" of the probability distribution, and has welldefined lower and upper bounds. The minimum entropy of 0 occurs when there is only one outcome with probability 1. For n possible outcomes, the maximum entropy is equal to log_{2} (n) and occurs for the uniform distribution. 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.
In the case of RNA secondary structure prediction based on a semantically unambiguous SCFG, the information entropy of the probability distribution over RNA secondary structures can be computed as the derivational entropy of the SCFG that generates the distribution. We restrict ourselves to semantically unambiguous SCFGs, in order to maintain a onetoone correspondence between SCFG derivations and secondary structures. Thus, throughout this paper we use "information entropy" and "derivational entropy" interchangeably.
Notation
Consider RNA alignments of k sequences (k ≥ 1), with the i'th column denoted c_{ i } ∈ Σ^{ k } = {"A", "C", "G", "U", ""}^{ k } \ {""}^{ k }. A stochastic contextfree phylogrammar (phyloSCFG) on such alignments is a tuple G = ((Σ^{ k }, N, S, R), P), where:

Σ^{ k }forms the (finite) set of terminal symbols

N is a finite set of nonterminal symbols, such that ${\sum}^{k}\cap N=\mathrm{0\u0338}$

S is the start symbol, S ∈ N

R is a finite set of production rules, each rule of the form A → α, A ∈ N and α ∈ (∑^{ k }∪ N)*

P is a function from R to real numbers in the interval [0,1]
In the case of a phylogrammar, P can be interpreted as Bayesian probabilities equal to the product of prior probabilities that only depend on the type of rule being used, and a likelihood factor that is typically derived from a phylogenetic model and is a function of the alignment columns. We will return to this more formally later. Furthermore, we assume the grammar is proper, that is ∀A ∈ N: ∑_{π = (A → α)}P(π) = 1.
where f_{ d } (A → α) is the number of times rule A → α occurs in derivation d.
The grammar is consistent if ∑_{ d } p(d) = 1, where the sum is over all possible derivations of the grammar. In the case of phylogrammars, consistency implies that the total probability of the grammar emitting all alignments of k sequences (of all lengths) is 1.
The I(α) and O(A) variables can also be computed for a particular string C_{ l } of length l, using the insideoutside algorithm in O(l^{3}) time.
for A, B, C ∈ N, c, c^{ ' } ∈ Σ^{ k }. Apart from generating empty strings, all SCFGs modelling nested RNA secondary structures can be written in double emission normal form, so the methods presented here can be adapted to RNA secondary structure grammars of all types. Additionally, they can also be adapted for specific mildly contextsensitive RNA grammars that generate specific types of pseudoknots.
where P_{ G }(A → c) only depends on A, and ${P}_{T}\left(c\mathit{\u010b}\right)$ is the likelihood of observing column c under the phylogenetic model, assuming that it is unpaired in the consensus structure (denoted by $\mathit{\u010b}$).
where P_{ G }(A → cBc') only depends on A and B, and ${P}_{T}\left(c,\phantom{\rule{2.77695pt}{0ex}}{c}^{\prime}\hat{c,{c}^{\prime}}\right)$ is the likelihood of observing column pair c, c^{ ' } under the phylogenetic model, assuming that they are paired with each other in the consensus structure (denoted by $\hat{c,{c}^{\prime}}$).
where P_{ G }(A → BC) only depends on A, B and C.
for r_{ a } ∈ R of Type 1, r_{ b } ∈ R of Type 2, r_{ c } ∈ R of Type 3.
Given a single RNA alignment, there are typically a large number of possible derivations, each corresponding to a possible secondary structure for the alignment. In the rest of this work, we restrict ourselves to this space of derivations, which we characterize by its derivational entropy as described below.
Results and discussion
Algorithm
which can be computed using the expected rule frequencies obtained from the insideoutside algorithm (cf. equations 7 and 8). This can be seen by noting that Equations 4 and 7 still hold when the set of derivations is restricted to a subset, as opposed to the entire space.
We observe that ${\sum}_{d\in \Phi}{\mathbb{1}}^{d}\left(i,\phantom{\rule{2.77695pt}{0ex}}j\right)p\left(d\right)$ is just the total probability under the model that positions i and j form a pair (i.e. they are emitted from rule type 2), and ${\sum}_{d\in \Phi}{\mathbb{1}}^{s}\left(i\right)p\left(d\right)$ is just the total probability under the model that position i is unpaired (i.e. it is emitted from rule type 1). The quantity ∑_{d∈Φ}p(d) log_{2}(p_{ T }(d)) can therefore also be computed using the expected rule frequencies.
Once the values of the insideoutside variables have been calculated for an input string of length n, the expected rule frequencies can be computed in O(n) time for rules of Type 1, O(n^{2}) time for rules of Type 2, and O(n^{3}) time for rules of Type 3. As the time complexity of the insideoutside algorithm is also O(n^{3}), the computation of the entropy over the possible derivations of the input string does not increase the time complexity of RNA secondary structure prediction.
Interpretation of the derivational entropy
H_{ max }(l) provides an upper bound on the value of the information entropy for an RNA alignment of length l, and can aid the user in the interpretation of the entropy corresponding to a particular input alignment.
Using information entropy to interpret low reliability scores
The probability of column i remaining unpaired is $P\left(\dot{i}\right)=1{\sum}_{\begin{array}{c}\hfill i,j\hfill \\ \hfill i\ne j\hfill \end{array}}P\left(\hat{i,j}\right)$. These probabilities function as "reliability scores" for every predicted basepair or unpaired nucleotide in the structure. The overall reliability score for the secondary structure can be computed as the average of the reliability scores of all positions.
Importantly, while the reliability scores depend both on the structures and their probabilities, the derivational entropy is only a function of probabilities, and does not depend on the similarity of the structures to each other. Derivational entropy therefore provides complementary information to reliability scores. For example, if the reliability scores predict a low accuracy, the entropy can help reveal the underlying reasons. A low reliability score can be observed in different situations, for example (a) if there is insufficient structure signal, so there are no structures of high probabilities and the probability distribution is "spread", or (b) if there are two or more possible (topologically different) structures of high probabilities, so the probability distribution has several "peaks". Entropy will be high in the first case, but low in the second case, and can therefore be used to distinguish the two situations from each other.
To illustrate this with a practical example, a PPfold prediction of the secondary structure of the random nucleotide sequence:
GACCAAACGCAGCCAGCGTCACTGTAGGATTTAAA
ACCGAGGGAATGCCGTCAGTAGGGTCGGGTTTAAC
reveals that the underlying probability distribution has an entropy of 28.55 bits, with an average reliability score 0.65 for the final predicted structure. In comparison, the following combinatorial sequence of the same length:
GGGGAAACCCCAAAGGGGAAACCCCAAAGGGGAAA
CCCCAAAGGGGAAACCCCAAAGGGGAAACCCCAAA
shows a significantly lower entropy of 18.82 bits, at the same time as a low average reliability score of 0.56. Despite the low reliability scores in both cases, the entropy value correctly reveals that the combinatorial sequence has a relatively small number of very different secondary structures dominating the probability space, whereas the probabilities are more uniformly distributed over a large number of possible structures in the case of the random sequence.
Relationship to prediction accuracy
An important question with respect to comparative RNA secondary structure prediction is how the accuracy of predictions varies with the quality of input alignment, and to what extent the accuracy of a predicted structure can be predicted. Reliability scores and entropy both measure variation in the secondary structure space, so both are expected to be correlated with prediction accuracy. We note, however, that a natural limitation of both entropy and reliability scores is that they are computed under the model, which effectively assumes that the model itself is an accurate description of the biological folding processes. If this is not the case, high confidence values computed under the model can still correspond to low prediction accuracies in reality.
We observe that despite the PPfold reliability scores generally suggest somewhat higher prediction accuracies than what was actually observed, they convey absolute information about the accuracy of predictions. Instead, entropy functions as a relative accuracy measure, when comparing several structure predictions for the same sequence. Entropy can, therefore, complement the currently existing measures for an increased understanding of RNA structure variability.
Lastly, we note that there is large variation in the reliability measures reported by different RNA structure prediction programs. By contrast, information entropy is a methodindependent measure: for the same alignment, it will have the same, welldefined interpretation, regardless of the method producing the underlying probability distributions. Additionally, all classic information theoretic results can be applied to it. Information entropy can therefore be used to compare the behaviour of any secondary structure prediction methods that produce a probability distribution over RNA secondary structures, in a way that is not possible with currently existing methods.
Future perspectives include the computation of measures similar to entropy, such as the selfinformation of particular outcomes, or the KullbackLeibner divergence of distributions to compare constrained and unconstrained models for RNA folding. We also expect that entropy may provide a number of other possible future applications in RNA secondary structure prediction.
Conclusions
The information entropy of the probability distribution generated by phylogrammars can be computed efficiently from the insideoutside variables, and has been implemented as part of PPfold. Information entropy is a welldefined characteristic of the underlying probability distribution, which complements the reliability values already reported by algorithms for an increased understanding of RNA structure variability. It is also a methodindependent measure of prediction certainty, providing theoretical advantages over existing methods.
Methods
Implementation
We note that the KH99 grammar is not in double emission normal form, as it includes the rule S → L. Even though the grammar can be expressed in the desired form, adapting the algorithm to include this additional rule is both straightforward and highly efficient; this is therefore what we have implemented.
In the case of the KH99 grammar, there is a maximum of one bifurcating rule originating from each nonterminal symbol. It is also known that the expected frequency of a nonterminal symbol can be computed in O(n^{2}) time for any SCFG, and the expected frequencies of rules from the same nonterminal symbol sum to the expected frequency of the nonterminal symbol. Hence, in the case of the KH99 grammar, the time complexity of the computation of the derivational entropy (given the insideoutside variables) could further be reduced to O(n^{2}).
As described in the Results and Discussion section, the value of the entropy depends on the length of the alignments. Hence, small adaptations were made in the algorithm to be able to compare the entropies of alignments that include a particular sequence. The default option in PPfold is to remove columns where fewer than 75% of the sequences have nucleotides. In the case of entropy computation, this is replaced with removing only the columns where the selected sequence has gaps. This ensures that all alignments that include the selected sequence have an equal prediction length.
Test sequences and alignments
Random sequences of different nucleotide compositions were generated using the online FaBox tool [19].
For the generation of alignments, among families with an experimentally verified secondary structure in the Rfam database, we chose two families (RF01960  Eukaryotic small subunit ribosomal RNA, and RF00177  Bacterial small subunit ribosomal RNA), for which our initial tests indicated that a particularly wide range of structure prediction accuracies may be achieved depending on the choice on the sequences included in the alignment. Starting with a sequence of interest in both families, we constructed our datasets by randomly adding sequences from the family alignment one by one, up to a final size of maximum 15 sequences. For each alignment size (between 1 and 15 sequences in the alignment), the process was repeated 50 times. This way we obtained 1+14 × 50, not necessarily distinct cases per family, with 115 sequences each (including the starting sequence as a standalone case). The alignments of the selected sequences and the reference secondary structures thereof were adapted from those in Rfam, by deleting gaponly columns and any base pairs involved with those columns.
Comparing accuracies
where TP is the number of correctly predicted basepairs.
Declarations
Declarations
This research has been funded by the Danish Council for Strategic Research under the project #09061856. JWJA would like to thank the EPSRC for funding. ÁN would like to thank BBSRC for continued funding through OCISB. CNSP is supported by Centre for Membrane Pumps in Cells and Disease (PUMPKIN), a Centre of Excellence funded by the Danish National Research Foundation. The publication costs for this article were also funded by PUMPKIN.
This article has been published as part of BMC Bioinformatics Volume 14 Supplement 2, 2013: Selected articles from the Eleventh Asia Pacific Bioinformatics Conference (APBC 2013): Bioinformatics. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/14/S2.
Authors’ Affiliations
References
 GrithsJones S, Moxon S, Marshall M, Khanna A, Eddy SR, Bateman A: Rfam: annotating noncoding RNAs in complete genomes. Nucleic acids research. 2005, 33 (suppl 1): D121D124.Google Scholar
 Andronescu M, Bereg V, Hoos H, Condon A: RNA STRAND: The RNA Secondary Structure and Statistical Analysis Database. BMC Bioinformatics. 2008, 9: 34010.1186/147121059340.PubMed CentralView ArticlePubMedGoogle Scholar
 Zuker M, Stiegler P: Optimal computer folding of large RNA sequences using thermodynamics and auxiliary information. Nucleic Acids Research. 1981, 9: 133148. 10.1093/nar/9.1.133.PubMed CentralView ArticlePubMedGoogle Scholar
 Markham NR, Zuker M, Keith JM, Walker JM: UNAFold. 2008, Humana Press. Bioinformatics, 331. [Methods in Molecular Biology; SP: 3]Google Scholar
 Hofacker I, Fontana W, Stadler P, Bonhoeer L, Tacker M, Schuster P: Fast folding and comparison of RNA secondary structures (The Vienna RNA Package). Chemical Monthly. 1994, 125 (2): 167188. 10.1007/BF00818163. [SP: 167]View ArticleGoogle Scholar
 Knudsen B, Hein J: RNA secondary structure prediction using stochastic contextfree grammars and evolutionary history. Bioinformatics. 1999, 15: 446454. 10.1093/bioinformatics/15.6.446.View ArticlePubMedGoogle Scholar
 Knudsen B, Hein J: Pfold: RNA secondary structure prediction using stochastic contextfree grammars. Nucleic Acids Res. 2003, 31 (13): 34233428. 10.1093/nar/gkg614.PubMed CentralView ArticlePubMedGoogle Scholar
 Gardner P, Giegerich R: A comprehensive comparison of comparative RNA structure prediction approaches. BMC Bioinformatics. 2004, 5: 14010.1186/147121055140.PubMed CentralView ArticlePubMedGoogle Scholar
 Sukosd Z, Knudsen B, Vaerum M, Kjems J, Andersen E: Multithreaded comparative RNA secondary structure prediction using stochastic contextfree grammars. BMC Bioinformatics. 2011, 12: 10310.1186/1471210512103.PubMed CentralView ArticlePubMedGoogle Scholar
 Mathews D, Sabina J, Zuker M, Turne D: Expanded Sequence Dependence of Thermodynamic Parameters Improves Prediction of RNA Secondary Structure. J Mol Biol. 1999, 288: 911940. 10.1006/jmbi.1999.2700.View ArticlePubMedGoogle Scholar
 Markham NR, Zuker M: UNAFold: software for nucleic acid folding and hybridization. Methods Mol Biol. 2008, 453: 331. 10.1007/9781603274296_1.View ArticlePubMedGoogle Scholar
 Liu L, Chen SJ: Computing the conformational entropy for RNA folds. The Journal of chemical physics. 2010, 132 (23): 23510410.1063/1.3447385.PubMed CentralView ArticlePubMedGoogle Scholar
 Kiryu H, Asai K: Rchange: Algorithms for computing energy changes of RNA secondary structures in response to base mutations. Bioinformatics. 2012, 28: 10931101. 10.1093/bioinformatics/bts097.View ArticlePubMedGoogle Scholar
 Wang Y, Manzour A, Shareghi P, Shaw TI, Li YW, Malmberg RL, Cai L: Stable stem enabled Shannon entropies distinguish noncoding RNAs from random backgrounds. Computational Advances in Bio and Medical Sciences (ICCABS), 2011 IEEE 1st International Conference on. 2011, 184View ArticleGoogle Scholar
 Lari K, Young SJ: The estimation of stochastic contextfree grammars using the InsideOutside algorithm. Computer Speech & Language. 1990, 4: 3556.View ArticleGoogle Scholar
 Nederhof M, Satta G: KullbackLeibler distance between probabilistic contextfree grammars and probabilistic finite automata. Computer Speech & Language, Association for Computational Linguistics. 2004, 71es, [http://dl.acm.org/citation.cfm?id=1220366]Google Scholar
 Anderson JWJ, Tataru P, Staines J, Hein J, Lyngso R: Evolving Stochastic ContextFree Grammars for RNA Secondary Structure Prediction. BMC Bioinformatics. 2012, 13: 7810.1186/147121051378.View ArticleGoogle Scholar
 Stein PR, Waterman MS: On some new sequences generalizing the Catalan and Motzkin numbers. Discrete Mathematics. 1978, 26: 261272.View ArticleGoogle Scholar
 Villesen P: FaBox: an online toolbox for fasta sequences. Molecular Ecology Notes. 2007, 7 (6): 965968. 10.1111/j.14718286.2007.01821.x.View ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.