Coarse-Grained Prediction of RNA Loop Structures

One of the key issues in the theoretical prediction of RNA folding is the prediction of loop structure from the sequence. RNA loop free energies are dependent on the loop sequence content. However, most current models account only for the loop length-dependence. The previously developed “Vfold” model (a coarse-grained RNA folding model) provides an effective method to generate the complete ensemble of coarse-grained RNA loop and junction conformations. However, due to the lack of sequence-dependent scoring parameters, the method is unable to identify the native and near-native structures from the sequence. In this study, using a previously developed iterative method for extracting the knowledge-based potential parameters from the known structures, we derive a set of dinucleotide-based statistical potentials for RNA loops and junctions. A unique advantage of the approach is its ability to go beyond the the (known) native structures by accounting for the full free energy landscape, including all the nonnative folds. The benchmark tests indicate that for given loop/junction sequences, the statistical potentials enable successful predictions for the coarse-grained 3D structures from the complete conformational ensemble generated by the Vfold model. The predicted coarse-grained structures can provide useful initial folds for further detailed structural refinement.


Introduction
The ability to predict RNA 3D structure is critical for understanding RNA functions. Recent developments in de novo prediction of RNA 3D structures have led to highly promising results [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18] (for review, see [19][20][21][22][23][24]). In particular, several de novo structure prediction methods have been developed based on the knowledge-based energy function. For example, Dima and coworkers [25] extracted the base-pair stacking parameters from RNA native structures and the extracted parameters agree with the experimental data [26]. Wu et al. [27] explored the correlation between RNA secondary structural motifs and their thermodynamic stability to derive energy parameters for base-pair stackings, and free segments such as hairpin loops, internal loops and bulge loops. Bernauer and co-workers [15] further extracted a set of distance-dependent energy parameters between any two bases, irrespective of the locations of the bases (in a helix or a loop). Das and Baker [5,9] obtained energy function -made available in the Rosetta software package -based on the base orientations and interactions. Parisien and Major [7] predicted 3D structures with a pipeline of two computer programs: MC-Fold and MC-Sym, developed based on the nucleotide cyclic motifs (NCMs). All of these methods have provided valuable insights into the correlation between loop sequence and their stability. Especially, these methods are particularly useful for selecting the most probable conformation from an ensemble of near-native structures.
A key issue in the predictions of RNA stability is how to compute the loop free energy. For hairpins and RNA secondary structures in general, the nearest neighbor model, which assumes that the total free energy is an additive sum of the free energy of each elements (base-pair stacking, loop), and other models [6,26,[28][29][30][31][32] has enabled successful predictions for RNA structures and stabilities. In most of the existing models, loop stability is often assumed to depend on loop size, the identity of the closing base pair, the interaction of the first mismatch with the closing base pair, and an additional stabilization term for loops with GA or UU first mismatches [33][34][35]. Further detailed sequencedependence of the loop stability has been ignored. Experimental results suggest that loop stability may be sensitive to the sequence context inside the loop. For example, for unusually stable RNA hairpin loops, Dale et al. [36] performed optical melting studies for a series of hairpins. The study led to a set of different stability parameters for different loop sequences, such as GNRA and UUCG (where N is any nucleotides and R is a purine) tetraloops, hexaloops with UU first mismatches, and hairpin loop of iron responsive element, GAGUGC, all of which are significantly more stable than other hairpin loops of the same length.
The prediction of sequence-dependent loop free energy requires a model that goes beyond simple fitting with the experimentally measured empirical parameters. The formation of the intraloop base pairs and stacks would cause significant restriction of the loop conformational space and the loop entropy. It is practically impossible to exhaustively measure the loop free energy for all the different possible intraloop contacts for the different sequences and loop lengths. Therefore, evaluation of loop free energy with an ensemble of possible intraloop base pairing/stacking interactions cannot be achieved by experiment alone. We also need a computational model. The main purpose of the present study is to develop a statistical potential model that enables predictions of loop and junction three-dimensional structures from the sequence.
In general, there are two classes of physics-based models for RNA structure prediction: molecular dynamics (For review, see [37]) and Monte Carlo simulation methods and polymer statistical mechanics methods. Molecular dynamics simulations have provided much insights into the atomic details of intraloop interactions and their contributions to the loop stabilities [1,[38][39][40]. The polymer statistical mechanical models often employ lowresolution (coarse-grained) conformational models in order to capture the complete conformational ensemble. Along this line there have been different ways to construct the low-resolution RNA structures. For example, with a knowledge-based potential, Jonikas et al. [8] developed a structure filter model (N AST) where nucleotides are represented by the C Ã 3 atoms. In another coarse-grained model where nucleotide are represented by the P and C Ã 3 atoms, Keating and Pyle [41] developed a semi-automated approach to build RNA structures with a directed rotameric search strategy. Furthermore, in an attempt to develop a highresolution RNA model (HiRE-RNA), Pasquali and Derreumaux [11] used six to seven beads for each nucleotide (one bead for the phosphate P, four beads for the sugar O Ã 5 , C Ã 5 , C Ã 4 , C 1 , respectively, and one bead for a pyrimidine base and two beads for a purine base).
Our RNA folding model is based on a virtual bond-based RNA conformational model (called ''Vfold'' model; Fig. 5) [42]. The Vfold model uses two virtual bonds (C Ã 4 {P{C Ã 4 ) for each nucleotide and samples RNA conformations through self-avoiding walks in a diamond lattice. It provides an effective tool to sample RNA conformations and to evaluate the conformational entropy. The model has shown a high promise in predicting the 2D and 3D structures and the folding stabilities from the sequence [14,[42][43][44][45][46][47][48]. However, the Vfold model does not account for the sequencedependent conformational propensity of the loop, namely, the model assumes that for any given loop, all the loop conformations generated in the model have the same energy. Such a simplification could cause inaccuracy in the prediction of loop stability and structure [36,49,50]. Physically, the sequence-dependence of loop stability arises from the local interactions, which affect the loop flexibility and hence conformational propensity [50], as well as the nonlocal interactions between the different nucleotides. In this study, we develop a method to extract a set of virtual bondbased (coarse grained) statistical potentials (scoring functions) from a set of non-redundant RNA structures such as the RNA09 database [51] and the Leontis database (http://rna.bgsu.edu/ nrlist/). Specifically, we aim to derive a set of dinucleotide statistical potentials u ij (h, g) as a function of the (4|4) types of the dinucleotides base i and j and the backbone pseudo-torsion angles (h, g) of the dinucleotide conformation. We use the RNA09 database [51], the Leontis database, the Capriotti's database [52] and the PDB database [53,54], respectively, to test the extracted potential functions. We note that the dinucleotide in this work refers to the two continuous nucleotides within the same loop. The goal is for a given RNA loop/junction sequence, to identify the lowest-RMSD structures from the Vfold-generated complete conformational ensemble.

Pseudo-torsion angles (h, g)
The h and g values of dinucleotides in the 152 RNA loops/ junctions are plotted as a 2D scatter plot (Fig. 6), where each point represents the h-g coordinates for a dinucleotide. In contrast to the analysis in the previous studies [55,56], here we find a distinct category of dinucleotides conformation around (h = 150 0 , g = 225 0 ). This region was once considered as the helical region, because the dinucleotides with h-g coordinates located in this group are most likely found within the RNA helix. However, as we only count the dinucleotides in RNA loops and junctions, the plot shows that the loop/junction residues can also have the tendency to have the helix-like conformation. This observation provides a rational in the next step for building the 3D all-atom structures by adding the helical residues back to the coarse-grained backbone model.   Quasi-chemical approximation-based potentials The 152 RNA loops/junctions within our training dataset consist of a total number of N total = 572 nucleotides, with the numbers of each type of the nucleotide (r A ,r C ,r G ,r U )( 217,105,121,129) and the corresponding mole fraction of each nucleotide (x A ,x C ,x G ,x U )~(0:38,0:18,0:21,0:23).
From the dataset of the coarse-grained ''correct'' structures (Table 1), for each pair of the nucleotides i~(1,2,3,4) and j~(1,2,3,4) and pseudo-torsion angles (h, g), we calculate the numbers r obs ij (h, g) of the dinucleotides conformations and the total observed number N(h, g) = P i,j r obs ij (h, g). According to Eq. 8, we then calculate the expected number r exp ij (h, g) of each type of dinucleotides (i,j) with pseudo-torsion angles (h, g). If no pseudo-torsion angles (h,g) is observed for dinucleotide (i,j), we assign an unfavorable potential u ij (h, g)~z2:0 k B T kcal/mol. From Eq. 7, we compute the potentials as a 4|4|3|3 tensor.

Iterative approach-derived potentials
The convergence speed of the iteration depends on the selection of convergence criteria. For instance, if the convergence threshold parameter in Eq. 13 is set to 10 23 , the iterative process would converge after around 3000 iterations. In contrast, if the convergence threshold parameter is set to 10 22 , the iterative process would converge after 340 iterations. On an Intel(R) Xeon(R) CPU 5150 @ 2.66 GHz on Dell EM64T cluster system, the 3000-cycle iteration process took about 30 minutes and a 300step iteration took less than 6 minutes.
Comparison between the two sets of the derived potential parameters shows that the ''true'' potentials SP fold are less uniform than the ''extracted'' potentials ( Fig. 2), suggesting that the ''true'' potentials SP fold are more sensitive to the sequences and conformations of the dinucleotide and thus have the ability to discriminate the ''correct'' conformation from an ensemble of conformations for a given RNA loop/junction sequence.

Tests on training datasets
We test the statistical potentials for the accuracy in loop/ junction structure prediction for a large number of sequences. For a given sequence, we generate the full ensemble of the virtual bond loop/junction conformations using our Vfold model. Each conformation is then scored by the total statistical potential, which is evaluated as the sum of the statistical potentials (SP native or SP fold ) for the dinucleotide pairs in the conformation. The conformation with the lowest value of the total statistical potential is identified as the predicted native structure. If the predicted native structure is the same as the coarse-grained ''correct'' structure (i.e., the virtual bond structure that is closest to the PDB Table 4. The numbers of sequences of the successfully predicted loop/junction structures for (I) all 8452 RNA loops/junctions in TEST-I (II) the 7459 RNA loops and junctions in TEST-II and (III) the 1119 RNA loops and junctions in TEST-III.   For each dataset, the numbers in columns ''A'' are calculated from SP fold (''# SPfold '', ''true'' potential parameters) and SP native (''# SPnative '', ''extracted'' potential parameters), obtained from the RNA09 database, and the numbers in columns ''B'' are calculated from SP fold (''# SPfold '') and SP native (''# SPnative ''), obtained from the Leontis' database respectively. The ''Top-#'' in the first column means that the ''correct'' structure is in the top # lowest-potential conformations. doi:10.1371/journal.pone.0048460.t004 structure, evaluated by RMSD), the prediction is successful for the sequence; otherwise, the prediction fails. We first apply the two sets of statistical potentials to the RNA loops and junctions in the training dataset. Such a test for structure prediction is nontrivial because the SP native and SP fold are derived based on the frequency r ij (h, g) (see Eqs. 7 and 10) instead of the structure. As described above, 152 loops and junctions are constructed from the 262 RNA structures in the RNA09 dataset, including 72 RNA hairpin loops, 25 internal/bulge loops, 14 pseudoknot loops, 35 multibranched loops and 6 junctions (the free segments other than the above 4 types) (Table 3A). Here we note that the two loops in an internal loops are counted separately in our calculation, as our model does not specify specific types of loops or junctions. This rule is also applied to the pseudoknot loops and multibranched loops.
We found that SP native succeeded in finding 130 coarse-grained ''correct'' (lowest potential) structures out of 152 loops and junctions. In contrast, SP fold can give successful predictions for all the 152 the coarse-grained ''correct'' conformations (Table 3A). SP fold is more reliable than SP native in structure prediction with success rate 100% vs 85.5%.

Test on the Capriotti's dataset
To rigorously test the reliability of the two sets of potential parameters, we need to perform the test on loops and junctions outside the training dataset. We will perform tests for several such test sets. We first choose a dataset collected by Capriotti, et al. [52], consisting of 85 structures with length §20 nucleotides and solved at resolution better than 3.5 Å . The 3DNA software determined a total of 72 loops/junctions with lengths ranging from 3 nt to 8 nt, excluding the 59/39-terminal dangling regions. The 72 loops/junctions can be classified into 37 hairpin loops, 14 internal/bulge loops, 5 pseudoknot loops, 1 junction and 15 multibranched loops (Table 3B).
Our results indicate that SP native and SP fold potentials can successfully find out 62 and 70 coarse-grained ''correct'' confor-mations, respectively, out of the 72 test cases (Table 3B). The result indicates that the ''true'' potential parameters SP fold are more reliable than the directly extracted potentials SP native (success rate 97% vs 86%). One of the two failed predictions for SP native is for a pseudoknot loop, which has special loop structure due to the tertiary interactions (base triplets) with the helices. Another failed prediction is a junction located in a large RNA molecule and the junction structure is determined by not only its sequence content but also the surrounding structural environment.

Test on the PDB dataset
The January 2012 version of PDB database contains 2227 structures that contain at least one strand of RNA sequence. These structures range from hairpin-loop structures to RNA-protein complexes or RNA-DNA hybrids. We found 8452 loops and junctions in the 2227 PDB structures (TEST-I). All the loops and junctions (excluding the 39/59-terminal dangling regions) have lengths from 3 nt to 8 nt. Within these 2227 RNA molecules, 1609 have structures determined by X-ray crystallography, which contain 7459 RNA loops and junctions (TEST-II), and 934 of these 1609 RNAs have high-resolution structures (ƒ3.0 Å ), which contain 1119 RNA loops and junctions (TEST-III).
Our test results show that the numbers of the correct predictions with the statistical potentials SP fold and SP native are 7364 (success rate = 87%) v.s. 6992 (success rate = 83%) for TEST-I, 6553 (success rate = 88%) v.s. 6222 (success rate = 83%) for TEST-II, and 1070 (success rate = 95%) v.s. 1001 (success rate = 89%) for TEST-III, respectively. Fig. 10 shows illustrations for two of the results (a hairpin loop from PDB structure 1IVS and an internal loop from PDB structure 1JJ2).
If we include the top five lowest-potential conformations, the numbers of successfully determined loops and junctions are increased to 8006 with SP fold v.s. 7857 with SP native for TEST-I, 7089 v.s. 6949 for TEST-II and 1102 v.s. 1065 for TEST-III (Table 2). Fig. 11 shows the minimal-RMSD structure, the 9-th potential structure computed with SP native , and the predicted structure with the lowest-potential, for multibranched loop from PDB structure 3CCM.
Moreover, we test the predictions of the first 20 lowest-potential conformations for the loops and junctions in all the three databases TEST-I, TEST-II and TEST-III as well as the influence of the convergence threshold parameter on the accuracy of the structure prediction. Fig. 12 shows the numbers of successfully determined loops and junctions in the first 20 lowest-potential conformations for the three databases. The comparisons between the predicted results with SP fold and SP native supports our conclusions in the previously two benchmark tests that SP fold is more reliable than SP native . The loops/junctions that we failed to predict mostly involve tertiary interactions with helices or other cofactors such as protein and DNA.
The predicted results using the statistical potentials SP fold and SP native , extracted from the Leontis dataset also support the above conclusions ( Fig. 9 and Table 4). Fig. 13 shows the sensitivity of the convergence threshold parameter (Du n ij (h, g) = 0.001 v.s. 0.01) to the RNA structural predictions. The comparisons of the numbers of the successfully determined loops and junctions in the three databases show that the convergence threshold parameters Du n ij (h, g) do not have strong influence on the accuracy of the structure predictions.
respectively. In each figure, the red bars represent the statistical potentials SP fold , the green bars represent the statistical potentials SP native , both of which are obtained from the RNA09 dataset, and the x-axis stands for the dinucleotides with different nucleotides (i, j = ) AA, AC, AG, AU, CA, CC, CG, CU, GA, GC, GG, GU, UA, UC, UG and UU from 1 to 16. doi:10.1371/journal.pone.0048460.g002 Furthermore, we categorize the sequences in the dataset TEST-III and the predicted results according to the types of RNA loops and junctions. The 1119 RNA loops and junctions can be grouped into 408 hairpin loops, 166 internal/bulge loops, 148 pseudoknot loops, 329 multibranched loops and 68 junctions. The correct predictions with SP fold for each type of RNA loops and junctions are 395, 165, 147, 297 and 66, respectively (Table 3C). For junctions, multibranched loops, pseudoknots, internal/bulge loops and five of the hairpin loops, the failed predictions are due to the interactions beyond the dinucleotide context, such as the loophelix tertiary interactions and RNA-protein interactions. The possible reason for other six failed predictions for hairpin loops is that these hairpin loops are not closed by canonical base pairs.   (success rate = 92.4%), respectively. The total number of the successful predictions is 2969 (success rate = 91.9%).

Statistical potentials and the Leontis dataset
The March 17, 2012 version of the Leontis dataset contains 642 RNA sequences with structures of resolution better than 4.0 Å . The 3DNA software identified a total of 435 RNA loops and junctions with lengths ranging from 3 nt to 8 nt, excluding the 59/ 39-terminal dangling regions. Our calculations show that the difference between the diamond lattice-represented structures and the PDB structures varies for the different loop lengths and sequence contents with RMSD from 0.74 Å to 3.93 Å (Fig. 8), and the mean and standard deviation of RMSD values are 1.37 Å and 0.29 Å , respectively. We use such coarse-grained structures to calculate the observed dinucleotide frequencies, to extract SP native and to search for the ''true'' potential functions SP fold . Figs. 3 and 4 show the comparison between the potentials SP native /SP fold extracted from the dataset RNA09 and from the Leontis dataset. We apply the ''true'' potential SP fold extracted from the Leontis dataset to predict the loop/junction structures in three testing datasets: TEST-I, TEST-II and TEST-III, constructed from the January 2012 version of PDB database. Our test results (Table 4) show a success rates of 87% (7369 out of 8452) for TEST-I, 88% (6567 out of 7459) for TEST-II, and 96% (1077 out of 1119) for TEST-III, respectively.  Moreover, the predictions of the top 20 lowest-potential conformations for the loops and junctions in all the three test datasets, TEST-I, TEST-II and TEST-III, are shown in Fig. 9, with comparisons with the predictions based on the SP fold from the RNA09 dataset. The comparisons for the three test datasets (Table 4 and Fig. 9) show that the two ''true'' statistical potentials SP fold extracted from the RNA09 dataset and the Leontis dataset lead to similar success rates in structure prediction and the predictions are not sensitive to the choice of the specific training set.

Discussion
Motivated by the biological significance to predict sequencedependent loop and junction structures, we have developed a knowledge-based scoring functions/potentials to predict the structures of RNA loops and junctions. We use a coarse-grained conformational model (the virtual bond model) to sample RNA loop/junction conformations. From the known RNA structures, we extract a set of sequence-dependent dinucleotide-based statistical potentials using two methods. In the first method, the statistical potentials (SP native ) are derived from the distributions of the dinucleotide conformations in the known (native) structures. In the second method, the statistical potentials (SP fold ) are derived based on the folding stability of the native structures (against all the other nonnative folds).
For a given sequence, the extracted statistical potentials enable ranking of the different conformations with the top ranked (the lowest-potential) structure as the predicted native structure. Extensive tests indicate that the statistical potentials can successfully predict the native structure for a large class of loop and junction sequences and that SP fold consistently outperforms SP native in structure prediction. Our test results also indicate that our results are not sensitive to the choice of the specific training dataset ( Fig. 9 and Table 4).
The present approach has several advantages. First, the extraction of the statistical potential SP fold is based on the sampling of the complete conformational ensemble, including the native and all the nonnative folds. Second, because we consider all the nonnative folds in the derivation of the statistical potential, our statistical potentials can be used to predict folding from the sequence. The build-up strategy for RNA loop structure in this study is a de novo approach. Compared with other 3D loop structure prediction models, such as ModeRNA [22] and RLooM [50], our loop/junction structure prediction method is based on the complete ensemble of the (coarse-grained) conformations and does not rely on the information of template structures or homologous RNA structures. The only input information for the prediction is the sequence. Therefore, the model can predict the low-resolution structure from the sequence if no known homol-  ogous conformations can be found in the PDB. With the predicted low-resolution scaffold, one may further predict the all-atom structures using the all-atom potential methods that are derived based on near-native structures [5,15]. The present coarse-grained model offers a useful complement to the other all-atom based models.
We have performed extensive benchmark tests using the different databases. However, a direct comparison between our model other methods is not very straightforward. This is because our model is a coarse-grained model while other models mainly focus on the all-atom structures. Furthermore, our model aims to fold a low-resolution structure from the sequence without using any input information such as homologous templates, while other models mainly focus on the prediction of the native structures from near-native folds. Future development of our method, which may give all-atom structures from the low-resolution folds, would make direct comparison between our model and other models possible.
Applications of the present statistical potentials to the Vfold structure prediction model [14] may provide an effective strategy for better structure prediction. Despite the success, there are several limitations of our model. First, the current study is based on the coarse-grained Vfold model, which only provides a lowresolution approximation for the all-atom structure of RNA loops and junctions. Ultimately an all atom-based model is required to treat the detailed tertiary interactions. Future development of the model should address the issue to build all-atom RNA structures based on the low-resolution (Vfold-generated) representations. Second, the statistical potentials are derived for dinucleotide conformations. In realistic loop and junction structures, the longrange sequence effects within the loop, such as intra-loop  interactions, can influence the loop structure. The present model cannot explicitly account for such long-range intraloop interactions, which occur frequently in loops and between the different loops. Following the same procedure outlined above, one may develop a trinucleotide or higher order many-body statistical potentials to address this issue. Third, in the current form of the model, we use the same set of statistical potentials for the different types of loops/junctions. More refined potentials according to RNA loop/junction types may lead to further improvement of the accuracy of model.

Vfold model
We use the Vfold model to generate the full conformation ensemble of a given RNA loop or junction sequence. The Vfold model is a virtual bond-based RNA folding model [42,45,48], which is developed based on the two observations: the C{O torsion in the nucleotide backbone of RNA tend to adopt the trans (t) rotational isometric state (Fig. 5A) [57]. Therefore, the nucleotide backbone conformations can be reduced into two effective virtual bonds P{C Ã 4 and C Ã 4 {P [58][59][60]. The length of each backbone virtual bond is about 3.9 Å . The virtual bonds show rotamer-like configurations gauche z (g z ), trans (t) and gauche { (g { ) (Fig. 5C) [55,56,61]. Such virtual bond conformations can be well represented by the bonds in a diamond lattice. Therefore, we can use self-avoiding walks on a diamond lattice to enumerate the conformations of RNA sequences. This approach promises proper treatment of the excluded volume effect between the different atoms and the complete sampling of the conformational ensemble.
The loop structures are enumerated on a diamond lattice through exhaustive self-avoiding walks. If the first atom (P) is fixed at position X 0 ! , then the coordinates (X N ! ) of N-th atom (P or C Ã 4 ) can be calculated with where, l i~3 :9A is the bond length of virtual bond b i and b i ! is the unit vector of virtual bond b i . Also, the relation between the (i{1)th virtual bond (b i{1 ) and the i-th virtual bond (b i ) is: where, b i and h i are the related bond angles b i = 120 0 and pseudotorsion angle h i = g z (60 0 ), g { (300 0 ) or t (180 0 ) (Fig. 5C). The matrix T is defined as By considering excluded volume effect, i.e., two atoms cannot occupy the same site on a diamond lattice, we can generate the atomic coordinates using equations 1, 2, 3 and 4 and the conformation ensemble of RNA loops; see Ref. [42] for the detailed calculations. Here we give an example for illustration. In Fig. 5, following the 59?39 direction, if the first atom P is fixed at X 0 ! = col (0, 0, 3.9) (Å ), the second atom C Ã 4 is fixed atx x 1 = col (0, 0, 0) (Å ), then, b 1 ! = col (0, 0, 21). According to Eq. 1, the coordinate of the third atom P can be computed: where l 1~l2~3 :9 Å are the bond lengths, b P (~120 0 ) is the bond angle and g is the torsion angle (we select 180 0 for illustration). According to Eq. 3, the matrix is computed as: Therefore, the coordinate of the third atom P is col (0, ffiffi ffi 3 p /2, 21/ 2) (3.9 Å ). Following the procedure and considering the excluded volume between any two atoms, we can compute the coordinates for other atoms and generate the conformations for a given RNA loop.

Statistical potential
Assuming an equilibrium Boltzmann distribution for the structures in the PDB [53,54] and NDB [62] database, one can extract the interaction potentials:  Here, k B is the Boltzmann constant, T is the absolute temperature, r(r) is the observed density, and r Ã (r) is the density in a ''reference'' state where no interactions occur.
The above approach can give continuous distance-dependent pairwise potentials [15,[63][64][65][66] or discrete base pairs/stacks potentials [27,67]. In this work, we extract the potentials for the different dinucleotide conformations described by the pseudotorsion angles h (P{C Ã 4 {P{C 4 ) and g (C Ã 4 {P{C 4 {P) and the different dinucleotide sequences ( Fig. 5). The potential u ij (h, g) is extracted based on the following formula: where i and j denote the types of the nucleotides (bases) of the two consecutive nucleotides within the same loop and (h, g) are the pseudo-torsion angles of the dinucleotide.

Reference state
The proper choice of the reference state is a key issue in deriving the statistical potentials from the known structures. As pointed out by Thomas and Dill [68,69], an accurate ideal reference state is not achievable. Different approximations such as the quasichemical approximation [64] have been used to model the reference state. Furthermore, to circumvent the reference state problem, an iterative method for successive refinement of the statistical potentials was developed and was shown to give reliable scoring functions for protein folding [68,69] and protein-ligand interactions [63]. Here, for comparison, we employ the above two approaches to extract two different sets of RNA knowledge-based potentials and apply the extracted potentials to RNA loop and junction structure prediction.
Quasi-chemical approximation-based approach. In the quasi-chemical approximation [64] we use an ''expected state'' as the ''reference state''. Such an approximation leads to the Boltzmann relation (Eq. 6): where r obs ij (h, g) and r exp ij (h, g) are the observed number and the expected number, respectively, of the dinucleotides (i, j) with pseudo-torsion angles (h, g) in the entire training dataset. r exp ij (h, g) is computed from the following formulas: where N(h,g) is the total observed number of dinucleotides with pseudo-torsion angles h and g (calculated with a compositionindependent scale, meaning each nucleotide is independent to other nucleotides in the state (Eq. 8)) [70], and x i and x j are the mole fraction of nucleotide i and j of the sequences in the entire training dataset, respectively [70]. The quasi-chemical approximation-based statistical potentials are the ''extracted'' energy-like parameters derived from the observed density in the database of native structures [68]. We denote such extracted statistical potentials as SP native . The SP native parameters are not the ''true'' potentials extracted from the criteria of identifying the ''correct''/native structure from an ensemble of the nonnative conformations. We denote the potentials that can account for the nonnative conformations as the ''true'' statistical potentials (SP fold ).
Convergent iterative approach. To extract the ''true'' energy (SP fold ), Thomas and Dill developed an iterative approach based on the folding stability of the native structure [68,69]. The method has the advantage of accounting for the effect of the distribution of the whole conformational ensemble, including both the native and the nonnative states, for a given sequence. The overall strategy of the iterative approach is to train a set of potential parameters iteratively until the collection of the native structures in the training dataset and the full conformational ensemble (native and nonnative) lead to the same frequencies of the different dinucleotide conformations [63,68,69].
In our calculation, the iterative process starts with a set of initial values for the potentials, u (0) ij (h, g). The superscript (n) denotes the n-th iterative step. We use the quasi-chemical approximationbased statistical potentials SP native as the initial input.
At each step, we compute the total potential/score for each conformation of loop/junction by summing over the potential of all the dinucleotides contained in the loop/junction: Then the predicted frequencies (numbers) of the different dinucleotide conformations are computed from the Boltzmann average over all the conformations: where E (n) (s, l) is the potential/score for the l-th conformation of the s-th RNA loop/junction sequence, computed with the potentials u (n) ij (h, g) at the n-th step (Eq. 9), n ij (h, g) is the number of dinucleotide (i, j) with pseudo-torsion angles (h, g) in the l-th conformation of the s-th RNA loop/junction, L is the number of conformations of the s-th RNA loop/junction in the training dataset, and S is the number of RNA loop/junction sequences in the training dataset.
In each step, we also calculate the difference between the observed and predicted dinucleotide frequencies: where r pred(n) ij (h, g) is calculated by Eq. 10 and r obs ij (h, g) is the dinucleotide frequency observed from the ''correct'' (native) structures in the training dataset.
In general, the initially guessed potential parameters are not equal to the ''true'' potentials and thus there are differences between the observed dinucleotide frequencies and the predicted dinucleotide frequencies (Du (n) ij (h, g)&0). For each n-th step, we refine the potentials by accounting for the differences between the ''observed'' and the ''predicted'' frequencies of the dinucleotide conformations (see the parameter Du (n) ij (h,g) in Eq. 11): We repeat the above iterative process until DDu (n) ij (h, g)D approaches 0 (smaller than a threshold number, e.g. 10 {3 ). The final set of potentials u (n) ij (h, g) is the ''true'' potentials SP fold .

Flowchart
The iterative method is summarized as follows ( Fig. 1 4. Calculate the potential/score for each generated conformation of each RNA loop/junction sequence by using Eq. 9 with iterative potentials u (n{1) ij (h, g). Then calculate the weighted predicted frequencies r pred(n) ij (h, g) of each dinucleotide (i, j) with pseudotorsion angles (h, g), by employing Eq. 10. 5. Calculate the convergence parameter Du (n{1) ij (h, g) according to Eq. 11. If the following convergence condition is satisfied: with a preset small value of the threshold parameter C 0 , say, 10 {3 , then skip to the final step and the ''true'' potentials SP fold are obtained; otherwise, continue to the next step. 6. Adjust the current potentials u (n{1) ij (h, g) by adding the correction Du (n{1) ij (h, g) (Eq. 12), and obtain a new set of potentials u (n) ij (h, g). Move to the next (n-th) iterative step and return to Step 4 for the next cycle.

Preparation of the training dataset
We use the RNA 2009 database (RNA09; http://Kinemage. biochem.duke.edu/databases/rnadb.phb, an updated version of previous 2005 database (RNA05) [51]) to extract the SP native potential and use the database as the training dataset to search for the ''true'' potential SP fold . The RNA09 dataset consists of 262 RNA sequences with the experimentally determined structures of resolution ƒ3.0 Å .
RNA loops and junctions are constructed by all the unpaired nucleotides in RNA structure. We detect the loops and junctions by removing the nucleotides involved in the base pairs. Furthermore, we remove all the loop structures with modified bases or non-RNA atoms and molecules. Advances in the knowledge of the structures of both the canonical (Watson-Crick and wobbles) base pairs and non-canonical base pairs (normally classified as tertiary interactions) [71,72] enable the development of several automated tools to detect the base pairs in RNA structures. In this work, we use the 3DNA software package (http://3dna.rutgers.edu/home) [73] to search for all the possible base pairs from the RNA structures.
In addition, the time consumption to enumerate the loop conformations on a diamond lattice increases exponentially with the length of the loop [45]. Therefore, we only choose loops/ junctions with length ƒ8 nt in the dataset due to the long computer time for the exhaustive enumeration of the loop conformations for longer loops [45]. In this study, we use the algorithm reported in Ref. [74] to search for the best fits on diamond lattice for the RNA loops/junctions. Larger loops often involve tertiary interactions with other subunits of RNA structures, which are not accounted for in this study. As a result, we found 152 loops/junctions with the length ranging from 3 nt to 8 nt in the RNA09 dataset. These loop and junction structures are used as the training dataset in our iterative approach.
To further test the influence of the use of the different training dataset on the statistical potentials and the predictive power of the statistical potentials, we also use another non-redundant highresolution RNA structure dataset, collected by Leontis' lab (http://rna.bgsu.edu/nrlist/).

Conformational ensembles of RNA loops/junctions
We use the Vfold model to generate the full conformation ensemble of a given RNA loop or junction sequence. In the Vfold-generated loop/junction conformational ensemble, because each of the two pseudo-torsion angles (h and g) of a dinucleotide occupies three torsional states in a diamond lattice, there exist 4|4|3|3 parameters for the potentials for the four types of each nucleotide (A, C, G and U) and the three possible rotamerlike states for each torsional angle (and hence 3|3 types of the pseudo-torsion angle pairs h and g).
For each PDB structure of the loop/junction, we find the minimum-RMSD fit of the virtual bond conformation on the diamond lattice. We call such a coarse-grained (correct) structure as the ''coarse-grained correct structure''; see Step 1 in Fig. 1. Our calculations show that the differences between the diamond latticerepresented structures and the PDB structures varies for the different loop lengths and sequence contents with RMSD in the range from 0.67 Å to 2.07 Å (Fig. 7), the mean and standard deviation of RMSD values are 1.35 A and 0.30 Å , respectively.. We will use the coarse-grained correct structure to calculate the observed dinucleotide frequencies, from which SP native and SP fold are extracted.