DIALIGN-TX: greedy and progressive approaches for segment-based multiple sequence alignment

DIALIGN is a widely used software for multiple alignment of nucleic acid and protein sequences [1, 2] that combines local and global alignment features. Pairwise or multiple alignments are composed by aligning local pairwise similarities. More precisely, pairwise local gap-free alignments called fragment alignments or fragments are used as building blocks to assemble multiple alignments. Each possible fragment is given a score that is related to the P values used by BLAST [3, 4], and the program then tries to find a consistent set of fragments from all possible sequence pairs, maximizing the total score of these fragments. Gaps are not penalized. Here, consistency means that a set of fragment alignments can be included into one single alignment without contradictions, see for example [5] for a more formal definition of our notion of consistency.

The main difference between DIALIGN and more traditional alignment approaches is the underlying scoring scheme or objective function. Instead of summing up substitution scores for aligned residues and subtracting gap penalties, the score of an alignment is based on P-values of local sequence similarities. Only those parts of the sequences are aligned that share some statistically significant similarity, unrelated parts of the sequences remain unaligned. This way, the method can produce global as well as local alignments of the input sequences, whatever seems more appropriate. Combining local and global alignment features is particularly important if genomic sequences are aligned where islands of conserved homologies may be separated by non-related parts of the sequences. Thus, DIALIGN has been used for comparative genomics [6], for example to find protein-coding genes in eukaryotes [7–9].

As with traditional objective functions for sequence alignment, numerically optimal pairwise alignments can be calculated efficiently in the segment-based approach. In DIALIGN, this is done by a space-efficient fragment-chaining algorithms [10, 11]. However, it is computationally not feasible to find mathematically optimal multiple alignments. Thus, heuristics must be used if more than two sequences are to be aligned. All previous versions of DIALIGN used a greedy algorithm for multiple alignment. In an initial step, optimal pairwise alignments are calculated for all possible pairs of input sequences. Since these pairwise alignments are completely independent of each other, they can be calculated on parallel processors [12]. Fragments from these pairwise alignments are then sorted by their scores, i.e. based on their P-values, and then included one-by-one into a growing consistent set of fragments involving all pairs of sequences – provided they are consistent with the previously included fragments.

The greedy algorithm used in DIALIGN is vulnerable to spurious random similarities. It has been shown, that the numerical score of alignments produced by this heuristic can be far below the optimum [5, 13]. Consequently, alternative optimization algorithms have been applied to the optimization problem defined by the DIALIGN approach, e.g. Integer Linear Programming [14].

A set of fragments is called consistent if all fragments can be included into one single alignment, see [15]. Fragments in A are allowed to overlap if different pairs of sequences are involved. That is, if two fragments f₁, f₂∈ A involve sequence pairs s _i , s _j and s _j , s _k , respectively, then f₁ and f₂ are allowed to overlap in sequence s _j . If two fragments involve the same pair of sequences, no overlap is allowed. It can be shown that the problem of finding an optimal consistent set A of fragments is NP-complete (Constructing multiple sequence alignments from pairwise data, Subramanian et al., in preparation). Therefore, we are motivated in finding intelligent approximations that deliver a good tradeoff between alignment quality and CPU time.

To decrease the computational complexity of this problem, we restrict ourselves to a reduced subset F' ⊂ F and we will first search for a consistent subset A ⊂ F' with maximum total score. As in previous versions of DIALIGN, we use pairwise optimal alignments as a filter. In other words, the set F' is defined as the set of all fragments contained in any of the optimal pairwise alignments of the sequences in our input data set. Here, we also restrict the length of fragments using some suitable constant.

For multiple alignment, previous versions of DIALIGN used the above outlined greedy approach. We call this approach a direct greedy approach, as opposed to the progressive greedy approach that we introduce in this paper. A modification of this 'direct greedy' approach was also used in our reimplementation DIALIGN-T. Here, we considered not only the weight scores of individual fragments (or their overlap weights [1]) but also took into account the overall degree of similarity between the two sequences involved in the fragment. The rationale behind this approach is that a fragment from a sequence pair with high overall similarity is less likely to be a random artefact than a fragment from an otherwise non-related sequence pair, see [16] for details.

2.3 The overall algorithm

In the previous section, we discussed all ingredients that are necessary to give a high-level description of our algorithm to compute a multiple sequence alignment. For clarity, we omit algorithmical details and data structures such as the consistency frontiers or consistency boundaries, respectively, that are used to check for consistency as these features have been described elsewhere [23]. We use a subroutine PAIRWISE_ALIGNMENT (s _i , s _j , A) that takes two sequences s _i and s _j and (optionally) an existing consistent set of fragments A as input and calculates an optimal alignment of s _i and s _j under the side constraint that this alignment is consistent with A and that only those positions in the sequences are aligned that are not yet aligned by a fragment from A. Note that in DIALIGN, an alignment is defined as an equivalence relation on the set of all sequence positions, so a consistent set of fragments corresponds to an alignment. Therefore, we do not formally distinguish between alignments and sets of fragments.

Next, a subroutine GREEDY_ALIGNMENT (A, F') takes an alignment A and a set of fragments F' as arguments and returns a new alignment A' ⊃ A by adding fragments from the set F' in a 'directly' greedy fashion. For details on these subroutines see also [16]. Furthermore we use a subroutine BUILD_UPGMA (F') that takes a set F' of fragments as arguments and returns a tree and a subroutine MERGE(p, F') that takes the parent node p and the set of fragments F' as argument and returns an alignment of the set of sequences represented by node p. Those two subroutines are described in the previous two subsections. A pseudo-code description of the complete algorithm for multiple alignment is given in Figure 1. As in the original version of DIALIGN [1], the process of pairwise alignment and consistency filtering is carried out iteratively. Once a valid alignment A has been constructed by removing inconsistent fragments from the set F' of the fragments that are part of the respective optimal pairwise alignments, this procedure is repeated until no new fragments can be found. In the second and subsequent iteration steps, only those parts of the sequences are considered that are not yet aligned and optimal pairwise alignments are calculated under the consistency constraints imposed by the existing alignment A.

Beside the above described improvements of the optimization algorithm, we incorporated new features into DIALIGN-TX that were already part of the original implementation of DIALIGN. DIALIGN-TX now supports anchor points the same way DIALIGN 2.2 does [5, 15]. Anchor points can be used for various purposes, e.g. to speed up alignment of large genomic sequences [25, 6], or to incorporate information about locally conserved motifs. This can been done, for example, using the N-local-decoding approach [26, 27] or other methods for motif finding.

DIALIGN-TX also now comes with an option to specify a threshold parameter T in order to exclude low-scoring fragments from the alignment. Following an approach proposed in [28], the alignment procedure can be iterated, starting with a high value of T and with lower values in subsequent iteration steps. By default, in the first iteration step of our algorithm, we use a value of T = -log₂(0.5) for the pairwise alignment phase, while in all subsequent iteration steps, a value of T = 0 is usedd. With a user-specified threshold of T = 2 for the first iteration step, the threshold value remains -log₂ (0.5) in all subsequent steps, and with a chosen threshold value of T = 1, the value for the subsequent iteration steps is set to -log₂ (0.75).

An optimal pairwise alignment in the sense of our fragment-based approach is a chain of fragments with maximum total weight score. Calculating such an optimal alignment takes O(l³) time if l is the (maximum) length of the two sequences since all possible fragments are to be considered. If the length of fragments is bounded by a constant L, the complexity is reduced to O(l² × L). In practice, however, it is not meaningful to consider all possible fragments. Our algorithm processes fragments starting at a pair of positions i and j, respectively, with increasing fragment length. To reduce the number of fragments considered, our algorithm stops processing longer fragments starting at i and j if the previously visited short fragments starting at the same positions have low scores. More precisely, we consider the average substitution score of aligned amino acids or the average number of matches for DNA or RNA alignment, respectively, to decide if further fragments starting at i and j are considered.

To reduce the run time for pairwise alignments, we implemented an option called fast mode. This option uses a lower threshold value for the average subsitution scores or number of matches. By default, during the pairwise alignment phase, fragments under consideration are extended until their average substitution score is at least 4 for amino acids (note that our BLOSUM62 matrix has 0 for the lowest score possible) and 0.25 for nucleotides, respectively. With the fast mode option, this threshold is increased by 0.25 which has the effect that the extension of fragments during the pairwise alignment phase is interrupted far more often than by default. This option, however, reduces the sensitivity of the program. We observed speed-ups up to factor 10 on various benchmark data when using this option while the alignment quality was still reasonably high, in the sense that the average sum-of-pair score and average column score on our benchmarks deteroriated around 5% – 10% only. We recommend to use this option for large input data containing sequences that are not too distantly related. Hence, this option is not advisable for strictly locally related sequences where we observed a reduction of the alignment quality almost down to a score of zero. However in the latter case this option is not necessary since the original similarity score thresholds of 4 and 0.25, respectively, are effective enough to prevent DIALIGN-TX of unnecessarily looking at too many spurious fragments.

In order to evaluate the improvements of the new heuristics we had several benchmarks on various reference sets and compared DIALIGN-TX with its predecessor DIALIGN-T 0.2.2 [16], DIALIGN 2.2 [29], CLUSTAL W2 [30], MUSCLE 3.7 [31], T-COFFEE 5.56 [32] POA V2 [33, 34], PROBCONS 1.12 [35] & PROBCONSRNA 1.10, MAFFT 6.240 L-INSi and E-INSi [21, 22]. We performed benchmarks for DNA as well as for protein alignment. As globally related benchmark sets we used BRAliBase II [36, 37] for RNA and BALIBASE 3 [38] for protein sequences.

The benchmarks on locally related sequence sets were run on IRMBASE 2 for proteins and DIRMBASE 1 for DNA sequences, which have been constructed in a very similar way as IRMBASE 1 [16] by implanting highly conserved motifs generated by ROSE [39] in long random sequences. IRMBASE 2 and DIRM-BASE 1 both consist of four reference sets ref1, ref2, ref3 and ref4 with one, two, three and four (respectively) randomly implanted ROSE motives. The major difference compared to the old IRMBASE 1 lies in the fact that in 1/s cases the occurrence of a motive in a sequence has been omitted randomly, whereby s is the number of sequences in the sequence family. The results on IRMBASE 2 and DIRMBASE 1 now tell us how the alignment programs perform in cases when it is unknown if every motive occurs in every sequence thus providing a more realistic basis for assessing the alignment quality on locally related sequences compared to the situation in the old IRMBASE 1 where every motive always occurred in every sequence.

Each reference set in IRMBASE 2 and DIRMBASE 1 consists of 48 sequence families, 24 of which contain ROSE motifs of length 30 while the remaining families contain motifs of length 60. 16 sequence families in each of the reference sets consist of 4 sequences each, another 16 families consist of 8 sequences while the remaining 16 families consist of 16 sequences. In ref1, random sequences of length 400 are added to the conserved ROSE motif while for ref2 and ref3, random seqences of length 500 are added. In ref4 random sequences of length 600 are added.

For both BAliBASE and IRMBASE, we used two different criteria to evaluate multi-alignment software tools. We used the sum-of-pair score where the percentage of correctly aligned pairs of residues is taken as a quality measure for alignments. In addition, we used the column score where the percentage of correct columns in an alignment is the criterion for alignment quality. Both scoring schemes were restricted to core blocks within the reference sequences where the 'true' alignment is known. For IRMBASE 2 and DIRMBASE 1, the core blocks are defined as the conserved ROSE motifs. To compare the output of different programs to the respective benchmark alignments, we used C. Notredame's program aln_compare [32].

In this paper, we introduced a new optimization algorithm for the segment-based multiple-alignment problem. Since the first release of the program DIALIGN in 1996, a 'direct' greedy approach has been used where local pairwise alignments (fragments) are checked for consistency one-by-one to see if they can be included into a valid multiple alignment. In this approach, the order in which fragments are checked for consistency is basically determined by their individual weight scores. Some modifications have been introduced, such as overlap weights [1] and a more context-sensitive approach that takes into account the overal significance of the pairwise alignment to which a fragment belongs [16]. Nevertheless, a 'direct' greedy approach is always sensitive to spurious pairwise random similarities and may lead to alignments with scores far below the possible optimal score (e.g. [13, 5], Pöhler and Morgenstern, unpublished data).

The optimization method that we introduced herein is inspired by the so-called progressive approach to multiple alignment introduced in the 1980s for the classical multiple-alignment problem [17]. We adapted this alignment strategy to our segment-based approach using an existing graph-theoretical optimization algorithm and combined it with our previous 'direct' greedy approach. As a result, we obtain a new version of our program that achieves significantly better results than the previous versions of the program, DIALIGN 2 and DIALIGN-T.

To test our method, we used standard benchmark databases for multiple alignment of protein and nucleic-acid sequences. Since these databases are heavily biased towards global alignment, we also used a benchmark database with simulated local homologies. The test results on these data confirm some of the known results on the performance of multiple-alignment programs. On the globally related sequence sets from BAliBASE and BRALIBASE, the segment-based approach is outperformed by classical, strictly global alignment methods. However, even on these data, we could achieve a considerable improvement with the new optimization algorithm used in DIALIGN-TX. On the simulated local homologies, our method clearly outperforms other alignment approaches, and again the new algorithm introduced in this paper achieved significantly better results than older versions of DIALIGN. Among the methods for global multiple alignment, the program MAFFT [21, 22] performed remarkably well, not only on globally, but also on locally related sequences.

We will conduct further studies to investigate to what extent the optimization algorithms used in the different versions of DIALIGN can be improved and which alternative algorithms can be applied to the optimization problem given by the segment-based alignment approach.

DIALIGN-TX is online available at Göttingen Bioinformatics Compute Server ( GOBICS ) at [41] The program source code and our benchmark databases IRMBASE 2 and DIRMBASE 1 are downloadable from the same web site.