Predicting pseudoknotted structures across two RNA sequences.

MOTIVATION: Laboratory RNA structure determination is demanding and costly and thus, computational structure prediction is an important task. Single sequence methods for RNA secondary structure prediction are limited by the accuracy of the underlying folding model, if a structure is supported by a family of evolutionarily related sequences, one can be more confident that the prediction is accurate. RNA pseudoknots are functional elements, which have highly conserved structures. However, few comparative structure prediction methods can handle pseudoknots due to the computational complexity.
RESULTS: A comparative pseudoknot prediction method called DotKnot-PW is introduced based on structural comparison of secondary structure elements and H-type pseudoknot candidates. DotKnot-PW outperforms other methods from the literature on a hand-curated test set of RNA structures with experimental support. AVAILABILITY: DotKnot-PW and the RNA structure test set are available at the web site http://dotknot.csse.uwa.edu.au/pw. CONTACT: janaspe@csse.uwa.edu.au SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

1 INTRODUCTION

Macromolecules such as DNA, RNA and proteins have the ability to form diverse tertiary structures, which enable functionality and thus, life. For many decades, proteins were deemed the global players in the cell until RNA entered the spotlight. For example, RNA structures have been found to be catalytically active, which was assumed to be the privilege of proteins. Furthermore, small RNAs are known to regulate gene expression and RNA viruses employ a plethora of structure elements to invade the host cell.

To gain insight into macromolecule function, one must investigate the structure. The first step in RNA folding is stable base pairing that leads to a secondary structure. As RNA structure formation is of hierarchical nature, secondary structure is the basis for the tertiary fold that produces the functional structure. Especially for RNAs, structure determination by experimental means is an intricate and expensive task. Computational RNA structure prediction is therefore an invaluable tool for biologists. Comparative structure prediction is considered the most reliable approach for computational RNA structure prediction. Single sequence structure prediction is always limited by the accuracy of the underlying folding model. Three main streams have been identified for comparative RNA secondary structure prediction: (i) predict a structure from a pre-computed sequence alignment; (ii) simultaneously compute an alignment and a structure and (iii) alignment-free methods (Gardner and Giegerich, 2004).

Tools for multiple sequence alignments such as ClustalW (Thompson ) are readily available and thus, structure prediction from an alignment is a tempting approach [e.g. RNAalifold (Hofacker )]. Such methods heavily depend on the sequence conservation and quality of the underlying alignment. However, ncRNAs are conserved rather on the structure level than on the sequence level. The gold standard of RNA comparative structure prediction is the Sankoff approach as it does not rely on a high-quality sequence alignment and captures the structural conservation of ncRNAs. Sankoff (1985) introduced a theoretical dynamic programming algorithm for simultaneous folding and aligning for a set of N sequences that takes O(n3) time and O(n2) space. Practical variants have been derived which more or less retain the Sankoff principle by sacrificing optimality. Alignment-free methods aim to avoid the pragmatic restrictions made in a practical Sankoff approach as well as the reliance on a high-quality alignment [e.g. CARNAC (Perriquet )]. Note that all of these comparative structure prediction methods exclude the prediction of RNA pseudoknots.

RNA pseudoknots are crossing structure elements with diverse functions. The principle of pseudoknot formation is that bases within a loop region pair with complementary unpaired bases outside the loop. From an algorithmic point of view, even the simplest type of pseudoknot adds considerable computational demands due to crossing base pairs. In fact, the majority of comparative RNA structure prediction methods exclude pseudoknots. Biologists have delivered a wealth of studies, which show that pseudoknots have an astonishing number of diverse functions and occur in most classes of RNA (Staple and Butcher, 2005). RNA viruses use pseudoknots for hijacking the replication apparatus of the host (Brierley ).

A limited number of RNA comparative structure prediction methods can handle pseudoknots due to the computational complexity. Several of these methods take a sequence alignment as an input. ILM is an algorithm that takes as an input either individual sequences or a sequence alignment (Ruan ). A base pair score matrix is prepared initially and helices are added to the structure in an iterative fashion. In the approach hxmatch, a maximum weighted matching algorithm with combined thermodynamic and covariance scores is used (Witwer ). This program gives the option to be combined with RNAalifold. KNetFold is a machine learning method, which takes a sequence alignment as an input and outputs a consensus structure allowing pseudoknots (Bindewald and Shapiro, 2006). Simulfold takes an alignment as an input and simultaneously calculates a structure including pseudoknots, a multiple-sequence alignment and an evolutionary tree by sampling from the joint posterior distributions (Meyer and Miklos, 2007). Tfold combines stem stability, covariation and conservation to search for compatible stems and subsequently for pseudoknots for a set of aligned homologous sequences (Engelen and Tahi, 2010). Several comparative structure prediction methods including pseudoknots do not rely on an initial sequence alignment. The graph-theoretical approach comRNA computes stem similarity scores and uses a maximum clique finding algorithm to find pseudoknotted structures (Ji ). SCARNA performs pairwise structural alignment of stem fragments with fixed lengths derived from the probability dot plot (Tabei ).

In the following, a novel comparative approach for predicting structures including H-type pseudoknots called DotKnot-PW will be introduced. The input consists of two unaligned, evolutionarily related RNA sequences. Similarity scores between structure elements will be calculated. Statistically significant pairs will be used to find the set of conserved structure elements common to two sequences, which maximize a combined thermodynamic and similarity score. Using a hand-curated test set of pseudoknotted structures with experimental support, the prediction accuracy of DotKnot-PW will be compared with methods from the literature.

2 APPROACH

Pseudoknots are functional elements in RNA structures and therefore, the most promising approach for comparative prediction is a structure comparison with less focus on exact sequence matching. In fact, perfect conservation on the sequence level can be more of a curse than a blessing. Especially ncRNAs are known to evolve quickly and so-called consistent and compensatory base pairs in both sequences will give much more confidence for structure conservation than a sequence alignment. One strong point of the DotKnot method for single sequence pseudoknot prediction (Sperschneider and Datta, 2010; Sperschneider ) is that the set of possible H-type pseudoknot candidates (and secondary structure elements) is explicitly computed and thus readily available for further investigation. The main steps in the pairwise pseudoknot prediction approach DotKnot-PW are as follows (Fig. 1):

Fig. 1.

For two unaligned RNA sequences Seq and Seq, DotKnot-PW produces structure element dictionaries derived from the probability dot plot. Similarity scores and P-values are computed to detect conserved elements

Run DotKnot for two unaligned sequences Seq and Seq. This returns secondary structure element and H-type pseudoknot candidate dictionaries.

Calculate pairwise base pair similarity scores for secondary structure elements and H-type pseudoknot candidates. Keep significant pairs that have a low estimated P-value.

Use significant pairs to calculate the set of conserved structure elements and pseudoknots for the two sequences that maximizes a combined free energy and similarity score.

The key point of the DotKnot-PW approach is how to score the similarity of stems, secondary structure elements and H-type pseudoknot candidates derived from sequences Seq and Seq. Related work has been done for stem finding in unaligned sequences, where stem candidates are assigned a matching score across unaligned sequences, e.g. in SCARNA. Another point is how to assess the significance of a similarity score using P-values. These points will be explained in detail in the following section.

3 MATERIALS AND METHODS

For two unaligned RNA sequences Seq and Seq, the single sequence prediction method DotKnot (Sperschneider and Datta, 2010; Sperschneider ) returns two stem dictionaries D(x) and D(y) derived from the probability dot plot. It also returns secondary structure element dictionaries and and H-type pseudoknot candidate dictionaries D(x) and D(y) (Fig. 1). To detect conserved structure elements for the two sequences, a pairwise structural comparison is performed. Instead of a full structure-to-structure alignment, which takes O(n4) time and O(n3) space, pairwise base pair similarity scores are calculated using the RIBOSUM85-60 matrix for base pair substitutions (Klein and Eddy, 2003).

3.1 Base pair similarity scores of stems

For two given stems s(x) and s(y) with fixed lengths in sequences Seq and Seq, respectively, the base pair similarity score sim[s(x), s(y)] is calculated using an ungapped local structure alignment of the base pairs with the RIBOSUM85-60 matrix. As an example, consider the following optimal ungapped local structure alignment of the two stems with base pair similarity score of sim[s1(x), s2(y)] = 22.04 using the RIBOSUM85-60 matrix.

s 1(x)

UCUCUAUC … … .GAUAGAGA

((((((((… … .))))))))

s 1(y)

–UUGUAC… … .GUACAA–

–((((((… … .))))))–

To evaluate the significance of base pair similarity scores instead of the raw score, one has to find out what the underlying probability distribution is. Similar to the case of ungapped local sequence alignments (Karlin and Altschul, 1990), it is assumed here that the base pair similarity scores follow an extreme value distribution. However, the main difference is that a comparison between fixed-length stem fragments is made. It is important to remember that parameters λ and K describe the extreme value distribution of optimal local alignment scores in the asymptotic limit of long sequences (Altschul ). Here, the parameters for the generalized extreme value distribution are pre-calculated using maximum likelihood fitting of a distribution to the histogram of a large sample of random base pair similarity scores. The maximum likelihood fitting was performed using the ismev package of the R statistical language for a range of stem lengths (see Supplementary Material). The P-value is defined as the probability to obtain a score greater than or equal to the observed score strictly by chance. A stem s(x) in sequence Seq and a stem s(x) in sequence Seq are a significant pair if the score sim[s(x), s(y)] has an estimated P-value less than α. Stem pairs with a P-value larger than α are not considered in the following.

3.2 Base pair similarity scores of interrupted stems

For two interrupted stems, the base pair similarity score is calculated by deleting bulges and internal loops and scoring stems as consecutive base pairs. Base pair similarity scores for regular and interrupted stems are also calculated if the difference in number of base pairs is less than 5. For example, a stem with one bulge might be a conserved match with a regular stem. A stem s(x) in sequence Seq and a stem s(y) in sequence Seq are a significant pair if the score sim[s(x), s(y)] has an estimated P-value less than α.

3.3 Base pair similarity scores of multiloops

Calculating the base pair similarity score for two multiloop structures is complex due to the variety of inner loop elements, which may be regular or interrupted stems. A multiloop can be decomposed into an outer stem and a set of inner structure elements S(x) = [s1(x),…,s(x)], where k ≥ 2. The base pair similarity score for the outer stems of two multi-loops can be easily obtained from the previously calculated base pair similarity scores. If the outer stem is a conserved match, a local alignment on the set of inner structure elements is used to find the base pair similarity score. Here, gaps are allowed in the local alignment of inner structure elements; however, no gap penalty is used. Let two sets of inner structure elements S(x) = [s1(x),…, s(x)] and S(y) = [s1(y),…, s(y)] be given. Let H(i, j) be the maximum similarity score between a suffix of S(x) and a suffix of S(y). The optimal local alignment is calculated as follows:

A multiloop in sequence Seq and a multiloop in sequence Seq are a significant pair if the similarity score has an estimated P-value less than α.

3.4 Base pair similarity scores of H-type pseudoknots

A H-type pseudoknot has two pseudoknot stems S1 and S2. The prerequisite for a conserved pseudoknot pair is that both core H-type pseudoknot stem pairs [S1(x), S1(y)] and [S2(x), S2(y)] are significant. The base pair similarity score for two H-type pseudoknots p(x) and p(y) in sequences Seq and Seq, respectively, is the sum of base pair similarity scores for the core pseudoknot stems as well as the base pair similarity score from a gapped local alignment of the recursive secondary structure elements in the loops (as described for multiloops). A pseudoknot p(x) in sequence Seq and a pseudoknot p(y) in sequence Seq are a significant pseudoknot pair if the similarity score sim[p(x), p(y)] has an estimated P-value less than α.

3.5 Dissimilarity and weight of significant structure elements pairs

The base pair similarity score calculated in the previous sections might not be powerful enough to distinguish true positive conserved structure element pairs from false-positive structure element pairs due to the finite lengths of stems and exclusion of loop sequences in the alignment. Therefore, a dissimilarity score is also used to confirm whether a pair is significant. The dissimilarity for two given structure elements s(x) and s(y) in sequences Seq and Seq is defined as: where dissim1 is the difference in the stem lengths and dissim2 is the difference in the number of loop lengths. As an example, consider the pseudoknot pair p1(x) and p1(y) in sequences Seq and Seq, respectively, with stems S1, S2 and loops L1, L2, L3. The pseudoknot pair has dissimilarity of 6.

p 1(x)

UCUCUAUCAGAAUGGAUGUCUUGCUGCUAUAAUAGAUAGAGAAGGUUAUAGCAG

((((((((… … … … …[[[[[[[[[[.)))))))).]]]]]]]]]]

p 1(y)

UUGUACAGAAUGGUAAGCCAAGUGUCAAUAGGAGGUACAAGCAACCUAUUGCAU

((((((… … … … …[[[.[[[[[[[.))))))… .]]]]]]]]]]

A weight is assigned to a significant pair, which is a combination of the free energy, covariation and dissimilarity. The overall weight s of a significant structure element pair [s(x), s(y)] in sequences Seq and Seq is a combination of the free energy weights w[s(x)] and w[s(y)], base pair similarity score sim[si(x), s(y)] and dissimilarity dissim[s(x), s(y)]:

Only structure element pairs with positive score s are allowed in the following dynamic programming algorithm. Here, α and γ are set to 0.5 and β is set to 1.

3.6 Finding best set of significant structure elements

Let be the number of structure elements in the first sequence Seq and be the number of structure elements in the second sequence Seq. Each structure element has a left and right endpoint in the sequence and is a stem, interrupted stem, multiloop or H-type pseudoknot. Structure elements can also be represented as nodes in a graph. In each sequence, the structure elements are ordered by their right endpoints. An edge is drawn between two structure elements in the first and the second sequence if their base pair similarity score has a P-value less than α. Given the set of edges between nodes and , the goal is to find the set of edges with maximum weight that are non-crossing. This relates to finding the set of non-overlapping structure elements in the two sequences that maximize the score under the requirement that the interval ordering is preserved. A set of structure elements in the first and second sequence, which preserves the interval ordering is called a feasible structure element alignment and must satisfy the following two requirements. Each structure element can be aligned with at most one other structure element in the other sequence. The order of structure elements must be preserved with respect to the alignment. That is, if structure elements and in the first sequence are aligned with and in the second sequence, respectively, the pairs may never overlap: (Fig. 2).

Fig. 2.

A set of edges with positive scores is given between nodes p1,…,p7 in the first sequence and p1,…,p6 in the second sequence. The goal is to find the best set of non-overlapping structure elements in the two sequences such that the interval ordering is preserved. The optimal structure element alignment, which preserves the interval ordering includes structure elements p1, p4, p7 in the first sequence and p1, p4, p6 in the second sequence

Given nodes in the first sequence Seq and in the second sequence Seq, let f(i, a) be the maximum sum of edge weights for nodes between 1 and i in the first sequence and 1 and a in the second sequence such that the edges are non-crossing (i ≤ n and a ≤ m). The nodes that maximize the sum of edge weights are called an optimal structure element alignment for the two sequences. The optimal structure element alignment is calculated using dynamic programming. For a given structure element p with start point a and end point b, let pre(i) be the non-overlapping predecessor. For each structure element, its predecessor is pre-computed using the sorted list of structure elements. The recursion for calculating the optimal structure element alignment is as follows:

Furthermore, nested structures are taken into account for significant outer stem pairs, which have estimated P-value less than α. For each significant outer hairpin loop pair, the optimal structure element alignment of inner elements is computed.

3.7 Time requirements

For two unaligned RNA sequences Seq and Seq, the single sequence prediction method DotKnot returns structure element dictionaries derived from the probability dot plot. Let n and m be the number of structure elements in sequences Seq and Seq, respectively. Calculating the similarity scores and the optimal structure element alignment takes O(nm) time. Furthermore, nested structures are taken into account for significant outer stem pairs, which have estimated P-value less than α. Let a be the number of significant stem pairs, where both stems are hairpin loops. For each significant outer hairpin loop pair, the optimal structure alignment of inner elements is computed. In the worst case, this increases time requirements to O(a × nm). The number of structure elements depends on the base composition of the sequence. Empirically, n and m can be observed to grow linearly with the length of the sequence for uniform base distribution (see Supplementary Material). In practice, DotKnot-PW can be expected to run in the order of minutes for sequences shorter than 500 nt.

4 RESULTS

Many pseudoknot prediction programs have been evaluated using all the entries in the PseudoBase database (van Batenburg ). There are several caveats in this approach. First, the sequences given in PseudoBase are those which exactly harbor the pseudoknot. However, in practice structure prediction algorithms will be applied to longer sequences without prior knowledge of the pseudoknot location. Second, long-range pseudoknot entries appear in a truncated version in the database. Third, some classes of pseudoknots have a large number of entries (such as short H-type pseudoknots in the 3′-untranslated regions of plant viruses), whereas more complex types of pseudoknots only have one representative (such as long-range rRNA pseudoknots). Therefore, a hand-curated dataset of pseudoknot structures will be used here.

When it comes to pseudoknots, many structures have been published based on a secondary structure predicted by free energy minimization. These predicted secondary structures are used as a working model and refined using experimental techniques such as chemical and enzymatic probing. However, the native structure remains unsolved unless tertiary structure determination methods such as X-ray crystallography are used. Testing structures that are based on computer predictions with no experimental support creates a bias in the benchmark and will be avoided in this evaluation.

A total of 16 pseudoknotted reference structures from different RNA types were collected, which have strong experimental support. For each reference structure, a supporting set of 10 evolutionarily related sequences was obtained from the database (Gardner ). Note that for the vast majority of supporting sequences, no experimentally determined structures are available. The average pairwise sequence identities vary from 55% to 99%. Given a reference structure, the performance of prediction algorithms is evaluated in terms of sensitivity (S), i.e. the percentage of base pairs in the reference structure, which are predicted correctly, as well as positive predictive value (PPV), i.e. the percentage of predicted pairs, which are in the reference structure. The Matthews correlation coefficient (MCC) is also reported and is in the range from −1 to 1, where 1 corresponds to a perfect prediction and −1 to a prediction that is in total disagreement with the reference structure. The performance of each method for predicting the reference structure was evaluated as described in Gardner and Giegerich (2004).

DotKnot-PW was compared with methods that are freely available and use standard input and output formats. The comparative methods are CARNAC, Tfold and hxmatch (with the -A option using RNAalifold). All of these methods return structure predictions for only the reference structure with regards to the support set of evolutionarily related sequences. Tfold and hxmatch take a sequence alignment as the input. ClustalW with the default parameters was used to produce the initial sequence alignment. DotKnot-PW and CARNAC take a set of unaligned sequences as the input. Furthermore, prediction results for the reference sequence (not the supporting sequences) were obtained from the single sequence methods DotKnot (Sperschneider and Datta, 2010; Sperschneider ), ProbKnot (Bellaousov and Mathews, 2010), IPknot (Sato ) and RNAfold (Hofacker ). Note that all methods except CARNAC and RNAfold allow pseudoknot prediction.

The results are shown in Table 1. DotKnot-PW has the highest average MCC of 0.75 for the test sequences. For each reference structure with the 10 support sequences from the corresponding RFAM family, 10 predictions are returned ordered by the combined free energy and similarity score. If only the pairwise prediction with highest combined free energy and similarity score is taken, DotKnot-PW has an improved average MCC of 0.81. Tfold and hxmatch have average MCC of 0.6 and 0.59, respectively. CARNAC has average MCC of 0.45 with much higher average specificity than sensitivity.

Table 1.

Prediction results using a test set of different RNA classes

Type	ID	Info	Family		Prediction of common structure				Reference structure prediction
					DotKnot-PW, First (Average)	CARNAC	Tfold	hxmatch	DotKnot	ProbKnot	IPknot	RNA fold
Frameshift	SARS-CoV (Plant et al., 2005)	NMR	RF00507	S	92.3 (79.2)	38.5	57.7	38.5	92.3	69.2	73.1	73.1
		82 nt	82 nt	PPV	100 (90)	90.9	93.8	58.8	100	72	79.2	73.1
		1 PK	62%	MCC	0.96 (0.84)	0.59	0.73	0.47	0.96	0.7	0.76	0.73
	VMV (Pennell et al., 2008)	NMR	RF01840	S	100 (100)	0	50	42.9	100	50	50	50
		68 nt	55 nt	PPV	100 (100)	0	50	54.5	100	41.2	41.2	41.2
		1 PK	89%	MCC	1 (1)	0	0.49	0.48	1	0.44	0.44	0.44
Ribozyme	HDV (Ferre-D’Amare et al., 1998)	X-ray	RF00094	S	93.8 (90.7)	21.9	21.9	12.5	93.8	40.6	62.5	37.5
		87 nt	90 nt	PPV	100 (97.1)	100	30.4	57.1	100	48.1	80	42.9
		1 PK	74%	MCC	0.97 (0.94)	0.46	0.24	0.26	0.97	0.43	0.7	0.39
	glmS-Ba (Klein and Ferre-D’Amare, 2006)	SC, MG	RF00234	S	76.4 (71.7)	7.3	50.9	43.6	63.6	76.4	72.7	69.1
		151 nt	178 nt	PPV	95.5 (92.9)	44.4	87.5	80	67.3	82.4	85.1	82.6
		2 PKs	55%	MCC	0.85 (0.81)	0.18	0.66	0.59	0.65	0.79	0.78	0.75
	EC-RNaseP (Harris et al., 2001)	X-ray	RF00010	S	55.3 (51)	26.8	29.3	46.3	53.7	74	69.9	60.2
		377 nt	380 nt	PPV	68.7 (64.4)	86.8	65.5	83.8	56.4	77.8	86	64.9
		2 PKs	64%	MCC	0.61 (0.57)	0.48	0.44	0.62	0.55	0.76	0.77	0.62
Untranslated regions	BCV (Williams et al., 1999)	SP, MG	RF00165	S	100 (95.6)	55.6	61.1	0	100	55.6	94.4	55.6
		63 nt	63 nt	PPV	100 (96.7)	100	84.6	0	81.8	66.7	100	66.7
		1 PK	78%	MCC	1 (0.96)	0.74	0.71	-0.01	0.9	0.6	0.97	0.6
	BaMV (Lin et al., 2007)	SP, MG	RF00290	S	59.5 (53.8)	0	*	50	40.5	64.3	50	40.5
		170 nt	145 nt	PPV	69.4 (62.4)	0	*	65.6	41.5	69.2	67.7	45.9
		1 PK	99%	MCC	0.64 (0.58)	0	*	0.57	0.4	0.66	0.58	0.43
	HPeV1 (Nateri et al., 2002)	MG	RF00499	S	95.3 (86)	11.6	37.2	79.1	83.7	88.4	86	79.1
		116 nt	112 nt	PPV	95.3 (88.2)	100	100	91.9	83.7	92.7	97.4	94.4
		1 PK	88%	MCC	0.95 (0.87)	0.34	0.61	0.85	0.83	0.9	0.91	0.86
Telomerase	Tthe-telo (Theimer and Feigon, 2006)	SP	RF00025	S	81.6 (61.1)	68.4	60.5	57.9	73.7	65.8	71.1	60.5
		159 nt	160 nt	PPV	91.2 (70.7)	86.7	79.3	84.6	70	55.6	67.5	54.8
		1 PK	72%	MCC	0.86 (0.65)	0.77	0.69	0.7	0.72	0.6	0.69	0.57
Riboswitch	preQ1 (Rieder et al., 2010)	NMR	RF00522	S	100 (73.4)	55.6	*	55.6	55.6	55.6	55.6	55.6
		34 nt	44 nt	PPV	100 (100)	100	*	100	100	55.6	100	71.4
		1 PK	80%	MCC	1 (0.84)	0.74	*	0.74	0.74	0.53	0.74	0.61
	SAM-I (Montange and Batey, 2006)	X-ray	RF00162	S	73.5 (53.8)	55.9	55.9	64.7	47.1	73.5	64.7	64.7
		94 nt	102 nt	PPV	100 (83.7)	79.2	100	91.7	57.1	100	100	88
		1 PK	73%	MCC	0.85 (0.66)	0.66	0.74	0.77	0.51	0.85	0.8	0.75
IRES	PSIV (Pfingsten et al., 2006)	X-ray	RF00458	S	65.5 (63.8)	39.7	94.8	77.6	70.7	70.7	63.8	72.4
		194 nt	198 nt	PPV	66.7 (73)	100	96.5	91.8	69.5	67.2	78.7	71.2
		2 PKs	57%	MCC	0.66 (0.68)	0.63	0.96	0.84	0.7	0.69	0.71	0.72
	CSFV (Kolupaeva et al., 2000)	SP, MG	RF00209	S	39 (44.9)	45.1	68.3	74.4	62.2	75.6	65.9	68.3
		244 nt	272 nt	PPV	42.1 (53.1)	97.4	80	85.9	68	80.5	68.4	71.8
		1 PK	83%	MCC	0.4 (0.49)	0.66	0.74	0.8	0.65	0.78	0.67	0.7
mRNA	S15 (Philippe et al., 1995)	SP, MG	RF00114	S	100 (81.8)	52.9	52.9	41.2	100	58.8	58.8	58.8
		74 nt	112 nt	PPV	100 (83.4)	52.9	40.9	46.7	100	47.6	52.6	52.6
		1 PK	77%	MCC	1 (0.82)	0.52	0.46	0.43	1	0.52	0.55	0.55
	repZ (Asano and Mizobuchi, 1998)	SP, MG	RF01087	S	73.8 (71.4)	9.5	31	83.3	71.4	57.1	66.7	57.1
		149 nt	149 nt	PPV	86.1 (80.2)	44.4	61.9	85.4	85.7	70.6	82.4	64.9
		1 PK	90%	MCC	0.8 (0.76)	0.2	0.43	0.84	0.78	0.63	0.74	0.6
tmRNA	Ec-tmRNA (Nameki et al., 1999)	NMR, SC	RF00023	S	37.5 (50)	6.7	34.6	27.9	75	56.7	66.3	50
		363 nt	377 nt	PPV	45.9 (58.4)	100	67.9	61.7	79.6	60.8	76.7	47.7
		4 PKs	57%	MCC	0.41 (0.54)	0.26	0.48	0.41	0.77	0.59	0.71	0.49
	Average			S	77.7 (70.5)	31	50.4	49.7	74	64.5	67	59.5
	Average			PPV	85.1 (80.9)	73.9	74.2	71.2	78.8	68	78.9	64.6
	Average			MCC	0.81 (0.75)	0.45	0.6	0.59	0.76	0.66	0.72	0.61

Note: Each reference structure is given by its ID (see Supplementary Material for dot-bracket notation). The following column gives the method of experimental support (NMR, NMR spectroscopy; X-ray, X-ray crystallography; SC, sequence comparison; MG, mutagenesis; SP, structure probing), length of the sequence and number of pseudoknots. For each reference structure, the corresponding RFAM family ID, average sequence length and average pairwise sequence identity is shown. The * symbol means that the method failed to run. The ‘first’ prediction for DotKnot-PW is the pairwise prediction with highest combined free energy and similarity score.

The prediction results for single sequence structure prediction for each of the reference sequences with experimentally determined structures are also shown in Table 1. Note that this does not include the prediction for the support sequences from RFAM, as no experimentally determined structures are available. All single sequence pseudoknot prediction methods show improved results over using RNAfold. DotKnot has the highest average MCC of 0.76, followed by IPKnot and ProbKnot.

As an example, consider the S15 mRNA pseudoknot that binds to specific proteins in the autoregulation mechanism of ribosomal protein S15 synthesis (Philippe ). For the reference sequence S15 and 10 support sequences from the corresponding RFAM family, DotKnot-PW returns pairwise predictions ordered by the combined free energy and similarity score. The top two pairwise predictions with the highest scores are shown in Figure 3.

Fig. 3.

Pairwise prediction results for the S15 mRNA pseudoknot (RFAM family RF00114) with the top two combined free energy and similarity scores. The reference structure is shown at the top and folds into two conformations in dynamic equilibrium: a H-type pseudoknot or a series of hairpins. For the pairwise prediction with highest score, the pseudoknot structure is returned. For the second-best pairwise prediction, the alternative hairpin loop structure is returned

5 DISCUSSION

We presented DotKnot-PW for prediction of structures common to two RNA sequences, including H-type pseudoknots. Both DotKnot and DotKnot-PW have been designed as dedicated pseudoknot prediction tools. In the following, important aspects of pseudoknot prediction will be discussed.

5.1 The underlying folding model

Single sequence prediction methods are always limited by the underlying RNA folding model. This may be the set of free energy parameters used by free energy minimization methods or the underlying methodological framework such as maximum expected accuracy methods. DotKnot-PW shows excellent results on H-type pseudoknots with short interhelix loops. For this type of pseudoknots, DotKnot-PW uses free energy pseudoknot parameters by Cao and Chen (2006, 2009) based on polymer statistical mechanics. Improvements of the accuracy of free energy parameters, both for secondary structures and pseudoknots, will lead to more accurate prediction methods. However, one has to keep in mind that the algorithms themselves must be designed in such a fashion that novel parameters can be efficiently incorporated. The heuristic framework of DotKnot-PW has been designed such that it can incorporate sophisticated free energy parameters for pseudoknots, secondary structures and coaxial stacking. In the future, DotKnot-PW could also use contributions from basic tertiary structure elements such as base triples around the pseudoknot junction or stem-loop interactions.

5.2 The type of pseudoknot

Pseudoknot prediction algorithms come in two flavors: either they can predict a certain, restricted class of pseudoknots or they do not have a restriction on the type of pseudoknot that can be predicted. For methods using free energy parameters, the inclusion of general types of pseudoknots might be more of a curse than a blessing, as no reliable free energy parameters for complex pseudoknots are available. DotKnot-PW has restrictions on the type of pseudoknot that can be predicted. However, this does not always lead to poor prediction results in practice. For example, DotKnot-PW shows the best result for the HDV ribozyme, which is a complex double nested pseudoknot.

5.3 The trouble with benchmarking

The results from the benchmark for structure prediction in Table 1 must be interpreted with care. First, the tested methods can be run with different parameters, possibly producing better results. However, as a typical user has no prior knowledge about the structure, the default parameters for each method are used. Of course, a comprehensive benchmark should include a larger number of structures to obtain a more reliable evaluation. However, in this study, the focus has been on a test set where the structures are supported experimentally. Many structures have been published, which were determined using computational tools and this will inevitably create a bias in a benchmark, and thus they were excluded here.

5.4 Gaining confidence with multiple sequences

Here, an extension of the single sequence prediction method DotKnot was presented based on the pairwise comparison of structure elements. This approach called DotKnot-PW is designed as an algorithm for finding the structure including H-type pseudoknots common to two sequences. As shown in Table 1, DotKnot-PW can greatly improve structure predictions for RNA families when compared with the single sequence prediction using DotKnot. In some cases, a comparative approach might have lower sensitivity than a single sequence prediction; however, this should not generally be judged as ‘inferior’. For example, ncRNAs might preserve some integral base pairs throughout evolution and only these will be detected by a comparative approach, which returns the set of base pairs common to a set of evolutionarily related sequences. DotKnot-PW uses a set of unaligned sequences as the input; therefore, no expert user intervention is required. In the future, DotKnot-PW will be extended to include intramolecular kissing hairpins. Furthermore, constrained folding will be implemented to predict a structure subject to constraints, e.g. enforce certain base pairs or regions, which must remain unpaired.

6 CONCLUSION

DotKnot-PW has been designed as a dedicated pseudoknot prediction tool and should be applied to RNA sequences where pseudoknotted interactions are suspected in the structures. Prediction accuracy will inevitably decrease for sequences, which are longer than say 400 nt for any single sequence structure prediction method (Reeder ). To achieve reliable results, short sequences should be folded using DotKnot and predictions should be compared with results from other methods from the literature. To gain confidence in predictions, subsequent comparative prediction using DotKnot-PW and other comparative methods is highly recommended. Ideally, experimental verification of computationally predicted pseudoknots should be sought.

Click here for additional data file.