MUSCA: An Algorithm for Constrained Alignment ofMultiple Data Sequences

Laxmi Parida

parida@us.ibm.com

Aris Floratos

floratos@us.ibm.com

Isidore Rigoutsos

rigoutso@us.ibm.com

Bioinformatics and Pattern Discovery Group
Computational Biology Center
IBM Thomas J. Watson Research Center
Yorktown Heights, NY10598, USA

Abstract

Given a set of N sequences, the Multiple Sequence Alignment problem is to align these N
sequences, possibly with gaps, that brings out the best commonality of the N sequences. MUSCA1
is a two-stage approach to the alignment problem by identifying two relatively simpler sub-problems
whose solutions are used to obtain the alignment of the sequences. We rst discover motifs in the
N sequences and then extract an appropriate subset of compatible motifs to obtain a \good"alignment. The motifs of interest to us are the irredundant motifs which are only polynomial in
the input size. In practice, however, the number is much smaller (sub-linear). Notice that thisstep aids in a direct N -wise alignment, as opposed to composing the alignments from lower order(say pairwise) alignments and the solution is also independent of the order of the input sequences;hence the algorithm works very well while dealing with a large number of sequences. The second
part of the problem that deals with obtaining a good alignment is solved using a graph-theoretic
approach that computes an induced subgraph satisfying certain simple constraints. We reduce aversion of this problem to that of solving an instance of a set covering problem, thus o er the best
possible approximate solution to the problem (provided P6=NP). Our experimental results, while
being preliminary, indicate that this approach is ecient, particularly on large numbers of longsequences, and, gives good alignments when tested on biological data such as DNA and protein
sequences. We introduce the the notion of an alignment number K (2  K  N ), a user-controlledparameter, that lends a useful exibility to the aligning program: this additional requirement
constrains the alignment to have at least K sequences agree on a character, whenever possible, in
the alignment. The usefulness of the alignment number is corroborated by the users who view thisas a natural constraint while dealing with a large number of sequences.

1

Introduction

Given a set of sequences, the Multiple Sequence Alignment problem is to align these sequences,possibly with gaps, that brings out the best commonality of the sequences. Various alignment costfunctions [2, 3, 4, 6, 8, 7, 14, 15, 12, 9], have been used in literature. The general approach to solving
the pairwise ( = 2) sequence alignment problem has been a dynamic programming technique using
di erent mechanisms of scores which is a function of the edit distance, along with gap penalties, toevaluate the similarity of the sequences. In [16, 13] the case of
2 has been handled by rst doinga pairwise alignment for some or all possible pairs in some order and then building a -wise alignment
from these.
MUSCA2 uses a two-stage approach to the alignment problem by identifying two relatively simpler
sub-problems which deal separately with the two issues, one of identifying the \local similarities" and
N

N

N

N

N >

N

1

Musca is a constellation in the polar region of the Southern Hemisphere near Apus and Carina. Also, MUSCA is ananagram of the salient characters in Constrained Multiple Sequence Alignment.

2

Musca is a constellation in the polar region of the Southern Hemisphere near Apus and Carina. Also, MUSCA is ananagram of the salient characters in Constrained Multiple Sequence Alignment.

112

the other of aligning the similarities appropriately. We rst discover motifs in the sequences, andthen use these motifs to obtain a \good" alignment. Informally, a motif is a repeated pattern thatappears more than once in a sequence. In the alignment context a motif is a pattern that appears in
two or more input sequences A major point of criticism regarding using motifs is that they are usually
very large in number (exponential in input size); however, we show that the motifs that are relevantto the alignment problem are the irredundant motifs, and the number of these motifs is polynomialin the input size [10]. Moreover, in practice, this number is much smaller (sub-linear). Thus, using
motifs for the alignment helps in at least two ways: (1) it aids in a direct -wise alignment, as
opposed to composing the alignments from lower order (say pairwise) alignments and (2) the solution
is independent of the order of the input sequences. We believe that, in practice, these have importantconsequences.
The second sub-problem of the alignment problem is that of obtaining a good alignment. Noticethat any arbitrary set of motifs need not necessarily give rise to an alignment, under the premise thatthe alignment that uses a motif does not introduce gaps in the motif. Having obtained all possible
motifs in the rst stage, this stage involves pruning this set to obtain a (sub)set that gives an alignment.
We solve this problem by mapping the motifs of the rst stage to a suitable directed graph. Next weshow that obtaining an alignment of the motifs is equivalent to solving a set-covering problem. Thus
we present a very systematic way of aligning sequences based on motifs.
It is well known that the multiple sequence alignment problem, in addition to being a hard-tosolve problem, is also very hard to model to the satisfaction of evolutionary biologists, geneticists
and other users. Does our approach have any theoretical contributions to the multiple sequencealignment problem in general? We introduce the the notion of an alignment number (2   ),a user-controlled parameter, that lends a useful exibility to the aligning program: this additional
requirement constrains the alignment to have at least sequences agree on a character, wheneverpossible, in the alignment. This is particularly of interest when a large number of sequences are beingaligned. The utility of the alignment number is corroborated by the users who view this as a natural
constraint while dealing with a large number of sequences.
N

N

K

K

N

K

Roadmap. We describe our two-stage approach of motif discovery in Section 2 and aligning sequences in Section 3. We discuss the issues involved in using motifs for alignment and present a
simple graph theoretic formulation in Section 3.1. We present heuristic algorithm for this problem bymapping it to a set covering problem in Section 3.2.
2

Motif Discovery (Stage 1)

We begin by giving a rigorous de nition of a motif.
De nition 1 (K -motif m, location list Lm ) Given a string s on alphabet 6 and an integer K , 2 
K < jsj, a string m on 6 [ `:` is a K -motif with location list Lm = (l1 ; l2 ; : : : ; lp ), if all of the following
hold:
1. m[0]; m[jmj 0 1] 2 6,

[First and last characters of the motif are solid characters; if \dont care" (the `.' character)characters are allowed at the ends, the motifs can be made arbitrarily long in size without
conveying any extra information.]

2. p  K ,3. there does not exist a location l, l 6= li , 1  i  p such that m occurs at l on s (the location list
is of maximal size), and,[This ensures that any two distinct location lists must have distinct motifs associated with each.]113

4. for every \dont care" character at position j in m, there exist at least two distinct occurrences
li1 and li2 , 1  i1 ; i2  p, such that s[li1 + j ] 6= s[li2 + j ].

[No \dont care" character can be replaced by a solid character, otherwise an arbitrary numberof \dont care" characters could be introduced without conveying any extra information.]
If is a string of 6[ f`.'g, is called a rigid motif. In the rest of the discussion a -motif is referred
to as a motif.
We give the the de nitions of maximality and irredundancy below.
mmK

De nition 2 (Maximal Motif) Let p1 , p2 , : : :, pk be all the motifs in the sequence s. A motif pi ismaximal if and only if there is no pj (j 6= i) such that pi is a substring of pj , or, if pi is a substringof pj , then there exists at least one occurrence of pi in s that is not covered by pj in s.
De nition 3 (Redundant motif) A maximal motif m, with location list Lm , is redundant if there exist
maximal motifs mi , 1  i  p, such that Lm = Lm1 [ Lm2 : : : [ Lmp .
A motif that is not redundant is called an irredundant motif.

The motifs of interest to the sequence alignment problem are the irredundant motifs (see lemma 4)
which are only quadratic in the input size and there exists a polynomial time algorithm to extract
them from the input [10].
Similar de nitions extend to discovering common motifs from multiple sequences. Teiresias isan ecient implementation of the motif discovery problem [11] and we use this in our experiments.
Notice that a strip is a restricted version of the motif, in the sense that it is a motif of size one. In
the following discussion we discuss motifs while bearing in mind that the results would hold also forstrips (the results that are not critically dependent on the motif size being greater than one).
3

Sequence Alignment

We skip the process of discovering motifs from input sequences. Similar de nitions extend to discovering common motifs from multiple sequences. TEIRESIAS is an ecient implementation of the motif
discovery problem [11] and we use this in our experiments.
We obtain the irredundant motifs and the position of the motif in each sequence it appears in.
The o set list, associated with each motif i , is a list of two-tuples ( i i ), where i is a pointer tothe sequence and i is the o set in the sequence. We assume that every o set list has exactly one
occurrence of each sequence. At the end of the motif discovery step we have motifs, 1 2N;is the number of the motifs (not necessarily distinct), each with an o set list consisting of at least 2 (distinct) sequences and the position in the sequence where the motif appears.
Can all the motifs be used in an alignment? If the answer is yes, we form an alignment that respects
all the motifs. If the answer is no, we remove the \o ending" motifs in a manner that optimizes acost function and gives an alignment. To investigate this further, we now explore the conditions underwhich two motifs can be used simultaneously in an alignment.
p

s ;l

s

l

p ;p ;:::;pNk

De nition 4 (Motif Overlap) Two irredundant motifs pi and pj overlap if there exists a sequence s
containing both these motifs and the following holds. Let ni and nj be the sizes of the motifs pi andpj and let li and lj be the locations (o sets) in a sequence s respectively, then the motifs overlap if theintervals [li ; li + ni ] and [lj ; lj + nj ] have a non-empty intersection.
De nition 5 (Pairwise Compatible Motifs) Two motifs, p1 and p2 , are pairwise compatible if there
exists an alignment of the sequences that does not introduce gaps in the motifs p1 and p2 .
Lemma 1 Two irredundant motifs pi and pj are pairwise compatible if and only if none of the following holds:

114

1. If pi and pj do not overlap in all the sequences, then pi is to the left of pj , without loss ofgenerality (otherwise a domain crossing mismatch is said to have occurred).
2. If pi and pj overlap in any sequence, then pi is at some xed distance d to the left of pj , without
loss of generality (otherwise an overlap mismatch is said to have occurred).

We de ne the alignment using a set of motifs as follows.
De nition 6 (sequence alignment, compatible set) Given a set S of motifs, vi1 ; vi2 ; : : : ; vin , a motifalignment of the sequences, s1 ; s2 ; : : : ; sm , is the alignment such that in all the sequences, with no gaps
in the motifs, the motifs vi1 ; vi2 ; : : : ; vin , are aligned (in all the sequences they appear). If such analignment exists, the set s1 ; s2 ; : : : ; sm , is called a compatible set.
De nition 7 (linear ordering of motifs) Given a set of compatible motifs, a consistent ordering of
the motifs such that, in every sequence, the set of motifs that are present in the sequence appear in
the left to right order and this ordering is called the linear ordering.

Is it sucient to just check for pairwise incompatibility of motifs while seeking an alignment?
De nition 8 (domain crossing error) Given a set of motifs, m1 ; m2 ; : : : ; mn , a domain crossing error
is said to occur if there exists a linear ordering of the motifs mi1 ; mi2 ; : : : ; min , yet there exists noalignment that respects all the n motifs.
Lemma 2 A set of irredundant motifs p1 ; p2 ; : : : ; pn is feasible if and only if none of the following
holds:
1. There exist distinct motifs pi and pj such that pi and pj are pairwise incompatible.
2. There exists a non-empty subset of the motifs without a linear ordering.
3. There exists a non-empty subset of the motifs that demonstrate domain crossing error.

3.1 The Graph-theoretic Formulation

Next we wish to capture these conditions in a graph as follows. Construct a directed graph = ( )where every motif i corresponds to a vertex i , thus = j j. The directed edges are introduced asfollows:
1. There is no edge between two vertices where the two corresponding motifs do not occur simultaneously in any sequence.
2. If i is to the left of j in every sequence that the two motifs are present, then a directed edge
is placed from i to j . This is to indicate that in the alignment the motif i appears to the leftof j . The edges are labeled as follows:(a) Label forbidden, if the motifs corresponding to 1 and 2 are not pairwise compatible.
(b) Label overlap, if the motifs corresponding to 1 and 2 overlap.
(c) Label nonoverlap, if the motifs corresponding to 1 and 2 are pairwise compatible but do
not overlap.
The linear ordering of motifs is captured by checking for cycles in the graph. However the domaincrossing mismatch requires a more careful handling as described below.
G

p

v

p

N

V

p

v

v

p

p

v

v

v

v

v

115

v

V; E

Handling domain crossing mismatches. We associate a distance Dv1 v2 with every edge that is
not labeled forbidden. This is used to compute the feasibility of a collection of motifs corresponding toa solution; this does not contribute to the cost of the alignment. (We discuss the weight corresponding
to the cost in the next section.) Let pi and pj be the two motifs corresponding to the vertices. Thenif d is the minimum of the distance between the occurrences of the two motifs in every sequence thatboth of them appear in, Dvivj = d.
To detect the domain crossing mismatches of motifs (that are pairwise compatible), we de ne thenotion of a consistent graph w.r.t. a vertex.
De nition 9 Let G = (V; E ) be a labeled, weighted, directed, graph with weights on the edge uv givenby D(u; v) and a label 2 ff orbidden; overlap; nonoverlapg. A path, P , is valid if it has no edges
labeled forbidden. Further, a valid path, P , is called an overlap-path if all the edges in the path arelabeled overlap. The weight of the valid path P , D , is the sum of the weights of its constituent edges.
Let p 2 V . The graph is consistent w.r.t p if 8 q 2 V , for all pairs of vertex-disjoint valid paths
from p to q, P 1 and P 2 ,P

1. D 1 = D 2 , if P 1 and P 2 are both overlap-paths, or,PP

2. D 1  D 2 , if P 1 is an overlap-path and P 2 is not.
PP

We now present the straightforward observation that relates the set of compatible motifs to afeasible subgraph.
Lemma 3 The following two statements are equivalent: Given a subset of motifs 1 2

n from the set of all motifs from m sequences of input, thesubset is compatible, if the following holds:
p ;p ;:::;p1. the motifs are not pairwise incompatible,
2. there exists a linear ordering of p1 ; p2 ; : : : ; pn , and,3. there is no domain crossing mismatch in p1 ; p2 ; : : : ; pn .

 Given a subset of vertices 1 2v ; v ; : : : ; vn , constructed as de ned in this section. The induced
subgraph on v1 ; v2 ; : : : ; vn is feasible, if the following holds:

1. there is no edge labeled forbidden in the induced subgraph,
2. the induced subgraph is acyclic, and,3. the induced subgraph is consistent w.r.t. every vertex vi , 1  i  n.

Lemma 4 If p is a redundant motif, then using the motif p does not improve the cost of the optimization problem.
Proof Sketch. Let p be rendered redundant by motifs p1 ; p2 ; : : : ; pn , n  1. By de nition, motif p
has less number of solid-characters than each of pi , 1  i  n. Thus if an alignment can use motifp, it can certainly use all the motifs p1 ; p2 ; : : : ; pn , giving a larger number of solid-characters; thus a2higher cost for the optimization problem.

3.2 Algorithm to compute the \best" alignment

Given a set of incompatible motifs, the set can be grouped into sets (not necessarily disjoint) suchthat each set violates exactly one of the three conditions of Lemma (2). These sets are called basic
incompatible sets. Next, it can be easily shown that we can remove exactly one motif from a basic
incompatible set to make it compatible.
The algorithm proceeds in the following three steps.
116

1. Detect the basic infeasible (sub)sets.
2. Eliminate motifs to obtain a feasible set that maximizes the cost.
3. Render the alignment.
Step 1. In this step we form subset of vertices that lead to incompatibility of the motifs. Usinglemma (3), we compute the following sets:
1. Construct the sets 1 , 2, , nf where each set consists of two vertices which are the end
points of an edge labeled Forbidden. This is done by simply scanning all the edges and collecting
the end points of the edges labeled Forbidden.
2. Construct the sets 1 , 2 , , nc where each set consists of vertices that form a directed cyclein the graph. All the cycles are captured by carrying out a depth rst search (DFS) of the graph.
3. Construct the sets 1 , 2 , , np where each set consists of vertices that form a closed pathin the graph. These are captured by carrying out a breadth rst search (BFS) rooted at eachvertex.
It is easy to see that the basic incompatible sets are 1 2
nf 1 2
nc 1 2np .
Step 2. Set-covering problem. An instance (Y ) of the set-covering problem consists of a nite set
and a family Y of subsets of , such that every element of belongs to at least one subset in Y := [S . We say that a subset 2 Y covers its elements. The problem is to nd a minimum-size
subset A  Y whose members cover all of : = [S . Any A satisfying this condition covers .
See [5] for details on this problem.
We construct an instance of a set cover problem (the dual of our problem) as follows. Let
F

F

:::

F

C

C

:::

C

P

P

:::

P

F ;F

:::;F

;C ;C

:::;C

;P ;P

:::;P

X;

X

X

X

2Y S

X

S

X

f1 2

X

2A S

X

n g = F1 [ F2 : : : [ Fnf [ C1 [ C2 : : : [ Cnc [ P1 [ P2 : : : [ Pnpv ;v ;:::;v

The elements of the the set ( of the set cover problem) areX

F1 ; F2 : : : ; F

nf ; C1 ; C2 : : : ; Cnc ; P1 ; P2 : : : ; Pnp :The subset i corresponds to each i , 1   (where Y = f 1
S

S

v

i

S ; S2 : : : ; Snn g), and is de ned as

i = fFl jvi 2 Fl ; 1  l  nf g [ fCl jvi 2 Cl ; 1  l  nc g [ fPl jvi 2 Pl ; 1  l  np g:

Thus i denotes all the basic incompatible set that has a common element i , and removing the motif
corresponding to i suces to make all the corresponding basic sets compatible. Now, it is easy to
see that a solution to the set-covering problem gives the minimum number of motifs that need to be
removed so that the remaining set of motifs is compatible.
Associating weights to i which re ect the weight of each motif (depending on the cost function),gives a weighted set cover problem. The set-covering problem is known to be MAX SNP hard and thegreedy algorithm is the best known approximation algorithm for the problem, assuming 6=[1].
To re ect the underlying cost function a weight is associated with every motif in the following manner.
Let be the number of solid characters in the corresponding motif, i , and let the number of sequences
containing i be , then for the -MSA problem the associated weight is and for -MSAmax it is .
We ignore the change in cost due to the common cost of a set of motifs 3.
Step 3. This step consists of the following substeps:
S

v

v

S

P

c

p

p

3

NP

l

k

cl

k

c

For a more accurate computation of the cost an appropriate common cost due to a set of overlapping motifs must

be used.

117

3.1 We compute a linear ordering (see De nition 7) using the graph G of Step 2, of the motifs ij1 ,
i , , i for each sequence . Such a linear ordering of the motifs exists, since the set of motifs
j2
ji
is feasible.
3.2 From the original sequence i , we obtain the llers (if any) between two consecutive motifs asi i i, jii . 0i is the leftmost portion, possibly empty. For example, let sequence
0 , j1 , j2 ,and let two motifs be as follows: i1 =and i2 = . Then 0i = ,i=
i1 = , 2i = .
3.3 We obtain an alignment of the sequences by appropriately aligning each ( ijl + jil ) and 0i , =12
i , lling the gaps with `-'. The symbol `+' denotes a string concatenation operation.
The alignment is such that each motif of a sequence is perfectly aligned with the corresponding
motif in all the other sequences.
For example let 1 =and 2 =. Then 11 = 21 =and 12 =
2 =
1
1
1
2
2
2
and 0 = , 1 = , 2 = , 0 =empty, 1 =, 2 = . Then the alignment
2of the sequences are as follows (the motifs are shown in bold):
0

p

:::

p

f

f

f

p

i

s

:::

f

s

abcdef ghij kl

f

f

f

f

p

cde

p

ghi

p

;

ghi

f

(1)
(2)
4

ab

f

f

l

;:::;j

s

p

f

jkl

abcdef ghijkl

ab

f

a

b

f

f

c d e
c d e

s

jkl

f
x

cdexyzpqrghitu

f

f

p

xyzpqrf

0 0 0 0 0 g h i
y

z

p

q

r

g h i

p

cde

p

tu

j

k

t

u

l

Summary

We have proposed a two-stage approach to the alignment problem by handling two relatively simpler
sub-problems which deal separately with the two issues, one of identifying the \local similarities" andthe other of aligning the similarities appropriately. In the rst stage we identify all possible -wise
motifs, i.e., all motifs that appear simultaneously in at least of the input sequences (2   ).
In the second stage, we give plausible alignments of a carefully chosen subset of these motifs (thatoptimize certain cost functions). Using this approach for the alignment helps in at least two ways:(1) it aids in a direct -wise alignment, as opposed to composing the alignments from lower order(say pairwise) alignments and (2) the resulting alignment is independent of the order of the input
sequences. is an input parameter and is called the alignment number. In practice, our approachworks particularly well for alignment of a large set of (long) sequences. We have presented the resultof running our alignment algorithm on biological data and the results look very promising. In the fullversion of the paper we discuss the computational complexity of the underlying optimization problem
along with an analysis of the approximation factor of any alignment that results from the algorithm.
K

K

N

K

N

N

K

References[1] Arora, S., Karger, D., Karpinski, M., Polynomial time approximation schemes for dense instancesof NP-hard problems, Proceedings of the Twenty-Seventh Annual ACM Symposium on the Theoryof Computing, ACM, 284{93, 1995.
[2] Altschul, S.F., Gap costs for multiple sequence alignment, J. Theor. Biol., 138:297{309, 1989.
[3] Chao, K.M., Hardison, R., and Millter, W., Recent developments in linear-space alignment methods: a survey, J. Comput. Biol., 3:271{291, 1994.
[4] Carrillo, H. and Lipman, D., The multiple sequence alignment problem in biology, SIAM J. Appl.
Math., 481073{1082, 1988.
[5] Cormen, T.H., Leiserson, C.E., and Rivest, R.L., Introduction to Algorithms, The MIT Press,
Cambridge, Massachusetts, 1990.
118

[6] Gupta, S.K., Kececioglu, J.D., and Scha er, A.A., Improving the practical space and time eciency of the shortest-paths approach to sum-of-pairs multiple sequence alignment, J. Comput.
Biol., 2(3):459{472, 1995.
[7] Higgins, D.G. and Sharp, P.M., CLUSTAL: a package for performing multiple sequence alignment
on a microcomputer, Gene, 73:237{244, 1988.
[8] Hirosawa, M., Totoki, Y., Hoshida, M., and Ishikawa, M., Comprehensive study on iterativealgorithms of multiple sequence alignment, Comput. Appl. Biosci., 11(1):13{18, 1995.
[9] Lawrence, C.E., Altschul, S.F., Boguski, M.S., Liu, J.S., Neuwald, A.F., and Wootton, J.C.,
Detecting subtle sequence signals: a Gibbs sampling strategy for multiple alignment, Science,
262(5131):208{214, 1993.
[10] Parida, L., Algorithmic Techniques in Computational Genomics, Ph.D. thesis, Courant Instituteof Mathematical Sciences, New York University, September, 1998.
[11] Rigoutsos, I. and Floratos, A., Motif discovery in biological sequences without alignment orenumeration, In Proceedings of the Annual Conference on Computational Molecular Biology (RECOMB98), 221{227, ACM Press, 1998.
[12] Rocke, E. and Tompa, M., An algorithm for nding novel gapped motifs in DNA sequences,
In Proceedings of the Annual Conference on Computational Molecular Biology (RECOMB98),
228{233, ACM Press, 1998.
[13] Vihinen, M., An algorithm for simultaneous comparison of several sequences, Comput. Appl.
Biosci., 4(1):89{92, 1988.
[14] Waterman, M.S., Parametric and ensemble alignment algorithms, Bull. Math. Biol., 56(4):743{
767, 1994.
[15] Waterman, M.S., An Introduction to Computational Biology: Maps, Sequences and Genomes,
Chapman Hall, 1995.
[16] Miller, W. and Zhang, Z., and He, B., Local multiple alignment vis subgraph enumeration,
Discrete Appl. Math., 71:337{365, 1996.

119