Accurate multiple sequence-structure alignment of RNA sequences using combinatorial optimization

Tree-based models

Tree-based structural alignment algorithms view an RNA secondary structure as a tree. Depending on the particular model (either tree-editing [19] or tree alignment [20]), one either searches for the minimal number of operations (node inserting, node deletion, and node substitution) to transform one tree into the other, or into a common supertree. Algorithms employing the model from [20] have time complexities in O(n⁴), thus making the computation expensive. Here and in the following, n denotes the size of the longest sequence. Tree-alignment algorithms have complexities that are on average only slightly worse than conventional sequence alignment. More precisely, their running time is in O(n²·Δ²), where Δ denotes the maximum number of branches of a multiloop in the input structures.

A tool that builds upon the tree paradigm is RNAFORESTER [21]. It computes multiple structure-to-structure alignments of RNA sequences by performing tree-alignment in a progressive fashion.

Annotated sequences

We call a sequence that is augmented by structural information an annotated sequence. Classical dynamic programming (DP) algorithms can be extended to annotated sequences. The DP solution for the structure-to-structure and structure-to-unknown problem then typically requires O(n⁴) and O(n³) in time and space, respectively. Bafna, Muthukrishnan, and Ravi describe an algorithm that simultaneously aligns the sequence and secondary structure of two RNA sequences [22]. Their method runs in time O(n⁴), which still does not make it applicable to instances of realistic size. Eddy [23] proposes an algorithm that reduces the memory consumption to O(n² log n). The STRAL tool [24] uses the values of the base pair probability matrices, as given by the partition function [18], to compute the maximal pairing probability of a single nucleotide and to align the sequences in a CLUSTALW-like fashion.

In the restricted structure-to-structure scenario, one can resort to more sophisticated edit-models like the one proposed by Jiang in [25] where the authors specify operations both on the sequence and the structure level. The dynamic programming algorithm is in O(n⁴), making the computation rather tedious for longer sequences. A program that implements the Jiang model is MARNA [26]: it computes pairwise sequence-structure alignments, but is additionally able to compute multiple alignments. To this end, MARNA computes all pairwise structural alignment and uses T-COFFEE to compute the actual multiple alignment incorporating the structural information of the pairwise alignments.

The unknown-to-unknown scenario requires the simultaneous computation of the alignment and consensus structure. The computational problem of simultaneously considering sequence and structure of an RNA molecule was initially addressed by Sankoff in [27], where the author proposed a DP algorithm to align and fold a set of RNA sequences at the same time. The CPU and memory requirements of the original algorithm are O(n^3k) and O(n^2k), respectively, where k is the number of sequences and n is their maximal length. Current implementations modify Sankoff's algorithm by imposing limits on the size or shape of substructures, e.g., DYNALIGN [28, 29], or FOLDALIGN [30] that combine a sliding window and banded alignment approach. Hofacker, Bernhart, and Stadler [31] have presented the PMMULTI software to align base pair probability matrices. Their recursions are essentially the same as the ones given by Sankoff in [27] and subsequently used for sequence-structure alignment by Bafna et al. in [22] with the only difference that they consider probabilities instead of fixed structures. By banding the range of possible alignment positions they bring the time and space complexity of the pairwise case down to O(n⁴) and O(n³), respectively. For the multiple case, they align consensus base pair probability matrices in a progressive fashion. Similar in spirit are FOLDALIGNM [32] or LOCARNA [33], two recent reimplementations of the PMMULTI approach. FOLDALIGNM provides both several restrictions on the alignment and a two-stage procedure to fill the DP matrix: this further reduces the running time to O(n²δ²) where n is the length of the longer sequence and d is the maximal length difference of the alignment of two subsequences. LOCARNA on the other hand takes advantage of the sparse base pair probabilities matrices to reduce the running time.

Probabilistic models

Eddy and Durbin [34] describe covariance models for measuring the secondary structure and primary sequence consensus of RNA sequence families. They present algorithms for analyzing and comparing RNA sequences as well as database search techniques. Since the basic operation in their approach is an expensive dynamic programming algorithm, their algorithms cannot analyze sequences longer than 150–200 nucleotides. Therefore, recent approaches reduce the running time by incorporating additional information, e.g. Holmes et al.'s STEMLOC [35, 36] where the authors propose the concept of alignment/fold envelopes that constrain possible alignments. Along these lines, in [37] the authors keep a set of probabilistically derived alignment positions fixed: these alignment positions serve subsequently as anchors for the structural alignment which prune away large parts of the search space. The authors of [38] describe a method based on conditional random fields to align an RNA sequence with known structure to one with unknown structure. They estimate their parameters using conditional random fields and compute the alignment using the recursions from [39].

Graph-based models

Kececioglu [40] has introduced a graph-theoretical model for the classical primary sequence alignment problem. In [41] the authors present a first extension of this model to RNA structures and propose a branch-and-cut approach based on an integer linear programming formulation. Based on this formulation and inspired by the successful application of Lagrangian relaxation by Lancia and Caprara [42] to the related contact map overlap problem, in [43] the authors switch from branch-and-cut to the Lagrangian relaxation technique. They are able to solve instances a magnitude larger by simultaneously reducing the running time significantly. In [44] the authors give a graph-theoretic model for the computation of multiple sequence alignments with arbitrary gap costs. In the next section we will combine the formulations given in [43] and [44], resulting in a novel graph-based formulation for sequence-structure alignment with arbitrary gap costs.

Note that the graph-based model naturally deals with all three alignment scenarios. In addition, unlike other algorithmic approaches, the graph-based algorithms do not restrict the input in any way and hence can handle arbitrary pseudoknots: Pseudoknots have been shown to play important roles in a variety of biological processes, see [45] for a recent review. Most DP-based algorithms assume nested secondary structures to compute subproblems efficiently. Few exceptions exist, for example [46], but these algorithms are always restricted to certain classes of pseudoknots (like H-type pseudoknots) and do not handle the general case.

Problem definition

Figure 4 gives a toy example showing two annotated sequences and two possible alignments, one maximizing the score of sequence and structure, and the other one just the sequence score alone. This problem definition comprises the one addressed by Bafna et al. in [22]: Our model, however, also allows tertiary elements, which is not covered by their recursions.

Basic model

Let s = s₁, ..., s _n be a sequence of length n over the alphabet Σ = {A, C, G, U}. A pair (s _i , s _j ) is called an interaction if i < j, and nucleotide i pairs with j. In most cases, these pairs will be Watson-Crick or wobble base pairs. The set p of interactions is called the annotation of sequence s. Two interactions (s _k , s _l ) and (s _m , s _o ) are said to be inconsistent, if they share one base; they form a pseudoknot if they "cross" each other, that is, if k <m <l <o or m <k <o <l. A pair (s, p) is called an annotated sequence. Note that a structure where no pair of interactions is inconsistent with each other forms a valid secondary structure of an RNA sequence, possibly with pseudoknots.

We are given two annotated sequences (s ^A , p ^A ) and (s ^B , p ^B ) and model the input as a structural graph G _S = (V, L). The set V denotes the vertices of the graph, in this case the bases of the sequences, and we write $v_i^A$ and $v_i^B$ for the i th base in sequence A and B, respectively. The set L contains undirected alignment edges between vertices of sequences A and B, for sake of better distinction called lines. A line l ∈ L with l = ($v_k^A,v_l^B$) represents the alignment of the k-th character in sequence A with the l-th character in sequence B. By s(l) and t(l) we refer to the adjacent vertices of line l in sequence A and B, respectively. A subset $ℒ$ ⊂ L represents a valid sequence alignment of sequence A and B, if there are no two lines k, l ∈ $ℒ$ such that k and l cross or touch each other [40]. Crossing or touching lines induce ordering conflicts in the alignment (see Fig. 5 for an illustration). We denote with the set ${\mathcal{C}}_L$ the collection of all maximal sets of mutually conflicting lines.

We extend the original graph G _S = (V, L) by the edge set I to model the annotation of the input sequences in our graph. Consequently, we have interaction edges between vertices of the same sequence, i.e., an edge ($v_i^A,v_j^A$) representing the interaction between nucleotides i and j in sequence A. Figure 6 illustrates these definitions by means of an example. Note that at this stage gaps are not modelled in our formulation. Hence, we have to extend our model to incorporate gap penalties in our model.

Gap edges

The initial model containing only lines (the set L) and interaction edges (the set I) is augmented by a set of gap edges G, which represents gaps in the alignment. For sake of compactness, we just describe the gap edges of sequence A, the gap edges of sequence B are defined analogously: We have an edge $e_{kl}^A$ from $v_k^A$ to $v_l^A$ with k, l ∈ 1, ..., |s ^A | representing the fact that no character of the substring $s_k^A\mathrm{...}{s}_l^A$ is aligned to any character of the sequence B, whereas $s_{k-1}^A$ (if k - 1 > 0) and $s_{l+1}^A$ (if l + 1 ≤ |s ^B |) are aligned with some characters in sequence B. We say that $v_k^A,\mathrm{...},v_l^A$ are spanned by the gap edge $e_{kl}^A$. Figure 7 shows the graph extended by gap edges.

Two gap edges $e_{kl}^A$ and $e_{mn}^A$ ∈ G are in conflict with each other if {k, ..., l + 1} ∩ {m, ..., n} ≠ ∅, that is, if they overlap or touch. This is intuitively clear, because we do not want to split a longer gap into two separate gaps: Consequently, there has to be at least one aligned character between any two realized gap edges. See Fig. 8 for an example. We denote with the set ${\mathcal{C}}_G$ the collection of all maximal sets of mutually conflicting gap edges. Finally, we define $G_{v_k^A\leftrightarrow v_l^A}$ as the set of gap edges that span the nodes $v_k^A,\mathrm{...},v_l^A$.

Interaction match

We call two interactions $(s_k^A,s_l^A)\in p^A$ and $(s_m^B,s_n^B)\in p^B$ an interaction match if there exist two alignment edges $a=(v_k^A,v_m^B)$ and $b=(v_l^A,v_n^B)$ that do not cross each other. We say that a subset $\mathcal{S}$ ⊆ L realizes the interaction match if {a, b} ⊆ $\mathcal{S}$. Interaction matches realized by a set $\mathcal{S}$ represent common interactions that are preserved by aligning the begin and end nucleotides of the interaction. Figure 9 illustrates the definitions.

Gapped structural trace

Figure 10 visualizes these properties by showing a toy example for a gapped structural trace.

We assign weights w _l and w _kl for each line l and interaction match (k, l) that represents the benefit of realizing l or (k, l). By default, we set these scores along the lines of standard scoring methods, e.g., BLOSUM matrices for the weight of the lines, base pair probabilities [18] for the interaction match scores, or by using the RIBOSUM scoring matrices derived from alignments of ribosomal RNAs [47]. Our model, however, is not limited to standard scoring schemes. Since we can set each (sequence or structure) weight separately, the user can assign completely arbitrary scores to each line or interaction match which makes the incorporation of expert knowledge into the computation of structural alignments easy. Furthermore, we assign negative weights to gap edges $a_{kl}^A$ with representing the gap penalty for aligning substring $s_k^A,\mathrm{...},s_l^A$ with gap characters. Note that the model allows for arbitrary, position-dependent gap scoring.

Approaches for traditional sequence alignment aim at maximizing the score of edges in an alignment $ℒ$. Structural alignments, however, must also take the structural information encoded in the interaction edges into account. The problem of structurally aligning two annotated sequences (s ^A , p ^A ) and (s ^B , p ^B ) corresponds to finding an alignment such that the weight of the sequence part (i.e., the weight of selected lines plus gap penalties) plus the weight of the realized interaction matches is maximal. More formally, we seek to maximize ${\displaystyle {\sum }_{l\in ℒ}{w}_l}+{\displaystyle {\sum }_{g\in \mathcal{G}}{w}_g}+{\displaystyle {\sum }_{(i,j)\in ℐ}{w}_{ij}}$, where ($ℒ$, $\mathcal{G}$) represents an alignment with arbitrary gap costs, and $ℐ$ contains the interaction matches realized by $ℒ$. Observe that this graph-theoretical reformulation matches the problem statement given at the beginning of this section.

Biological aspects

The basic entities of our model are the alignment, interaction, and gap edges in the structural graph, which contribute to the objective function rather independently. Hence, one could argue that the model does not capture important features of RNA structures, like the incorporation of stacking energies or loop scores that depend on the actual size of the loop. We are aware of these limitations.

Nevertheless, the results of our computational experiments presented in Sect. 3 show that this approach yields high-quality structural alignments. In the pairwise case, our graph-based model is competitive with state-of-the-art approaches and develops its strength with an increasing number of sequences, outperforming all other programs that we tested (for details see Sect. 3). Additionally, the authors of [48] showed that models that do not capture stacking energies and loops are still competitive.

Beyond, our graph-based approach offers the possibility to change the model from nucleotides as the working entities to stems: Instead of taking single nucleotides as the vertices of the structural graph, we could search for candidate stems in the sequences and introduce a vertex for each half-stem. This would allow us to incorporate energy-based scoring into our model, which then, however, will have to be adapted to take into account overlapping stem candidates.

Lemma 2.1 (Proof in [16])

A feasible solution to the ILP (1)–(8) corresponds to a valid gapped structural trace of weight equal to the objective function and vice versa.

Observe that constraints (2)–(6) exactly correspond to the properties of a gapped structural trace as described in Sect. 2.1.

In [49] the authors show that the problem of computing an optimal gapped structural trace is already NP-hard, even without considering gap costs. Hence, we cannot hope to find an optimal solution to the problem in polynomial time.

Commonly used mathematical programming techniques for NP-hard problems therefore resort to various relaxation techniques that are the basis for further processing. A relaxation results from the removal of constraints from the original ILP formulation, and is often solvable in polynomial time. A popular relaxation is the so called LP relaxation where the integrality constraints on the variables are dropped, yielding a standard linear program, for which solutions can be found efficiently.

Another possible relaxation technique is Lagrangian relaxation: Instead of just dropping certain inequalities, we move them to the objective function, associated with a penalty term that becomes active if the dropped constraint is violated. By iteratively adapting those penalty terms using, for instance, subgradient optimization, we get better solutions with each iteration. A crucial parameter is therefore the number of iterations that we perform: the higher the number, the more likely it is to end up with an optimal or near-optimal solution.

Inspired by the successful approach of Lancia and Caprara for the contact map overlap problem, we consider the relaxation resulting from moving constraint (7) into the objective function.

Lemma 2.2 (Proof in [16])

The relaxed problem is equivalent to the pairwise sequence alignment problem with arbitrary gap costs.

Score vs. alignment accuracy

We were interested to what extent the accuracy of our alignments correlates with the actual BPP score that we computed. Since the score depends on the length of the input sequences, we normalized the score with respect to the number of paired bases in the minimum free energy structure. Note that we did not use the actual structure, but the number of base pairs in the structure to get a rough estimate of how many pairings we expect in the structure. Then, let $\widehat{p}$ and n be the average score and the number of base pairs in the MFE structure, then the base-pair normalized score is given by $\widehat{p}$/n. The left side of Fig. 13 shows the results for all 189 k 10 instances with an average pairwise sequence identity less than 50%. The great majority of instances behaves as expected: the higher the bp-score is, the better is the corresponding Compalign score: There is, however, a group of 10 outliers (represented by the red boxes). Although they have a high bp-score (greater than 10.0), the alignment accuracy is bad: it turned out that these 10 instances are all SECIS-elements, indicating that the BPP scoring scheme is not appropriate for this group.

Furthermore, we assumed that there should a correlation between the actual performance of our algorithm and, again, the quality of our alignments: Remember that each Lagrange iteration results in a new valid solution and a new upper bound for the problem instance. Dividing the value of the highest lower bound by the value of the lowest upper bound gives an optimality ratio, i.e., a measure of how close the best solution is to an optimal one. Assuming an inverse correlation between the gap between lower and upper bound and the quality of the alignment, we again took all k 10 BRALIBASE instances of low pairwise sequence identity and computed the arithmetic mean of the optimality ratios of all pairwise alignments. The right side of Fig. 13 shows the plot for all 189 k 10 instances with a sequence similarity lower than 50%. Most of the instances behave as expected: the higher the average optimality ratio is, the closer is the computed alignment to the reference alignment (and vice versa). There is, however, a group of 19 instances that behave differently (marked as red boxes in Fig. 13): Although their average optimality ratio is high (> 0.7), the corresponding Compalign value is rather low compared to instances of a similar average optimality ratio. A closer inspection revealed that all instances of the upper left corner (that is instances having a Compalign value lower than 0.65 and an average optimality ratio of greater than 0.7, represented by red boxes in Fig. 13) comprises almost all instances of either bacterial SRP RNAs or SECIS elements (just one SRP RNA instance is not among the 19 instances). We therefore increased the number of iterations for one SECIS instance to see whether this would positively influence the quality of the alignment. By setting the number of iterations to 500, 1000, and 2000 we got average optimality ratios of 0.83, 0.85, and 0.87, by simultaneously yielding Compalign values of 0.39, 0.38, and 0.36, respectively. Obviously, the better the computed alignments in terms of the optimality ratio are, the worse they got with respect to the reference alignment.

Consequently, for the outlier instances described above, we changed the scoring from BPP to RIBOSUM scores. Figure 14 shows the change in terms of the Compalign score and optimality ratio for the 19 outlier instances: 16 instances had better Compalign scores by using the RIBOSUM scoring, whereas the optimality ratio decreased in the majority of instances.

In general, however, our experiments showed that RIBOSUM scoring is not superior to BPP scoring (at least for the BRALIBASE benchmark and LARA): Figure 15 shows a comparison of all low homology k 5 instances using either base pair probability matrices or RIBOSUM scoring, and it is obvious that base pair probability scoring yields better results on these input instances.

Comparison of running times

We were especially interested in how the running times of the programs that use structure information scaled with respect to the number of the input sequences: FOLDALIGNM is a progressive approach which computes (n - 1) pairwise alignments given n input sequences. MARNA and LARA, however, compute all $\frac{n(n-1)}{2}$ pairwise alignments. Figure 18 shows the execution time of all five programs on all k 2, k 3, k 5, k 7, k 10, and k 15 instances. As one can see, with an increasing number of input sequences, a progressive alignment strategy pays off compared to the computation of all pairwise alignments.