Nucleic acid structure prediction

Nucleic acid structure prediction is a computational method to determine secondary and tertiary nucleic acid structure from its sequence. Secondary structure can be predicted from one or several nucleic acid sequences. Tertiary structure can be predicted from the sequence, or by comparative modeling (when the structure of a homologous sequence is known).

The problem of predicting nucleic acid secondary structure is dependent mainly on base pairing and base stacking interactions; many molecules have several possible three-dimensional structures, so predicting these structures remains out of reach unless obvious sequence and functional similarity to a known class of nucleic acid molecules, such as transfer RNA (tRNA) or microRNA (miRNA), is observed. Many secondary structure prediction methods rely on variations of dynamic programming and therefore are unable to efficiently identify pseudoknots.

While the methods are similar, there are slight differences in the approaches to RNA and DNA structure prediction. In vivo, DNA structures are more likely to be duplexes with full complementarity between two strands, while RNA structures are more likely to fold into complex secondary and tertiary structures such as in the ribosome, spliceosome, or transfer RNA. This is partly because the extra oxygen in RNA increases the propensity for hydrogen bonding in the nucleic acid backbone. The energy parameters are also different for the two nucleic acids. The structure prediction methods can follow a completely theoretical approach, or a hybrid one incorporating experimental data.^[1][2]

A common problem for researchers working with RNA is to determine the three-dimensional structure of the molecule given only a nucleic acid sequence. However, in the case of RNA much of the final structure is determined by the secondary structure or intra-molecular base pairing interactions of the molecule. This is shown by the high conservation of base pairings across diverse species.

Secondary structure of small RNA molecules is largely determined by strong, local interactions such as hydrogen bonds and base stacking. Summing the free energy for such interactions should provide an approximation for the stability of a given structure. To predict the folding free energy of a given secondary structure, an empirical nearest-neighbor model is used. In the nearest neighbor model the free energy change for each motif depends on the sequence of the motif and of its closest base-pairs. (The "neighbor" aspect allows for describing the base stacking energy.) The energy of fragments are added together to get the total free energy. The structure with the minimal free energy (MFE) is the most stable.^[3]

Examples of RNA models include:

• Turner 1999^[4] and Turner 2004.^[5]
• Andronescu 2007 parameters for RNA with any of the four canonical bases. This is derived from earlier wet lab experiments and from other forms of RNA structure data.^[6][7]
• Langdon 2018 parameters for RNA with any of the four canonical bases.^[8]

As of 2020, most software packages default to Turner 2004. A 2020 benchmark shows that Turner 2004 is the worst of the four for RNA-RNA interaction tasks, with the best (Andronescu 2007) outperforming it by 5 pairs. It is known that Langdon 2018 is the best for structure prediction.^[9]

The simplest way to find the MFE structure would be to generate all possible structures and calculate the free energy for it, but the number of possible structures for a sequence increases exponentially with the length of RNA: number of secondary structures = (1,8)N, N- number of nucleotides.^[10] For longer molecules, the number of possible secondary structures is huge: a sequence of 100 nucleotides has more than 10²⁵ possible secondary structures. This calls for smarter methods.^[3]

Most popular methods for predicting RNA and DNA's secondary structure involve dynamic programming.^[11][12] One of the early attempts at predicting RNA secondary structure was made by Ruth Nussinov and co-workers who developed a dynamic programming-based algorithm that maximized the length and number of a series of "blocks" (polynucleotide chains).^[11] Each "block" required at least two nucleotides, which reduced the algorithm's storage requirements over single base-matching approaches.^[11] Nussinov et al. later published an adapted approach with improved performance that increased the RNA size limit to ~1,000 bases by folding increasingly sized subsections while storing the results of prior folds, now known as the Nussinov algorithm.^[12] In 1981, Michael Zuker and Patrick Stiegler proposed a refined approach with performance comparable to Nussinov et al.'s solution but with the additional ability to find also find "suboptimal" secondary structures.^[13]

Dynamic programming algorithms provide a means to implicitly check all variants of possible RNA secondary structures without explicitly generating the structures. First, the lowest conformational free energy is determined for each possible sequence fragment starting with the shortest fragments and then for longer fragments. For longer fragments, recursion on the optimal free energy changes determined for shorter sequences speeds the determination of the lowest folding free energy. Once the lowest free energy of the complete sequence is calculated, the exact structure of RNA molecule is determined.^[3]

Dynamic programming algorithms are commonly used to detect base pairing patterns that are "well-nested", that is, form hydrogen bonds only to bases that do not overlap one another in sequence position. Secondary structures that fall into this category include double helices, stem-loops, and variants of the "cloverleaf" pattern found in transfer RNA molecules. These methods rely on pre-calculated parameters which estimate the free energy associated with certain types of base-pairing interactions, including Watson-Crick and Hoogsteen base pairs. Depending on the complexity of the method, single base pairs may be considered, and short two- or three-base segments, to incorporate the effects of base stacking. This method cannot identify pseudoknots, which are not well nested, without substantial algorithmic modifications that are computationally very costly.^[14]

One of the issues when predicting RNA secondary structure is that the standard free energy minimization and statistical sampling methods can not find pseudoknots.^[4] The major problem is that the usual dynamic programing algorithms, when predicting secondary structure, consider only the interactions between the closest nucleotides, while pseudoknotted structures are formed due to interactions between distant nucleotides. Rivas and Eddy published a dynamic programming algorithm for predicting pseudoknots.^[14] However, this dynamic programming algorithm is very slow. The standard dynamic programming algorithm for free energy minimization scales O(N³) in time (N is the number of nucleotides in the sequence), while the Rivas and Eddy algorithm scales O(N⁶) in time. This has prompted several researchers to implement versions of the algorithm that restrict classes of pseudoknots, resulting in performance gains. For example, pknotsRG tool includes only the class of simple recursive pseudoknots and scales O(N4) in time.^[15]

The general problem of pseudoknot prediction has been shown to be NP-complete.^[16]

The accuracy of RNA secondary structure prediction from one sequence by free energy minimization is limited by several factors:

1. The free energy value list in a nearest neighbor model may be incomplete or inaccurate.
2. Not all known RNA folds in such a way as to conform with the thermodynamic minimum.
3. Some RNA sequences have more than one biologically active conformation (i.e., riboswitches).
4. Modified bases are often assumed to fold in the same way as the unmodified version, but they actually have differences ranging from subtle (free energies changes on the order of 0.5 kcal/mol) to major (modification physcally blocks hydrogen bond from forming).^[17]

For this reason, the ability to predict structures which have similar low free energy can provide significant information. Such structures are termed suboptimal structures. MFOLD is one program that generates suboptimal structures.^[18]

A collection of structures can be represented by a "dot plot", a visualization of a square matrix containing the frequency of two bases being paired. This is essentially a type of heat map.

A collection of possible structures forms an ensemble. By sampling from the ensemble according to the Boltzmann distribution (as exemplified by the program SFOLD), many possible structures can be returned.^[19][20]

The states in an Boltzmann ensemble can also be represented by a "dot plot", with magnitudes being the probability of finding any two base being paired in the ensemble. These magnitudes implicitly represent thermodynamic data as well.

Multiple evolutionarily related RNA sequences share patterns of covariation between sites. If two sites are often seen changing in sync, there is a high possibility that there is a structurally required hydrogen bond between those positions. This offers an additional source of information for what the true structure could look like.^[16]

In general, the problem of alignment and consensus structure prediction are closely related. A few different approaches to the prediction of consensus structures can be distinguished:^[21][22]

1. Folding of alignment
2. Alignment of singular predicted structures, in some cases jointly with sequence alignment
3. Simultaneous alignment of sequence and ensembles of predicted structures

Each of these steps can also be mixed and matched in an iterative fashion.^[22]

A practical heuristic approach is to use multiple sequence alignment tools to produce an alignment of several RNA sequences, to find consensus sequence and then fold it. Covariation can be extracted from a multiple sequence alignment of multiple homologous RNA sequences with related but dissimilar sequences. The quality of the alignment determines the accuracy of the consensus structure model. Consensus sequences are folded using various approaches similarly as in individual structure prediction problem:

• Pfold implements SCFGs. Covariation is converted into a "dot plot" of base pair probabilities, which is considered on top of the usual single-sequence version of Pfold.^[23]
• RNAalifold from the Vienna suite is essentially an alignment-based variant of RNAfold from the same suite. Both use a thermodynamic folding approach and can run in two modes: a MFE (minimum free energy) mode, which returns a single structure, and a partitioning mode, which gives several possible structures.^[24] The score of a consensus structure is a combination of average thermodynamic energies and extra points derived from co-variation.^[22]
  • A linear-time variant is found as LinearAlifold.^[25]
• ILM (iterated loop matching) uses combination of thermodynamics and mutual information (covariance) content scores. Like suggested in its name, it is based on Nussinov style loop-matching between predicted structures. In each iteration, a "best" helix/stem is found in the predicted secondary structure and the rest of the structure is folded again. This allows the prediction of pseudoknoted structures.^[26]

Evolution frequently preserves functional RNA structure better than RNA sequence.^[24] Hence, a common biological problem is to infer a common structure for two or more highly diverged but homologous RNA sequences. Sequence alignments become unsuitable and do not help to improve the accuracy of structure prediction when sequence similarity of two sequences is less than 50%. This calls for an approach that combines alignment and secondary structure prediction.^[22]

The obvious approach to combine alignment and secondary structure prediction would be to fold the sequences using single-sequence structure prediction methods and align the resulting structures. Alignment of two structures is usually done by parsing the two structures into trees and then trying to find the minimum edit distance between the two trees.^[27] This method runs fairly quickly even when a joint sequence-structure alignment is performed and gives results better than starting from a fixed, unchanging, potentially incorrect alignment.^[28] It

The use of singular structures and tree comparison still has it weaknesses. It is well-known that the true (physiological) structure is often not the minimum free energy structure. It is also found that tree-editing has limited ability to repair mispairings. These problems call for yet another type of approach, one that compares whole ensembles of possible structures.

There is an existing representation of ensembles of possible structures already: the dot plot, which remains useful here. To align two sequences while considering ensemble structures, each sequence first has its base-pair probability dot plot predicted. The problem then becomes joint alignment of sequences and dot plots; the prototypic algorithm for this problem is the Sankoff algorithm,^[29] basically a merger of sequence alignment and Nussinov (maximal-pairing)^[11] folding dynamic programming method.^[30] The original Sankoff algorithm is a theoretical exercise because it requires extreme computational resources (O^(n3m) in time, and O^(n2m) in space, where n is the sequence length and m is the number of sequences).

Some notable attempts at implementing restricted versions of Sankoff's algorithm are Foldalign,^[31][32] Dynalign,^[33][34] PMmulti/PMcomp,^[30] Stemloc,^[35] Murlet,^[36] and LocARNA-P.^[37] These implementations are practical to run, but impose additional limits on properties such as the maximal length of alignment or variants of possible consensus structures. For example, Foldalign focuses on local alignments and restricts the possible length of the sequences alignment.

A different approach to joint alignment is used in CARNA to allow for more possible arrangements of base pairs such as pseudoknots or multiple stable structures. This is a MAX-SNP-hard problem, but constraint programming and hard limits keep the runtime practical.^[38]

No matter which joint alignment method is used, they only deal with two sequences (or two pre-aligned groups of sequences) at a time. This is solved by ordinary iterative alignment with a guide tree.

Once secondary structure of RNA is known, the next challenge is to predict tertiary structure. The biggest problem is to determine the structure of regions between double stranded helical regions. Also RNA molecules often contain posttranscriptionally modified nucleosides, which because of new possible non-canonical interactions, cause a lot of troubles for tertiary structure prediction.^{[39][40][41][42]}

In the comparative approach, a related known structure is extracted from a database (e.g. the PDB for 3D structures). This structure is used as a template, onto which the sequence in question is grafted and adjusted.^[43]

In the de novo or ab initio approach, no known 3D structures are used. The first step would be to predict the secondary structure using one of the approaches described above.^[44] The secondary structure is then converted into fragments of 3D structures and joined together to form starting points for further refinement:

• The more obvious method is to represent the temporary 3D structure using an all-atom representation containing the coordinates of every atom in 3D space.^[45] The free energy of such a structure can be improved using existing physics-based methods such as molecular dynamics^[46] or an RNA-specific random sampling of the conformational landscape^[47] followed by screening with a statistical potential for scoring.^[48]
• All-atom structures take a lot of memory to store and a lot of computation power to improve. A coarse-grained representation can be used first to bring the structure closer to the true shape^[49] before entering the high-resolution (all-atom) refinement stage, reducing the amount of work needed to be done in the high-resolution step.^[50]

As with the problem of protein folding, machine learning (ML) has also been used to find and make use of higher-order relationships between the RNA sequence and the 2D or 3D structure. They outperform traditional methods in predicting the global RNA fold, but traditional (non-ML) methods still have an advantage in modeling intramolecular interactions and ligand binding sites as of 2024.^[51]

• RNA
• RNA structure
• Non-coding RNA
• List of RNA structure prediction software
• Comparison of nucleic acid simulation software
• Comparison of software for molecular mechanics modeling

• Baker D, Sali A (2001). "Protein structure prediction and structural genomics". Science. 294 (5540): 93–6. Bibcode:2001Sci...294...93B. doi:10.1126/science.1065659. PMID 11588250. S2CID 7193705.
• Chiu D.K.; Kolodziejczak T. (1991). "Inferring consensus structure from nucleic acid sequences". Comput. Appl. Biosci. 7 (3): 347–352. doi:10.1093/bioinformatics/7.3.347. PMID 1913217.
• Do CB, Woods DA, Batzoglou S (2006). "CONTRAfold: RNA secondary structure prediction without physics-based models". Bioinformatics. 22 (14): e90–8. doi:10.1093/bioinformatics/btl246. PMID 16873527.
• Gutell R.R.; et al. (1992). "Identifying constraints on the higher-order structure of RNA: continued development and application of comparative sequence analysis methods". Nucleic Acids Res. 20 (21): 5785–5795. doi:10.1093/nar/20.21.5785. PMC 334417. PMID 1454539.
• Leontis NB, Lescoute A, Westhof E (2006). "The building blocks and motifs of RNA architecture". Curr Opin Struct Biol. 16 (3): 279–87. doi:10.1016/j.sbi.2006.05.009. PMC 4857889. PMID 16713707.
• Lindgreen S, Gardner PP, Krogh A (2006). "Measuring covariation in RNA alignments: physical realism improves information measures". Bioinformatics. 22 (24): 2988–95. doi:10.1093/bioinformatics/btl514. PMID 17038338.
• Lorenz, Ronny (2014). RNA secondary structure thermodynamics and kinetics. Vienna, Austria: University of Vienna, Dissertation.
• Macke T, Case D (1998). "Modeling Unusual Nucleic Acid Structures". Modeling unusual nucleic acid structures. In Molecular Modeling of Nucleic Acids. Edited by Leontes N, SantaLucia JJ. Washington, DC. ACS Symposium Series. Vol. 682. American Chemical Society. pp. 379–393. doi:10.1021/bk-1998-0682.ch024. ISBN 978-0-8412-3541-0.
• Major F (2003). "Building three-dimensional ribonucleic acid structures". Computing in Science & Engineering. 2003 (5): 44–53. Bibcode:2003CSE.....5e..44M. doi:10.1109/MCISE.2003.1225860. S2CID 17627934.
• Massire C, Westhof E. "MANIP: an interactive tool for modelling RNA". J Mol Graph Model. 1998 (16): 197–205, 255–257.
• Parisien M.; Major F. (2008). "The MC-Fold and MC-Sym pipeline infers RNA structure from sequence data". Nature. 452 (7183): 51–55. Bibcode:2008Natur.452...51P. doi:10.1038/nature06684. PMID 18322526. S2CID 4415777.
• Tuzet, H. & Perriquet, O., 2004. CARNAC: folding families of related RNAs. Nucleic Acids Research, 32(Web Server issue), W142-145.
• Touzet H (2007). "Comparative Analysis of RNA Genes". Comparative Genomics. Methods in Molecular Biology. Vol. 395. pp. 465–474. doi:10.1007/978-1-59745-514-5_29. ISBN 978-1-58829-693-1. PMID 17993692. S2CID 244726.
• Yingling YG, Shapiro BA (2006). "The prediction of the wild-type telomerase RNA pseudoknot structure and the pivotal role of the bulge in its formation". J Mol Graph Model. 25 (2): 261–274. Bibcode:2006JMGM...25..261Y. doi:10.1016/j.jmgm.2006.01.003. PMID 16481205.
• Zwieb C, Muller F (1997). "Three-dimensional comparative modeling of RNA". Nucleic Acids Symp Ser. 36 (36): 69–71. PMID 9478210.
• ModeRNA: A program for comparative RNA modeling