Pseudoknots, RNA structures in which the base pairs are not fully nested, are biologically important but cannot be handled by many computational procedures. Pseudoknots began as a theoretical prediction (Studnicka et al. 1978; Waterman 1978) but were found in viral RNAs a few years later (Rietveld et al. 1982; Pleij et al. 1985). They are now known to be critical for the structure and function of many RNAs, and evolutionarily conserved pseudoknots are involved in processes that include ribosomal frameshifting, self-cleavage, and self-splicing (Pleij 1994; Hilbers et al. 1998; Staple and Butcher 2005; Rødland 2006). The pseudoknotted region can be substantial: for example, at least 36% (8 out of 22) of the base pairs in the hepatitis delta virus (HDV) ribozyme structure must be removed to eliminate pseudoknots. In longer structures, pseudoknots often make up a relatively small but critical part of the molecule. For example, in the Escherichia coli 16S rRNA structure (Cannone et al. 2002), as few as 1.88% (9 out of 478) of the base pairs can be removed to produce a nested structure, but these regions include two pseudoknots that are universally conserved and essential for translation (Vila et al. 1994; Poot et al. 1996).

Removing pseudoknots from structural models is a relevant problem for a growing group of RNA researchers using computational tools. Because of limitations in software or algorithms, it is often necessary to work with nested structures. Pseudoknot elimination can be necessary, for example, to use the growing repository of RNA crystal structures (Ponty 2006; Tyagi and Mathews 2007), to compare different RNA structures (Andronescu et al. 2007), to classify structures into structural components (Smit et al. 2006), to create RNA covariance models (Eddy 2002), or to search for RNA homologs (Chang et al. 2006). Pseudoknots also increase the computational complexity of RNA structure prediction (Rivas and Eddy 1999; Lyngsø and Pedersen 2000) and visualization (Han and Byun 2003; Jossinet and Westhof 2005).

To make a knotted structural model pseudoknot-free, one or more base pairs must be “broken,” treating the corresponding part of the sequence as unpaired. The subtleties of this process are often underestimated. Which criteria should be used to designate one group of base pairs as “the helix” and another group as “the pseudoknot”? In the RNA molecule there may be no physical distinction between the conflicting helices that make up a pseudoknot at a structural level. Consequently, the decision about which helix to leave out can be arbitrary, undocumented, and difficult to reproduce. Pseudoknots are often removed manually or with unspecified computational methods (for example, Xayaphoummine et al. 2003; Smit et al. 2006; Andronescu et al. 2007; Metzler and Nebel 2008; Tyagi and Mathews 2007). An algorithm to calculate a nested structure containing the maximum number of base pairs has been described (Ponty 2006). In addition, related work has been done on computing the number of locally optimal secondary structures with respect to the Nussinov-Jacobson energy model (Clote 2005, 2006). Some information about pseudoknots is also available in the SSTRAND database (http://www.rnasoft.ca/sstrand). Progress is also being made in prediction and visualization of pseudoknotted structures, including a recent method for choosing pseudoknotted helices that minimize the sum of the stacking energies of the individual stems (Huang and Ali 2007).

The existing methods for pseudoknot removal are limited in their underlying assumptions: keeping the maximum number of base pairs is just one possibility and is not necessarily related to the RNA folding process. Other properties of the RNA that have been used to designate pseudoknots include the strength of a helix (either length or free energy), the folding pathway (which helix forms or melts first?), historical considerations (which helix was discovered first?), or ease of visualization. The manual approach is not feasible for large-scale analyses, and the computational approaches are often unspecified or not available to the wider community. The static pseudoknot information in a database is useful, but researchers need to be able to apply methods to their own structures as they are determined.

Current practices lead to duplicated efforts in automating the process of pseudoknot removal, inconsistent nomenclature, and difficulties reproducing analyses that require pseudoknot removal. Because the decision about how to unknot an RNA structure can have significant effects on structure-based downstream analyses, there is a need for explicit descriptions and user-friendly implementations of the possible algorithms.

We present a variety of algorithms for pseudoknot removal. We explore the use of different criteria to define pseudoknots: given a specific goal, each method points out the most critical base pairs that have to be removed. We demonstrate that different methods applied to the same initial structure lead to different results, indicating that pseudoknot removal affects structure-based analyses and that researchers should document their methods to improve reproducibility of their studies.

Our effort is in line with the goals of the RNA Ontology Consortium (Leontis et al. 2006), which aims to construct an ontology of RNA-related concepts to facilitate the integration of data flows in bioinformatics analyses (Jossinet et al. 2007). We intend to begin building a common vocabulary and set of reference software implementations for the removal process, such that pseudoknots can be removed in a consistent matter when the initial set of base pairs and removal method are specified. To make our software available to a wide audience, we provide not only the source code under the PyCogent project (Knight et al. 2007) but also a web interface (see Supplemental Data) and a standalone implementation offering command line control over the methods.

The main result of this study is the development and implementation of several methods to make knotted RNA structures pseudoknot-free. When pseudoknots are present in a structural model, different criteria can be used to remove them. Keeping the maximum number of base pairs ensures the least amount of information is lost. From a biological perspective, however, other criteria might come into play, such as the length of a helix or the distance between its upstream and downstream regions. We have implemented several heuristics with different underlying assumptions about what is important in the structure and a formal optimization approach that calculates all optimal solutions given some scoring function. All methods find saturated structures (Clote 2006), in which no base pair can be added without introducing a pseudoknot. This section begins with a short technical introduction outlining necessary concepts and definitions. We then describe the methods and their performance.

EC (elimination, conflicts) The EC method tries to reach a nested structure as fast as possible. It removes paired regions from the whole set beginning with the one with the most conflicts. If two regions have an equal number of conflicts, the region's gain and starting position are used to determine which region is processed first. This simple strategy might remove too many regions, resulting in an unsaturated structure, so nonconflicting regions are optionally added back (see the supplemental materials for details).

EG (elimination, gain) The EG heuristic eliminates conflicting paired regions on the basis of their “gain” (see RNA Structure and Pseudoknots). The algorithm processes the regions from the one with the smallest gain to the one with the largest gain, which means that the most unfavorable regions are removed first. It uses the number of conflicts and the starting position in case of equal gain. To prevent finding unsaturated structures, this method adds back nonconflicting regions (see the EC method).

IO (incremental, order) The IO method starts with an empty list and adds paired regions that do not conflict with any region that has already been added. It takes the most simplistic approach, which is to add paired regions one by one, from the 5′ end to the 3′ end or the other way around. The order of the regions is controlled by a parameter: the default order is from 5′ to 3′. This method is currently used by INFERNAL (E. Nawrocki, pers. comm.).

IL (incremental, length) The IL method operates under the idea that longer regions are more important than shorter regions. Thus, it adds nonconflicting paired regions one by one, starting with the longest region and working toward the shortest region. In case of equal lengths, the region starting closest to the 5′ end is added first (preferred end is controlled by a parameter).

IR (incremental, range) The IR method prefers short-range interactions over long-range interactions, thus in this scenario the structure is built-up starting with the formation of the hairpin loops. Paired regions are added to the list from short to long-range. If ranges are equal, the region starting closest to the 5′ end is added first (preferred end is controlled by a parameter).

Our optimization approach (OA) differs from the known DP algorithm in two respects. First, it calculates all optimal solutions rather than a single one. Second, it can handle arbitrary scoring functions, incorporating sequence or alignment information in addition to the structure. We have implemented the traditional scoring function to maximize the number of base pairs in the nested structure and a sequence-dependent function that maximizes the number of hydrogen bonds in Watson–Crick base-paired regions, where each GC base pair scores 3 points and each AU or GU pair scores 2 points (Mathews et al. 1999). A detailed description of the optimization approach can be found in the supplemental materials (see Supplemental Data).

A dynamic programming method is preferred to an exhaustive approach despite the relatively low complexity of many natural RNA structures in terms of pseudoknots. Although most collections of canonical base pairs in biological structures belong to a class called bisecondary structures (Haslinger and Stadler 1999), in which the topology of the pseudoknots is restricted, the number and size of the knot components increases rapidly in more exotic base pair collections. Even among sets of canonical base pairs, the complexity varies widely. For example, there are only 14 possibilities for the HDV ribozyme structure (Protein Data Bank [PDB] ID: 1DRZ) (Ferré-D'Amaré et al. 1998), but 3.8 × 10²⁹ possible nested structures for the E. coli large-subunit (LSU) rRNA structure (PDB: 2AW4) (Schuwirth et al. 2005).

Several factors influence the pseudoknot removal process. First, the chosen method is of major importance, because different criteria could all produce different solutions for the same knotted structure (Fig. 1). Second, the selection of base pairs forming the initial structure can lead to different conclusions under the same criterion. For example, in the set of canonical base pairs from the RNaseP structure (PDB ID: 2A64) (Kazantsev et al. 2005), the EC and IR methods found a nested structure containing 82 base pairs (helix P2 broken), while the EG, IO, and IL methods reached a solution of 81 base pairs (helix P4 broken). When adding the immediate helix extensions, there were two optimal solutions containing 95 base pairs each, indicating that P2 and P4 were of equal importance in terms of the number of base pairs. The fact that multiple optimal solutions for the same knotted structure might exist is also illustrated in Figure 2. Finally, we show (in Fig. 2) that applying a different scoring function will change the outcome of the optimization approach. Detailed information about the performance of the methods on biological and artificial structures can be found in Supplemental Tables 1 and 2.

FIGURE 1. Different behavior of pseudoknot removal methods. All algorithms can find different nested structures for a single knotted structure. This figure shows two different representations of a set of paired regions as found in a randomly generated artificial (more ...)

FIGURE 2. Different criteria for pseudoknot removal lead to different results for a pseudoknot in the group I intron structure (PDB ID: 1ZZN) (Stahley and Strobel 2005). Helices P3 and P7 are involved in a pseudoknot. There are two possible solutions to create (more ...)

We have presented a collection of automated methods for pseudoknot removal from RNA structures. This collection will be useful to the RNA community for two reasons. First, it allows pseudoknots to be removed in a consistent manner across different studies, so that analyses can be reproduced exactly. This consistency is especially important for benchmarking of new software that predicts or otherwise uses structural information. Second, it makes explicit the implicit assumptions that underlie the different methods, facilitating discussion of the pros and cons of the different methods for specific biological situations. In their role as operational definitions of pseudoknots, the presented methods can also spark a more informed discussion of which base pairs must be broken to form a nested structure.

The results demonstrate that there are many ways to remove pseudoknots and that, in general, these different methods give different results. For many pseudoknotted structures, all of the methods produce distinct solutions. The optimization method provides a formally optimal solution in terms of optimization of some score function, but this solution is often not unique. Therefore, we recommend that analyses are performed either by using each optimal solution and averaging the results or by using additional biologically informed criteria that are applied in a well-described and consistent manner to choose an optimal solution.

One important, but unresolved, question is how to decide which method is best in a given situation. Both the helices that contribute to a pseudoknot are usually standard A-form RNA helices: because designation of one group of base pairs as “the helix” and the other as “the pseudoknot” has often been performed manually in different families of sequences with no explicit justification—leading to the comment that “pseudoknots are pseudointeresting” (N. Pace, pers. comm.), comparison with existing nomenclature is likely to produce inconsistent results. This will be an important area of investigation for automatic use of structural information in databases such as Rfam (Griffiths-Jones et al. 2003) but can only be addressed with additional data. In particular, which nested structure is the best may depend on the goal: we can ask, for example, which of the methods is most likely to preserve conserved interactions, or which of the methods is most likely to produce covariance models that maximize the number of additional homologs found in the databases. However, these are empirical questions that require empirical studies to resolve.

In conclusion, the availability of a variety of pseudoknot resolution algorithms, along with reference implementations of these methods, fills an important and previously unappreciated need in the field. Because the procedure for pseudoknot removal often seems intuitively obvious, but because intuition differs between different researchers, the explicit rules and motivations for the different procedures have often been hidden. By making this information explicit and by providing a common vocabulary for describing the different methods, we now have a starting point for determining which of the various methods are optimal in specific cases and provide a platform on which more advanced methods can be constructed.

Supplemental material can be found at http://www.ibi.vu.nl/programs/k2nwww.

This work has been strongly motivated by discussions at the 2007 meeting of the RNA Ontology Consortium (NSF Grant No. 0443508) in Madison, WI. We are also grateful to Robin Dowell, Sébastien Lemieux, Michael Yarus, and Norman Pace for their contributions to the discussion. The comments of two anonymous reviewers have greatly helped us to improve this manuscript. This work is supported by the Centre for Medical Systems Biology. The Centre for Medical Systems Biology is supported by a Genomics Centre of Excellence grant from the Netherlands Genomics Initiative, which is funded by the Dutch government.