RNA-Puzzles Round IV: 3D structure predictions of four ribozymes and two aptamers.

RNA-Puzzles is a collective endeavor dedicated to the advancement and improvement of RNA 3D structure prediction. With agreement from crystallographers, the RNA structures are predicted by various groups before the publication of the crystal structures. We now report the prediction of 3D structures for six RNA sequences: four nucleolytic ribozymes and two riboswitches. Systematic protocols for comparing models and crystal structures are described and analyzed. In these six puzzles, we discuss (i) the comparison between the automated web servers and human experts; (ii) the prediction of coaxial stacking; (iii) the prediction of structural details and ligand binding; (iv) the development of novel prediction methods; and (v) the potential improvements to be made. We show that correct prediction of coaxial stacking and tertiary contacts is essential for the prediction of RNA architecture, while ligand binding modes can only be predicted with low resolution and simultaneous prediction of RNA structure with accurate ligand binding still remains out of reach. All the predicted models are available for the future development of force field parameters and the improvement of comparison and assessment tools.

INTRODUCTION

The RNA-Puzzles competition was initiated in 2011 following discussions and exchanges at the annual meetings of the RNA Society. The blind assessment of RNA predictions was started along broadly similar procedures for protein predictions organized within the Critical Assessment of Protein Structure Prediction (CASP) experiment (Moult et al. 2014). While CASP releases protein sequences publicly, the RNA Puzzles experiment follows strict rules of confidentiality. With the agreement of a structural biology group, an RNA sequence, with the X-ray structure recently determined in that laboratory, is submitted to modeling groups that all agree to keep the information confidential. Because the structures are unknown to the modelers and program developers, and absent from databases, the procedure minimizes bias in the prediction algorithms. The modeling groups have a deadline for submitting the results of the prediction of the puzzle. The deadline depends on the delay between submission and publication of the structure (with the public release of the coordinates) as well as on the type of modeling exercise that is being tested (with or without manual intervention). The prediction time for the automated web server is set to 48 h, while human experts have 3–4 wk. The results of the comparisons between the experimental and theoretical structures are then published together with all groups involved in the enterprise. The puzzles are grouped into rounds based on the similarity of function and approximate cooccurrence in prediction time. Three rounds of RNA-Puzzles have been previously reported in the literature (Cruz et al. 2012; Miao et al. 2015, 2017). Over these rounds, the RNA-Puzzles consortium, through numerous exchanges and discussions, has fostered a better understanding between the experimental and theoretical groups, as well as between the groups involved in the development of the modeling tools. Further, based on the feedback from the evaluation, the modeling and assessment tools have improved in quality and ease of use during the whole process.

Ribozymes (Doherty and Doudna 2000; Cech 2009) are RNA molecules that can catalyze specific biochemical reactions, in a way similar to the action of protein enzymes, while riboswitches (Breaker 2011; Garst et al. 2011) are regulatory segments of protein-coding RNAs that bind a small molecule, resulting in a change in the protein translation. Both ribozymes and riboswitches have biological functions achieved through ligand/small molecule binding (Zhang et al. 2010). In particular, some ribozymes and riboswitches may share similar structural folds, such as the c-di-GMP riboswitch, the guanine riboswitch and the hammerhead ribozyme (Porter et al. 2017). According to the Rfam database (Kalvari et al. 2018), 27 out of the 42 riboswitch families have members with experimentally determined structures, while for the 31 ribozyme families only six have structures experimentally solved. Predicting these functionally important RNA structures can provide significant structural insights to the RNA functions and help in RNA engineering. In particular, the self-cleaving nucleolytic ribozymes became a hot topic in the field of structural biology of the ribozyme (Ren et al. 2016). The fair assessment of the computational prediction of these newly discovered ribozymes underscores a key aspect of RNA structure prediction and function inference.

Here, we report our community effort to predict three-dimensional structures of ribozymes and riboswitches. The six puzzles in this round of prediction include four ribozymes: the pistol ribozyme, the hammerhead ribozyme, two twister-sister ribozymes, and two riboswitches: the 5-hydroxytryptophan and the Guanidine-III riboswitch. A few additional puzzles were carried out concomitantly with these six puzzles corresponding to non-coding RNA domains with functions distinct from chemical catalysis or small-molecule recognition and will be discussed in a separate manuscript. According to Round IV prediction results, we highlight: (i) the progress in automated web server approaches and the still existing gap between the web servers and human experts; (ii) key rules that play important roles in prediction; (iii) the existing bottleneck in ligand/ion binding prediction. The prediction methods developed within our community will provide a practical resource to RNA biologists.

RESULTS

Availability

All the prediction results are available at www.rnapuzzles.org. The structure manipulation tools and all the assessment metrics together (Magnus et al. 2020) with the predicted data are now available in an open-source repository at https://github.com/RNA-Puzzles.

The six prediction targets

Puzzle 9: The 5-hydroxytryptophan aptamer

5-hydroxytryptophan (5-HTP) is a naturally occurring amino acid and chemical precursor of serotonin, which is also marketed in many European countries for the treatment of major depression. Porter et al. (2017) introduced the aptamer as a scaffold for genetically encodable small-molecule biosensor, which selectively binds to 5-HTP. The 71-nucleotide (nt)-long sequence is:

5′-GGACACUGAUGAUCGCGUGGAUAUGGCACGCAUUGAAUUGUUGGACACCGUAAAUGUCCUAACACGUGUCC-3′

This puzzle structure was only predicted by human experts. The 2 Å crystal structure was deposited in the Protein Data Bank (PDB) (Berman 2000) (5kpy) after the prediction. According to the scaffold analysis (Porter et al. 2017), the secondary structure of 5-hydroxytryptophan riboswitch is similar to the guanine riboswitch (Batey et al. 2004), di-GMP riboswitch (Kulshina et al. 2009), and hammerhead ribozyme (Anderson et al. 2013), which have been reported. The 5-hydroxytryptophan aptamer adopts a topology identical to that of the guanine riboswitch, with a backbone root mean square deviation (RMSD) of 2.5 Å, Supplemental Figure S1a.

Puzzle 15: The hammerhead ribozyme

Although the hammerhead ribozyme was the first of the small nucleolytic ribozymes to be described (Uhlenbeck 1987), the initial structures were trapped in an inactive conformation (Pley et al. 1994; Scott et al. 1995). The binding to the pentavalent transition state analog vanadate results in significant core rearrangements in the structure (Mir and Golden 2016) compared to the initial structures. Considering the availability of the solved hammerhead ribozyme structure, this puzzle was predicted by both fully automated web servers and human experts to benchmark the difference between these two types of approaches. Two rounds of predictions were organized. In the first round, only the RNA sequences were given for prediction; in the second round, the secondary structure, based on the crystal structure, was given as additional information. Ligand or ion binding information was not given in the prediction. The structure was predicted as two chains, 48 nt + 20 nt. The sequences are:

5′-GGGUACUUAAGCCCACUGAUGAGUCGCUGGGAUGCGACGAAACGCCCA-3′

5′-GGGCGUCUGGGCAGUACCCA-3′

The crystal structure was solved at resolution 2.95 Å and deposited in PDB (5di4). In the crystal, the 5′ end of the first chain and the 3′ end of the second chain form intermolecular base pairs between symmetry-related molecules, forming, therefore, a dimer structure. Thus, the known base pairs between the 5′ end of the first chain and 3′ end of the second chain, verified through functional studies are not formed between the expected strands in the crystal structure. Unless explicitly stated, crystal contacts are not considered for comparison. Most models rightly contained a regular stem between the 5′- and 3′-ends. Automatic comparisons, therefore, lead to larger (and penalizing) deviations for those nucleotides (about 6). The catalytically relevant part of the molecule is in the core of the structure.

Puzzle 17: The Pistol ribozyme

The recent discovery of pistol (this puzzle), twister-sister (Puzzles 19 and 20) and hatchet (Puzzle 22) ribozymes (Weinberg et al. 2015) have invigorated the field of small self-cleaving nucleolytic ribozymes. Pistol ribozyme studies have yielded significant insights into the G-U cleavage mechanism, which is facilitated by its pseudoknot fold (Ren et al. 2016). The 62-nucleotide (nt)-long sequence is:

5′-CGUGGUUAGGGCCACGUUAAAUAGUUGCUUAAGCCCUAAGCGUUGAUAAAUAUCAGGUGCAA-3′

This puzzle was predicted by both web servers and human experts. After the prediction deadline, the crystal structure (Ren et al. 2016) of resolution 2.73 Å was deposited in the PDB (5k7c). The G53-U54 cleavage site is the key functional site and from a biological point of view, the accuracy of predicted G53 interactions is most relevant.

Puzzle 19: The Twister sister ribozyme 1

Another self-cleaving ribozyme is twister-sister, which includes two coaxial helical stacks organized by a three-way junction and two tertiary contacts (Liu et al. 2017). It functions predominantly as a metalloenzyme. Both human expert prediction and automated web server prediction were performed. The fold was predicted as a two-chain (40 nt + 22 nt) structure with the following sequences:

5′-GCAGGGCAAGGCCCAGUCCCGUGCAAGCCGGGACCGCCCC-3′

5′-GGGGCGCGGCGCUCAUUCCUGC-3′

The crystal structure of resolution 2 Å in PDB (5t5a) was deposited after the prediction. This puzzle was a good case to evaluate the prediction of coaxial stacking (Lescoute and Westhof 2006) as well as the long-distance tertiary contacts.

Puzzle 20: The Twister sister ribozyme 2

This is a variant of a twister-sister ribozyme that includes a four-way junction. As in Puzzle 19, the structure of Puzzle 20 contains several long-distance tertiary contacts and a number of RNA-bound magnesium ions. These two structures adopt a similar global topology, but they exhibit significant conformational differences (Supplemental Fig. S2), with a backbone RMSD > 10 Å. Although Puzzle 19 prediction was organized before Puzzle 20, Puzzle 19 structure was not made available as a template for homology modeling of Puzzle 20. The structure of Puzzle 20 was predicted both by automated web servers and by human experts. It was also predicted as a two-chain structure (18 nt + 50 nt), with the longer sequences below:

5′-ACCCGCAAGGCCGACGGC-3′

5′-GCCGCCGCUGGUGCAAGUCCAGCCACGCUUCGGCGUGGGCGCUCAUGGGU-3′

Structures coming from the native and MnCl2-soaked crystals were solved and deposited in the PDB as 5y85 (resolution 2.00 Å) and 5y87 (resolution 2.13 Å), respectively. The backbone RMSD between these two structures is only 0.4 Å. Since both the solved structures (5y85 and 5y87) can be considered as “ground truth” of the sequence, both of them were used as reference structures in the later assessment.

Puzzle 21: The Guanidine-III Riboswitch

Barrick et al. (2004) identified the ykkC family as probable riboswitches more than a decade ago. It was eventually found that the ykkC family RNAs respond to guanidine (Nelson et al. 2017), and the ykkC family was later classified into three subclasses by Breaker and colleagues, termed as the guanidine-I (Nelson et al. 2017), -II (Sherlock et al. 2017), and -III (Sherlock and Breaker 2017) riboswitches. This puzzle involves the prediction of the structure of the guanidine-III riboswitch (Huang et al. 2017b), which corresponds to the Rfam family (RF01763) (Kalvari et al. 2018). The 41-nucleotide (nt)-long sequence is:

5′-CCGGACGAGGUGCGCCGUACCCGGUCACGACAAGACGGCGC-3′

After both web server and expert predictions were performed, two crystal structures were deposited in the PDB: one structure solved in space group P3121 at 1.91 Å resolution (PDB: 5nwq) and another structure in space group P3212 at 2.94 Å resolution (PDB: 5nz6). The two structures share a backbone RMSD of 0.58 Å and both were considered as reference structures in the assessment, Supplemental Figure S3. The structure is dominated by a pseudoknot fold and the guanidine binding site is formed by a triplex left-handed helical core. This is the first puzzle for which the identity of the binding ligand was provided before prediction.

Overall prediction results

The distributions of all the assessment metrics in the sets of all models submitted for every puzzle are displayed in Figure 1. One important aspect to understand the predictions is based on comparisons with the best model. From the RMSD distribution, the best-predicted models start at RMSD values of 3.65 Å (PZ21) and 6.45 Å (PZ15), Figure 1A. But with consideration of the sequence length as in the Deformation Index, the best prediction models for all the six puzzles are comparable, Figure 1B, ranging from 4.75 (PZ21) to 8.64 (PZ9). The worst predicted models peak at about 25.55 Å RMSD (PZ20). In Figure 1B, roughly two ensembles of models can be distinguished for Puzzles 19 and 21. Although the RMSD range spreads, Puzzles 20 and 21 display lower RMSD predictions, while Puzzles 17 and 19 yielded predicted models centered around 17 Å. Although most of the predicted structures show very low Clash scores, still some outliers exist in the distribution, Figure 1C. Supplemental Figure S4 displays the various metrics as a function of RMSD. Deformation Index shows a distribution proportional to RMSD. For the Clash score, the range is independent of RMSD, except for large outliers, which is not easy to explain (Supplemental Fig. S4a).

FIGURE 1.

The distributions of structure comparison metrics of the six puzzles. The distribution of structure assessment metrics is plotted as the violin plots where each prediction model is shown as a solid line in the violin plot. Thus, the vertical spread displays the range of values for a given metric in all predicted models and the horizontal breadth reflects the number of predicted models at a given value. The assessment metrics include (A) RMSD, (B) Deformation Index (Parisien et al. 2009), (C) Clash score, (D) the Mean of Circular Quantities (Zok et al. 2014) and (E–H)Interaction network fidelity (Parisien et al. 2009) for all parameters together, Watson–Crick pairs, non-Watson–Crick pairs and stacking.

The correct prediction of the interactions is the premise of accuracy in structure prediction. From the distribution of INF (Interaction Network Fidelity for Watson–Crick base pairs), prediction of Watson–Crick interactions is near perfect in the best model for most of the puzzles, ranging from 84% (PZ15) to 100% (PZ20), Figure 1E–H. However, the INF ranges are rather large in several puzzles. The recurrent structural modules are normally formed by non-Watson–Crick interactions. For example, guanidine, as a ligand, in Puzzle 21 is hydrogen-bonded to the Hoogsteen edges of two guanine bases (Huang et al. 2017b). The prediction of non-Watson–Crick is more difficult than Watson–Crick base pairs, thus, showing a much wider range in the INF distribution, except for Puzzles 17 and 21, for which the accuracy of INF is very low. The prediction of base stacking is highly related to the Watson–Crick interactions and the distribution of INF generally follows the same pattern as the accuracy of Watson–Crick interactions INF. Interestingly, in Figure 1D, two groups of models appear with very long tails in the MCQ distribution (PZ17 and PZ20).

Detailed prediction results

Below, the results specific for each RNA-Puzzle are described. Figures 2–7 contain a superposition between the crystal structure and the model with the lowest RMSD, a secondary structure diagram, and a heat map of the Deformation Index (DI) (Parisien et al. 2009). In all figures, the secondary structure elements are colored and numbered (helix P1, P2,…, pseudoknot PK). The corresponding colored numberings are indicated on the 3D drawings. In the Deformation Index plots, the secondary structure elements are delimited by colored rectangles within the heat maps. The RMSD is used as the primary ranking parameter because of its simplicity. However, RMSDs spread errors over the whole structure and one cannot tell which parts of the molecule were correctly modeled and which were not. The Deformation Index attempts to alleviate this drawback in the following way. The structures to be compared are superimposed locally three nucleotides at a time (following the sequence) and the average deviations between the central nucleotide and all the other nucleotides are calculated and plotted as a heat map. Thus, all domains globally superimposable appear in white while those departing from that global core appear in red.

FIGURE 2.

Puzzle 9: 5-hydroxytryptophan riboswitch aptamer. (A) 3D structure comparison between the predicted model with the lowest RMSD (Chen group model 7, shown in blue) and the reference structure in green. (B) The heatmap shows the deformation profile, poorly predicted regions are shown in red; the red regions concentrate on the loop escaping P2 and P3 with the linking region between P2 and P3. (C) The numbering and secondary structure of the reference structure.

FIGURE 3.

Puzzle15: Hammerhead ribozyme. (A) 3D structure comparison between the predicted model with the lowest RMSD (Adamiak model 8, shown in blue) and the reference structure in green. (B) The heatmap shows the deformation profile, in which the poorly predicted regions are shown in red. The largest deviations happen between the 5′- and the 3′-ends. The large deviations, between the 5′- and 3′-ends in the blue box, are artifactual as discussed in the text (see also (D)). (C) When the 5′- and the 3′-ends are removed from evaluation, RNAComposer model automated server model 1 in the round 1 prediction ranks the best. The heatmap shows the deformation profile, in which the poorly predicted regions are shown in red. calculations of the heatmap. (D) The secondary structure of the reference structure. Nucleotides 1–6 and 63–68 (the gray dashes) within the blue rectangle do not form intramolecular base-pairing but instead, they form base pairs with symmetry-related molecules in the crystal (see text for discussion).

FIGURE 4.

Puzzle17: Pistol ribozyme. (A) 3D structure comparison between the predicted model with the lowest RMSD (SimRNA model 1 in the second round prediction, shown in blue) and the reference structure in green. (B) The heatmap shows the deformation profile, in which the poorly predicted regions are shown in red. The deviations are limited to the last two nucleotides of the long strand and the 5′-end of the short second strand (extremities of P3). (C) The secondary structure of the reference structure.

FIGURE 5.

Puzzle19: Twister sister ribozyme 1. (A) 3D structure comparison between the predicted model with the lowest RMSD (Chen group model 1, shown in blue) and the reference structure in green. (B) The heatmap shows the deformation profile, in which the poorly predicted regions are shown in red. (C) The secondary structure of the reference structure.

FIGURE 6.

Puzzle 20: Twister sister ribozyme 2. (A) 3D structure comparison between the predicted model with the lowest RMSD (Bujnicki group model 4, shown in blue) and the reference structure in green. (B) The heatmap shows the deformation profile, in which the poorly predicted regions are shown in red. (C) The secondary structure of the reference structure.

FIGURE 7.

Puzzle 21: Guanidine III Riboswitch. (A) 3D structure comparison between the predicted model with the lowest RMSD (LORES approach from Das group model 1, shown in blue) and the reference structure in green. (B) The heatmap shows the deformation profile, in which the poorly predicted regions are shown in red. The largest deviations occur at the single strand leaving P1 to join P2 (UCAG). (C) The secondary structure of the reference structure. (D) The guanidine binding site in the crystal structure, where guanidine is shown as magenta and proximal contacts are shown as cyan dashed lines.

Puzzle 9: The 5-hydroxytryptophan aptamer (Fig. 2)

Despite differences in local conformation, the 5-HTP riboswitch adopts the same topology as the guanine riboswitch. Therefore, the guanine riboswitch structure was used correctly by many groups as a global or local template in the prediction. Among the models, model 7 from the Chen group demonstrates the lowest backbone RMSD as 3.74 Å (6.06 Å all-atom RMSD), Figure 2. According to the structure comparison in Supplemental Figure S1b (guanine riboswitch-based similarity), the best model (Chen model 7) has a backbone RMSD of 1.54 Å from the previously reported guanine riboswitch (Batey et al. 2004). Homology modeling using templates may well capture the main topology but the local structural details in the loops of P2 and P3 regions show a difference from the crystal structure. It is not easy to evaluate how significant those changes are, since they occur in loop regions allowing dynamical movements and often with high crystallographic temperature factors as discussed previously (Miao and Westhof 2017; Miao et al. 2017; Magnus et al. 2020). The correct modeling of such local structure variability continues to be a major challenge for the modeling methods. In this puzzle, the crystal structure forms more turns around the ligand-binding region compared to the best-predicted model. As the ligand information was not disclosed during the prediction process, the structure was predicted as ligand-free. Knowing the ligand information may still be beneficial to the prediction process. A similar case was detected in Puzzle 10 (Round II): prediction groups pointed out that the kink-turn module in RNA can be easily predicted if the binding protein (YbxF) information is given. Model 5 from the Das group achieved the best INF score; Das group modeling used FARFAR (Das et al. 2010) with two starting templates: residues A:30–40 and A:56–70 came from 2XNW, while a tRNA T-loop from 3L0U furnished residues A:45–49. The RMSD values of the 34 submitted predictions range from 6.06 to 25.43 Å with an average value of 9.93 Å, as shown in Supplemental Table S1.

Puzzle 15: The hammerhead ribozyme (Fig. 3)

The prediction of the hammerhead ribozyme structure in the state of binding to a transition state analog is a difficult case for automatic prediction because the binding ligand was not given in the prediction as prior knowledge. RMSD of the 74 submitted models ranges from 6.45 Å to 24.78 Å. The predicted model of the lowest RMSD (6.45 Å) shows a topology similar to the experimentally determined structure. Figure 3 and Supplemental Table S2 contain the automatic comparisons between crystal structure and models. As discussed above, the strands forming stem I in hammerhead ribozymes involve intermolecular pairs between molecules in the crystal while in models stem I forms intramolecularly. Thus, automatic comparisons in that region are artifactual. This example shows the advantage of the Deformation Profile in which local errors do not spread throughout the structure (see Fig. 3B). In Supplemental Table S3, the last six nucleotides on each of the strands are removed from the comparison. Most of the resulting assessment parameters are much better. The corresponding Deformation Profile is shown in Figure 3C. As shown in Supplemental Figure S5, the RNAComposer (Popenda et al. 2012; Antczak et al. 2016) prediction shows better deformation profile and RMSD in the core structure (without stem I) than the SimRNA (Boniecki et al. 2016) prediction. SimRNA has a better global RMSD because of the orientation of the stem I region.

Puzzle 17: The Pistol ribozyme (Fig. 4)

The Pistol ribozyme structure is dominated by a pseudoknot structure and Watson–Crick interactions. According to Rfam, P1, P2 and PK regions of the Pistol ribozyme have significant covariations in the sequence alignment. The only uncertainty is the length of the pseudoknot: Rfam (Kalvari et al. 2018) structure includes five base pairs while R-scape (Rivas et al. 2017) predicts 7 bp, Supplemental Figure S6. And, in the crystal structure, the pseudoknot is made of six base pairs. Despite this minor difference in the number of base pairs between Rfam prediction and the crystal structure, the best prediction demonstrated a good accuracy with 5.15 Å RMSD (SimRNAweb server model 1 in the second round prediction), Figure 4. But only seven out of the total 107 predicted structures have RMSDs better than 10 Å, Supplemental Table S3. Models from each prediction group are ranked according to their estimated reliability. In this puzzle, the best RMSD model is SimRNA model 1 and the model with the best INF score (most interactions having been correctly predicted) is model 1 from the Das group (8.60 Å RMSD).

Puzzle 19: The Twister sister ribozyme 1 (Fig. 5)

The RMSD values of the 55 submitted predictions range from 5.52 to 22.44 Å with an average value of 14.33 Å, as shown in Supplemental Table S4. The challenges in the prediction lie in: a) the correct coaxial stacking of the P4 and P5 regions, b) the correct detection of long-distance contacts. These two issues are linked to each other. Some models predict the stacking between P4 and P1 rather than P4 and P5 (for example, model 3 from the Adamiak-Szachniuk group and model 5 from the Bujnicki group), Supplemental Figure S7. These models significantly deviate from the reference structure but also do not capture the key tertiary interactions between P4 and P5. As there is no gap between P4 and P5, the more likely stacking would have been predictable from the rules in Lescoute and Westhof (2006) and was predicted by other groups. As shown in Figure 5, the best model (Chen group model 1) correctly predicts the global topology as well as most of the contacts (seen from Deformation Profile). The only difference is that the reference structure is more compact. The question of assessing structure compactness has been brought up during the second RNA-Puzzles community meeting in 2018.

Puzzle 20: The Twister sister ribozyme 2 (Fig. 6)

Unlike Puzzle 19, Puzzle 20 includes four helices that form two coaxial stackings: between P2 and P3, between P4 and P5 (even though P5 has no gap between P2). Although Puzzle 20 has a longer sequence than Puzzle 19, it shows better prediction accuracy: 10 out of 50 predicted models have <10 Å RMSD, which comes from the Bujnicki, Das and Adamiak-Szachniuk groups, Supplemental Table S5. According to the 3D and 2D comparisons in Figure 6, the best-predicted model (Bujnicki model 4) demonstrates great consistency with the experimentally determined structure (RMSD 4.71 Å). However, in the low RMSD models, only model 5 from the Bujnicki group has a low Clash score, which indicates that there is still room for improving atomic clashes.

Puzzle 21: The Guanidine-III Riboswitch Aptamer (Fig. 7)

The guanidine-III riboswitch aptamer domain is a small structure forming a pseudoknot. The RMSD values of the 54 submitted predictions range from 3.65 to 21.57 Å with 20 models of <6 Å RMSD, as shown in Supplemental Table S6, while the lowest RMSD model (LORES model 1 from the Das laboratory) is compared with the reference structure in Figure 7. However, all of the high-quality models are from human expert prediction, while the best server predicted model (SimRNA model 4) has an RMSD of 10.6 Å. According to the RMSD distribution in Supplemental Figure S8, all human predicted models (except LORES model 3) perform consistently better than server-generated models. The majority of server predictions have low accuracy of INF, indicating the problem in getting the correct secondary structure. Predicting the pseudoknot structure, which can be manually achieved through inspection of phylogenetic alignments for covarying Watson–Crick pairs, remains a major challenge for web servers. Similar to the case in Puzzle 15, RNAComposer model 3 shows worse RMSD (12.12 Å) than SimRNA model 4 (10.57 Å), but RNAComposer's base-pair prediction (INF = 0.619) is still better than SimRNA (INF = 0.511). Another critical issue in this puzzle is the correct prediction of the guanidine (ligand) binding core structure formed by the three nucleotides. Despite the low-RMSD structures predicted by many groups, correctly predicting all contacts in the guanidine binding pocket is still difficult, Supplemental Figures S8h, S13. At the same time, there is a reason for optimism; in some models, contacts of guanidine with one RNA nucleotide, G7, were correctly predicted.

Server predictions compared to human expert predictions

In the present RNA-Puzzles Round IV, the server category was open for Puzzles 15, 17, 19, 20 and 21. The resulting distributions divided into human and server categories are shown in Supplemental Figure S12. Overall, server predictions were worse than human predictions in terms of accuracy, but compare favorably in terms of stereochemical quality. We compare server and human predictions for the five puzzles in turn.

In Puzzle 15, four servers, 3dRNA (Zhao et al. 2012), RNAComposer (Popenda et al. 2012), RW3D (Das et al. 2010), and SimRNAweb (Magnus et al. 2016), submitted a total of 54 predicted models. In the human category, there were only two participants who submitted 20 models in total. As far as RMSD is concerned, its value for the best model in the server category is 7.12 Å (best human = 6.45 Å), and for the worst one, it is 24.78 Å (worst human = 16.06 Å). The median RMSD value for server prediction is 17.03 Å and for human prediction 8.80 Å. Five server models (one SimRNAweb model and four RNAComposer models) have better RMSDs than the median score obtained for human predictions. Interestingly, there is no such discrepancy between human and automatic predictions in the INF scores, especially when considering Watson–Crick interactions only. When it comes to the Clash score, 3D models delivered by the servers actually outperform the predictions submitted in the human category. Nevertheless, even in the category of stereochemical quality, both server and human expert predictions could be improved. An analysis made with the use of the RCSB Maxit program (Berman 2000) shows that only 8 server-generated models (four submitted by RNAComposer and four by RW3D) and 10 models of the Adamiak-Szachniuk group have correct stereochemical parameters.

Puzzle 17 ended up with 38 models (in the two rounds of prediction: in the first round of prediction, no information on the 2D structure was given, while in the second round the 2D structure was given) submitted by two servers (RNAComposer and SimRNAweb) and 69 predictions by eight human groups. In this Puzzle, automatic predictions are qualitatively very close to the human ones, as far as the RMSD, INF and DI measures are concerned. For example, the lowest RMSD for server models is 5.15 Å (and for human models it is 7.13 Å) and the highest is 21.52 Å for server models (and for human models it is 21.75 Å). Thirteen server models obtained better scores than the human RMSD median value that reached 14.87 Å. On the other hand, 21 human models were worse than the automatic RMSD median (16.09 Å). 3D models in the human category have on average better MCQ and Clash score values. Only 23 models submitted in this Puzzle have correct stereochemistry. This set includes 12 automatic predictions submitted by RNAComposer and 11 human predictions (three of the Adamiak-Szachniuk group and eight by the Das group).

The collection of RNA structure 3D models submitted for Puzzle 19 includes 28 automatic models (from five servers: 3dRNA, LeeServer, RNAComposer, RW3D, and SimRNAweb) and 27 human models (from seven expert groups). Here, the median values of all considered measures for both categories are very similar. Although human predictions are more diverse—the ranges of their RMSD, INF and DI values are wider than for automatic models. For example, the latter have RMSDs in the range 9.59–19.58 Å, while the human models are in the range 5.51–22.44 Å. The DI computed for automatic models takes values between 15.48 and 38.28 Å and that of human expert models between 7.48 and 53.18 Å. The INF scores obtained for non-canonical and stacking interactions are very similar in both sets. Again, 11 server models are closer to the reference structure than the median human prediction (its RMSD value is 14.75 Å). The subset of human models contains 12 submissions with RMSD worse than the median score computed for server predictions which equal 15.25 Å. It should be emphasized that RCSB Maxit reports no errors in 14 cases only which come from nine automatic predictions (five from RNAComposer and four from RW3D) and five human predictions (all by the Adamiak-Szachniuk group). It clearly shows that in all subsequent predictions the check of stereochemical correctness should be performed. The RCSB Maxit program is a reliable option.

In Puzzle 20, there were three servers and five human participants. Servers (RNAComposer, RW3D, SimRNAweb) submitted 20 models, and human experts 35 models, respectively. In the case of this challenge, automatic models’ qualities visibly differ from what human experts predicted. The best RMSD for server submissions is 10.99 Å, while for human experts it is 4.70 Å. A significant difference can also be seen in the DI (14.01 Å vs. 5.68 Å) and MCQ (9.97° vs. 15.45°) values between the best automatic and the best human predictions. The surprise is the best automatic prediction (RNAComposer, model 1) with INF = 1, while the best human model (Adamiak, model 5) has 0.98 for this measure. Only four server models obtained lower RMSD scores than the median human prediction that has RMSD = 16.05 Å. On the other hand, five human submissions have higher RMSD than the median score computed for automatic models which equals 16.64 Å. In this Puzzle, only eight models pass the RCSB Maxit test with no stereochemical errors reported. These are three automatic submissions by RNAComposer and five human predictions by the Adamiak-Szachniuk group.

Finally, Puzzle 21 was closed with 25 automatic predictions from 4 servers (3dRNA, RNAComposer, RW3D, SimRNAweb) and 34 human predictions from five expert groups. Server models submitted in this challenge differ significantly in quality from the majority of human expert models. RMSD of automatic predictions is in the range 10.56-18.02 Å and the median value reaches 15.99 Å; Deformation Index is between 19.58 Å and 39.46 Å with the median equal to 31.46 Å. For human models, RMSDs are in 3.64–21.57 Å with the median equal to 5.31 Å and DIs are in 4.74–54.79 Å with median 8.07 Å. INF computed for Watson–Crick base pairs is the same for both the best server and the best human model. But when it comes to the evaluation of non-Watson–Crick interactions’ prediction, 3D models submitted by human experts are the winners of this challenge. Only when the Clash score is considered, automatic models perform better with their scores between 0 and 118.53, as compared to 7.52–218.94 range for human submissions. All models submitted in this Puzzle have stereochemical errors. For the best models, as far as the stereochemistry is concerned (13 automatic predictions from RNAComposer and RW3D, 10 human predictions from the Das and Adamiak-Szachniuk groups), RCSB Maxit reported problems with phospho-sugar bonds stereochemistry.

DISCUSSION

Overall, the results of RNA-Puzzles Round IV confirm observations made since earlier rounds of the puzzles (for reviews, see Miao and Westhof 2017; Magnus et al. 2020). In every test case, the current community of RNA structural modelers is able to produce blind predictions of 3D structures that are in agreement with crystal structures (Figs. 2–7). An important question then is whether these methods are ready to give biologically useful models for the many RNAs that have not been experimentally characterized but are suspected to have well-defined 3D structures. Only a very small fraction of the Rfam families (Kalvari et al. 2018) have structure(s) solved, which is 99 families of the 3016 Rfam families. Compared to the solved protein structures, the known RNA structure space is still limited. Computational prediction of RNA 3D structure may therefore exert great potential in understanding the unknown RNA structures, especially in generating potential structure models to explain the molecular functions. Automatic web servers will likely be critical to enable the broader RNA biology community to make such predictions, and so assessment of these servers compared to expert models has been an important theme of recent rounds of the RNA-Puzzles.

The prediction results in Round IV are compared with all the previously reported results in Supplemental Figure S13. One success achieved in Round IV is that all the best-predicted models have a RMSD < 7 Å. The high atomic clash of otherwise accurate models has been pointed out since the Round II prediction (Miao et al. 2015) and new structure relaxation approaches have been developed. Still, it appears to be difficult to achieve both accurate structure prediction and relaxation within the limited prediction time.

Since the first challenge of RNA-Puzzles, we have observed a growing interest in the development and use of automated methods supporting selected stages of the modeling of 3D RNA structures, including web servers that predict RNA 3D models. In RNA-Puzzles Round I (Cruz et al. 2012), the Major group tested MC-Fold/MC-Sym (Parisien and Major 2008), while the Das group applied their Rosetta SWA (Das 2013) algorithm. In Round II (Miao et al. 2015), the spectrum of computing methods used in RNA 3D structure prediction expanded to six programs (3dRNA, RNAComposer, DMD [Dokholyan et al. 1998; Ding and Dokholyan 2005; Ding et al. 2008; Proctor et al. 2011; Krokhotin et al. 2015], ModeRNA [Rother et al. 2011], Rosetta FARFAR, MC-Fold), some of them operating as fully-automated servers, while in Round III (Miao et al. 2017) there were five programs used (3dRNA, RNAComposer, RW3D, SimRNA, DMD). These methods are still mostly used as the elements of the modeling procedure and are complemented by expert knowledge, experimental data and manual or semi-automatic editing of the initial model. It was only in 2015 that rules for fully automated prediction were defined and, consequently, Puzzle 15 was opened in two categories, human and server. Since then, both the structures predicted by human experts over weeks and models simulated fully automatically within 2 d have been submitted for every puzzle. In the current Round IV, a large gap remains between the accuracy of server models compared to human expert models. Nevertheless, stereochemical quality (e.g., atom-atom clashes) of the server models was similar to and in some puzzles, better than human expert models, suggesting that further automation of expert procedures could help ameliorate these pervasive stereochemical errors. In addition, this round of puzzles has inspired several participating groups to improve automation of their methods, including the development of new servers that are being tested in more recent RNA-Puzzles (see Materials and Methods).

Although ligand binding prediction and ligand pocket prediction have been discussed since the previous round of RNA-Puzzles, the prediction of the RNA-ligand complex structure is difficult when no structural template is available to generate a RNA model that can be regarded as very close to the native RNA structure, and serve as a receptor to dock the ligand. Thus, in Puzzles 9 and 21, corresponding to aptamers for 5-hydroxytryptophan and guanidine, respectively (Supplemental Figs. S1, S8, S14), the accurate prediction of ligand contacts was still difficult. In Puzzle 21, the predicted structures within 5 Å RMSD from the crystal structure demonstrate ligand binding positions that are close to the correct binding pocket, and the model 1 of LORES prediction from the Das group correctly predicts the contacts with G7 nucleotide (Supplemental Fig. S14). In terms of a more incisive evaluation of accuracy, we still do not have a systematic program to analyze the quality of the ligand-binding pocket and the ligand-binding interactions for RNA. A major difficulty lies in the analysis of various interaction types between ligand and RNA. RNA-ligand binding remains an important aspect of both prediction and evaluation, and will be further discussed in studies summarizing results on later puzzles (e.g., PZ23, a Mango fluorogenic aptamer).

In addition to testing the predictive ability of new automated servers, Round IV provided information on the ability of expert and automated server predictors to identify homologous templates and ligand binding sites. To better understand the prediction performance, the human vs. server results are compared in Supplemental Figure S12, while radar plots discussing each puzzle are summarized in Supplemental Figures S15–S20. The diagrams generally show that polygon shapes are larger on diagrams with human predictions (the larger is the polygon, the better is the corresponding model). This is especially evident in Supplemental Figure S20 (Puzzle 21) which shows automatically predicted 3D models scored considerably worse than those predicted manually by the Das, Chen, and Adamiak-Szachniuk groups. In some cases, radar diagrams show clear winners for example, the Bujnicki model 1 in Puzzle 20. In other cases, radar plots highlight the relative weakness of each prediction approach.

In the Round IV prediction, we have included two twister sister ribozyme puzzles (PZ19, PZ20). Although belonging to the same sequence family, the two puzzles demonstrate significant structural differences. The Guanidine-III riboswitch in Puzzle 21 is also different from the previously solved Guanidine-II riboswitch structures (Huang et al. 2017a; Reiss and Strobel 2017) and Guanidine-I riboswitch (Reiss et al. 2017). Contrarily, some newly solved RNA structures could adopt the same structural topology as structures present in the database, for example, 5HTP riboswitch (5kpy) is similar to guanine riboswitch (4fe5), which was in fact its “parent” in an in vitro selection. It is crucial to understand how different sequences in the same family result in structural differences. And RNA-Puzzles aims to understand if such structural differences can be captured by existing prediction approaches. According to the results of PZ19 and PZ20, current prediction can already predict the difference between the two twister sister ribozymes. In our next RNA-Puzzles, manuscript, we will further discuss such structural variances in other RNA families, including the T-box (Puzzles 10, 26, 27, and 28).

Ions are important for stabilizing RNA structures by screening the electrostatic repulsions among the backbone phosphate groups, and divalent ions are known experimentally to be more efficient than monovalent ions (Woodson 2005). Post-blind-prediction studies by the Ding group showed that explicit modeling of divalent ions in all-atom DMD simulations could significantly enhance the sampling of near-native structures and improve the accuracy of predicted 3D structures (Supplemental Fig. S9C). Additional studies to incorporate explicit modeling of Mg2+ into the CG RNA simulations could further improve the prediction accuracy and efficiency of multiscale RNA modeling approaches.

Conclusions

The aims of the RNA-Puzzles collective work is to provide a snapshot and evaluation of the state-of-the-art in RNA 3D structure modeling. In Round IV, blind tests on ribozymes and riboswitch aptamers indicate that accurate prediction still requires an appropriate combination of different prediction strategies, including the fragment/structure library-based methods and the de novo prediction/force field-based simulation methods. Automated predictions from web servers still remain worse than human expert models for many cases. In future rounds, we speculate that further improvements may be achieved through machine-learning-based contacts prediction (Senior et al. 2020), integration of experimental data or better ranking of the predicted models. As discussed before (Miao and Westhof 2017), secondary structure (the knowledge of the helical Watson–Crick pairs) alone is far from sufficient to model RNA structures without prediction of non-WC base pairs. Non-WC base pairs, atomic clashes, and ligand binding site prediction still have room for improvement. Alternative structure comparison metrics, which are not sequence-length-dependent, may benefit the RNA-Puzzles assessment in later puzzles.

MATERIALS AND METHODS

Results assessment

At the first RNA-Puzzles community meeting, organized in 2016, it was decided to allow up to five structural models, which should be ordered according to the expected reliability for each prediction approach. Still, each research group may submit the results for multiple prediction approaches, while automated web servers are separated from human expert prediction.

We assess the prediction accuracy by comparing each prediction model with the reference structure of the lowest RMSD. Assessment metrics also include Deformation Index (Parisien et al. 2009), which is a metric to normalize RMSD with the sequence length; Clash score (Williams et al. 2018) shows the credibility of the structure in terms of atomic distance; Interaction Network Fidelity (Parisien et al. 2009) (INF) is an evaluation of the prediction accuracy of different types of interactions, which include Watson–Crick interaction, non-Watson–Crick interaction and base stacking. Considering some new prediction approaches may abstract RNA structure (such as the network model Kim et al. 2013), we use the Mean of Circular Quantities (Zok et al. 2014) to assess the accuracy in the torsion angle space.

Summary of structure prediction methods

In the Supplemental Notes section, we briefly introduce (1) the updates and improvement of the existing prediction approaches; (2) the development of novel prediction approaches; (3) the development of automated web servers; (4) the methods used and focus on the special treatments that were used in the prediction of the present ribozymes and riboswitches.

SUPPLEMENTAL MATERIAL

Supplemental material is available for this article.