Fast and accurate clustering of noncoding RNAs using ensembles of sequence alignments and secondary structures

Ensemble of all possible sequence alignments

To measure the similarity between two RNAs, one common approach is to perform pairwise alignment, and to calculate its alignment score. The Sankoff algorithm simultaneously models sequence alignments and secondary structures, and is extremely time-consuming. Therefore, we first approximate the Sankoff algorithm by the SW algorithm that only models sequence alignments apart from secondary structures. Although this is a strong approximation, we attempt to improve the reliability by utilizing all possible sequence alignments rather than one optimal sequence alignment.

For an RNA sequence x, we denote its length by |x|. For each position 1 ≤ i ≤ |x| in x, we denote the nucleotide by x _i ∈ {A, C, G, U}.

For two sequences, x and y, let Π _xy be the set of all possible sequence alignments in the SW algorithm. Let π _xy denote one particular sequence alignment in Π _xy .

We calculate the similarity between x and y by accumulating the alignment score of π _xy over Π _xy . For this purpose, we employ local alignment (LA) kernel [16] defined as follows:

(1)

where β ≥ 0 is a parameter, and Score(π _xy ) is the alignment score of π _xy under a given scoring scheme (gap penalties and match scores). In practice, we take the logarithm of LA kernel, and similarity values are normalized to range from 0 to 1:

(2)

LA kernel (1) can be computed by the variant of the SW algorithm as follows:

Initialization:

for i ∈ {0,…, |x|} and j ∈ {0, …, |y|} do

M(i, 0) = I _X (i, 0) = I _Y (i, 0) = T _X (i, 0) = T _Y (i, 0) = 0

M(0, j) = I _X (0, j) = I _Y (0, j) = T _X (0, j) = T _Y (0, j) = 0

end for

Iteration:

for i ∈ {1,…,|x|} and j ∈ {1,…,|y|} do

M(i, j) = e ^βS _xy^(i,j)(1 + I _X (i – 1, j – 1) + I _Y (i – 1, j – 1) + M(i – 1, j – 1))

I _X (i, j) = e ^βg M(i – 1, j) + e ^βd I _X (i – 1, j)

I _Y (i, j) = e ^βg (M(i, j – 1) + I _X (i, j – 1)) + e ^βd I _Y (i, j – 1)

T _X (i, j) = M(i – 1, j) + T _X (i – 1, j)

T _Y (i, j) = M(i, j – 1) + T _X (i, j – 1) + T _Y (i, j – 1)

end for

Termination:

K(x, y) = 1 + T _X (|x|, |y|) + T _Y (|x|, |y|) + M(|x|, |y|)

where the parameters g and d are the penalties for gap opening and gap extension, respectively, and S _xy (i, j) is a scoring function for matching the i-th position in x and the j-th position in y. The design of S _xy (i, j) impacts the performance of the resulting similarity measure, and will be described later.

At this point, we note that our method can take into account all possible sequence alignments in O(|x||y|) time. If we use the exact Sankoff algorithm instead, it takes prohibitive O(|x|³|y|³) time, which is not practical. In the case of the approximate Sankoff-style algorithms employed in the existing methods, all possible sequence alignments cannot be incorporated to the reduced dynamic programming matrix. Therefore, LA kernel based on the SW algorithm is an efficient way to deal with the ensemble of all possible sequence alignments.

Ensemble of all possible secondary structures

To design a scoring function S _xy (i, j) for LA kernel, we need secondary structures of x and y. As mentioned above, the Sankoff algorithm models secondary structures simultaneously with sequence alignments which we have already modeled by the SW algorithm. Therefore, we next employ the McCaskill algorithm that only models secondary structures apart from sequence alignments. Although this is an additional approximation, we attempt to improve the reliability by utilizing all possible secondary structures rather than one optimal secondary structure.

For an RNA sequence x, let Θ _x be the set of all possible secondary structures. Let θ _x denote one particular secondary structure in Θ _x . We represent a secondary structure as a set of binary variables θ _x = {θ _x (i, j)}_{1≤ i < j ≤| x |}, where θ _x (i, j) = 1 means that the i-th position and the j-th position in x form a base pair. For each position 1 ≤ i ≤ |x| in x, we represent the state of base-pairing using three kinds of binary variable: L _x (i) = ∑ _{j:j > i} θ _x (i, j) = 1 means that a base pair is formed with one of the downstream positions; R _x (i) = ∑ _{j : j < i} θ _x (j, i) = 1 means that a base pair is formed with one of the upstream positions; and U _x (i) = 1 – L _x (i) – R _x (i) = 1 means that the position is unpaired. Given a fixed pair of secondary structures, θ _x and θ _y , we can measure the similarity between the i-th position in x and the j-th position in y using their state of base pairing:

W _xy (i, j|θ _x , θ _y ) = α (L _x (i)L _y (j) + R _x (i)R _y (j)) + s(x _i , y _j ) U _x (i)U _y (j), (3)

where α ≥ 0 is a weight parameter for structural similarity, and s(x _i , y _j ) is a substitution matrix for RNA sequences like the RIBOSUM 85-60 matrix [17]. This scoring function takes a non-zero value in three different cases: it takes α when both of the two positions form a base pair with one of their downstream positions, respectively; it takes α when both of the two positions form a base pair with one of their upstream positions, respectively; and it takes s(x _i , y _j ) when both of the two positions are unpaired.

The McCaskill algorithm defines a probability distribution P(θ _x |x) over Θ _x . The binary variables θ _x (i, j) and {L _x (i), R _x (i), U _x (i)} are converted to the probabilities by taking the expectation over Θ _x . For θ _x (i, j), we obtain a base-pairing probability P _x (i, j) that the i-th and the j-th positions form a base pair:

For {L _x (i), R _x (i), U _x (i)}, we obtain three kinds of probability that the i-th position is paired with one of the downstream/upstream positions, or unpaired, respectively:

We design a scoring function S _xy (i, j) by taking the expectation of (3) over Θ _x and Θ _y :

(4)

The proposed method is obtained by combining the normalized LA kernel (2) with the scoring function (4).

It should be noted that our method can take into account all possible secondary structures in O(|x|³ + |y|³) time, thanks to the McCaskill algorithm. Just as in all possible sequence alignments, the exact Sankoff algorithm results in O(|x|³|y|³) time, and the existing methods cannot incorporate all possible secondary structures. Our method requires O(|x||y|) + O(|x|³ + |y|³) time in total, which is more efficient than the exact Sankoff algorithm. Therefore, our strategy that combines the SW algorithm and the McCaskil algorithm allows to utilize the ensemble information with the reasonable computational cost.

Variations of the proposed method

The scoring function (4) proposed in this study is similar to the scoring function used in BPLA kernel [18, 19]. BPLA kernel is a prediction method that we previously developed for detecting new members of known ncRNA families. Although BPLA kernel was not applied to clustering problems in our previous study, we here clarify its relation to the proposed method. The scoring function used in BPLA kernel is defined as follows:

(5)

where , , and . Therefore, the scoring function (5) can be regarded as a variation of the proposed scoring function (4) with the additional coefficients C ^L , C ^R , and C ^U . These coefficients take large values when the probabilities and are small. That is, BPLA kernel emphasizes the contribution of low-probability (unsure) secondary structures compared to the proposed method. In the next section, we experimentally verify this theoretical implication; the proposed method outperforms BPLA kernel.

Because of the resemblance between the scoring functions, (4) and (5), we set the parameters of the proposed method as used in BPLA kernel: α = 1.0, β = 0.1, g = –27, and d = –0.1

Dataset and experimental system

We compared our method with the state-of-the-art methods developed for the hierarchical clustering of ncRNAs: LocARNA v1.5.2 [8], FOLDALIGN v2.1.1 [11], and Stem kernel v216c [10]. We also performed the experiments with CLUSTALW v1.83 [20] and LA kernel by setting in our method (4).

We can summarize our method and the existing methods as follows. Our method utilizes all possible sequence alignments and all possible secondary structures. LocARNA and FOLDALIGN only use one optimal sequence alignment and one optimal secondary structures. Stem kernel utilizes a subset of all possible sequence alignments and a subset of all possible secondary structures. CLUSTALW and LA kernel ignore secondary structures; CLUSTALW only uses one optimal sequence alignment, while LA kernel utilizes all possible sequence alignments.

We produced three versions of dataset. First, we used ncRNA sequences without modification, and named them the “normal” dataset. Second, we concatenated random sequences to both ends of ncRNA sequences, and named them the “plus flanking regions” dataset. This dataset was intended to simulate the situation where we do not know the exact boundaries of unannotated transcripts. A random sequence was generated from a ncRNA sequence so that it had the quarter length and the same dinucleotide contents. Third, we added false reference clusters, each of which contains one random sequence, and named them the “plus unrelated sequences” dataset. This dataset was intended to simulate the situation where non-functional ncRNAs arise from transcriptional noises. Therefore, we evaluated whether a false reference cluster could be a resultant cluster with a single member. We used a quarter number of false reference clusters compared to true reference clusters. A random sequence was generated from a ncRNA sequence so that it had the same length and the same dinucleotide contents.

We evaluated the overall quality of the cluster tree by the ROC analysis proposed in [8]. (Note that we can obtain different resultant clusters from a cluster tree depending on a distance threshold to cut the branches.) Given a distance threshold, the number of true positives (TP) was defined as the number of sequence pairs that belong to the same reference cluster and are correctly assigned to the same resultant cluster. Analogously, the numbers of false positives (FP), true negatives (TN), and false negatives (FN) are defined, respectively, by counting the pairs from different reference clusters but the same resultant cluster, the pairs from different reference clusters and different resultant clusters, and the pairs from the same reference cluster but different resultant clusters. The ROC analysis was performed by plotting true positive rates TP/(TP + FN) versus false positive rates FP/(TN + FP) for different distance thresholds. The quality of the clustering was measured by the area under the ROC curve (AUC). We measured the total time for computing similarity matrices on a 2.53 GHz Intel Xeon processor.

Quality of the clustering

We first examined the quality of the clustering for the “normal” dataset (Figure 1). Our method achieved the better or comparable AUC to the existing methods in all the range of sequence identity. The accuracy of our method was especially remarkable in the sequence identity range below 60%, where the existing methods resulted in low AUC. This means that our method successfully grouped diverse member sequences in each reference cluster by detecting their remote homology.

Our results can be attributed to the design of each method. The AUC of CLUSTALW and LA kernel, which ignore secondary structures and only use sequence alignments, drastically fell down as the sequence identity decreased. LocARNA, FOLDALIGN, and Stem kernel, which consider secondary structures, kept the AUC relatively moderate in the low sequence identity range. However, their accuracy was still limited when the sequence identity was extremely low (20–39%) because these methods only use one optimal secondary structure or a subset of secondary structures. Our method, which utilizes all possible sequence alignments and all possible secondary structures, achieved the sufficiently high AUC in this region. These results suggest that our design of the similarity measure is effective for identifying a broad range of ncRNA families.

Figure 2 compares an example of the cluster tree between our method and FOLDALIGN in the sequence identity rage of 20–39%. As indicated by AUC, our method produced the more accurate cluster tree than FOLDALIGN, and reconstructed ncRNA families as compact clusters. Although the cluster tree of FOLDALIGN was largely consistent with the references in terms of its topology, boundaries of resultant clusters were quite unclear. In the actual application of hierarchical clustering, we need to choose a proper distance threshold for extracting clusters from a given tree. In this sense, the cluster tree of FOLDALIGN was not sufficient for the practical use. In fact, the previous studies that employed clustering approaches required manual inspection to compensate for ambiguous cluster trees [6, 8, 12]. The cluster tree of our method was much more clear and easier to interpret than the existing methods. These results suggest that our method can reduce human labor costs of clustering approaches, and help to identify novel ncRNAs families.

Next, we evaluated the quality of the clustering for the “plus flanking regions” dataset (Figure 3), and the “plus unrelated sequences” dataset (Figure 4). In both cases, we observed the same tendency as in the results for the “normal” dataset (Figure 1). Our method kept high accuracy in all the range of sequence identity, and achieved the best AUC in the sequence identity range below 60%. These results further support the effectiveness of our method in the practical situations that involve flanking regions and unrelated sequences.

Differences in the variations of the proposed method

As described in Methods, the proposed method has the theoretical advantage compared to BPLA kernel, which can be regarded as a variation of our method. To verify this point experimentally, we compare the proposed method and BPLA kernel using the scoring functions (4) and (5), respectively.

Figure 5 presents the experimental results. The proposed method achieved the slightly better AUC in the sequence identity range below 60%. These results are consistent with the fact that BPLA kernel emphasizes the contribution of unsure secondary structures compared to the proposed method. The proposed scoring function (4) has the theoretical justification as the expectation of the primitive scoring function (3) over all possible secondary structures. Our results provide an experimental verification of the superiority of the proposed scoring function.

Computational cost

In the design of the proposed method, our idea was to improve the reliability of approximate algorithms by the information of suboptimal solutions in their dynamic programming frameworks. Among LocARNA and FOLDALIGN, which only use one optimal solution in their approximate Sankoff-style algorithms, there was a trade-off that LocARNA was faster but less accurate than FOLDALIGN (Figure 1, and Table 2). Stem kernel, which utilizes a subset of solutions in the more approximate Sankoff-style algorithm, partly improved this problem, being faster and more accurate than LocARNA. Our method, which utilizes all possible solutions in the combination of the Smith-Waterman algorithm and the McCaskill algorithm, successfully overcome the trade-off. These results suggest that our strategy is essential to enable fast and accurate clustering of ncRNAs.