Efficient Dominant Point Algorithms for the Multiple Longest Common Subsequence (MLCS) Problem

1 Introduction

The multiple longest common subsequence problem (MLCS) is to find the longest subsequence shared between two or more sequences. It is a classical computer science problem with important applications in many fields such as information retrieval and computational biology [Masek and Paterson, 1980; Smith and Waterman, 1981]. For over 30 years, significant efforts have been made to find an efficient algorithm for the MLCS problem. The most significant contribution has been done to study the simplest case of MLCS of two or three sequences [Hirschberg, 1977; Hakata and Imai, 1998]. However, while several attempts towards finding an efficient algorithm for a general case of more than 3 sequences [Hakata and Imai, 1998; Chen et al., 2006], it is yet to be developed. A general case of MLCS is of a tremendous value to computational biology and computational genomics that deal with biological sequences [Korkin et al., 2008]. With the increasing volume of biological data and prevalent usage of computational sequence analysis tools, an efficient MLCS algorithm applicable to many sequences will have a significant impact on computational biology and applications.

In this paper, we present an efficient algorithm for the MLCS problem of three and more sequences. The new method is based on the dominant point approach. Dominant points are minimal points in a multidimensional search space. Knowing those points allows to reduce the search space size by orders of magnitude, hence significantly the computation time. Our algorithm performs a new divideand-conquer technique to construct dominant point sets efficiently. Unlike FAST-LCS [Chen et al., 2006], a MLCS algorithm that works with the whole dominant point set, our method takes advantages of the structure relationships among dominant points and partitions them into independent subsets, where the divide-and-conquer technique is applied. Compared to existing state-of-the-art MLCS algorithms, our dominant-point algorithm is significantly faster on multiple sequences longer than 1000. We have also developed an efficient parallel version of the algorithm. By dividing the problem into smaller sub-problems and solving the sub-problems in parallel, we have achieved a near linear speedup.

The paper is organized as follows. In the next section, we briefly review state-of-the-art methods for MLCS. Then, we present the basics of the dominant point method in Section 3 and the new sequential algorithm in Section 4. The new parallel algorithm is presented in Section 5. In Section 6, we show the experimental results. Finally, in Section 7, we summarize the paper.

2 Related work

Classical methods for the MLCS problem are based on dynamic programming [Sankoff, 1972; Smith and Waterman, 1981]. In its simplest case, given two sequences a1 and a2 of length n1 and n2 respectively, a dynamic programming algorithm iteratively builds a n1 × n2 score matrix L, in which L[i, j], 0 ≤ i ≤ n1, 0 ≤ j ≤ n2 is the length of a LCS between two prefixes a1[1, . . . , i] and a2[1, . . . , j].

In a straightforward implementation of dynamic programming, we calculates all entries in L. The resulting algorithm has time and space complexity of O(nd ) for d sequences of length n. Various approaches have been introduced to reduce the complexity of dynamic programming [Hirschberg, 1977; Masek and Paterson, 1980; Hsu and Du, 1984; Apostolico et al., 1992; Rick, 1994]. Unfortunately, these approaches primarily address the special case of 2 sequences.

In contrast to dynamic programming, the dominant point approach limits its search to exploring a smaller set of dominant points rather than the whole set of positions in L. The initial idea of dominant points as special points in a matrix was introduced by Hirschberg [1977]. The dominantpoint approach has been successfully applied to the twosequences cases [Hirschberg, 1977; Chin and Poon, 1990; Apostolico et al., 1992]. In [Hakata and Imai, 1998], several dominant-point algorithms for more than 2 sequences were proposed. One of the algorithms, called AlgorithmA, which was designed specifically for finding a LCS of 3 strings, is overwhelmingly faster than dynamic programming for 3 sequences. However, Algorithm A minimizes dominant point sets by enumerating points of the same coordinate values in each dimension. As a result, its complexity increases rapidly for more sequences. The other algorithm, Hakata and Imai's C algorithm [1998], works for arbitrary number of strings. It is similar to another MLCS algorithm published recently, FAST-LCS [Chen et al., 2006], in that they both use pairwise comparison algorithm to compute dominant points.

Parallel processing of the MLCS algorithms have been developed to further speed up the computation of LCS. Similar to sequential algorithms, early parallel algorithms were designed mostly for the special case of 2 sequences and cannot be easily generalized for three and more sequences [Babu and Saxena, 1997; Yap et al., 1998; Xu et al., 2005]. Korkin [2001] attempted to tackle the general MLCS problem but failed to achieve a near-linear speedup. In the most recent approaches, FAST-LCS [Chen et al., 2006] and parMLCS [Korkin et al., 2008], a near-linear speedup was reached for a large number of sequences. parMLCS is a parallel version of Hakata and Imai's C algorithm.

3 Fundamentals of Dominant-Point Method

In this section, we introduce the basics and main ideas of the new dominant-point method.

Assume sequences are strings of characters defined over a finite alphabet Σ. Let a and b be two sequences of lengths n and k correspondingly. For sequence a = s1s2 . . . sn, sequence b = si1 si2 . . . sik is called a subsequence of a if 1 ≤ ij ≤ n, for 1 ≤ j ≤ k, and is < it, for 1 ≤ s < t ≤ k. Let S = {a1, a2, . . . , ad} be a set of sequences over alphabet Σ. A multiple longest common subsequence (MLCS) for set S is a sequence b such that (i) b is a subsequence of ai, 1 ≤ i ≤ d, and (ii) b is the longest one satisfying (i).

Let L be the score matrix for a set of d sequences, a1, a2, . . . , ad, as defined in Eq. (1). A point p in matrix L is denoted as p = [p1, p2, . . . , pd] , where each pi is a coordinate of p for the corresponding string ai. The value at position p of the matrix L is denoted as L[p].

A point p = [p1, p2, . . . , pd] in L is called a match if a1[p1] = a2[p2] = . . . = ad[pd]. If a match p corresponds to character s ∈ Σ, i.e., ai[pi] = s, it is denoted as p(s).

For two points p = [p1, p2, . . . , pd] and q = [q1, q2, . . . , qd], we say that p dominates q if pi ≤ qi, for 1 ≤ i ≤ d. If p dominates q, we denote this relation as p ≤ q. Similarly, p strongly dominates q if pi < qi, 1 ≤ i ≤ d. We denote this relation as p < q. p does not dominate q (and denote this as p ̸≤ q), if there is an i, 1 ≤ i ≤ d, such that qi < pi. Note that p ̸≤ q does not necessarily imply q ≤ p. For some points p and q, both p ̸≤ q and q ̸≤ p can be true. A match p(s) is called a k-dominant point, or k-dominant, or dominant at level k, if (i) L[p] = k and (ii) there is no other match q(s) satisfying (i) and q ≤ p, i.e., p is not dominated by another point q with the same value k and for the same character s. The set of all k-dominants is denoted as Dk . The set of all dominant points is denoted as D.

An example of the score matrix L, the set of dominant points and rest of the matches for two sequences are shown in Figure 1.

A match p(s) is called an s-parent for a point q, if q < p and there is no other match r(s) such that q < r < p. The set of all s-parents for q is denoted as Par(q, s). The set of all s-parents for a set of points A is denoted as Par(A, s). The set of all parents, ∪s∈ΣPar(A, s), for A is denoted as P ar(A, Σ). A point p in a set of points A is called a minimal element of A, if q ̸≤ p for all q ∈ A- {p}. If |A| = 1, then its single element is defined to be a minimal element of A. The set of minimal elements is called minima of A.

It has been proven that (k+1)-dominants, Dk+1 , is exactly the minima of the parent set, Par(Dk ,Σ), of k-dominants, Dk [Hakata and Imai, 1998]. It is easy to find the minima of Par(Dk , Σ) directly [Chen et al., 2006; Korkin et al., 2008]. However, Theorem 1 below provides a more efficient way.

Theorem 1 Let Minima() be an algorithm that returns the minima of a set of dominant points. Then, Minima(Par(Dk ,Σ)) = ∪s∈ΣMinima({p(s) | p(s) ∈ Minima(Par(p,Σ)), p ∈ Dk }), for k = 1, . . . ,K, where K is the length of the longest common subsequence.

Proof of the theorem is omitted due to the size limitations of the paper. This theorem leads to a two-step procedure to compute the minima of the parent set Par(Dk , Σ) as follows,

1. Minimize Par(p, Σ) for each point p ∈ Dk ,
2. Minimize each s-parent set of Dk , s ∈ Σ.

Both Par(p, Σ) of p, p ∈ Dk , and s-parent set of Dk , s ∈ Σ, are smaller than Par(Dk ,Σ). By avoiding minimizing entire Par(Dk , Σ), computation time is saved.

We have developed a fast divide-and-conquer algorithm for the procedureMinima() mentioned above. It is based on the following Theorem 2.

Theorem 2 For d ≥ 3, the minima of N dominant points in the d-dimensional space can be computed in O(dNlogd-2N) time by a divide-and-conquer algorithm. The computation time is O(dNlogd-2n) if the sequence length n ≤ N.

Proof of the theorem can be derived from [Kung et al., 1975; Bentley, 1980; Hakata and Imai, 1998]. Here we just give the central idea. Consider a divide-and-conquer algorithm on a set of N dominant points in d-dimensional space. The algorithm first evenly partitions these points into two subsets R and Q with respect to the d-dimensional coordinate so that the d-dimensional coordinates of points in R are greater than those of points in Q. Then, the algorithm recursively minimizes R and Q, respectively. Finally, after getting the minima of both R and Q, the algorithm removes points in R that are dominated by points in Q.

4 A new sequential algorithm for MLCS

In this section, we present a new dominant-point algorithm for MLCS of any number of sequences. For convenience, we assume below that a1, a2, . . . , ad are sequences over alphabet Σ, and that the lengths of all sequences are equal to n. Our algorithm is based upon the following ideas:

1. Each position in a MLCS corresponds to a match in the score matrix L. Therefore, the search can be restricted to the set of all matches in L.
2. The number of dominant levels is equal to the size of MLCS. Moreover, there exists at least one dominant point at each level k, k = 1, 2, . . . , K, which corresponds to the k-th position in theMLCS. Therefore, the search can be further restricted to all dominant points.
3. The set of dominant points in L can be computed recursively, where D(k+1) is computed, based solely on the dominant points of the previous level, Dk . The next two ideas explain the computation in more detail.
4. The set of (k + 1)-dominants, D(k+1) , is the minima of the set of all parents for k-dominants, Dk , i.e. D(k+1) = Minima(P ar(Dk , Σ)).
5. Minima( Par(Dk ,Σ)) can be computed in two steps as indicated in Theorem 1: (i) For each dominant point p ∈ Dk , compute Minima(Par(p,Σ)) and

take a union of all such sets, Pars = {p(s) | p(s) ∈ Minima(Par(p,Σ)), p ∈ Dk }, s ∈ Σ; (ii) Dk+1 is a union of non-overlapping sets, Minima(Pars), i.e., Dk+1 = ∪s∈ΣMinima(Pars), for each s ∈ Σ.

1. Finally, Minima() can be implemented using a fast divide-and-conquer approach as indicated in Theorem 2.

Based on these ideas, the sequential dominant-point and divide-and-conquer (Quick-DP) algorithm is as follows:

Algorithm Quick-DP ({a1, a2, . . . , ad}, Σ ) Calculation of dominant points 01 Preprocessing; D 0 = {[0, 0, . . . , 0]}; k = 0; 02 while D k not empty do { 03 for p ∈ D k do { 04 B = Minima(Par(p, Σ)); 05 for s ∈ Σ do{ 06 Pars = Pars ∪ {p(s)|p(s) ∈ B}; }} 07 D k+1 = ∪s∈Σ Minima (Pars) ; 08 k = k+1; } Calculation of MLCS-optimal path 09 pick a point p = [p1, p2, . . . , pd] ∈ D k-1 ; 10 while k-1 > 0 do { 11 current LCS position = a1[p1]; 12 pick a point q such that p ∈ Par(q, Σ) ; 13 p = q; 14 k = k-1; }

Quick-DP consists of two parts. In the first part, the set of all dominants is calculated iteratively, starting from a 0dominant set (containing one element). The set of (k + 1)- dominants, D(k+1) , is obtained, based on the set of kdominants, Dk . In the second part, a MLCS-optimal path, corresponding to a MLCS, is calculated, tracing back through set of dominant points obtained in the first part of the algorithm, and starting with an element from the last dominant set. All MLCS can be enumerated systematically as well.

To efficiently enumerate all parents of each dominant point, we calculate a preprocessing matrix T = {T [s, j, i]} , s ∈ Σ, 0 ≤ j ≤ max1≤k≤d{|ak|}, 1 ≤ i ≤ d, where each element T[s, j, i] specifies the position of the first occurrence of character s in the i-th sequence, starting from the (j + 1)-st position in that sequence. If s does not occur any more in the i-th sequence, the value of T[s, j, i] is equal to 1 +max1≤k≤d{|ak|}. The calculation of this preprocessing matrix T takes O( n |Σ| d) time.

Let |Σ| be the size of alphabet Σ and |D| the size of dominant point set D. From Theorem 2, we can derive that it takes O(|Σ| d logd-2 |Σ|) time to compute the minima of each parent set Par(p,Σ), p ∈ D, and O(|D| d logd-2 n) time to compute the minima of each s-parent set Pars, s ∈ Σ. Hence the time complexity of Quick-DP is

Based on the experimental evaluation of |D| and the estimated complexities of Hakata and Imai's C algorithm [1998] and FAST-LCS [2006], we expect our approach to be significantly faster than C algorithm and FAST-LCS for large n.

5 A new parallel algorithm forMLCS

As shown in the previous section, calculating the set of (k + 1)-dominants, Dk+1 , requires computing the minima of parent set Par(p, Σ), for p ∈ Dk , and the minima of s-parent Pars, s ∈ Σ, of Dk . These sets can be calculated independently. Based on this observation, we propose the following parallelization of the sequential algorithm.

Given Np+1 processors, the parallel algorithm uses one as the master and Np as slaves and performs the following steps:

1. The master processor computes D0 .
2. Every time the master processor computes a new set Dk of k-dominants ( k = 1, 2, 3, . . . ), it distributes them evenly among all slave processors.
3. Each slave processor computes the set of parents and the corresponding minima of k-dominants that it has and then sends the result back to the master processor.
4. The master processor collects each s-parent set Pars, s ∈ Σ, as the union of the parents from slave processors and distributes the resulting s-parent set among slaves.
5. Each slave processor i is assigned to find the minimal elements only of one s-parent set, Pars.
6. Each slave processor i computes the set Dk+1 i of (k+1)- dominants of Pars and sends it to the master processor.
7. The master processor computes Dk+1 = Dk+1 1 ∪ Dk+1 2 ∪ . . . ∪ Dk+1 Np , and goes to step 2.

The pseudocode of the parallel algorithm Quick-DPPAR is as follows.

Algorithm Quick-DPPAR ({a1, a2, . . . , ad}, Σ, Np ) 01 Proc0 : Preprocessing; D 0 = {[0, 0, . . . , 0]}; k = 0; 02 while D k not empty do { 03 Proc0 : distribute elements of D k Each processor Proci, 1 ≤ i ≤ Np, performs : 04 get D k i from Proc0; 05 for p ∈ D k i do { 06 B = Minima(Parents(p)); 07 for s ∈ Σ do{ 08 Parsi = Parsi ∪ {p(s)|p(s) ∈ B}; }} 09 Send Parsi, s ∈ Σ, to Proc0 ; 10 Proc0 : calculate Pars = ∪1≤i≤Np Parsi, s ∈ Σ; 11 Proc0 : distribute Pars, s ∈ Σ; Each processor, Proci, 1 ≤ i ≤ Np, performs : 12 get Pars, s ∈ Σ ; 13 D k+1 i = Minima(Pars); 14 send D k+1 i to Proc0; 15 Proc0 : D k+1 = ∪1≤i≤Np D k+1 i ; 16 k = k+1; }

Quick-DPPAR assigns each s-parent set Pars, s ∈ Σ, of Dk to a slave processor to compute the minima of Pars using our divide-and-conquer method. So as many as |Σ| slave processors can work simultaneously in this step of minimization. To utilize more than |Σ| processors, we further parallelize the divide-and-conquer algorithm. Our parallel algorithm generates a binary tree during the execution of Minima(Pars).

The root node of the tree contains the entire data set Pars. The data set of each internal node of the tree is equally split into two subsets which correspond to the children of the node. The subsets are distributed among processors and the tree is built in parallel. Finally, the minima of each subset is solved recursively and the results are combined in the parent node.

6 Experimental Results

In our experiments, the algorithms were run on a Linux Server with 8 Intel Xeon(R) CPUs (2.826GHz) and 16GB memory. The programming environment is GNU C++. The algorithms were tested on random strings of length 100 to 4000 over alphabets of size 4 (e.g., nucleotide sequences) and 20 (e.g., protein sequences).

6.1 Sequential algorithmQuick-DP

The sequential method Quick-DP is compared with Hakata and Imai's A and C algorithms [Hakata and Imai, 1998]. The A algorithm is specifically designed for 3 strings and is the most efficient one so far for 3 strings. The C algorithm can work with any number of strings. We implemented Hakata and Imai's A algorithm and C algorithm according to their paper. The data structures of them were carefully designed to reduce their computation time.

For each string length of an alphabet, we generated 10 sets of 3 random strings. Quick-DP and Hakata and Imai's A and C algorithms were tested on the same data sets and their average running times for each length of sequence are shown in Fig. 2. Fig. 2 shows that Quick-DP is comparable and slightly faster than Hakata and Imai's A algorithm on 3 strings. We note that A was designed specially for the case of 3 strings, while Quick-DP works for an arbitrary number of sequences. Fig. 2 also shows that Hakata and Imai's C algorithm is significantly slower.

Next, we compared Quick-DP with another state-of-theart sequential algorithm FAST-LCS [Chen et al., 2006] and Hakata and Imai's C algorithm on test sets consisting of more than 3 sequences. The results in Fig. 3 shows that Quick-DP is several orders of magnitude faster than FAST-LCS and C algorithm. For instance, for randomDNA sequences, QuickDP is over 1,000 times faster than FAST-LCS on sequences of length 140 and over 1,000 times faster than C algorithm on problems of length 200.

6.2 Parallel algorithmQuick-DPPAR

The parallel algorithm Quick-DPPAR was implemented using multithreading. Following the parallelization scheme in Section 5, our implementation consists of 1 master thread and multiple slave threads. The master thread allocates dominant points to slaves to perform time-consuming computation.

We first evaluated the speedup of our parallel algorithm Quick-DPPAR over sequential algorithm Quick-DP. Speedup is defined as the ratio of sequential time and parallel time. Similar to previous test data, we generated 10 sets of 5 random strings each for alphabet 4 and 20, respectively. We ran Quick-DPPAR using different number of slave threads. Fig. 4 shows the results for sequence lengths 1000, 1500, and 2000. Near-linear speedups were achieved for all these cases.

Then, we measured the running time of Quick-DPPAR on strings of various lengths. We ran Quick-DPPAR using 8 slave threads on test sets of 5 random DNA sequences. The results in Fig. 5 demonstrate the efficiency of Quick-DPPAR. We have also applied Quick-DPPAR to determine a MLCS for 11 proteins from the family of melanin-concentrating hormone receptors (MCHR). The lengths of the protein sequences range from 296 to 423 amino acids. As a result, it took Quick-DP 939 seconds to detect a MLCS. For QuickDPPAR, the time was reduced to 185 seconds.

7 Conclusions

In this paper, we have presented fast sequential and parallel algorithms for MLCS problems of arbitrary number of strings. The algorithms improve existing dominant point methods through new efficient methods for computing dominant points. For problems of 3 strings, the sequential algorithm is as fast as the best existing MLCS algorithm designed for the 3-string case, Hakata and Imai's A algorithm. For problems of more than 3 strings, the sequential algorithm is much faster than the best existing algorithms, achieving more than a thousand times faster performance on larger size problems. The parallel implementation is efficient and achieves a near-linear speedup over the sequential algorithm.

Acknowledgement

This research was supported in part by the Shumaker Endowment in Bioinformatics and NIH Grant R33GM078601.