 Research
 Open Access
 Published:
Maximum expected accuracy structural neighbors of an RNA secondary structure
BMC Bioinformatics volume 13, Article number: S6 (2012)
Abstract
Background
Since RNA molecules regulate genes and control alternative splicing by allostery, it is important to develop algorithms to predict RNA conformational switches. Some tools, such as paRNAss, RNAshapes and RNAbor, can be used to predict potential conformational switches; nevertheless, no existent tool can detect general (i.e., not family specific) entire riboswitches (both aptamer and expression platform) with accuracy. Thus, the development of additional algorithms to detect conformational switches seems important, especially since the difference in free energy between the two metastable secondary structures may be as large as 1520 kcal/mol. It has recently emerged that RNA secondary structure can be more accurately predicted by computing the maximum expected accuracy (MEA) structure, rather than the minimum free energy (MFE) structure.
Results
Given an arbitrary RNA secondary structure S_{0} for an RNA nucleotide sequence a = a_{1},..., a_{ n }, we say that another secondary structure S of a is a kneighbor of S_{0}, if the base pair distance between S_{0} and S is k. In this paper, we prove that the Boltzmann probability of all kneighbors of the minimum free energy structure S_{0} can be approximated with accuracy ε and confidence 1  p, simultaneously for all 0 ≤ k < K, by a relative frequency count over N sampled structures, provided that $N>N\left(\epsilon ,p,K\right)=\frac{{\Phi}^{1}{\left(\frac{p}{2K}\right)}^{2}}{4{\epsilon}^{2}}$, where Φ(z) is the cumulative distribution function (CDF) for the standard normal distribution. We go on to describe the algorithm RNAborMEA, which for an arbitrary initial structure S_{0} and for all values 0 ≤ k < K, computes the secondary structure MEA(k), having maximum expected accuracy over all kneighbors of S_{0}. Computation time is O(n^{3} · K^{2}), and memory requirements are O(n^{2} · K). We analyze a sample TPP riboswitch, and apply our algorithm to the class of purine riboswitches.
Conclusions
The approximation of RNAbor by sampling, with rigorous bound on accuracy, together with the computation of maximum expected accuracy kneighbors by RNAborMEA, provide additional tools toward conformational switch detection. Results from RNAborMEA are quite distinct from other tools, such as RNAbor, RNAshapes and paRNAss, hence may provide orthogonal information when looking for suboptimal structures or conformational switches. Source code for RNAborMEA can be downloaded from http://sourceforge.net/projects/rnabormea/ or http://bioinformatics.bc.edu/clotelab/RNAborMEA/.
Background
RNA secondary structure conformational switches play an essential role in a number of biological processes, such as regulation of viral replication [1] and of viroid replication [2], regulation of R1 plasmid copy number in E. coli by hok/sok system [3], transcriptional and translational gene regulation in prokaryotes by riboswitches [4], regulation of alternative splicing in eukaryotes [5], and stressresponsive gene regulation in humans [6], etc. Due to the biological importance of conformational switches, several groups have developed algorithms that attempt to recognize switches  in particular, thermodynamicsbased methods such as paRNAss[7], RNAshapes[8], RNAbor[9], as well as an approach using the second eigenvalue of the Laplacian matrix [10].
Riboswitches are portions of the 5' untranslated region (UTR) of messenger RNAs, experimentally known to regulate genes in bacteria by allostery [4], and to regulate alternative splicing of the gene NMT1 in the eukaryote Neurospora crassa [5]. Riboswitches are composed of a 5' aptamer and a 3' expression platform. Since the aptamer binds to a specific ligand with high affinity (K_{ D } ≈ 5 nM), thus triggering the conformational change of the expression platform upon ligand binding [11], its sequence and secondary structure tend to be highly conserved. In contrast, there is lower sequence for the expression platform, which forms a bistable switch, effecting gene regulation by premature abortion of transcription (as in guanine riboswitches [12]), or by sequestering the ShineDalgarno sequence (as in thiamine pyrophosphate riboswitches [13]). Due to the conserved sequence and secondary structure within the aptamer, all existent algorithms (to the best of our knowledge), such as [14–16], attempt only to detect riboswitch aptamers, without the expression platform. In addition to these specific algorithmic approaches, more general computational tools that rely on stochastic context free grammars, such as Infernal[17] and CMFinder[18], have been trained to recognize riboswitch aptamers; in particular, Infernal was used to create the Rfam database [19], which includes 14 families of riboswitch aptamers.
Since current riboswitch detection algorithms do not attempt to predict the location of the expression platform, we have developed tools, RNAborSample and RNAborMEA, described in this paper, which yield information concerning alternative, or suboptimal, structures of a given RNA sequence. These tools can suggest the presence of a conformational switch; however, much more work must be done to actually produce a riboswitch gene finder, part of the difficulty due to the fact that riboswitch aptamers contain pseudoknots that cannot be captured by secondary structure.
In previous work [20, 21], we described a novel program RNAbor to predict RNA conformational switches. For a given secondary structure S of a given RNA sequence s, the secondary structure T of s is said to be a kneighbor of S, if the base pair distance between S and T is k. (Base pair distance is the minimum number of base pairs that must be either added or removed, in order to transform the structure S into T.) Given an arbitrary initial structure S_{0}, for all values 0 ≤ k < K, the program RNAbor[20], computes the secondary structure MFE(k), having minimum free energy over all kneighbors of S_{0}. (Note that K ≤ 2 · n, since the base pair distance between any two secondary structures of a length n RNA sequence is at most 2 · n.) As well, RNAbor computes for each value 0 ≤ k ≤ K, the Boltzmann probability ${p}_{k}=\frac{Z\left(k\right)}{Z}$, where Z(k) is the sum of all Boltzmann factors exp(E(S)/RT) of all structures S having base pair distance k from an initially given structure S_{0}, and where the partition function Z is the sum of all Boltzmann factors of all secondary structures of the given RNA sequence. Here E(S) is the free energy of secondary structure S, with respect to the Turner energy model [22, 23], R = 0.001987 kcal mol^{1} K^{1} is the universal gas constant, and T is absolute temperature. In the case that S_{0} is the minimum free energy structure, the existence of one or more 'peaks', or values k ≫ 0, where p_{ k } is relatively large, suggests that there are two or more low energy structures having large base pair distance k from S_{0}  i.e., a potentially distinct metastable structure, as shown in Figure 1.
In [24], Do et al. introduced the notion of maximum expected accuracy (MEA) secondary structure, determined as follows: (i) compute base pairing probabilities p(i, j) using a trained stochastic context free grammar; (ii) compute probabilities $q\left(i\right)=1{\sum}_{i<j}p\left(i,j\right){\sum}_{j<i}p\left(j,i\right)$ that position i does not pair; (iii) using a dynamic programming algorithm similar to that of Nussinov and Jacobson [25], determine that secondary structure S having maximum score ${\sum}_{\left(i,j\right)\in S}2\alpha \cdot p\left(i,j\right)+{\sum}_{i\text{unpaired}}\beta {q}_{i}$, where the first sum is over paired positions (i, j) of S and the second sum is over positions i located in loop regions of S, and where α, β > 0 are parameters with default values 1. Subsequently Kiryu et al. [26] computed the MEA structure by replacing the stochastic context free grammar computation of base pairs in (i) by using McCaskill's algorithm [27], which computes the Boltzmann base pairing probabilities
The sum in the numerator is taken over all secondary structures of the given RNA sequence, that contain base pair (i, j), while the sum in the denominator is taken over all secondary structures of the given RNA sequence. Thus p(i, j) is the sum of the Boltzmann factors of all secondary structures that contain the fixed base pair (i, j), divided by the partition function, which latter is the sum of Boltzmann factors of all secondary structures. In fact, Kiryu et al. [26] describe an algorithm to compute the MEA structure common to all RNAs in a given alignment. Later, Lu et al. [28] rediscovered Kiryu's method; in addition, Lu et al. computed suboptimal MEA structures by implementing an analogue of Zuker's method [29].
Our motivation in developing both RNAborSample and RNAborMEA, was to simplify and improve our previous software, RNAbor, in detecting conformational switches. Since RNAbor computes the minimum free energy structure, MFE(k), over all structures having base pair distance k from an initially given structure S_{0}, a complex portion of the code in RNAbor concerns the retrieval of free energy parameters from the Turner model [22, 23]. The idea of RNAborMEA was to compute the base pairing probabilities p(i, j)  see equation (1)  by McCaskill's algorithm using RNAfold, then to compute the maximum expected accuracy structure, MEA(k), which needs no retrieval of energy parameters, and which we hoped would be very similar to the MFE(k) structure, in light of previously mentioned results [26, 28]. Surprisingly, it turns out that MEA(k) structures are quite different from MFE(k) structures, as shown later in one of the figures.
In this paper, we begin by showing rigorously how to approximate the output of RNAbor by frequency counts from sampling, using Sfold[30]. We then extend the MEA technique to compute the maximum expected accuracy kneighbor of a given RNA secondary structure S_{0}; i.e., that secondary structure which has maximum expected accuracy over all structures that differ from S_{0} by exactly k base pairs. By analyzing the family of purine riboswitches, obtained by retrieving full riboswitch sequences (aptamer and expression platform) from corresponding EMBL genomic data, by extending the aptamers from the seed alignment of Rfam family RF00167 [31], we show that our software RNAborMEA produces strikingly different results from other software that produce suboptimal structures (RNAbor, RNAborSample, RNAlocopt, RNAshapes, UNAFold).
Since the detection of computational switches remains an open problem, despite the success of some tools such as RNAshapes and RNAbor, we feel the addition of the tool RNAborMEA could prove useful, since it appears to be orthogonal to all other methods of generating suboptimal secondary structures.
Results and discussion
In this paper, we describe the following new results, discussed in the 'Methods' section in greater detail with attendant definitions of unexplained concepts.

1.
We describe a Python script RNAborSample that approximates the output ${p}_{k}=\frac{{Z}_{k}}{Z}$ of RNAbor by frequency counts ${\widehat{p}}_{k}$ from sampled structures, for all 0 ≤ k ≤ 2n, using Sfold[30], or RNAsubopt p[32].

2.
We prove that for any desired accuracy 0 < ε and probability 0 < α < 1, if at least
$$N\left(\epsilon ,p,K\right)=\frac{{\Phi}^{1}{\left(\frac{p}{2K}\right)}^{2}}{4{\epsilon}^{2}}$$(2)
structures are sampled, then
for all 0 ≤ k < K; i.e., RNAborSample furnishes estimates ${\widehat{p}}_{k}$ of p_{ k }, for all 0 ≤ k < K, which with confidence 1  α are within ε of the actual values p_{ k }. Here, Φ(z) is the cumulative distribution function (CDF) for the standard normal distribution.

3.
We develop an algorithm, RNAborMEA, running in time O(n^{3} · K^{2}) and space O(n^{2} · K), which computes simultaneously for all 0 ≤ k ≤ K, the maximum expected accuracy kneighbors of a given RNA secondary structure S_{0}; i.e., for each 0 ≤ k ≤ K, RNAborMEA computes that structure S_{ k } which has maximum expected accuracy over all structures that differ from S_{0} by exactly k base pairs. The algorithm RNAborMEA additionally computes, for each 0 ≤ k ≤ K, the pseudo partition function
$$\stackrel{\u0303}{{Z}_{k}}=\sum _{\left\{S:{d}_{\text{BP}}\left(S,{S}_{0}\right)=k\right\}}\text{exp}\left(MEA\left(S\right)/RT\right).$$
Moreover, RNAborMEA allows the user to stipulate (partial) hard constraints, that stipulate whether particular nucleotides are unpaired, or basepair with certain other nucleotides. The implementation of hard constraints follows ideas from Mathews [33], albeit suitably modified to simultaneously consider all kneighbors, for 0 ≤ k ≤ K.
We now describe the 13 figures and 4 tables, corresponding to computational experiments performed with RNAborSample and RNAborMEA. These tables and figures are only briefly described, and we refer the reader to the captions of the figures and tables, which explain the results in greater detail.
Figure 1 illustrates the presence of two peaks, corresponding to the Boltzmann probability of each of the metastable structures for a 27 nt bistable switch previously considered by Hofacker et al. Figure 2 displays the Boltzmann probabilities p_{ k } from RNAbor, Boltzmann probabilities estimates ${\widehat{p}}_{k}$ from RNAborSample for the SAM riboswitch aptamer with GenBank accession code AP004597.1/1189411904. Clearly, probability estimates ${\widehat{p}}_{k}$ are close to actual values p_{ k }. The figure additionally shows probabilities r_{ k } from our software RNAlocopt[34], computed by ${r}_{k}=\frac{{Z}_{k}\left(LO\right)}{Z\left(LO\right)}$, where Z(LO) is the sum of Boltzmann factors of all locally optimal secondary structures, and Z_{ k }(LO) is the sum of all locally optimal kneighbors of S_{0}. A secondary structure S is said to be locally optimal, if its energy does not decrease by the addition or removal of a single (valid) base pair; i.e., E(S ∪ {(x, y)}) ≥ E(S), and E(S  {(x, y)}) ≥ E(S). Figure 3 displays the experimentally determined GENE ON and GENE OFF structures of an XPT guanine riboswitch from B. subtilis, taken from [35]. Figure 4 shows the outputs of RNAborMEA, RNAbor, and RNAshapes, which are most similar to the GENE ON structure from the previous Figure 3. Figures 5 and 6 determine the structural simlarity, as measured by the program NestedAlign[36], between that structure output by RNAborMEA (as well as structures output by RNAbor, RNAborSample, RNAlocopt, RNAshapes, and UNAFold), which are most similar to the XPT purine riboswitch, displayed in Figure 3. Figure 5 determines the structural similarity to the GENE ON structure (left panel of Figure 3), while Figure 6 determines the structural similarity to the GENE OFF structure (right panel of Figure 3). None of the structural neighbors, or sampled structures, are identical to the GENE ON or GENE OFF structures; however, there are some candidates that bear some resemblance to those structures. At this point, we can say that RNAborSample and RNAborMEA are methods that generate suboptimal structures, some of which may be similar to the metastable structures of a conformational switch; however, much additional work is necessary before a robust method can be developed to detect conformational switches.
Figure 7 shows that the MEA(k) structural neighbors, as computed by RNAborMEA, are very different than the MFE(k) structural neighbors, as computed by RNAbor. At present, such computational experiments show RNAborMEA computes suboptimal structures, which seem to share (chimeric) similarities between parts of low energy structures, but which themselves do not have very low energies. Such suboptimal structures appear to be 'orthogonal' to those output by all other methods, such as Sfold, RNAbor, RNAborSample, RNAlocopt, RNAshapes, UNAFold). Figure 8 displays the output of RNAborMEA, given the sequence of a TPP riboswitch with EMBL accession code AF269819/18111669. In this instance, RNAborMEA found two low energy structures having large base pair distance from each other. (Other computational experiments did not yield such a good example.) Figure 9 displays the free energy and maximum expected accuracy scores, for each of the kneighbors of the given TPP riboswitch sequence, just described in Figure 8. Figures 10 and 11 present the pseudocode for the RNAborMEA algorithm, which given an RNA sequence a_{1}, .. .,a_{ n } and initial structure S_{0}, computes the MEA(k) structure and pseudo partition function $\stackrel{\u0303}{{Z}_{k}}$, for each 0 ≤ k ≤ K in time O(n^{3} · K^{2}) and space O(n^{2} · K). Figure 12 presents pseudocode for the O(n^{2}) algorithm to sample structures from the ensemble of structures having high MEA scores  a maximum expected accuracy analogue of the sampling algorithm Sfold[30]. Figure 13 displays the pseudoBoltzmann probabilities $\stackrel{\u0303}{{p}_{k}}=\frac{{\stackrel{\u0303}{Z}}_{k}}{Z}$ for two small RNA sequences. While temperature T has a natural significance, when computing Boltzmann probabilities ${p}_{k}=\frac{{Z}_{k}}{Z}$, there is no natural meaning of temperature T, when computing pseudo Boltzmann factors exp(MEA(S)/RT), and indeed very different curves can be obtained with different temperatures.
We now briefly describe Tables 1, 2, 3, 4. Table 1 provides some sample sizes N, computed by the formula from equation (2), for an ε approximation of Boltzmann probabilities p_{ k }, 0 ≤ k < K, with 1  α confidence level. Tables 2 and 3 provide the numerical values for the earlier described Figures 5 and 6, where the NestedAlign structural similarity is computed for the most similar kneighbor, determined by RNAborMEA, RNAborSample and RNAlocopt. Table 4 presents the number of times that each of the methods RNAborMEA, RNAbor, RNAborSample, RNAlocopt, RNAshapes, UNAFold output the most similar structure to the GENE ON resp. GENE OFF structure for the XPT purine riboswitch described in Figure 3. This computational experiment was performed for all RNA sequences in the seed alignment of the Rfam purine riboswitch family RF00167 [31]. This table shows that RNAborMEA and RNAbor both outperform any other method in determining structures similar to the GENE OFF XPT structure; however, RNAborMEA uniquely outperforms all methods, including RNAbor, in determining structures similar to the GENE ON XPT structure. One of the reasons for this excellent result is that unlike other methods, RNAborMEA does not look for low energy structures, but rather for maximum expected accuracy structures.
The figures and tables show, in summary, that RNAborMEA provides useful suboptimal structures, which may be closer to metastable structures of a conformational switch than more traditional methods, which rely on searching for low energy structures.
Conclusions
We have applied the notion of maximum expected accuracy within the context of structural neighbors of a given RNA sequence a_{1}, .. ., a_{ n } and structure S_{0}. Our software RNAborMEA not only computes the structures MEA(k) having maximum expected accuracy over all structures S, whose base pair distance d_{BP}(S_{0}, S) is equal to k. In addition, RNAborMEA allows the user to enter structural constraints, which specify partial secondary structures required of all MEA(k) structures, if so desired. Additionally, RNAborMEA computes an analogue of the temperaturedependent partition function, defined by
and
Here, the expected accuracy score σ is defined by
where first sum is taken over all base pairs (i, j) belonging to S, and the second sum is taken over all unpaired positions in S, and where p_{i,j}[resp. q_{ i }] is the probability that i, j are paired [resp. i is unpaired] in the ensemble of low energy structures, as computed by McCaskill's algorithm [27]. Finally, RNAborMEA allows the user to sample structures from the maximum expected accuracy ensemble, in a fashion analogous to DingLawrence sampling from the low energy Boltzmann ensemble, as in Sfold[30].
Our preliminary investigations have not indicated a clear application of the partition function analogue, though it may be construed to provide a representation of the temperaturedependent mixing of various structures having large score σ. On the other hand, in computational experiments reported in the Results Section, it appears that RNAborMEA produces nearoptimal structures that are closer to the biologically functional structures, in the case of conformational switches that are difficult to predict by any method.
Indeed, in 18 [resp. 11] out of 34 instances, RNAborMEA produced the secondary structure most structurally similar to the experimentally determined XPT GENE ON [resp. GENE OFF] structure, as measured by NestedAlign[36]. See Table 4. Since there appears to be little to no correlation between the structures MFE(k) output by RNAbor[20] and the structures MEA(k) output by our current program RNAborMEA, it appears that RNAborMEA yields a signal that is orthogonal and complementary to that provided by stateoftheart thermodynamics software, such as UNAFold, RNAfold, RNAstructure, Sfold, RNAshapes, RNAbor, etc. For these reasons, we feel that RNAborMEA has a certain value, along with the programs UNAFold, RNAfold, RNAstructure, Sfold, RNAshapes, RNAbor, etc. when producing suboptimal structures. RNAborMEA is written in C and available at http://sourceforge.net/projects/rnabormea/ and http://bioinformatics.bc.edu/clotelab/RNAborMEA/.
Methods
Preliminaries
Recall the definition of RNA secondary structure.
Definition 1 A secondary structure S on RNA sequence a_{1}, .. ., a_{ n } is defined to be a set of ordered pairs (i, j), such that 1 ≤ i < j ≤ n and the following are satisfied.
1. WatsonCrick or GU wobble pairs: If (i, j) belongs to S, then pair (a_{ i }, a_{ j }) must be one of the following canonical base pairs: (A, U), (U, A), (G, C), (C, G), (G, U), (U, G).
2. Threshold requirement: If (i, j) belongs to S, then j  i > θ, where θ, generally taken to be equal to 3, is the minimum number of unpaired bases in a hairpin loop; i.e., there must be at least θ unpaired bases in a hairpin loop.
3. Nonexistence of pseudoknots: If (i, j) and (k, ℓ) belong to S, then it is not the case that i < k < j < ℓ.
4. No base triples: If (i, j) and (i, k) belong to S, then j = k; if (i, j) and (k, j) belong to S, then i = k.
The preceding definition provides for an inductive construction of the set of all secondary structures for a given RNA sequence a_{1}, .. ., a_{ n }. For all values of d = 0, .. ., n and all values of i = 1, .. ., n  d, the collection ${S}_{i,i+d}$of all secondary structures for a_{ i }, .. ., a_{i+d}is defined as follows. If 0 ≤ d ≤ θ, then ${S}_{i,i+d}=\left\{\varnothing \right\}$; i.e., the only secondary structure for a_{ i }, .. ., a_{i+d}is the empty structure containing no base pairs (due to the requirement that all hairpins contain at least θ unpaired bases). If d > θ and ${S}_{i,j}$ has been defined by recursion for all i ≤ j < i + d, then

Any secondary structure of a_{ i }, .. ., a_{i+d1}is a secondary structure for a_{ i }, .. ., a_{i+d}, in which a_{i+d}is unpaired.

If a_{ i }, a_{ j }can form a WatsonCrick or wobble base pair, then for any secondary structure S for a_{i+1}, .. ., a_{i+d1}, the structure S ∪ {(i, j)} is a secondary structure for a_{ i }, ..., a_{i+d}.

For any intermediate value i + 1 ≤ r ≤ j  θ  1, if a_{ r }, a_{ j }can form a WatsonCrick or wobble base pair, then for any secondary structure S for a_{ i }, .. .,a_{r1}and any secondary structure T for a_{r+1}, ..., a_{j1}, the structure S ∪ T ∪ {(r, j)} is a secondary structure for a_{ i }, .. ., a_{i+d}.
Given two secondary structures S, T, we define the base pair distance between S, T, denoted by d_{ BP } (S, T), to be the cardinality of the symmetric difference of S, T; i.e., d_{ BP } (S, T) = (S  T) ∪ (T  S).
RNAborSample
In this section, we describe how sampling from the Boltzmann ensemble, using Sfold[30], leads to a simple and fast approximation of RNAbor computations, provided that the input consists of an RNA sequence and the minimum free energy (MFE) secondary structure for that sequence. The work of this section is inspired by sampling approximations of the number of structural motifs, such as hairpins, multiloops, etc. of Ding and Lawrence [30], as well as the sampling approximation used in RNAshapes[8] for large sequences. A novelty of our work is that we provide a rigorous justification for the accuracy of the approximation, depending on sample size.
Let RNAborSample denote the protocol, where we apply Sfold[30] to sample N secondary structures S of an input RNA sequence a_{1}, .. .,a_{ n }, then subsequently compute the relative frequencies f_{ k } for 0 ≤ k < K, where ${f}_{k}=\frac{{N}_{k}}{N}$ is defined to be the number N_{ k } of sampled structures S, whose base pair distance with S_{0} is k, divided by N. Since Sfold appears to be only available as a web server, where the user can not stipulate the number of sampled structures, we instead use the Vienna RNA Package implementation of Sfold, given in RNAsubopt p[32].
Let a_{1}, .. .,a_{ n } be an arbitrary RNA sequence having MFE structure of S_{0}. Following [9], let Z_{ k } denote the sum of Boltzmann factors of all kneighbors of S_{0}; i.e.,
As usual, let Z denote the partition function, representing the sum of Boltzmann factors of all secondary structures of a_{1}, .. ., a_{ n }; i.e.,
and let ${p}_{k}=\frac{{Z}_{k}}{Z}$ denote the Boltzmann probability of all kneighbors.
Given a desired approximation accuracy of ε, a probability p, and an upper bound K, we wish to compute a value N = N(ε, p, K), such that whenever we sample at least N secondary structures from the Boltzmann ensemble using Sfold, the relative frequency f_{ k } of kneighbors sampled is within ε of the probability p_{ k } of kneighbors, for all 0 ≤ k < K, with confidence level of (1  p). Formally, this means that for each 0 ≤ k < K,
Consider the value k as bin number. For a fixed bin k, let us denote the exact value of $\frac{{Z}_{k}}{Z}$ by p_{ k }. If we sample N structures, each falling in the k th bin with probability p_{ k }, then the number of structures in the k th bin is given by the random variable X_{ k } having binomial distribution with mean N · p_{ k } and variance N · p_{ k }(1  p_{ k }). It follows that the proportion $\frac{{X}_{k}}{N}$ of structures in the k th bin has mean µ = p_{ k } and standard deviation $\sigma =\sqrt{\frac{{p}_{k}\left(1{p}_{k}\right)}{N}}<\frac{1}{2\sqrt{N}}$. To determine minimum sample size sufficient to ensure a certain approximation accuracy with certain confidence interval, we adapt a standard argument from statistics [37] (see equation (24.35) on p. 529), by approximating the binomial distribution by the standard normal distribution using Zscores.
Before starting, we mention that it will suffice for our intended application of RNAborSample to have a precise approximation of those p_{ k } which exceed some modest lower bound, such as δ = 0.01 or δ = 0.0001. Thus we intend to prove that for all 0 ≤ k < K, if p_{ k } ≥ δ, then Equation (4) holds.
Temporarily, we fix k. Let X be a Bernouilli random variable with success probability p_{ k }, corresponding to the indicator random variable that returns 1 if a single sampled secondary structure is a kneighbor of S_{0}. Provided that we sample a number N of structures, which satisfies N · p_{ k } ≥ 30, then the standard normal distribution can be used to approximate the left and right tail of the distribution of Zscores of sampled proportions ${f}_{k}=\frac{{\sum}_{i=1}^{N}{X}_{k}}{N}$, defined by
Let $\Phi \left(z\right)=\frac{1}{\sqrt{2\pi}}{\int}_{\infty}^{z}\text{exp}\left({x}^{2}/2\right)dx$ denote the cumulative distribution function (CDF) for the standard normal distribution. Given desired confidence interval of C = 1  α, recall that the critical value z_{α/2}is defined by
If ε is the margin of error in the left tail (∞, z_{α/2}) and right tail (z_{α/2}, +∞) for the normal approximation of the binomial distribution, then by a wellknown argument (e.g. equation (24.35) on p. 529 of [37]), we have
It follows that
provides a sufficient lower bound on number of samples necessary to guarantee margin of error ε. Let $\alpha =\frac{p}{K}$ and define
We have just shown that for N ≥ N(ε, p, K), $P\left(\leftz\right>\left{\Phi}^{1}\left(\frac{P}{2K}\right)\right\right)<\frac{p}{K}$, hence
The following is now a key step. If we have K bins, and we desire to have a small probability p that we are off by more than ε in our estimate of the probability of any bin (in our intended application, the k th bin, for 0 ≤ k < K, is the collection of kneighbors of S_{0}, i.e., those structures S, whose base pair distance with S_{0} is k) then it suffices that we have a probability $\frac{p}{K}$ that we are off by more than ε in any single bin. Indeed, let Y_{ k } denote the indicator random variable, with value 1 provided that f_{ k }  p_{ k } > ε, where f_{ k } is the relative frequency of sampling a kneighbor of S_{0}, after sampling N secondary structures, where by Equation (5), N is chosen so that
then
Putting everything together, we have shown that for given ε, p, K, we can define by defining N
we have
We have completed a more rigorous argument using Chernoff bounds, but prefer the exposition given here for simplicity. Note that the same argument, given verbatim, can be used for any binning procedure. In particular, this approach provides information on sufficient number of samples in order to approximate the result of RNAshapes [8, 38, 39] by means of sampling.
We can make some basic conclusions from the above analysis. The number of samples sufficient to ensure that f_{ k }  p_{ k } < ε for 0 ≤ k < K with confidence 1  p is reasonable, and only slightly increases for a higher number K of bins or to ensure greater confidence 1  p. However, the number of samples increases greatly when higher precision estimates (smaller ε values) are needed, even for one bin.
In the case of one bin, it is important to remember that the value N is a conservative estimate, sufficient to ensure our conclusion. This estimate uses the worst case of 50% probability of being in a bin, which maximizes the standard deviation. For bins with small probability, one can reestimate what is needed. However, for bins with smaller probability, higher precision is needed as well, unless all that is required is to verify that the bin has small probability. Also, clearly if a bin has probability of 10^{6}, then at least on the order of one million samples are needed, just for a reasonable probability of sampling the bin once.
Algorithm description
Given an RNA sequence a = a_{1}, .. ., a_{ n }, a secondary structure S_{0} of a, and a maximum desired value Kmax ≤ n, the RNAborMEA algorithm computes, for each 1 ≤ i < j ≤ n and each 0 ≤ k <Kmax ≤ n, the maximum score M(i, j, k)
where the first sum is taken over all base pairs (i, j) belonging to S, the second sum is taken over all unpaired positions in S, and where p_{ i,j } [resp. q_{ i }] is the probability that i, j are paired [resp. i is unpaired] in the ensemble of low energy structures, and α, β > 0 are weights. Our computational experiments, as in Figure 9, were carried out with default values of 1 for α, β. (See Equation 1 for the formal definition of Boltzmann base pairing probability p_{ i,j }.)
The dynamic programming computation of M(i, j, k) is performed by recursion on increasing values of j  i for all values 1 ≤ i ≤ j ≤ n and 0 ≤ k ≤ Kmax. The value of M(i, j, k), stored in the upper triangular portion of matrix M, will involve taking the maximum over three cases, which correspond to the inductive construction of all secondary structures on a_{ i }, .. ., a_{ j }, as described in the previous section. At the same time, the value M(j, i, k), stored in the lower triangular portion of matrix M, will consist of a triple r, k_{0}, k_{1} of numbers, such that the following approximately holds (in this section, we provide the motivating idea; the actual algorithm description, which deviates slightly from the description here, is given in the next section and in Figures 10 and 11). (i) If r = 0 then M(i, j, k) is maximized by a kneighbor S of S_{0}[i, j] for the subsequence a_{ i }, .. ., a_{ j } in which a_{ j } is unpaired. In this case, k_{0} = k and k_{1} = 0. (ii) If r = i, then M(i, j, k) is maximized by a kneighbor S of S_{0}[i, j] for the subsequence a_{ i }, ...,a_{ j } in which base pair (i, j) ∈ S. In this case, k_{0} = 0 and k_{1} = k  1. (i) If i < r ≤ j  θ  1 then M(i, j, k) is maximized by a kneighbor S of S_{0}[i, j] for the subsequence a_{ i }, .. .,a_{ j } in which base pair (r, j) ∈ S. The left portion of S, which is S[i, r  1] will be a k_{0} neighbor of S[i, r  1], while the right portion of S, which is S[r, j] must contain the base pair (r, j) and itself be a k_{1} neighbor of S[r, j]. In summary, the values r, k_{0}, k_{1} will be used in computing the traceback, where the maximum expected accurate structure that is a kneighbor of S[i, j] will be constructed by one of the following: (i) MEA kneighbor of S[i, j  1], in the event that a_{ j } is unpaired in [i, j]; (ii) MEA k  1neighbor of S[i + 1, j  1], in the event that a_{ i }, a_{ j } form a base pair; (iii) MEA k_{0}neighbor of S[i, r  1] and the MEA k_{1}neighbor of S[r, j], where k_{0} + k_{1} = k, in the event that a_{ r }, a_{ j } form a base pair.
Pseudocode for the algorithm RNAborMEA is given in Figures 10 and 11. An array M of size n × n × Kmax is required to store the MEA scores in M(i, j, k) for all subsequences [i, j] and all base pair distances 0 ≤ k ≤ Kmax between structures S[i, j] and initially given structure S_{0}[i, j]. For 1 ≤ i ≤ j ≤ n and all 0 ≤ k ≤ Kmax, the pseudocode in Figure 11 stores a value of the form (x, y, z) in the lower triangular portion, M(j, i, k), of the array. Here, x = 0 indicates that the optimal structure on [i, j], i.e., having maximum MEA score over all kneighbors of S_{0}[i, j], is obtained by not pairing j with any nucleotide in [i, j]; for values x > 0, hence x ∈ [i, j  θ  1], the optimal kneighbor of S_{0}[i, j] is obtained by pairing x with j. The values y, z correspond to the values k_{0}, k_{1}, such that: (i) if x = 0, then the optimal kneighbor of S_{0}[i, j] is obtained by first computing the optimal k_{0}neighbor of S_{0}[i, j  1], where k_{0} = k  b_{0}, then leaving j unpaired; (ii) if x = i, then the optimal kneighbor of S_{0}[i, j] is obtained by first computing the optimal k_{1}neighbor of S_{0}[i + 1, j  1], where k_{1} = k  b_{1}, then adding the enclosing base pair (i, j); (iii) if x = r ∈ [i + 1, j  θ  1], then the optimal kneighbor of S_{0}[i, j] is obtained by first computing the optimal k_{0}neighbor of S_{0}[i, r  1] as well as the optimal k_{1}neighbor of S_{0}[r + 1, j  1], then adding the base pair (r, j). This last calculation must be done over all values k_{0}, k_{1} such that k_{0} + k_{1} = k. Using the values M(j, i, k) = (x, y, z), the traceback can be easily computed by recursion; see Figure 12 for pseudocode of traceback.
In a manner similar to the pseudocode of Figures 10 and 11 (essentially, one replaces the operation of taking the maximum by the a summation, and one replaces the MEA score by the pseudoBoltzmann factor exp(MEA(S)/RT)), RNAborMEA also computes the pseudoBoltzmann partition function values
We have graphed the Boltzmann probabilities $\frac{{Z}_{1,n}^{\left(k\right)}}{{Z}_{1,n}}$ as well as the uniform probabilities $\frac{{N}_{1,n}^{\left(k\right)}}{{N}_{1,n}}$, where ${N}_{1,n}^{\left(k\right)}$ is the number of kneighbors, and N_{1,n}is the total number of secondary structures. When RT = n, which normalizes the MEA score to a maximum of 1, it appears that the Boltzmann distribution is the same as the uniform distribution, as shown in Figure 13.
References
 1.
Olsthoorn R, Mertens S, Brederode F, Bol J: A conformational switch at the 3' end of a plant virus RNA regulates viral replication. EMBO J 1999, 18: 4856–4864. 10.1093/emboj/18.17.4856
 2.
Repsilber D, Wiese S, Rachen M, Schroder A, Riesner D, Steger G: Formation of metastable RNA structures by sequential folding during transcription: timeresolved structural analysis of potato spindle tuber viroid ()stranded RNA by temperaturegradient gel. RNA 1999, 5: 574–584. 10.1017/S1355838299982018
 3.
Franch T, Gultyaev AP, Gerdes K: Programmed Cell Death by hok/sok of Plasmid R1: Processing at the hok mRNA 3Hend Triggers Structural Rearrangements that Allow Translation and Antisense RNA Binding. J Mol Biol 1997, 273: 38–51. 10.1006/jmbi.1997.1294
 4.
Mandal M, Boese B, Barrick J, Winkler W, Breaker R: Riboswitches control fundamental biochemical pathways in Bacillus subtilis and other bacteria. Cell 2003, 113(5):577–586. 10.1016/S00928674(03)00391X
 5.
Cheah MT, Wachter A, Sudarsan N, Breaker RR: Control of alternative RNA splicing and gene expression by eukaryotic riboswitches. Nature 2007, 447(7143):497–500. 10.1038/nature05769
 6.
Ray PS, Jia J, Yao P, Majumder M, Hatzoglou M, Fox PL: A stressresponsive RNA switch regulates VEGFA expression. Nature 2009, 457(7231):915–919. 10.1038/nature07598
 7.
Voss B, Meyer C, Giegerich R: Evaluating the predictability of conformational switching in RNA. Bioinformatics 2004, 20(10):1573–1582. 10.1093/bioinformatics/bth129
 8.
Voss B, Giegerich R, Rehmsmeier M: Complete probabilistic analysis of RNA shapes. BMC Biol 2006., 4(5):
 9.
Freyhult E, Moulton V, Clote P: Boltzmann probability of RNA structural neighbors and riboswitch detection. Bioinformatics 2007, 23(16):2054–2062. Doi: 10.1093/bioinformatics/btm314 Doi: 10.1093/bioinformatics/btm314 10.1093/bioinformatics/btm314
 10.
Barash D: Second eigenvalue of the Laplacian matrix for predicting RNA conformational switch by mutation. Bioinformatics 2004, 20(12):1861–1869. 10.1093/bioinformatics/bth157
 11.
Mandal M, Breaker RR: Adenine riboswitches and gene activation by disruption of a transcription terminator. Nat Struct Mol Biol 2004, 11: 29–35. 10.1038/nsmb710
 12.
Serganov A, Yuan Y, Pikovskaya O, Polonskaia A, Malinina L, Phan A, Hobartner C, Micura R, Breaker R, Patel D: Structural Basis for Discriminative Regulation of Gene Expression by Adenine and GuanineSensing mRNAs. Chem Biol 2004, 11(12):1729–1741. 10.1016/j.chembiol.2004.11.018
 13.
Serganov A, Polonskaia A, Phan AT, Breaker RR, Patel DJ: Structural basis for gene regulation by a thiamine pyrophosphatesensing riboswitch. Nature 2006, 441(7097):1167–1171. 10.1038/nature04740
 14.
AbreuGoodger C, Merino E: RibEx: A web server for locating riboswitches and other conserved bacterial regulatory elements. Nucleic Acids Res 2005, 33: W690W692. 10.1093/nar/gki445
 15.
Bengert P, Dandekar T: Riboswitch finder  A tool for identification of riboswitch RNAs. Nucleic Acids Res 2004, 32: W154W159. 10.1093/nar/gkh352
 16.
Chang TH, Huang HD, Wu LC, Yeh CT, Liu BJ, Horng JT: Computational identification of riboswitches based on RNA conserved functional sequences and conformations. RNA 2009., 15(7):
 17.
Nawrocki EP, Kolbe DL, Eddy SR: Infernal 1.0: inference of RNA alignments. Bioinformatics 2009, 25(10):1335–1337. 10.1093/bioinformatics/btp157
 18.
Weinberg Z, Barrick JE, Yao Z, Roth A, Kim JN, Gore J, Wang JX, Lee ER, Block KF, Sudarsan N, Neph S, Tompa M, Ruzzo WL, Breaker RR: Identification of 22 candidate structured RNAs in bacteria using the CMfinder comparative genomics pipeline. Nucleic Acids Res 2007, 35(14):4809–4819. 10.1093/nar/gkm487
 19.
Gardner PP, Daub J, Tate JG, Nawrocki EP, Kolbe DL, Lindgreen S, Wilkinson AC, Finn RD, GriffithsJones S, Eddy SR, Bateman A: Rfam: updates to the RNA families database. Nucleic Acids Res 2009, 37(Database):D136D140. 10.1093/nar/gkn766
 20.
Freyhult E, Moulton V, Clote P: Boltzmann probability of RNA structural neighbors and riboswitch detection. Bioinformatics 2007, 23(16):2054–2062. 10.1093/bioinformatics/btm314
 21.
Freyhult E, Moulton V, Clote P: RNAbor: a web server for RNA structural neighbors. Nucleic Acids Res 2007, 35(Web):W305W309. 10.1093/nar/gkm255
 22.
Matthews D, Sabina J, Zuker M, Turner D: Expanded sequence dependence of thermodynamic parameters improves prediction of RNA secondary structure. J Mol Biol 1999, 288: 911–940. 10.1006/jmbi.1999.2700
 23.
Xia T, SantaLucia JJ, Burkard M, Kierzek R, Schroeder S, Jiao X, Cox C, Turner D: Thermodynamic parameters for an expanded nearestneighbor model for formation of RNA duplexes with WatsonCrick base pairs. Biochemistry 1999, 37: 14719–35.
 24.
Do CB, Mahabhashyam MS, Brudno M, Batzoglou S: ProbCons: Probabilistic consistencybased multiple sequence alignment. Genome Res 2005, 15(2):330–340. 10.1101/gr.2821705
 25.
Nussinov R, Jacobson AB: Fast Algorithm for Predicting the Secondary Structure of Single Stranded RNA. Proceedings of the National Academy of Sciences, USA 1980, 77(11):6309–6313. 10.1073/pnas.77.11.6309
 26.
Kiryu H, Kin T, Asai K: Robust prediction of consensus secondary structures using averaged base pairing probability matrices. Bioinformatics 2007, 23(4):434–441. 10.1093/bioinformatics/btl636
 27.
McCaskill J: The equilibrium partition function and base pair binding probabilities for RNA secondary structure. Biopolymers 1990, 29: 1105–1119. 10.1002/bip.360290621
 28.
Lu ZJ, Gloor JW, Mathews DH: Improved RNA secondary structure prediction by maximizing expected pair accuracy. RNA 2009, 15(10):1805–1813. 10.1261/rna.1643609
 29.
Zuker M: On finding all suboptimal foldings of an RNA molecule. Science 1989, 244(7):48–52.
 30.
Ding Y, Lawrence CE: A statistical sampling algorithm for RNA secondary structure prediction. Nucleic Acids Res 2003, 31: 7280–7301. 10.1093/nar/gkg938
 31.
Gardner PP, Daub J, Tate J, Moore BL, Osuch IH, GriffithsJones S, Finn RD, Nawrocki EP, Kolbe DL, Eddy SR, Bateman A: Rfam: Wikipedia, clans and the "decimal" release. Nucleic Acids Res 2011, 39(Database):D141D145. 10.1093/nar/gkq1129
 32.
Gruber A, Lorenz R, Bernhart S, Neubock R, Hofacker I: The Vienna RNA websuite. Nucleic Acids Research 2008, 36: 70–74.
 33.
Mathews DH, Disney MD, Childs JL, Schroeder SJ, Zuker M, Turner DH: Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of RNA secondary structure. Proc Natl Acad Sci USA 2004, 101(19):7287–7292. 10.1073/pnas.0401799101
 34.
Lorenz W, Clote P: Computing the partition function for kinetically trapped RNA secondary structures. Public Library of Science One (PLoS ONE) 2011, 6: 316178. Doi:10.1371/journal.pone.0016178 Doi:10.1371/journal.pone.0016178
 35.
Wakeman CA, Winkler WC, Dann C: Structural features of metabolitesensing riboswitches. Trends Biochem Sci 2007, 32(9):415–424. 10.1016/j.tibs.2007.08.005
 36.
Blin G, Denise A, Dulucq S, Herrbach C, Touz H: Alignments of RNA structures. IEEE/ACM Transactions on Computational Biology and Bioinformatics 2010.
 37.
Zar J: Biostatistical Analysis. PrenticeHall, Inc; 1999.
 38.
Giegerich R, Voss B, Rehmsmeier M: Abstract shapes of RNA. Nucleic Acids Res 2004, 32(16):4843–4851. 10.1093/nar/gkh779
 39.
Steffen P, Voss B, Rehmsmeier M, Reeder J, Giegerich R: RNAshapes: an integrated RNA analysis package based on abstract shapes. Bioinformatics 2006, 22(4):500–503. 10.1093/bioinformatics/btk010
 40.
Hofacker I: Vienna RNA secondary structure server. Nucleic Acids Res 2003, 31: 3429–3431. 10.1093/nar/gkg599
 41.
Markham NR, Zuker M: UNAFold: software for nucleic acid folding and hybridization. Methods Mol Biol 2008, 453: 3–31. 10.1007/9781603274296_1
 42.
Ponty Y: Efficient sampling of RNA secondary structures from the Boltzmann ensemble of lowenergy: The boustrophedon method. J Math Biol 2008, 56(1–2):107–127.
Acknowledgements
Funding for the research of P. Clote and F. Lou was provided by the Digiteo Foundation for the project RNAomics. Additional funding was provided to P. Clote by National Science Foundation grants DMS1016618 and DMS0817971. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.
This article has been published as part of BMC Bioinformatics Volume 13 Supplement 5, 2012: Selected articles from the First IEEE International Conference on Computational Advances in Bio and medical Sciences (ICCABS 2011): Bioinformatics. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/13/S5.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
RNAborSample was conceived by W.A.L., who provided a proof for sample size sufficient to ensure a desired degree of accuracy with a desired level of confidence. RNAborMEA was conceived by P.C., then designed and programmed by P.C. (prototype implementation in Python) and F.L. (implementation in C). Computational experiments were carried out by F.L., figures were created by F.L. and P.C. and the manuscript was written by P.C.
Rights and permissions
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.
About this article
Cite this article
Clote, P., Lou, F. & Lorenz, W.A. Maximum expected accuracy structural neighbors of an RNA secondary structure. BMC Bioinformatics 13, S6 (2012). https://doi.org/10.1186/1471210513S5S6
Published:
Keywords
 Secondary Structure
 Minimum Free Energy
 Expression Platform
 Conformational Switch
 Boltzmann Probability