Workflow and algorithms

RNAsalsa implements a workflow that makes use of several well-established algorithms for both RNA secondary structure prediction and structure enhanced alignment (Figure 1). It consists of three major stages: (i) the determination of structural constraints for each input sequence that emphasize the common features; (ii) the computation of a detailed structural model for each input sequence; and (iii) the construction of a final multiple sequence alignment together with a consensus structure.

The starting point for RNAsalsa is an initial alignment 𝔸⁰ of a collection {x¹, …, x^N} of homologous RNA sequences (produced e.g. simply by clustalW) and an a priori secondary structure constraint σ for a single sequence x⁰ which is contained in the alignment 𝔸⁰. The sequence x⁰ and its structural model are used only to initialize the structure prediction and alignment process.

In the first step, RNAsalsa checks the consistency of the initial alignment 𝔸⁰ and the initial constraint σ⁰: for each base pair in σ⁰, we evaluate whether the corresponding aligned positions of a sufficient number of sequences in 𝔸⁰ can also pair. If so, we retain the base pair, otherwise it is removed from the constraint. The resulting ‘relaxed’ constraint

can be seen as base pair-wise filtering of the initial constraint σ⁰ that removes pairs from σ⁰ that are largely inconsistent with the initial alignment. Projecting the relaxed constraint σ separately onto each aligned sequence (i.e. retaining only the canonical base pairs of σ that can be formed by the input sequence xⁱ ), then produces distinct initial structure constraints for each sequence: Up to this point, the result heavily depends upon the initial alignment. It may not cover the input sequences uniformly, in particular it will often be concentrated on the well-conserved (and therefore properly aligned) regions.

In the second step, RNAsalsa utilizes the improved accuracy of the predicted consensus structures. To this end, we construct a collection of pairwise sequence alignments 𝔸^ij from the input sequences xⁱ and x^j. These can be constructed in different ways, either by dynamic programming alignment, or by projecting the corresponding sub-alignment of 𝔸⁰. For details we refer to the RNAsalsa manual, which is part of the Supplementary Data. For each of the pairwise alignments 𝔸^ij we compute the consensus minimum free-energy structures

using the base pairs common to the projected structures σⁱ and σ^j, respespectively, as constraint. This step uses the Vienna RNA Package library functions underlying RNAalifold (²²) to perform the constrained folding computations.

For each sequence xⁱ, the collection of structures {τ^ij, i ≠ j} taken together defines a set of base pairs on xⁱ that are both thermodynamically plausible and conserved in at least one other sequence of the input set. From this set of pairs, we select a single secondary structure

for sequence xⁱ using a majority voting procedure. RNAsalsa currently implements a simple greedy procedure that selects the most frequent base pairs first and rejects pairs that would cross previously selected ones to avoid the formation of pseudo-knotted structures. Alternatively, one could also use Nussinov's Maximum Circular Matching algorithm (³⁴) to retrieve a maximum weight sub-set of non-intersecting pairs. The base pairs of τⁱ, which by construction typically contain most of the initial constraint-derived pairs σⁱ, are now used as a constraint for computing the final secondary structure prediction for each input sequence xⁱ. Again, we employ the Vienna RNA Package library here.

The purpose of the entire—rather complex and computationally expensive—procedure is to use as much information as possible in guiding the last step, the computation of the secondary structure models ψⁱ for each input sequence. This guiding information is derived from two sources: the initial structural constraint σ and the ensemble of plausible base pairs generated from all pairwise alignments.

In the next step, the sequence–structure pairs (xⁱ, ψⁱ) are realigned. To this end, RNAsalsa uses a hierarchical progressive alignment based on pairwise dynamic programming alignments with affine gap costs (³⁵). The scoring function explicitly incorporates the secondary structure annotation: the (mis)match score s(x_i, y_j) of position i from sequence x with position j from sequence y is defined as follows:

where x_π(i) and y_π(j) denote the pairing partners of x_i and y_j in their respective secondary structures. The coefficient b₀ = 1 if both x_i and y_j are paired nucleotides. The coefficient b₁ is set to 1 if x and y share sufficient structural conservation to a certain extent that overcame the precedent filtering steps and if x_π(i) and y_{π( j)} are located either both upstream or both downstream of x_i and y_j, respectively. Otherwise the structural contribution is ignored, b₁ = 0. Finally, if one x_i or y_j are unpaired, then b₀ = b₁ = 0 and c = 1. In regions without structural information we therefore use a pure nucleic acid sequence score s_p, where as in structured regions, the modified scoring functions s_m and s_n are used. For instance, within trusted structural regions A-G is scored as a match because both may pair with U, while it is not in regions without sufficiently trusted structural information. Default scoring tables are listed in the manual and Supplementary Data. The final result is a global re-alignment 𝔹 of the input sequences which integrates all the secondary structure information obtained in the previous steps.

The individual folds ψⁱ and the alignment 𝔹 are used to derive a consensus structure

Since we now have the trusted alignment 𝔹, we can again employ a simple voting strategy: we start from the set of all base pairs that appear sufficiently often in superposition of the ψⁱ. Again, we use a greedy strategy to avoid conflicting base pairs (note that no conflicts can arise if we consider only base pairs that occur at least N/2 times).

Several parameters can be adjusted in the process. In particular, the stringency of the initial filtering of base pairs, Equation (1), and the two majority voting procedures, Equations (4) and (7) can be adjusted by the user to the peculiarities of the data sets. In each case, a threshold for the minimum number of consistent pairs can be specified. Some guidelines for practical use can be found in the RNAsalsa manual.

Secondary structure prediction

Performance of RNAsalsa's structure predictions was evaluated in comparison with three other relevant methods: MXSCARNA (¹⁵) computes pairing probabilities and considers potential stem information in the subsequent alignment process; RNAfold (¹⁷) produces individual secondary structures of RNA sequences; and RNAalifold (²²) generates the consensus structure for a given input alignment. We compared RNAfold predictions with RNAsalsa's individual predictions ψⁱ (Equation 5), while the MXSCARNA and RNAalifold results are compared with RNAsalsa's final consensus structure ω. As a reference model we employed the mammalian 16S rRNA secondary structure adopted from (³⁶) (Figure 2).

The RNAsalsa secondary structure model for the mammalian 16S rRNA sequences is highly congruent to the Bos taurus reference model proposed by Burk et al. (³⁶), see Supplementary Data for an illustrating graphical representation. In particular, 44 of the 52 helices within the conserved core of the structure are correctly predicted. Much of the remaining discrepancy can be explained by differences in the stringency of rules (Figure 2A). While the literature reference was derived with very conservative rules, we allow more variability among the organisms, and accept stems even if there is no direct evidence from compensatory substitutions. MXSCARNA and RNAalifold capture only 27 and 23 helices, respectively. In contrast to RNAsalsa, they failed to detect in particular long-range interactions.

Furthermore, to demonstrate the capabilities of constraints in thermodynamic folding algorithms, we compared the results of unconstrained foldings by the RNAfold software (¹⁷) only with those generated by RNAsalsa. In both cases, we used the same thermodynamical parameter sets and algorithms as implemented in RNAfold. The only difference is the usage of folding constraints in the case of RNAsalsa. These structure constraints are automatically generated and optimized for each RNA sequence within the data set.

RNAsalsa's predictions of the individual structures always use an individual structure constraint that was optimized by RNAsalsa's procedure itself up to that point. Therefore, they can match the reference model much better than thermodynamic folds by RNAfold, which have been calculated without supporting constraints (Figure 2B). See Supplementary Data for a set of comparisons performed on 16S rRNA.

Structure-aware RNA sequence alignments

The impact of secondary structures on the alignments as well as the overall performance was investigated by comparison with two commonly used sequence alignment methods, the classical ClustalW (³⁸) and the more modern MAFFT (³⁹) approach, and with the structural alignment method MXSCARNA. For this purpose we employed both simulated and real data sets.

Manually curated alignments for rRNAs are a rare commodity. At present, only partial data sets are available. Here, we used the alignments of 26 mammalian mitochondrial rRNAs published in (⁴⁰,⁴¹). We only compared the subset of alignment positions that are marked as ‘trusted’ in these alignments. For the 16S rRNAs, we find that all four alignment programs reproduce 97.5–99% of the ‘trusted’ pairwise (mis)matches. The two structural alignment programs, MXSCARNA and RNAsalsa, have a slightly increased fraction of ‘trusted’ (mis)matches, Table 1. Very similar results were observed for the 12S mitochondrial rRNAs from the same source.

The relative quality of alignments was assessed following the procedures outlined in (⁴²,⁴³). In particular, we calculated the sum of pairs score (SPS) and the total column score (TC) as implemented in the baliscore software (⁴²), as well as the structure conservation index (SCI) (⁴⁴). Both SPS and TC are similarity metrics for multiple alignments. While SPS evaluates the number of equivalently aligned residue pairs summed over all pairs of sequences, TC counts the number of columns exactly shared by the two alignments. When one of the alignments is used as a benchmark reference, both SPS and TC can be interpreted as measures of alignment quality. The SCI, on the other hand, measures to what extent the aligned sequences form a common secondary structure in such a way that the sequence alignment also represents a structural alignment. Originally introduced as a measure of structural conservation given the alignment, it was used in (⁴⁵) to measure the ability of alignment algorithms to reconstruct conservative consensus folds. Of course, the SCI can only be employed to measure alignment quality when the existence of a common global consensus structure can safely be assumed.

Simulated data were generated as reference alignments using the new software rnasim (⁴⁶). The tool takes both a structure annotated sequence string and a pre-defined phylogenetic tree as input and simulates the evolution of RNA molecules with fitness constraints derived from thermodynamic stability along the given tree. The tool provides us with the true reference alignment of the leaf sequences at the end of each simulation run. Several sets of examples representing smaller and larger RNAs (tRNAs, mammalian 12S rRNAs, and the Saccharomyces cerevisiae 28S rRNA) were used as roots. We used various values between 1 and 100 000 for rnasim's ‘branch scaling’ option -s, whereas the default was used for all other options. Since the branch scaling parameter has a strong influence on both sequence and structure, we systematically analysed how various alignment algorithms behave as a function of divergence.

Table 2 summarizes the results of some of the RNAsim-derived data sets. For the tRNAs, all programs except clustalW perform similarly well. In case of the 28S rRNA data set, which constitutes a very good example with highly divergent sequences, RNAsalsa performs slightly better than its competitors, albeit on a low level. These examples illustrate a general pattern: most programs, including RNAsalsa, perform with comparable scores on a wide range of test cases.

As a function of the branch length scaling parameter, we observe that all alignment algorithms perform comparably well or poorly in terms of the SPS, TC and SCI values. As expected, when the alignment problem becomes harder, i.e. for large scaling values, the performance metrics decrease uniformly (Supplementary Data). The resulting alignment then becomes more and more divergent between different methods.

We also used the BRAliBase-II set of structural alignments (http://projects.binf.ku.dk/pgardner/bralibase/bralibase2.html) (⁴²), which covers group II introns, 5S rRNA, tRNA, U5 and SRP RNA [following (⁴²), we did not use the SRP RNA alignments]. Again, we observe no major differences in the performance metrics (Supplementary Data). As in the case of the simulated data, the alignments become more divergent between methods as the problems become harder.

It appears that they tend to deviate from the reference in different features. RNAsalsa avoids poorly supported (mis)matches and instead tends to introduce more gaps in regions where the alignment becomes ambiguous. While this feature does not positively influence the scores in usual performance metrics, it does have a desirable effect for phylogenetic reconstruction, because gap-rich regions tend to be down-weighted or even excluded in typical phylogenetic applications. In that sense, RNAsalsa apparently reduces conflicting information arising from spurious (mis)matches. The RNAsalsa alignments emphasize the well-supported regions, and—as we have seen in the previous section—lead to quite accurate assignments of the most conserved secondary structure elements.

Exemplary applications in phylogeny reconstruction

In order to demonstrate the usefulness of RNAsalsa in phylogenetic applications, we consider two distinct data sets in detail. For each alignment, we performed phylogenetic analyses using a likelihood-based approach and compared the results with the published analyses of the data set. After individually aligning the LSU and SSU sequences, the alignments were concatenated. Then the program Aliscore (⁴⁷), a method to identify ambiguously aligned regions in multiple sequence alignments, was used to extract the informative parts of the alignment. Aliscore identifies ambiguously aligned regions in multiple sequence alignments based on comparisons of aligned sequences in a sliding window and a Monte Carlo (MC) approach. The MC re-sampling compares the score of the originally aligned sequences in a given window position with scores of randomly drawn sequences of similar character composition. Maximum likelihood (ML) analysis was conducted using the RAxML software under a GTR+Γ model set up. For further details we refer to the Supplementary Data.