Maximum expected accuracy structural neighbors of an RNA secondary structure

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 15-20 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 S0 for an RNA nucleotide sequence a = a1,..., an, we say that another secondary structure S of a is a k-neighbor of S0, if the base pair distance between S0 and S is k. In this paper, we prove that the Boltzmann probability of all k-neighbors of the minimum free energy structure S0 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(ε,p,K)=Φ-1p2K24ε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 S0 and for all values 0 ≤ k < K, computes the secondary structure MEA(k), having maximum expected accuracy over all k-neighbors of S0. Computation time is O(n3 · K2), and memory requirements are O(n2 · 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 k-neighbors 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 stress-responsive 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, thermodynamics-based 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 Shine-Dalgarno 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][15][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, RNAbor-Sample 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 k-neighbor 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 k-neighbors 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 = Z(k) 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(i) = 1 − i<j p(i, j) − j<i p(j, i) 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 (i,j)∈S 2α · p(i, j) + iunpaired β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 a, b >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 RNAbor-Sample 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 k-neighbor 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, RNAbor-Sample, 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.  ((((....)))))))) with free energy -10.3 kcal/mol. The 16-neighbor of S 0 is the metastable structure ((((((((....))))))))....... with free energy -9.9 kcal/mol. The MFE structure appears above the leftmost peak, while the MFE(16) structure appears above the rightmost peak. The output of RNAbor includes a graph of the Boltzmann probabilities p k = Z k Z , and MFE(k) structures, for all 0 ≤ k ≤ 2n. The existence of distinct 'peaks' suggests the presence of a conformational switch.

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 RNAbor-Sample that approximates the output p k = Z k Z of RNAbor by frequency countsp 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 < a <1, if at least structures are sampled, then for all 0 ≤ k < K; i.e., RNAbor-Sample furnishes estimatesp k of p k , for all 0 ≤ k < K, which with confidence 1 -a 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 k-neighbors of a given RNA secondary structure S 0 ; i.e., for each 0 ≤ k ≤ K, RNA-borMEA 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 Moreover, RNAborMEA allows the user to stipulate (partial) hard constraints, that stipulate whether particular nucleotides are unpaired, or base-pair with certain other nucleotides. The implementation of hard constraints follows ideas from Mathews [33], albeit suitably modified to simultaneously consider all k-neighbors, for 0 ≤ k ≤ K.
We now describe the 13 figures and 4 tables, corresponding to computational experiments performed with RNAbor-Sample 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 estimatesp k from RNAbor-Sample for the SAM riboswitch aptamer with GenBank accession code AP004597.1/11894-11904. Clearly, probability estimatesp k are close to actual values p k . The figure additionally shows probabilities r k from our software RNAlocopt [34], computed by r k = Z k (LO) Z(LO) , 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 k-neighbors 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, RNAbor-Sample, 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 RNAbor-Sample 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 RNA-borMEA 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, RNAbor-Sample, RNAlocopt, RNAshapes, UNAFold). Figure 8 displays the output of RNAborMEA, given the sequence of a TPP riboswitch with EMBL accession code AF269819/1811-1669. 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 k-neighbors 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 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 pseudo-Boltzmann probabilitiesp k = Z k Z for two small RNA sequences. While temperature T has a natural significance, when computing Boltzmann probabilities p k = 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 -a 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 k-neighbor, determined by RNAborMEA, RNAbor-Sample and RNAlocopt. Table 4 presents the number of times that each of the methods RNAborMEA, RNAbor, RNAbor-Sample, RNAlocopt, RNAshapes, UNAFold output the most similar structure to the GENE ON resp. GENE OFF structure for the the initial structure (taken as the minimum free energy here). The script RNAbor-Sample calls Sfold on 1000 structures, in order to compute a relative frequence f k ≈ p k of all k-neighbors of S 0 . Finally, we compute relative frequency of RNAlocopt [34], a program that samples only locally optimal secondary structures, having the property that one cannot obtain a lower energy structure by adding or removing a single base pair.
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 RNA-borMEA 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  [35]. The free energy of GENE ON resp. GENE OFF secondary structrure is -16.46 kcal/mol resp. -22.6 kcal/mol. Free energies were determined using RNAeval and figures produced using RNAplot, both from the Vienna RNA Package [40].

Figure 4
Given riboswitch sequence X83878/168-267 and initial structure S 0 , the minimum free energy structure, a structure output by RNAborMEAis most structurally similar to the XPT GENE ONstructure, as measured by NestedAlign [36]. The NestedAlign score for RNAborMEA is 87.5, while optimal score for RNAbor is 60.0, and for RNAshapes is 64.0.
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 temperature-dependent partition function, defined by Here, the expected accuracy score s 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 Ding-Lawrence 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 Figure 5 For each RNA sequence in the seed alignment from Rfam family RF00167 of purine riboswitch aptamers, we retrieved downstream flanking residues from the appropriate EMBL files, in order to ensure likelihood that the expression platform was included. Then the following six programs were run: RNAbor, RNAborMEA, RNAbor-Sample, RNAlocopt, RNAshapes, UNAFold. Each program outputs a number of near-optimal secondary structures, each according to different criteria. Taking RNAbor and RNAborMEA as examples, the programs RNAbor and RNAborMEA were run, in order to compute the MFE(k) structure and the MEA(k) structure, which have minimum free energy resp. maximum expected accuracy among all k-neighbors of the intial minimum free energy structures S 0 . Subsequently, we applied the program NestedAlign described in [36] to compute the structural similarity between the experimentally determined GENE ON structure for XPT guanine riboswitch of B. subtilis; i.e. the left panel of Figure 3. (Similar structures have positive scores; dissimilar structures have negative scores.) For each RNA in the seed alignment of RF00167, we determined the value k 1 , such the MEA(k 1 ) structure for that RNA has the greatest structural similarity with the XPT GENE ON structure, as determined by NestedAlign. (See the left panel of Figure 3 for the experimentally determined GENE ON structure of XPT.) As earlier explained, we performed similar computations for the programs RNAshapes [39] and UNAFold [41], the programs RNAborMEA and RNAbor-Sample, described in this paper, and programs RNAbor [9] and RNAlocopt [34], developed by our lab. In 21 out of 34 instances, RNAborMEA produced the secondary structure most structurally similar to the experimentally determined XPT GENE OFF structure, as measured by NestedAlign.
of the temperature-dependent mixing of various structures having large score s. On the other hand, in computational experiments reported in the Results Section, it appears that RNAborMEA produces near-optimal 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, RNAbor-MEA produced the secondary structure most structurally similar to the experimentally determined XPT GENE ON [resp. GENE OFF] structure, as measured by Nes-tedAlign [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 state-of-the-art 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/.

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. Watson-Crick 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. Figure 6 For each RNA sequence in the seed alignment from Rfam family RF00167 of purine riboswitch aptamers, we retrieved downstream flanking residues from the appropriate EMBL files, in order to ensure likelihood that the expression platform was included. Then the following six programs were run: RNAbor, RNAborMEA, RNAbor-Sample, RNAlocopt, RNAshapes, UNAFold. Each program outputs a number of near-optimal secondary structures, each according to different criteria. Taking RNAbor and RNAborMEA as examples, the programs RNAbor and RNAborMEA were run, in order to compute the MFE(k) structure and the MEA(k) structure, which have minimum free energy resp. maximum expected accuracy among all k-neighbors of the intial minimum free energy structures S 0 . Subsequently, we applied the program NestedAlign described in [36] to compute the structural similarity between the experimentally determined GENE OFF structure for XPT guanine riboswitch of B. subtilis; i.e. the right panel of Figure 3. (Similar structures have positive scores; dissimilar structures have negative scores.) For each RNA of the seed alignment of RF00167, we determined the value k 1 , such the MEA(k 1 ) structure for that RNA has the greatest structural similarity with the XPT GENE OFF structure, as determined by NestedAlign. (See the right panel of Figure 3 for the experimentally determined GENE OFF structure of XPT.) As earlier explained, we performed similar computations for the programs RNAshapes [39] and UNAFold [41], the programs RNAborMEA and RNAbor-Sample, described in this paper, and RNAbor [9] and RNAlocopt [34]. In 22 out of 34 instances, RNAborMEA produced the secondary structure most structurally similar to the experimentally determined XPT GENE OFF structure, as measured by NestedAlign.
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, .. ., nd, 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 = {∅} ; 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 < For each RNA sequence in the seed alignment from Rfam family RF00066 of U7 small nuclear RNAs, both RNAbor and RNAborMEA were run, in order to compute the MFE(k) structure and the MEA(k) structure, which have minimum free energy resp. maximum expected accuracy among all k-neighbors of the intial minimum free energy structures S 0 . We computed the base pair distance between the MFE(k) structure and the MEA(k) structure over all sequences in the seed alignment of RF00066. The figure displays the average ± one standard deviation of base pair distance. Figure 8 Sample outputs from RNAborMEAon a 143 nt TPP-riboswitch, AF269819/1811-1669 with sequence CUACUAGGGG AGCCAAAAGG CUGAGAUGAA UGUAUUCAGA CCCUUAUAAC CUGAUUUGGU UAAUACCAAC GUAGGAAAGU AGUUAUUAAC UAUUCGUCAU UGAGAUGUCU UGGUCUAACU ACUUUCUUCG CUGGGAAGUA GUU. We took the TPP riboswitch aptamer from the Rfam database [19], then extracted right-flanking nucleotides from the corresponding EMBL file, in order to include the expression platform. Displayed from left to right are the structures MEA(0) and MEA(61) (the structure MEA(52) is similar to that of MEA(61) and corresponds to a free energy local minimum in the left figure.) The structure MEA(61) had the highest MEA score over all structural neighbors, including the original structure S 0 = MEA(0), and had free energy, -46.0 kcal/mol, that was equal to that of the initial structure S 0 = MEA(0), which is the minimum free energy structure for the given sequence.

Figure 12
Pseudocode for the O(n 2 ) traceback computed by our RNAborMEAalgorithm. Note that run time could be reduced to O(n ln n) by applying the boustrephedonic method described in [42].
i + d, then • Any secondary structure of a i , .. ., a i+d-1 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 Watson-Crick or wobble base pair, then for any secondary structure S for a i+1 , .. ., a i+d-1 , 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 Watson-Crick or wobble base pair, then for any secondary structure S for a i , .. .,a r-1 and any secondary structure T for a r+1 , ..., a j-1 , 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)|.

RNAbor-Sample
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 RNAbor-Sample 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 = N k N is defined to be the number N k of sampled structures S, whose base pair distance with S 0 is k, ).))).))))))). Here, the (formal) parameter RT taken to be 49 (length of sequence), in order to uniformize MEA scores to range  ((((....)))))))). The left curve is when RT = 0.6, the approximate value obtained by multiplying the universal gas constant 0.00198 kcal/mol times 310 Kelvin. In contrast, the right curve is when RT = 27 (length of sequence). Though not shown in this graph, the pseudo-Boltzmann distribution is identical with the uniform distribution, when RT = n, where n is sequence length. of samples sufficient to guarantee that |f k -p k | < ε with confidence 1 -p, for all 0 ≤ k < K, in the application RNAbor-Sample. Here p k = Z k Z , the Boltzmann probability, as computed exactly by RNAbor, for a k-neighbor of S 0 , and f k is the relative frequency of k-neighbors among N sampled structures. 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 k-neighbors 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., 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  NestedAlign similarity scores between the GENE ON structure of the XPT guanine riboswitch of B. subtilis, experimentally determined using in-line probing (see [35]), and the structurally most similar secondary structure among near-optimal structures generated by each of the following six methods: RNAbor, RNAborMEA, RNAbor-Sample, RNAlocopt, RNAshapes, UNAFold. These values are plotted in Figure 5, where more details on the computational experiment are given.  NestedAlign similarity scores between the GENE OFF structure of the XPT guanine riboswitch of B. subtilis, experimentally determined using in-line probing (see [35]), and the structurally most similar secondary structure among near-optimal structures generated by each of the following six methods: RNAbor, RNAborMEA, RNAbor-Sample, RNAlocopt, RNAshapes, UNAFold. These values are plotted in Figure 6, where more details on the computational experiment are given. Number of times that the most similar structure to the GENE ON resp. GENE OFF structure of the B. subtilis XPT riboswitch was produced by each of the six methods investigated. Although the test was made with 34 sequences from the seed alignment of Rfam family RF00167 [31], the sums of each column may exceed 34; this is because If two or more methods produced the maximum score, then each was counted. Structural similarity was measured using the NestedAlign structural alignment algorithm [36]. While the GENE OFF structure involves a terminator loop, that is often correctly found by thermodynamicsbased software, the GENE ON secondary structure, having higher free energy (hence less stable thermodynamically) is less likely to be found using thermodynamics-based approaches.
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 Z k Z by p k . If we sample N structures, each falling in the kth bin with probability p k , then the number of structures in the kth bin is given by the random variable X k having binomial distribution with mean N · p k and variance N · p k (1p k ). It follows that the proportion X k N of structures in the kth bin has mean µ = p k and standard deviation 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 Z-scores.
Before starting, we mention that it will suffice for our intended application of RNAbor-Sample 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 k-neighbor 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 Z-scores of sampled proportions z −∞ exp(−x 2 /2)dx denote the cumulative distribution function (CDF) for the standard normal distribution. Given desired confidence interval of C = 1 -a, recall that the critical value z a/2 is defined by If ε is the margin of error in the left tail (-∞, -z a/2 ) and right tail (z a/2 , +∞) for the normal approximation of the binomial distribution, then by a well-known 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 α = p K and define We have just shown that for N ≥ N(ε, p, K), 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 kth bin, for 0 ≤ k < K, is the collection of k-neighbors 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 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 kp k | > ε, where f k is the relative frequency of sampling a k-neighbor of S 0 , after sampling N secondary structures, where by Equation (5), N is chosen so that 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 kp k | < ε for 0 ≤ k < K with confidence 1p is reasonable, and only slightly increases for a higher number K of bins or to ensure greater confidence 1p. 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 re-estimate 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 a, b >0 are weights. Our computational experiments, as in Figure 9, were carried out with default values of 1 for a, b. (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 ji 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 k-neighbor 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 k-neighbor 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 k-neighbor 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 k-neighbor of S[i, j] will be constructed by one of the following: (i) MEA k-neighbor of S[i, j -1], in the event that a j is unpaired in [i, j]; (ii) MEA k -1-neighbor 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 k-neighbors 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 k-neighbor 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 k-neighbor of S 0 [i, j] is obtained by first computing the optimal k 0 -neighbor of S 0 [i, j -1], where k 0 = kb 0 , then leaving j unpaired; (ii) if x = i, then the optimal k-neighbor 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 = kb 1 , then adding the enclosing base pair (i, j); (iii) if x = r [i + 1, j -θ -1], then the optimal k-neighbor 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 pseudo-Boltzmann factor exp(MEA (S)/RT)), RNAborMEA also computes the pseudo-Boltzmann partition function values 1,n is the number of k-neighbors, 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.