RNA tertiary structure prediction with ModeRNA

Abstract

Noncoding RNAs perform important roles in the cell. As their function is tightly connected with structure, and as experimental methods are time-consuming and expensive, the field of RNA structure prediction is developing rapidly. Here, we present a detailed study on using the ModeRNA software. The tool uses the comparative modeling approach and can be applied when a structural template is available and an alignment of reasonable quality can be performed. We guide the reader through the entire process of modeling Escherichia coli tRNA^Thr in a conformation corresponding to the complex with an aminoacyl-tRNA synthetase (aaRS). We describe the choice of a template structure, preparation of input files, and explore three possible modeling strategies. In the end, we evaluate the resulting models using six alternative benchmarks. The ModeRNA software can be freely downloaded from http://iimcb.genesilico.pl/moderna/ under the conditions of the General Public License. It runs under LINUX, Windows and Mac OS. It is also available as a server at http://iimcb.genesilico.pl/modernaserver/. The models and the script to reproduce the study from this article are available at http://www.genesilico.pl/moderna/examples/.

INTRODUCTION

Noncoding RNAs play many important roles in the cell. Their functions are strongly connected with their three-dimensional structure [1]. There is a scarcity of detailed structural information compared to more than 10⁵ known ncRNA sequences. To close this gap, many tools and algorithms for RNA structure prediction from sequences have been developed [2]. Many of them use approaches similar to those used in protein structure prediction, a field with a longer history and with many success stories [3].

An example of such an RNA prediction tool is ModeRNA [4], inspired by SwissModel—one of the most popular programs for protein structure modeling [5]. The ModeRNA software is based on the comparative modeling approach, which requires a template: another molecule whose 3D structure is already known and which is expected to be similar to the one to be modeled. In particular, this condition usually holds for homologous molecules, therefore this approach is often called ‘homology modeling’. Additionally, a target–template sequence alignment indicating corresponding residues in both molecules must be provided. ModeRNA interprets the alignment nucleotide wise and assigns adequate operations e.g. copying residues that are invariant in both target and template, exchanging bases for mismatches, modeling nucleotide modifications, or inserting fragments for insertion–deletion (indel) regions. These operations can be performed in a single step; however, additional script-based editing of parts of the model may improve the quality of the resulting model. To facilitate such modeling cases, ModeRNA provides a robust set of commands, which can be used via a Python scripting interface.

In the original article describing ModeRNA, we have reported a large-scale benchmark of ModeRNA on a predefined set of targets, templates and alignments. In this article, we would like to highlight the entire modeling process starting from a target sequence. We present a detailed case study of using ModeRNA to model the tRNA^Thr molecule from Escherichia coli in a conformation corresponding to the one in complex with an aminoacyl-tRNA synthetase (aaRS), a protein enzyme that charges this tRNA with threonine (see Figure 1 for target sequence). The choice of this particular example has been motivated by the following rationale:

First, tRNA is an important family of RNAs that participate in protein biosynthesis as carriers of amino acid residues. As its biological function was identified already 50 years ago [6], it is now one of the most well-researched RNA families. Its first three-dimensional structures (tRNA^Phe from yeast) were solved in 1974 [7,8] and their number has grown to more than 170 as of 2011.

Second, the overall structure of tRNA is well conserved because they all need to fit into the same sites in the ribosome. Although tRNA sequences exhibit high divergence, the overall fold for all known tRNA structures is the same and thus the basic assumption for homology modeling (target and template must have a similar architecture) holds true.

Third, an experimental structure of E. coli tRNA^Thr is available so that the results of this case study can be assessed critically, i.e. the model can be compared to the experimentally solved structure.

Fourth, tRNAs are abundant in modified nucleotides that play an important role for both structural integrity and function. They modulate structural flexibility and rigidity, regulate thermostability and optimize codon–anticodon binding [9–11]. Hypomodified tRNAs appear to be less stable and thus prone to degradation. Modification of some positions is highly conserved in most tRNAs (Figure 2). For example, modifications occurring in the T-loop stabilize a tertiary interaction between the D-loop and the T-loop, which is crucial to form the L-shaped tRNA tertiary structure [10]. Dihydrouridine residues present in the D-loop destabilize the C3′-endo sugar conformation and promote conformational flexibility in the loop in order to facilitate the tertiary interaction between D-loop and T-loop. The impact of dihydrouridine is not limited to the modified residues themselves, but also conveyed to the 5′-neighboring residues [12]. Finally, highly conserved modifications in the anticodon loop restrict its dynamics and enhance the accuracy of codon binding [13]. Generally, the most frequent modifications are methylations and pseudouridylations. Since ModeRNA is to our knowledge the only method for comparative modeling of RNA that enables modeling of modified residues, we have put special attention on this aspect during modeling of E. coli tRNA^Thr.

Modeling of an RNA–protein complex requires also the structure of the protein component. The protein structure may be already known or would have to be modeled independently with protein-specific tools. Protein modeling is a mature field and has been reviewed extensively [14–19]. If models of both components of the complex are constructed in appropriate ‘bound’ conformations, rigid body docking can be used to construct the complex. Such analysis is beyond the scope of this article, therefore we focus only on modeling of the RNA.

Finding a template and preparing the alignment

We started the modeling experiment with only the 76-nt sequence of tRNA^Thr from E. coli (the target). As an input ModeRNA requires a template structure and a pair-wise alignment of the target and template sequences. The first step in modeling is to find an experimentally solved RNA structure, for which we have strong reason to believe that it shares a similar structure with the target molecule. The most obvious candidate with a potentially similar structure would be some evolutionary relative (a homolog), as molecules derived from a common ancestor are believed to preserve their structure despite accumulation of mutations [20]. To find a template for modeling the tRNA^Thr molecule from E. coli, we identified all members of the tRNA family in the Rfam database (RF00005) [21] that have experimentally determined structures in the PDB database [22] with a resolution better than 2.5 Å. This way, we limited the number of potential templates to 23. We used ModeRNA to extract the sequences of all these structures (using the get_sequence function) and examined PDB files for features that could cause problems during modeling (using the function examine_structure). Many chains contained water and ions (visible as a ‘_._._._._.’ tail in the sequence). These were removed by the clean_structure function.

In the next step, we wanted to find out which of these tRNAs align well with the target, i.e. have a high number of confidently identifiable homologous residues between the target and template. This typically requires high percent sequence identity and a low number of indels. We calculated alignments for each of the candidate templates using the Infernal program [23] and the covariance model for tRNAs from Rfam. Two templates (PDB codes 1EVV and 1EHZ) were found to produce an alignment with the highest sequence identity (0.59). As both structures represent the same tRNA^Phe from S. cerevisiae, we decided to use 1EHZ because it has a slightly better resolution (1.93 Å) than 1EVV (2.0 Å). According to the alignment, 16 out of 17 Watson–Crick and wobble base pairs in the four helical stems had to be replaced by other Watson–Crick/wobble pairs during modeling. It is worth noting that tRNA^Phe from S. cerevisiae includes posttranscriptionally modified nucleotides. Three base pairs exist in 1EHZ where a canonical residue is replaced by a modified residue in our modeling target or vice versa. In two cases, the modified base is m5C, whose additional methyl group does not disturb the Watson–Crick edge, and the third is m2G (residue 10) with the methyl group oriented in such a way that it does not disturb the formation of a canonical base pair with C25. The alignment contains no insertions or deletions.

We examined the template structure 1EHZ more closely. The dot–bracket secondary structure (function get_secstruc) showed that 1EHZ possesses the cloverleaf structure typical for tRNA (Figure 2). In the tertiary structure, some uncommon values of dihedral angles were found in residues 14–16 by the functions analyze_geometry and find_clashes. These residues are located in the DHU loop, which has a flexible conformation due to the presence of two dihydrouridines. According to observations, a few nonstandard dihedral angles can be found in many PDB structures, and especially in this position they are not very unusual. Finally, identity of the template structure and the template sequence from alignment was successfully verified with the match_template_with_alignment function. Thus, the structure can safely be used as a template.

One observation that required special attention was that the residues in the anticodon loop were folded compactly in a stacked conformation, typical for tRNA in an unbound state (Figure 2). However, we wanted to model a state corresponding to the tRNA complexed with the anticodon-binding domain of a cognate class II aaRS protein. To do so, we searched for template candidates for tRNA complexed with an aaRS protein. In the structure 1C0A, E. coli tRNA^Asp interacts with a type II aaRS [24]. However, the sequence identity of that tRNA with our modeling target was lower (0.46). A solution was to build a preliminary model of the whole tRNA^Thr using 1EHZ as the template and incorporate the anticodon loop from 1C0A. For comparison, we also wanted to prepare a model based on the template 1C0A alone.

Summarizing, we tested three modeling strategies: (i) we used the most similar template to tRNA^Phe from S. cerevisiae (1EHZ) in a different biological state, than desired; (ii) we used a more remotely related template tRNA^Asp from E. coli (1C0A) in the desired biological state; and (iii) we combined the two templates. The alignment for both templates is shown in Figure 1.

Model building from a template structure and alignment, with both the template and alignment prepared, can be done in a single step by ModeRNA. We used a Python script that loads the input data, creates the model and saves the model structure. It uses four simple commands: load_template, load_alignment, create_model and write_model. Alternatively, ModeRNA could have been run from a terminal providing the names of all input and output files (python moderna.py -t template_file.pdb -c template_chain -a alignment_file.fasta -o output_file.pdb). The input script for ModeRNA with all commands used for modeling along with the input files are available at the ModeRNA homepage (http://www.genesilico.pl/moderna/examples/).

During the modeling, ModeRNA analyzed the alignment and assigned particular operations to each position. A distinctive feature of ModeRNA is that it can perform these operations on modified nucleotides as well. In cases where residues were identical, the atomic coordinates were directly copied from the template. In case of substitutions, backbone and ribose atoms were copied and the proper base was added according to superposition of three atoms closest to the glycosidic bond. For each insertion and deletion, ModeRNA searched for a fragment in a library created from a nonredundant set of RNA structures, based on the geometry of residues to be connected (so called ‘anchor residues’). In case a fragment did not fit ideally between the anchors, the backbone coordinates were optimized using the FCCD Loop Closer algorithm [25], which closes small breaks in the backbone. A detailed description of these operations can be found in Rother M et al. [4].

In the modeling exercise based on the 1EHZ template, the alignment was 76 positions long, with 47 identical target–template residue pairs. The 29 remaining positions contained mismatches. In this modeling episode, we had no insertions or deletions. In some of these positions, modified nucleosides occurred and for these, additional chemical groups were added or removed.

MODELING OF MODIFIED NUCLEOSIDES

The target sequence and both template structures contained modified nucleosides mentioned above (many tRNA structures in the PDB are however modeled without modifications). For instance, the 1EHZ structure contained methylations (e.g. 2′-O-methylguanosine in position 34; Figure 2G) and more complicated chemical groups (e.g. wybutosine in position 37, Figure 2H). Dihydrouridine (residue 16 in 1EHZ, Figure 2F) is a unique modification, as it possesses a nonplanar ring. It is important to note that the modification pattern between the target and the template differs (Table 1).

In cases where modifications from the template matched those from the target, they were simply copied in the same manner as unmodified residues (e.g. the dihydrouridines in position 16 and 17). Some modifications needed to be changed into unmodified nucleosides (e.g. position 10). In such cases, the unmodified base was introduced by superposition of its three atoms nearest to the glycosidic bond onto the modified base to be replaced.

In the opposite situation, i.e. when an unmodified residue needed to be changed into a modified one (e.g. in position 20), one or more structural fragments containing the additional chemical groups were added to the base or ribose. When one modified residue needed to be replaced by another modified residue (e.g. in position 37), the first operation introduced a new unmodified residue, followed by an addition of the new modification. ModeRNA contains a set of 70 structural fragments that enable building 115 known modifications. Addition of the fragments is guided by a set of rules describing the atom triplets used for superposition and atoms to be added and removed.

In order to add dihydrouridine, the entire base needed to be exchanged due to the nonplanarity of the partially saturated ring. The same applied to pseudouridine, because the base ring needed to be rotated and connected with ribose via the C6 atom instead of N1. ModeRNA automatically identified and executed the proper rules for adding these modifications.

In addition to modeling modified residues, ModeRNA features a few other operations for working with modified nucleosides. First, modified nucleosides in an RNA structure can be detected (using the find_modifications function). The recognition is based on atomic coordinates and therefore independent of the nomenclature in the PDB file. Modifications can also be directly added to or removed from a specified residue (functions add_modification and remove_modification). It is also possible to remove all modifications in one step (function remove_all_modifications). ModeRNA contains a modification nomenclature implemented in the MODOMICS database [26]. For each modification, it stores a full name, a one-letter and a few-letter abbreviation, a common PDB residue name and a numerical code. In the alignment, the one-letter abbreviations are used, if available. Because there are more modifications than reasonably usable ASCII characters, the numerical code can be used alternatively, e.g. 001U for pseudouridine.

ADVANCED MODELING OF THE ANTICODON LOOP

Constructing a complete model from a template and alignment is the simplest way to obtain a model. However, in some cases it is not sufficient and some additional editing is required. ModeRNA provides many commands that enable changes on different levels: the entire molecule, a particular region and a single residue (see http://www.genesilico.pl/moderna/commands/).

In case of E. coli tRNA^Thr, we used the model built on the 1EHZ template and grafted the anticodon loop modeled on another template (1C0A). We built three models in addition to the one obtained previously. First, we added a fragment containing just the anticodon loop (9 residues), then one with almost the entire anticodon arm plus one terminal base pair (15 residues). Finally, we tried a fragment with the anticodon arm and a few additional nucleotides from an adjacent helical stem (19 residues).

In order to insert the three anticodon arm fragments, we cut out parts of 1C0A and saved them as separate PDB files. To extract residues from 1C0A, we used ModeRNA instead of manually copying residues between PDB files. We did that, because ModeRNA guarantees avoiding accidental formatting mistakes by using a unified nomenclature and numbering of atoms and residues. The structure 1C0A was loaded, unnecessary residues were deleted and the remaining structure was saved (the load_model, delete_residue, write_model commands were used). Thereby, three fragments were generated: U-625 to G-645 (UACCUGCCUQUCACGCAGGGG), C-627 to G-643 (CCUGCCUQUCACGCAGG) and G-630 to C-640 (GCCUQUCACGC). The two terminal residues of each fragment were used as anchor residues for superposition with the corresponding residues from the model.

The prepared PDB files were used to create fragment objects with ModeRNA (create_fragment). We specified two anchor residues from the model to guide the insertion and a new sequence for the anticodon fragment. During the insertion process, the fragments were superimposed using the anchor residues from fragment and model (the fragment was mobile and the model stayed in a fixed position). In particular the atoms O3′, C3′, C4′, C1′ and N1 or N9 on the 5′-end and C5′, C4′, C3′, C1′, N1 or N9, and O5′ on 3′-end were used for superposition. All residues present in the model between the anchors specified during fragment creation were removed and new residues from the fragment were added.

Eventually, the insertion of a fragment with a nonideal match with the anchor residues resulted in minor gaps in the backbone of the model. However, ModeRNA managed to rebuild backbone between such residues (function fix_backbone).

USAGE OF THE MODERNA SERVER

To make using ModeRNA simpler, we implemented a web interface available at http://iimcb.genesilico.pl/modernaserver/ [27], which provides many of the functions described above. First, a template and alignment can be obtained according to the procedure described above using the ‘Find template’ tab. Thus, it is possible to use just a sequence of E. coli tRNA^Thr to find possible templates (in this case 1EVV and 1EHZ in the first position and 1C0A in the fifth) and subsequently build a model. The server also provides alignments prepared using Infernal and a covariance model from Rfam. Moreover, the geometry and secondary structure of the chosen template and model can be analyzed by the server, and PDB structures can be cleaned from ions and water molecules using the ‘Analyze structure’ tab. Modeling based on a template and alignment can be carried out using the ‘Build model’ tab. In addition, a series of nucleotide exchanges including modified nucleotides can be performed in a straightforward way by providing a template structure and a target sequence.

More advanced manipulations of the RNA structure, e.g. extending a helix, adding a particular structural fragment, and searching for fragments with a given secondary structure require the standalone version of the software and are not yet possible with the server.

Summarizing, all steps described in this manuscript, except the replacement of the anticodon loop, can be done online by the ModeRNA server, without installing ModeRNA locally. Thus, the ModeRNA server is a convenient way to get started in many straightforward cases of homology modeling.

EVALUATION OF TRNA MODELS

As a result of the modeling procedure described above, we obtained five complete tRNA models (Figure 3). One was built without further editing on the 1EHZ template, one was built using the same procedure on the 1C0A template, and three were built based on the 1EHZ template, followed by replacement of an anticodon hairpin loop model based on 1C0A. Despite the fact that the 1EHZ template had the highest sequence identity, its anticodon loop was in a stacked conformation. Conversely, 1C0A had lower sequence identity but its structure was changed due to its interaction with class II tRNA synthetase, so it was in a similar state as we would like to obtain for our model.

In order to evaluate the models, we checked simple structural features like unusual geometry and interatomic clashes (functions find_clashes, analyze_geometry) (Table 2). The test revealed that all models exhibited minor geometrical problems. Nonetheless, after manual inspection in PyMOL [28], we saw that none of them were very serious, all except one originated from unusual lengths of the glycosidic bond, which were 0.01 Å shorter than the allowed boundary. Moreover, none of the models seemed to have accumulated a particularly high level of errors during modeling. The model built with just one command on the 1EHZ template contained the lowest number of geometric distortions and the one built solely on 1C0A had the most problems. These measures provided information on how far the chemical structure and details of stereochemistry were from the ideal values in the model, but not whether the macromolecular structure was modeled correctly.

As the structure of E. coli tRNA^Thr interacting with its cognate aaRS enzyme has been experimentally solved [29] and is available as a PDB entry 1QF6 (chain B), we were able to assess the real accuracy of our models. To do so, we applied six different benchmarks of structural similarity (Table 2).

The root mean square deviation (RMSD) is the most standard measure of similarity between two 3D structures, and is still the one most frequently found in the literature. Generally, the smaller the RMSD, the more similar the compared structures (with 0 corresponding to identity); however, this relationship may not hold for structures with very different sizes, or those exhibiting conformational changes. We calculated RMSDs between the experimentally solved structure and our five models in two ways: using all heavy atoms and only two backbone atoms (P and C4′). The template modeling (TM) score is a normalized measure of the overall structure similarity that does not depend on the structure size. It ranges from 1 (identity) to 0 (no similarity) and follows an extreme value distribution; for protein structures values above 0.5 typically indicate similar structures, while values below 0.5 are characteristic of dissimilar structures [30]. Global Distance Test—Total Score (GDT–TS) is another measure whose value ranges between 0 (no identity) to 1 (100% identity) and indicates the number of corresponding atom pairs between the compared structures found below four different distance thresholds: 1, 2, 4 and 8 Å divided by four times the total number of atoms [14]. We also applied the deformation index (DI) and deformation profile (DP) metrics [31]. The DI indicates how well the base pairing and stacking interactions were modeled, with 1 being the ideal value. On the contrary, DP highlights the dissimilarity between two structures and is calculated by superimposing each pair of corresponding residues and generating a matrix of per-residue RMSD values (hence, the lower the DP, the more similar the structures). Here, we present the average of all values from the DP matrix, which depends on the molecule size.

As indicated earlier, we intended to generate a model of E. coli tRNA^Thr in a bound conformation; therefore, we used the E. coli tRNA^Thr–aaRS complex as the benchmark, and in this case ‘model accuracy’ indicates the ability of the predicted structure to approximate this particular conformation solved under specific experimental conditions. The analysis using five of the six measures (RMSD, backbone RMSD, TM, GDT–TS and average DP) showed that insertions of the anticodon loop modeled in the bound conformation improved the accuracy of E. coli tRNA^Thr models over the model built on 1EHZ alone (anticodon loop in unbound conformation), and that the longest insertion gave the best result. The model based on the template 1C0A (the one with lower sequence similarity, but in the most relevant functional state) turned out to be the most accurate. The DI reached the best value for the model constructed with the 1EHZ template alone (0.82), the same DI score was achieved for the model with the smallest insertion, and slightly lower for the model built with the 1C0A template (0.79). Summarizing, the models based on 1EHZ (template with the highest sequence identity to the target sequence) had the best local quality and the worst global accuracy, and the models built on 1C0A (the template in the desired functional state) had the worst local quality, but were globally most accurate.

Having inspected the models superimposed on the crystal structure 1QF6 (Figure 3), we found that the models with the replaced anticodon loop were more similar to the crystal structure than the model built on the 1EHZ template alone. In the models with the anticodon loop replaced, the anticodon bases were displaced from the corresponding residues in the 1QF6 structure. In the model built on the 1C0A template, the anticodon bases were closer to those in 1QF6, caused by a different orientation of the anticodon stem in 1C0A compared to 1EHZ. Another region that had a large impact on the quality of models was the 3′-CCA tail, which is known to crystallize in different conformations [32]. Another problematic region was the D loop. It is also known as a flexible region in tRNA because of two dihydrouridine residues, as described above.

We also compared the backbone RMSD between the native structure 1QF6 and the templates. In both cases, the values were identical: 4.37 for the 1EHZ template and 2.60 for the 1C0A template. Thus, the modeling process did not affect the overall similarity of the backbone between the target and template molecules. The values hint at the structural variety among tRNA structures. In an analysis of 99 tRNA structures, we found that tRNAs with an identical sequence differ by up to 4.2 Å [4].

DISCUSSION

This article presents a case study of RNA 3D structure prediction using the comparative modeling approach and the ModeRNA software. We built five models of E. coli tRNA^Thr in a conformation bound to an aaRS, starting from a target sequence. The choice of the template was not obvious and a few alternatives have been compared.

Two alternative templates were selected based on the known fact that the structural architecture of tRNA is strongly conserved. The target and the two templates originate from different phylogenetic domains and are specific for three different amino acids. Their sequences thus diverged already around the time of fixation of the genetic code. The conservation of the L-shaped architecture in all experimentally determined tRNA structures suggests that the same structure is likely to be found in their homologs. Although it is known that there is considerable structural variability in tRNA [33], it occurs mainly in the anticodon loop and CCA tail, two regions for which modeling was problematic.

Further support for the feasibility of the chosen templates comes from the alignments. The target sequence had the same length as the 1EHZ template, and a single-residue deletion compared to 1C0A. Recently, Capriotti et al. [34] described the ‘twilight zone of sequence identity’ at which homology modeling of RNA becomes difficult and it was reported to be around 0.42–0.71. Sequence identities of the templates 1EHZ and 1C0A to our target were 0.59 and 0.46, respectively. These numbers have been calculated with modified nucleotide residues present (which were not considered by Capriotti et al.). By default, D and U were considered as nonidentical, even though dihydrouridine is a derivative of uridine. ‘Demodification’ of the sequences (replacement of all Ds by Us, etc.) makes their per cent identity cross the level of 0.72, considered as the ‘safe zone’ for RNA homology modeling [34].

As exemplified by the modeling of E. coli tRNA^Thr, the sequence identity to the target molecule should not be the only criterion for the selection of a template. In comparative modeling, the user is responsible for defining the functional state desired for the target molecule. Macromolecules such as proteins or RNAs exhibit local and/or global conformational changes that depend on their environment, in particular interactions with other molecules. The binding of a particular molecule type often induces a conformational transition common to many members of the same family of macromolecules. Therefore, the conformational variability should be kept in mind and only those templates, whose conformation is consistent with the functional state to be modeled, should be selected.

Five tRNA models obtained in the modeling described above had an all-atom RMSD to the experimentally solved structure in the range of 3.4–5.1 Å. When only RMSD is concerned as structure similarity measure, using the 1C0A template in that was in the protein-bound state gave the best result, despite its lower sequence identity to the target, compared to the 1EHZ template. Five of the six benchmarks showed that a model based on the 1C0A is the most similar to the crystal structure of E. coli tRNA^Thr complexed with aaRS. However, according to the DI, the model built on 1EHZ is more accurate than the one on 1C0A, i.e. the base interactions are more similar to those in the experimentally determined structure. ‘Grafting’ an anticodon loop from the 1C0A-based model onto the 1EHZ-based model improved the global accuracy, but led to a decrease in the local quality. It is clear that building a model accurate both on global and local scale using a comparative approach is difficult. Further, the combined use of different templates for different regions of the target molecule significantly increases the time and complexity of the modeling exercise, compared to the ‘one step’ modeling based on target sequence alone.

In order to check the performance of ModeRNA in its simplest, maximally automated mode, it was tested on a set of 99 known tRNA structures using each other as templates [4]. The conformational heterogeneity of these structures was considerable, with the average P and C4′ RMSD reaching 4.9 Å. Hence, a dataset of all-versus-all target–template pairs contained many very similar structures, some of them in a completely different conformation than the crystal structure 1QF6. Besides, there were many tRNAs with a large insertion (a variable loop), which could not be modeled accurately on templates lacking a corresponding loop. The resulting 9675 tRNA models exhibited average all-atom RMSD values of 5.6 Å, 5.2 Å for P & C4′-atom RMSD, 0.5 for GDT-TS, 0.62 for DI, and 13.82 for DP. Comparison of these results to the values obtained for E. coli tRNA^Thr showed that by carefully choosing the template and modeling strategy, it is possible to build reliable models. Searching templates in the Rfam database by using covariance models is applicable to other families of ncRNA as well.

Models for tRNAs have been built with other programs as well. Several models for yeast tRNA^Phe have been generated, clustered and ranked with Nucleic Acid Simulation Tool (NAST), which is based on a coarse-grained knowledge-based potential. The three resulting clusters had an average all-atom RMSD to the native structure of 8.0, 13.6 and 15.6 Å and a GDT-TS of 0.2, 0.08 and 0.07, respectively [35]. Flores and Altman have also built a model for tRNA^Phe using constraints from a limited number of tertiary contacts, stacking interactions and NMR data which resulted in a P-only RMSD of 9.6 Å compared to the native structure [36]. Lavender et al. [37] have built a tRNA^Asp model with a P-only RMSD of 6.2 Å. Ding and Dokholan have modeled tRNA^Phe using the DMD software and secondary structure restraints with a RMSD 7.2 Å compared to the native structure [38], and the same authors have built a tRNA model with a RMSD of 4.0 Å using Molecular Dynamics alone [39]. More recently, Cao and Chen reported modeling of a tRNA with an all-atom RMSD of 4.2 Å to the native structure using a combination of knowledge-based potential and Molecular Dynamics simulation [40]. These values require two comments: first, one needs to be very careful when comparing RMSD values. The above paragraph alone contains four different modes of calculation that cannot be compared directly. Second, most of the methods mentioned above start from an unfolded RNA sequence and simulate the folding process. ModeRNA uses a template, which of course facilitates the model building, but taking into account the template identification and preparation of a target–template alignment, comparative modeling presents its own challenges. While it has been demonstrated that methods that simulate folding with the use of experimental restraints are capable of building accurate tRNA models, the ones built by ModeRNA are at least competitive. One decisive advantage of the homology modeling approach is calculation time. Where knowledge-based potentials often require hours to fold a tRNA-sized structure on a single processor, and full-atom Molecular Dynamics much longer, ModeRNA calculates such a model within seconds.

ModeRNA allows editing models explicitly, e.g. by fragment insertion, and by adding or removing helical regions and single base pairs. The recombination of different structural parts may result in backbone breaks and ModeRNA can reconstruct the backbone in such regions (function fix_backbone). This functionality was successfully used to build a model of the Azoarcus group I intron that compares well (4.3 Å versus 4.4 Å) with the one generated using the RNABuilder software by Flores et al. [41]. An initial model for the Azoarcus intron has been created by Rangan et al. [42]. Homology modeling was also used to construct a model of the 30S ribosomal subunit from E. coli [43]. Compared to the crystal structure (PDB-ID 3R8N), the model reaches a RMSD of 3.3 Å.

Comparative models (and in fact any models built with any method) may contain geometrical distortions and other steric problems. The identification of such flaws (e.g. by the use of ModeRNA) indicates that the model should undergo further optimization, for instance local optimization by programs for energy minimization like MMTK [44] or OpenMM Zephyr [45]. In the absence of the experimentally determined structure to be used as a reference for model quality assessment (i.e. in ‘real life’ cases of comparative modeling), the accuracy of the theoretical model can be predicted using computational tools, such as the recently developed knowledge-based potentials RASP [20] or KB [46]. It is also advisable to check the validity of the model against the available experimental data, e.g. with FILTREST3D [47] or with NAST [35].

The further use of RNA 3D structure models, e.g. to model interactions with other molecules, requires the use of other bioinformatics tools. In the example described in this article, the next step would require the acquisition of a model of the protein partner in the appropriate functional state, which can be achieved by an analogous modeling protocol, with protein-specific tools such as SwissModel [5] or Modeller [48]. The assembly of a complex can be guided by homology (e.g. by superposition onto another related complex) or de novo, by protein–RNA docking, e.g. with HADDOCK [49]. Protein–RNA docking is an emerging field, and currently no standard protocols exist, especially for analyzing RNA molecules with modified residues. For instance, HADDOCK cannot automatically process modified residues, and such analysis would require ‘demodification’ of the RNA with ModeRNA prior to docking. A detailed description of docking is however beyond the scope of this article.

The ModeRNA script and input files to reproduce all models presented in this study are available on the ModeRNA website (http://genesilico.pl/moderna/examples/). We believe they are useful as a starting point for developing further RNA comparative modeling experiments, on tRNA or other targets.

FUNDING

The German Academic Exchange Service (D/09/42768 to K.R.); the Polish Ministry of Science and Higher Education (N N301 035 539 to T.P.); the European Research Council (RNA+P=123D to J.M.B.); Foundation for Polish Science (‘Ideas for Poland’ fellowship to J.M.B.).

References