Related work

Computational structure prediction can be divided into de novo structure prediction and comparative structure prediction, which derive structures from a single sequence and multiple homologous ncRNAs respectively. As our method is to derive the consensus shape and structure of homologous ncRNAs, we briefly introduce related work in comparative ncRNA structure derivation. There are three general approaches for structure derivation from multiple sequences: simultaneously align and fold, align-then-fold, and fold-then-align. It is computationally expensive to simultaneously align and fold pseudoknot structures. The performance of the align-then-fold pseudoknot prediction heavily depends on the quality of the alignment. Usually multiple sequence alignment (MSA) tools such as ClustalW [23] are used to generate the alignment as the input to the folding tool. However, common structures can be missed due to misalignment between sequences lacking significant similarity [24]. In this work, we design a pseudoknot prediction tool using the fold-then-align strategy that does not require an alignment as input. Tools based on fold-then-align use a de novo folding tool to construct a small but representative sample of the folding space, which consists of the predicted optimal and sub-optimal structures. Structures from the folding space are chosen to maximize the structural and sequence similarity.

A number of software packages exist to predict the abstract shape for a single sequence. The sum of energies or the accumulated Boltzmann probabilities of all structures within a shape have been used as main features for shape prediction. The latter is often referred to as the shape probability. Usually the shapes with small sum of energies or high shape probabilities are more likely to be the correct shapes. It is claimed in RapidShapes [22] that using shape probabilities has superior performance over free energy-based approach because of its independence on sequence length and base composition. However, exact computation of the shape probability incurs exponential computational cost to the sequence length [22]. Thus, various heuristics or restrictions [25, 26] have been adopted for fast shape probability computation.

RNAcast [12] derives the consensus shape from homologous pseudoknot-free sequences based on the fold-then-align strategy. Structures are grouped based on their shapes and shapes are ranked by sum of free energies of structures within the shape in ascending order. The first-ranked shape is presented as the consensus shape. The consensus structure is derived from the lowest free energy structures of each sequence within the shape.

RNA structures and their representations

RNA structures and pseudoknots

RNA molecules fold into complex three dimensional structures by nitrogenous bases that are connected via hydrogen bonds. The secondary structure of an ncRNA is defined by the interacting base pairs. Some RNA molecules fold into pseudoknot structures by paring bases in loop regions with bases outside the stem loop (Figure 1a).

In this work, two types of ncRNA secondary structure representations are used. The first type is the arc-based representation, where nucleotides and hydrogen bonds are represented by vertices and arcs, respectively (Figure 1b). For pseudoknot-free secondary structures, all arcs are either nested or in parallel. Crossover arcs indicate pseudoknots. The second type is based on dot-bracket notation, where '.' represents unpaired bases and matching parenthesis '(' and ')' indicate base-pairing nucleotides. Following the annotation of Rfam [3], we use an extended dot-bracket notation to represent pseudoknot structures. The base-pairing nucleotides forming pseudoknots are represented by upper-lower case character pairs, such as A..a or B..b, as shown in Figure 1c.

Abstract shapes

Abstract shapes were formally introduced by Giegerich et al. [9]. The folding space of a given RNA sequence is partitioned into different classes of structures, by means of abstracting from structural details. These classes are called abstract shapes, or shapes for short.

An RNA secondary structure can be mapped to an abstract shape with different levels of abstraction [12]. In the abstract shape, details about the lengths of the loop and stacking regions are removed (see Figure 1 for examples of stacking and loop regions). Stacking regions are represented by pairs of brackets and unpaired regions are represented by underscores.

Pseudoknots are represented by pairs of upper-lower case characters. Figure 2 presents examples of the abstract shapes of level 1, 3, and 5 of a pseudoknot-free structure and a pseudoknot. Level 5 represents the strongest abstraction and ignores all bulges, internal loops, and single-stranded regions. Level 3 adds the helix interruptions caused by bulges or internal loops. Level 1 only abstracts from loop and stack lengths while retains all single-stranded regions.

Shape prediction

In this section we describe KnotShape, a comparative shape prediction tool for homologous ncRNA sequences that allows pseudoknots. The task of shape prediction is to select the best representative shape for given homologous sequences. In order to identify the best shape, various features such as shape probability [22], sum of energies of all structures in this shape [12], and the rank sum score [12] are evaluated to rank shapes. It has not been systematically assessed whether combinations of multiple features can lead to better shape prediction. In this work, we incorporate Support Vector Machine (SVM) [27] and feature selection techniques to determine important features for shape identification. In addition, we adopted a machine learning-based scoring function to evaluate the qualities of shapes.

The method contains two important components. The first one is the consensus shape prediction (Knot-Shape) and the second one is structure prediction using predicted shape as input (KnotStructure). We will first describe KnotShape, focusing on the feature construction and selection strategy. Then we will describe how to derive the consensus structure given the consensus shape.

Notation

Figure 3 illustrates the mapping between sequences, structures, and shapes. The input is a set of homologous ncRNAs and the output is the predicted consensus shape. Notations used in this paper correspond to this mapping.

• The N homologous ncRNAs constitute the input sequence space. X _i represents the i th sequence.

• Each sequence can be folded into different secondary structures. Let S ⁱ represent the set of folded structures of the i th sequence X _i . The set of structures predicted from all N input sequences is the union of S ⁱ : S = S¹ ∪ S² ∪ . . . ∪ S ^N .

• $S_j^i$ is the j th structure in the folding space of X _i . Its free energy is denoted by $\mathrm{\Delta }G\left(S_j^i\right)$. For a sequence X _i , the minimum free energy MFE(X _i ) is the lowest free energy among the energies of all predicted structures of X _i , i.e. $\text{MFE}(X_i)={\min }_{1\le j\le \left|S^i\right|}\mathrm{\Delta }G(S_j^i)$.

• All structures in S can be classified into a set of abstract shapes. For a shape P, we record its associated sequences and structures. P.LX denotes the set of associated sequences, each of which can fold into a structure with shape P. P.LS denotes all structures with shape P.

• $\widehat{P}$ is the predicted shape of the given homologous sequences X₁, X₂, .., X _N .

In order to explore the large folding space of multiple homologous sequences, we use a de novo folding tool to output the optimal and sub-optimal structures within a given energy cutoff. This heuristic does not allow us to explore the complete folding space. Given the observation that the correct structure is usually close to the optimal structure, this heuristic works well in practice [28].

Consensus structure prediction given a shape

Once we determine the shape, we will predict the structure in the shape class for the given homologous ncRNAs. Structures corresponding to the same shape can differ significantly in the details of the loop and stacking regions. A strategy is needed to choose the correct structure inside the shape class for each input sequence. The simplest strategy is to output the MFE structure for the chosen shape, which has been used in previous work [12]. However, the MFE structure in a shape may not be the native structure. In particular, the accuracy of the MFE prediction for ncRNAs containing pseudoknots is low.

In this section we describe KnotStructure, a comparative structure prediction method for homologous sequences given the shape. The rationale behind comparative structure prediction is that the secondary structures and sequences are conserved during evolution. Thus, finding the structures to maximize both the sequence and the secondary structure similarity among homologous ncRNAs provides the basis for comparative structure prediction. Various methods for evaluating structural and sequence similarity exist. The major challenge is to efficiently select $\left|\widehat{P}.LX\right|$ representative structures to achieve the highest structural and sequence similarity.

As we already derived the consensus shape $\widehat{P}$ using KnotShape, only structures with shape $\widehat{P}$ will be allowed. In addition, for each sequence $X_i\in \widehat{P}.LX$, only one structure with shape $\widehat{P}$ can be selected. The total number of combinations of structures for measuring the similarity is thus ${\prod }_{i=1to\left|\widehat{P}.LX\right|}|\widehat{P}.LS\cap S^i|$, where $\widehat{P}.LS\cap S^i$ contains structures with shape $\widehat{P}$ for a sequence X _i . Although efficient heuristics exist to measure the similarity among multiple structures and sequences, the sheer amount of combinations poses a great computational challenge.

Procedure 1 Representative structures selection

Input: $\widehat{P}$, $\widehat{P}$.LX, $\widehat{P}$.LS

Output: The representative structures

  1. Initialization

  for Every two structures $S_i^x$ and $S_j^y$do

    //only evaluate similarity of structures from different sequences

    if x ≠ y then

      Evaluate the similarity of $S_i^x$ and $S_j^y$

    else

      $S_i^x$ and $S_j^y$ has similarity -∞

    end if

  end for

  2. Select the set of representative structures using hierarchical clustering

  //Each structure is a cluster by itself

  repeat

    Combine a pair of clusters with the highest similarity

    For any structure $S_i^x$ added to the cluster, remove all other structures $S_j^x$ for j ≠ i

    Re-evaluate the similarity between clusters

  until the cluster has size $\left|\widehat{P}.\mathsf{\text{LX}}\right|$

In order to efficiently select representative structures, we use a similar method to Unweighted Pair Group Method with Arithmetic Mean (UPGMA), an agglomerative hierarchical clustering technique [33]. Each object (i.e. secondary structure) starts in its own cluster. The closest pair of clusters is selected and merged into a single cluster as one moves up the hierarchy. The distance between clusters is measured using arithmetic mean defined in UPGMA. Compared to the standard clustering procedure, we have constraints on the objects that can be selected into the same cluster. Given the shape, only structures that have shape $\widehat{P}$ and come from different ncRNAs can be combined in the same cluster. The detailed clustering process is described in Procedure 1.

During clustering, the structural and sequence similarity is evaluated using grammar string-based approach [34, 35]. Grammar strings encode both secondary structure and sequence information for an ncRNA sequence. Grammar string alignment score can accurately quantify the structural and sequence similarity of two ncRNAs. In addition, grammar string can encode pseudoknot structures [34, 35]. For a sequence X _i and one structure $S_j^i$ in the folding space of X _i , X _i and $S_j^i$ are encoded in a grammar string $g{s}_j^i$. We measure the similarity between any two grammar strings using the normalized grammar string-based alignment score over the alignment length. The similarity between groups of grammar strings is measured by arithmetic mean in UPGMA.

Figure 4 sketches the representative structure selection based on clustering procedure. Let $g{s}_j^i$ be a grammar string converted from X _i and $S_j^i$, X _i ∈ P.LX. Once $g{s}_j^i$ is selected, all the other grammar strings derived from the folding space of X _i will be removed from further processing.

The progressive MSA is performed on the set of representative structures using the clustering path as a guide tree. We then derive the consensus secondary structure from the alignment. The consensus structure can be mapped to each aligned sequence to accomplish the predicted structure of an individual sequence.

Data sets

The training data set is the K10 from BraliBASE [36]. It contains 845 sequence sets, each of which has 10 homologous ncRNAs. There are two test data sets. The first one is the K15 from BraliBASE. K15 contains 503 sequence sets, each of which has 15 homologous ncRNAs. As existing shape prediction tools are not designed for handling pseudoknots, we use the pseudoknot-free sequence sets in K15 to compare the performance of shape prediction. After removing the sets containing pseudoknots, we have 452 sequence sets left. To test the performance of pseudoknot prediction, we constructed the second test set R15 from psuedoknot families of Rfam [3]. In Rfam 10.0, there are 71 families containing pseudoknots. 25 of them have published structures. Of the 25 families, only families with at least 15 seed sequences are used for testing our tools. For each chosen family, sets of 15 sequences are chosen randomly to construct the test sets. Finally R15 contains 160 test sets. The average pairwise sequence identities range from 60-93%. For all sequence sets, the reference shapes were derived from Rfam [3].

SVM training

For both the training and testing data sets, we need to apply de novo folding tools to the sequences. We choose a folding tool using the following criteria. First, this tool is able to output both the optimal and sub-optimal structures. Second, this tool has high accuracy and can be efficiently applied to a large number of ncRNAs. Finally, if the target sequences contain pseudoknots, this tool should be able to output pseudoknot structures. As a result, we chose TT2NE [20]. Different from many other pseudoknot prediction tools that have constraints on the type of the pseudoknot, TT2NE is more exible about the types of the target sequences. However, when it was applied to K10, TT2NE failed to output structures for some sequences due to the length limit (200 nt) and also existence of IUPAC characters in some sequences. Thus, for the training data set K10, we applied quikfold [37] because K10 rarely contains pseudoknots. Although it is ideal to use the same folding tool to the training and testing data set to achieve optimal classification performance, the complexity of the training and test data sets together with the performance of de novo folding tools lead to the current combination. In the Discussion section we will briefly discuss how de novo folding tools affect the performance.

We employed the SVM model implemented by LIBSVM tool [38] for classification. For each sequence in K10, we applied quikfold with the energy range 5% to obtain both optimal and sub-optimal structures of each sequence. The predicted structures were grouped based on their corresponding shapes. Associated features were extracted and enclosed with each shape. We normalized feature values to fit the different properties of test sets to the same scale.

To label shapes, we used the shapes extracted from the consensus structures in Rfam [3] as the reference. Shapes are labeled according to their correctness. We label a shape as 1 if it is as same as the reference shape. Otherwise, it is labeled as -1.

Shape prediction comparison

Structure prediction comparison

TP is the number of correctly predicted base pairs. FN is the number of base pairs that are in the reference structure but not in the predicted structure. FP is the number of base pairs that are in the predicted structure but not in the reference structure. Figure 5 is the boxplot of the sensitivity and PPV over all 160 test sets. KnotStructure has the best overall performance on the whole data set. The median values of sensitivity and PPV are 54.55% and 46.15% for KnotStructure. Hxmatch has the next highest sensitivity and PPV (42.11% and 42.86% respectively). The abstract shapes of these families are shown in Table 2. Three families contain simple H-type pseudoknots while the other three families contain more complicated pseudoknots. In order to show the effect of shape prediction in structure prediction, we predicted the structures of R15 using 10 randomly selected shapes. The average sensitivity and PPV of predicted structures with the predicted shapes are higher than those using random shapes as shown in Table 3. Table 4 shows the average sensitivity and PPV over all sequences of each family compared to other tools. The average running time of KnotStructure on each family compared to other tools is shown in Table 5.