Methods

LAGAN and Multi-LAGAN: Efficient Tools
for Large-Scale Multiple Alignmentof Genomic DNA
Michael Brudno,1 Chuong B. Do,1 Gregory M. Cooper,2 Michael F. Kim,1Eugene Davydov,1 NISC Comparative Sequencing Program,1 Eric D. Green,3Arend Sidow,2 and Serafim Batzoglou1,41Department of Computer Science, Stanford University, Stanford, California 94305-9010, USA; 2Department of Pathology
and Department of Genetics, Stanford University, Stanford, California 94305-5324, USA; 3Genome Technology Branch
and NIH Intramural Sequencing Center, National Human Genome Research Institute, National Institutes of Health,Bethesda, Maryland 20892, USA
To compare entire genomes from different species, biologists increasingly need alignment methods that are
efficient enough to handle long sequences, and accurate enough to correctly align the conserved biologicalfeatures between distant species. We present LAGAN, a system for rapid global alignment of two homologous
genomic sequences, and Multi-LAGAN, a system for multiple global alignment of genomic sequences. We tested
our systems on a data set consisting of greater than 12 Mb of high-quality sequence from 12 vertebrate species.
All the sequence was derived from the genomic region orthologous to an ∼1.5-Mb region on humanchromosome 7q31.3. We found that both LAGAN and Multi-LAGAN compare favorably with other leadingalignment methods in correctly aligning protein-coding exons, especially between distant homologs such as
human and chicken, or human and fugu. Multi-LAGAN produced the most accurate alignments, while requiring
just 75 minutes on a personal computer to obtain the multiple alignment of all 12 sequences. Multi-LAGAN is apractical method for generating multiple alignments of long genomic sequences at any evolutionary distance.
Our systems are publicly available at http://lagan.stanford.edu.
Comparing genomic sequences across related species is a fruitful source of biological insight, because functional elements
such as exons tend to exhibit significant sequence similarity,
whereas regions that are not functional tend to be less conserved. The first step in comparing genomic sequences is toalign them—that is, to map the letters of one sequence tothose of the others. There are several categories of alignments:
local alignments that identify local similarities between regionsof each sequence, and global alignments that find a monotonically increasing map between the letters of each sequence;
pairwise alignments that compare two sequences, and multiplealignments that compare several sequences.
Local pairwise alignment methods such as SmithWaterman (1981), BLAST (Altschul et al. 1990, 1997), BLASTZ(Schwartz et al. 2000), SSAHA (Ning et al. 2001), and BLAT
(Kent 2002) are able to pinpoint locations of rearrangements
between two sequences, and are suitable for aligning draft
sequences or individual reads. Global alignments are important because they reveal the shared order of biological featuresin the compared species, and produce a more accurate alignment at the base-pair level when the features are in the sameorder. The best-known global alignment algorithm is Needleman-Wunsch (1970), which requires time proportional to the
product of the lengths of the aligned sequences. Unfortunately this algorithm is too inefficient for comparing longgenomic sequences. Faster methods have been developed re4
Corresponding author.
E-MAIL serafim@cs.stanford.edu; FAX (650) 725-1449.
Article and publication are at http://www.genome.org/cgi/doi/10.1101/
gr.926603. Article published online before print in March 2003.

cently: DIALIGN (Morgenstern et al. 1998, Brudno and Morgenstern 2002), MUMmer (Delcher et al. 1999, 2002), GLASS(Batzoglou et al. 2000), WABA (Kent and Zahler 2000), and
AVID (Bray et al. 2003). Most of these methods have proveneffective in aligning genomic sequences from two closely related organisms, such as human and mouse or Caenorhabditis
elegans and C. briggsae, but have not been tested in alignments
between distant relatives such as human and fugu.
Multiple alignments, a natural extension of twosequence comparisons, are a powerful way to study biologicalsequences. Even weak similarity across several sequences usually reveals an important conserved biological feature (Dubchak et al. 2000; Göttgens et al. 2002). Moreover, multiplealignments enable the computation of local rates of evolution, giving a quantitative measure of the strength of evolutionary constraints and the functional importance of local
regions (Simon et al. 2002). Multiple alignments are considerably more difficult to compute than are pairwise alignments: the running time scales as the product of the lengthsof all the sequences. Formally, the problem is NP-complete(Wang and Jiang 1994; Bonizzoni and Vedova 2001). For thisreason heuristic approaches are usually applied, of which themost widely used is progressive alignment, which constructs a
multiple alignment by successive applications of a pairwisealignment algorithm. The best-known system based on progressive alignment is perhaps CLUSTALW (Thompson et al.
1994). Some other systems include MULTALIGN (Barton andSternberg 1987), MULTAL (Taylor 1988), YAMA (Hardison etal. 1993, 1994), and PRRP (Gotoh 1996). DIALIGN (Morgenstern 1999) does not use progressive alignment; instead it uses

13:721–731 ©2003 by Cold Spring Harbor Laboratory Press ISSN 1088-9051/03 $5.00; www.genome.org

Genome Researchwww.genome.org

721

Brudno et al.

another heuristic approach to chain local conserved blocks
between several sequences into a multiple alignment. Thesesystems can effectively align proteins and relatively short genomic regions, but are not efficient enough to align entire
genomes. MGA (Höhl et al. 2002) is a rapid multiple alignersuitable for comparing very close homologs, such as different
strains of a bacterium, but is not designed to align distanthomologs.
Here we describe novel systems for pairwise and multiplealignment of genomic sequences: LAGAN (Limited Area Global Alignment of Nucleotides), an efficient and reliable pairwise aligner that is suitable for genomic comparison of distantly related organisms, and MLAGAN (Multi-LAGAN), a
multiple aligner based on progressive alignment with LAGAN.
We tested our systems on sequence from 12 species generated
for the genomic segment harboring the cystic fibrosis transmembrane conductance regulator (CFTR) gene (J.W. Thomas,
J.W. Touchman, R.W. Blakesley, G.G. Bouffard, S.M. Beckstrom-Sternberg, E.H. Margulies, M. Blanchette, A.C. Siepel,
P.J. Thomas, J.C. McDowell, B. Maskeri, N.F. Hansen, M.S.
Schwartz, R.J. Weber, W.J. Kent, D. Karolchik, T.C. Bruen, R.
Bevan, D.J. Cutler, S. Schwartz, L. Elnitski, J.R. Idol, A.B. Prasad, S.-Q. Lee-Lin, V.V.B. Maduro, M.E. Portnoy, N.L. Dietrich, N. Akhter, K. Ayele, B. Benjamin, K. Cariaga, C.P.
Brinkley, S.Y. Brooks, S. Granite, X. Guan, J. Gupta, P.
Haghighi, S-L. Ho, M.C. Huang, E. Karlins, P.L. Laric, R. Legaspi, M.J. Lim, Q.L. Maduro, C.A. Masiello, S.D. Mastrian,
J.C. McCloskey, R. Pearson, S. Stantripop, E.E. Tiongson, J.T.
Tran, C. Tsurgeon, J.L. Vogt, M.A. Walker, K.D. Wetherby, L.S.
Wiggins, A.C. Young, L-H. Zhang, K. Osoegawa, B. Zhu,
B. Zhao, C.L. Shu, P.J. De Jong, C.E. Lawrence, A.F. Smit,
A. Chakravarti, D. Haussler, P. Green, W. Miller, and E.D.
Green, in prep.). Based on comparisons with other availablealignment programs and benchmarking on standard desktop
computer systems, we conclude that LAGAN and MLAGAN
are practical and reliable methods for large-scale pairwise andmultiple genomic alignment that should prove useful for obtaining alignments of the entire human, mouse, fugu, rat, andother genomes in the context of a whole-genome alignment
pipeline.

proportional to the product of the sequence lengths. Efficient
global aligners such as MUMmer (Delcher et al. 1999, 2002),
GLASS (Batzoglou et al. 2000), and AVID (Bray et al. 2003) relyon anchoring. First, they detect local similarities between thetwo sequences. Second, they select and fix an ordered set of
local similarities, the anchors. Finally, they align the interleaving regions.
Anchoring reduces computation time because the initiallarge alignment problem is subdivided into many smaller,manageable ones. To perform anchoring reliably, the aligner
needs a sensitive local-alignment detection method, an accurate procedure for selecting the anchors, and an efficientmethod for computing the global alignment based on theanchors. Anchoring is easier when the two sequences arehighly similar and contain frequent matching words. Existing
global aligners based on anchoring are primarily designed forcomparing highly similar sequences, such as human andmouse. MUMmer, GLASS, and AVID, for instance, use exact
matching words for detecting local alignments.
LAGAN performs anchoring with techniques designed towork well on distant, as well as close organisms. Specifically,
LAGAN uses the CHAOS algorithm (Brudno and Morgenstern2002), a highly sensitive method that detects local alignments
using multiple short inexact words instead of longer exact
words. LAGAN then constructs the global alignment by first
applying CHAOS recursively in areas with sparse anchors, sothat each consecutive pair of anchors is separated by a distance smaller than a given maximum, and then performing alimited-area dynamic programming algorithm on the areaaround the anchors.
LAGAN aligns two sequences in three steps outlined
next, and described in Methods.

Algorithm LAGAN
Given a pair of sequences X, Y: A global alignment between thetwo sequences is a path from the top-left to the bottom-rightcorner of the alignment matrix X ⳯ Y (Fig. 1A)

RESULTS
Outline of Algorithms
LAGAN is a global alignment system that aligns two genomic
sequences in three main steps: (1) generation of local alignments between the two sequences, (2) construction of a roughglobal map, by chaining an ordered subset of the local alignments, and (3) computation of the final global alignment, byfinding the best alignment that stays within a limited areaaround the rough global map. Based on LAGAN, we developed MLAGAN, a new multiple alignment system that aligns
genomic sequences in two main phases: (1) a progressive alignment phase that constructs a multiple alignment by successively aligning two sequences, or intermediate multiple alignments, with the LAGAN algorithm, and (2) an optional iterative improvement phase that successively removes eachsequence from the multiple alignment, and realigns it to therest of the alignment, until no significant improvements areobserved.

Global Sequence Alignment, LAGAN
In aligning two long genomic sequences, the NeedlemanWunsch algorithm is too inefficient because it requires time

722

Genome Researchwww.genome.org

Figure 1 The LAGAN algorithm. (A) A global alignment between
two sequences is a path between the top-left and the bottom-rightcorner of their alignment matrix. (B) LAGAN first finds all local alignments between the two sequences. (C) LAGAN computes a maximalscoring ordered subset of the alignments, the anchors, and puts together a rough global map. (D) LAGAN limits the search for an optimal alignment to the area included in the boxes and around theanchors, and computes the optimal Needleman-Wunsch alignment
limited to that area. LAGAN uses memory proportional to the area ofthe largest box plus the memory to hold the optimal alignment.

LAGAN: Multiple Alignment of Genomic DNA

Step 1. Generation of local alignments. Compute the localalignments between the two sequences (Fig. 1B). Assign a weight to each local alignment.
Step 2. Construction of rough global map. Two local alignments can be chained if the end of one precedes thestart of the other in both X and Y. Compute the highest-weight chain of local alignments, using the Longest Increasing Subsequence algorithm (Fig. 1C). Calleach local alignment in the chain an anchor. Apply
Steps 1, 2 recursively between every pair of anchorsthat are more than a threshold distance apart.
Step 3. Computation of global alignment. Define a subset ofthe cells in the alignment matrix:
1. Include all cells at distance at most r from cells of theanchor.
2. Between each consecutive pair of anchors, where thefirst ends at positions (i, j) and the second begins at
positions (k, l) in X, Y, respectively, include the fullalignment box (Xi…Xk) ⳯ (Yj…Yl) between the two
endpoints (Fig. 1D).
Compute the global alignment using NeedlemanWunsch-like dynamic programming within the subsetof the cells defined above.

Multiple Alignment, MLAGANMLAGAN is based on progressive alignment: A multiple alignment of K sequences is constructed in K–1 pairwise alignment
steps, where in each step two sequences, or intermediate multiple alignments, are aligned. MLAGAN uses LAGAN as the
pairwise-alignment subroutine, and introduces new methods
for scoring a multiple alignment with affine gaps, and forrefining a multiple alignment.
The multiple aligner CLUSTALW infers a phylogenetic
tree of the sequences, using hierarchical clustering. Inferring
such a tree is useful when aligning protein families with several members within one genome. On the other hand, whenaligning orthologous genomic regions from different species(e.g., human, baboon, mouse, rat, and chicken)—which is the
problem that we are targeting—the phylogenetic tree is usually known. Accordingly, MLAGAN assumes that the phylogenetic tree is given.
MLAGAN aligns K sequences in three main steps that are
outlined next. Step 1, the generation of rough global maps,
uses Steps 1 and 2 of LAGAN. Steps 2 and 3 are described in
Methods.

Algorithm MLAGAN
Given K sequences X1,…, XK, and a phylogenetic binary tree
between them.
Step 1. Generation of rough global maps. Find the rough global
map between each pair of sequences.
Step 2. Progressive multiple alignment with anchors.
2.1 Perform a global alignment between the two closest
sequences according to the phylogenetic tree, using
LAGAN.
2.2 Each alignment between two or more sequences produced in step 2 is a new multi-sequence. Find therough global maps of this multi-sequence to all othermulti-sequences. The details of this procedure and of
the scoring metric are described in Methods.
2.3 Iterate Steps 2.1 and 2.2, performing at every step a
global alignment between the two closest (multi)-

sequences, according to the phylogenetic tree. Repeat
until left with a multiple alignment of all sequences.
Step 3 (Optional). Iterative refinement with anchors. For eachsequence Xi in the multiple alignment:
3.1 Find segments of Xi that align better than a given
cutoff, in the existing multiple alignment. These segments are the anchors between Xi and the other sequences.
3.2 Align Xi to the multiple alignment of the other sequences with LAGAN.

Evaluation of Performance
We tested the performance of LAGAN and MLAGAN againstthat of leading global and local alignment programs, by measuring the ability of each aligner to correctly align proteincoding exons in orthologous sequences. We used the following two datasets: (1) the ROSETTA set (Batzoglou et al. 2000),
which contains 129 orthologous annotated genes with complete intron sequences between human and mouse of average
length 10 Kbp, and (2) the CFTR region (J.W. Thomas, J.W.
Touchman, R.W. Blakesley, G.G. Bouffard, S.M. BeckstromSternberg, E.H. Margulies, M. Blanchette, A.C. Siepel, P.J. Thomas, J.C. McDowell, B. Maskeri, N.F. Hansen, M.S. Schwartz,
R.J. Weber, W.J. Kent, D. Karolchik, T.C. Bruen, R. Bevan,
D.J. Cutler, S. Schwartz, L. Elnitski, J.R. Idol, A.B. Prasad, S.-Q.
Lee-Lin, V.V.B. Maduro, M.E. Portnoy, N.L. Dietrich, N.
Akhter, K. Ayele, B. Benjamin, K. Cariaga, C.P. Brinkley, S.Y.
Brooks, S. Granite, X. Guan, J. Gupta, P. Haghighi, S-L. Ho,
M.C. Huang, E. Karlins, P.L. Laric, R. Legaspi, M.J. Lim, Q.L.
Maduro, C.A. Masiello, S.D. Mastrian, J.C. McCloskey, R. Pearson, S. Stantripop, E.E. Tiongson, J.T. Tran, C. Tsurgeon, J.L.
Vogt, M.A. Walker, K.D. Wetherby, L.S. Wiggins, A.C. Young,
L-H. Zhang, K. Osoegawa, B. Zhu, B. Zhao, C.L. Shu, P.J.
De Jong, C.E. Lawrence, A.F. Smit, A. Chakravarti, D. Haussler,
P. Green, W. Miller, and E.D. Green, in prep.), which for thestudies described here consisted of 12 orthologous sequencesfrom human, chimpanzee, baboon, cat, dog, cow, pig, mouse,rat, chicken, fugu, and zebrafish; these ranged in length from160 Kbp to 1.8 Mbp, with an average of 1 Mbp. For the CFTR
region sequences, we used the human gene annotations to
identify the orthologous exons in the other 11 species’ sequences. This was done on the unaligned sequences before
any of the tests were conducted. We first aligned each known
human exon to the other sequences using TBLASTN and filtered out all hits that were <50% identical at the amino-acidlevel or covered <50% of the human exon. Then we designated the best match for each species as the orthologous
exon. Using this method, we annotated 75 or more exons inall species, except chicken (24) and zebrafish (34), whose sequences were grossly incomplete. We recorded the positionsof the exons for use in the tests. After producing the alignments of the entire region, we tested whether the exons werecorrectly aligned: For each exon, we calculated the fraction ofits length aligned to the corresponding exon in the orthologous sequences.
Tables 1 and 2 summarize the results for the ROSETTA
and CFTR datasets, respectively. In the ROSETTA set, we
tested the performance of MUMmer, GLASS, DIALIGN, AVID,BLASTZ, and LAGAN. In the CFTR set, we did not test GLASSor DIALIGN because these methods require too much computation time to align longer sequences. We tested the performance of MLAGAN by projecting the multiple alignment
to the 11 pairwise alignments between the human and the

Genome Researchwww.genome.org

723

Brudno et al.

alignment-finding steps (the calls to CHAOS during LAGAN
Step 1) are performed on the translated amino-acid level,
rather than on the nucleotide level (as explained in Methods).
We tested the performance of translated anchoring on the
CFTR dataset. Table 3 summarizes the performance of LAGAN
and MLAGAN, when applied under our default translated anchoring parameters (see Methods for details on the parameters). LAGAN’s performance changes in the more distantalignments: it worsens in human/chicken, and improves significantly in human/fugu, and human/zebrafish. MLAGAN’s
performance is only affected in the human/fugu alignment,
where it improves. Overall the performance improves whentranslated anchoring is used, and this option may be significant when aligning very distant species such as human/fugu.
However, translated anchoring in principle biases the aligners
to perform best in detecting homologies in protein-coding
regions, and consequently could perform worse in noncoding regions. For that reason, we disabled translated anchoring in the default parameters of the programs.
In anchoring-based global alignment, accuracy of the final alignment depends on the accuracy of the anchors. We
found that the majority of human exons are covered correctly
during the anchoring phase (calls to CHAOS, Steps 1 and 2 of
LAGAN; Table 4).
Pairwise alignments produced by LAGAN can be visual-

Table 1. Performance of Aligners on the ROSETTA Datasetof 1160 Total Exons in Human and Mouse
Aligner

100%
exons

90%
exons

70%
exons

Time
(sec)

DIALIGN
MUMmer
GLASS
AVID
BlastZ
LAGAN

89
0
91
90
94
94

96
1
97
95
97
97

98
3
98
97
98
98

388
17
154
19
17
48

Columns show the percentage of exons annotated in human thatare aligned to the orthologous mouse exon over at least 70%,
90%, and 100% of their length, and the time required to align the
129 sequences.

other species’ sequences. Of the pairwise aligners tested,
LAGAN, MLAGAN, and BLASTZ were the most accurate ataligning the close (mammalian) sequences. LAGAN and MLAGAN were the best in aligning distant homologs. MLAGAN
was slower than the other aligners, but was also the mostaccurate.
LAGAN and MLAGAN can be applied with a translated
anchoring option, according to which one or more of the local-

Table 2. Performance of Aligners on the CFTR Region

Number of exons
MUMmer100%
90%–100%
70%–100%
AVID
100%
90%–100%
70%–100%
BlastZ
100%
90%–100%
70%–100%
LAGAN
100%
90%–100%
70%–100%
MLAGAN100%
90%–100%
70%–100%
Time (sec)MUMmerAVID
BlastZLAGANMLAGANMax memory (MB)MUMmerAVID
BlastZLAGANMLAGAN

BaboonChimpanzee

Mouse

Rat

Cow

Pig

Cat

Dog

Chicken

Zebrafish

Fugu

Overall

232

176

230

230

224

174

176

182

68

48

150

1914

100
100
100

99
100
100

6
8
14

7
9
16

28
40
52

32
44
55

38
47
56

28
37
45

0
0
0

0
0
0

0
0
0

36
41
47

100
100
100

100
100
100

94
98
100

95
100
100

98
99
100

97
100
100

99
100
100

93
98
99

66
79
85

33
42
44

19
21
29

88
91
92

100
100
100

100
100
100

97
100
100

97
100
100

96
98
100

97
100
100

100
100
100

94
99
99

96
97
97

73
79
79

66
72
80

94
97
98

100
100
100

100
100
100

97
100
100

97
100
100

98
99
100

97
100
100

100
100
100

94
99
99

96
99
100

83
88
92

72
77
81

95
98
98

100
100
100

100
100
100

97
100
100

97
100
100

98
99
100

97
100
100

100
100
100

94
99
99

99
100
100

88
98
100

73
84
90

96
98
99

8
82
31
56
—

6
57
24
50
—

7
215
46
78
—

7
221
43
82
—

8
165
40
60
—

6
111
36
68
—

6
139
34
62
—

6
131
33
78
—

3
83
7
338
—

2
53
5
158
—

2
76
6
133
—

61
1775
305
1135
4547

40
598
239
90
—

39
551
276
90
—

40
581
202
90
—

40
584
212
89
—

40
578
204
90
—

39
498
200
87
—

39
522
208
88
—

39
502
206
87
—

39
387
188
88
—

38
340
185
87
—

38
360
185
89
—

40
598
276
90
670

The annotated human exons were aligned with TBLASTN to create pseudo-annotations of exons in the other organisms. The table reports thepercentage of exons covered by alignments over at least 70%, 90%, and 100% of their length, and the time and memory required to obtain
the alignments. The last column reports the percentages of exons aligned out of total number of exons, the total time required for the 11alignments (for the single 12-sequence multiple alignment in the case of MLAGAN), and the maximum memory required.

724

Genome Researchwww.genome.org

LAGAN: Multiple Alignment of Genomic DNA

evidence for an alternative annotation. For example, consider
exon 10 of the cMet gene, as annotated by RefSeq as well as
Ensembl and Twinscan gene predictions. According to ouralignment, the RefSeq-annotated start position is not wellconserved among the mammals analyzed (Fig. 3). There is alsoa large gap (∼250 positions) in the alignment after the first 20
nucleotides of this exon, generated by sequence in cat thathas no clear homolog in any of the other species aligned.
However, a consensus splice acceptor site and the start position of an Acembly gene prediction exon (exon 9 in Met.f andMet.e, exon 10 in Met.a and Met.h), beginning 54 nucleotides
downstream of the RefSeq exon 10, is neatly conserved acrossall the mammals studied. Additionally, no human ESTs correspond to the RefSeq annotation, whereas all ESTs (as annotated in the UCSC genome browser) align precisely with the
Acembly gene-prediction exon. The comparative information
provided by MLAGAN thus complements the experimental
evidence from the human ESTs in identifying this as a potential misannotation.
Finally, we investigated whether the availability of genomic sequence of intermediate distances can help align distant sequences correctly. We applied MLAGAN to the human,fugu, and mouse sequences of the CFTR region, and compared
the projected human/fugu alignment with the correspondingalignment produced by LAGAN. The MLAGAN alignment resulted in a 3% increase in the number of exons perfectlyaligned (from 72%–75%), and a 6% increase in the number ofexons aligned over at least 70% of their length (from 81%–87%), compared to the LAGAN human/fugu alignment. The
well-aligning regions of human and mouse effectively attract
the conserved regions in fugu, making it more likely thatthose will be correctly aligned to human.
In summary, our results demonstrate that LAGAN and
MLAGAN are capable of efficiently solving difficult multiple
sequence alignment problems.

Table 3. Changes in Performance of LAGAN and MLAGAN
on the CFTR Dataset When Translated Anchoring Is Enabled
Alignment

Aligner

100%

90%+

70%+

Time
diff.

Human/chicken

LAGAN
MLAGANLAGANMLAGANLAGANMLAGANLAGANMLAGAN

ⳮ2%
0%
+6%
+4%
+7%
0%
+1%
0%

ⳮ3%
0%
+8%
+3%
+6%
0%
0%
+1%

ⳮ3%
0%
+8%
+1%
+4%
0%
+1%
0%

+25%
—
+93%
—
+29%
—
+60%
+66%

Human/fugu
Human/zebrafish
Total for all 11alignments

We report the differences in the percentage of exons correctlyaligned with translated anchoring vs. with the default options
where translated anchoring is disabled. No changes were observed in the alignments between mammalian sequences.

ized using the VISTA (Mayor et al. 2000) program. Multiplealignments produced by MLAGAN can also be visualized withVISTA, by projecting them into the pairwise alignments withrespect to one reference sequence. For example, Figure 2A
shows the VISTA plot of a multiple alignment between human, chimpanzee, cow, mouse, chicken, and fugu, projected
into the human/chimpanzee, human/cow, human/mouse,
human/chicken, and human/fugu alignments.
We compared the alignment generated by LAGAN and
CLUSTALW for the first intron of the cMet gene in eight mammalian sequences (human, baboon, cat, dog, cow, pig, mouse,
and rat). Although the alignments between all of the species
besides rodents were similar, CLUSTALW misaligned themouse sequence around 4 kb and 10 kb, as is shown in Figure
2B. CLUSTALW required ∼40 minutes to compute the eightspecies alignment after it generated the guide tree, whereas
MLAGAN required ∼2 minutes. We did not compare the performance of MLAGAN and CLUSTALW in longer sequences,because the running time of CLUSTALW scales quadraticallywith sequence length.
MLAGAN’s accuracy in aligning known exons permits
identification of potential misannotations, and can provide

DISCUSSION
Sequence alignment is one of the oldest and most successful
applications of Computer Science to Biology. Despite the considerable advances achieved after several decades of researchin this area, many important challenges remain (Miller 2001).

Table 4. Percentage of Exons Correctly Aligned During the CHAOS-Based Anchoring Phase of LAGAN(Steps 1,2)
90%–100%

50%–89%10%–49%0%–9%Alignment
human vs.

defaulttransl.

defaulttransl.

defaulttransl.

defaulttransl.

Chimpanzee
Baboon
Cow
Pig
Dog
Cat
Mouse
Rat
Chicken
Fugu
Zebrafish

97
77
21
19
20
20
10
13
9
12
27

97
77
21
19
20
20
10
10
9
29
33

3
21
43
33
37
41
35
47
26
19
29

3
21
43
30
37
41
37
46
41
29
42

0
3
36
47
41
39
52
40
62
57
38

0
3
36
50
42
39
51
42
47
35
25

0
0
0
1
2
0
2
1
3
12
4

0
0
0
1
1
0
2
2
3
7
0

Percentage of human exons that have a given proportion of their lengths covered (90%–100%, 50%–89%, 10%–49%, and 0%–9%). Translated anchoring performs better at covering a greater proportion of exons. However, we do
not have measurements on whether translated anchoring performs worse in covering biologically significant noncoding features.

Genome Researchwww.genome.org

725

Brudno et al.

Figure 2 Visualization of a multiple alignment using VISTA. (A) MLAGAN alignments can be visualized using VISTA, if they are projected topairwise alignments with respect to one reference sequence. This plot shows the conservation between human and chimpanzee, cow, mouse, andfugu around the first intron of the cMet gene. The human/chimpanzee conservation is uniformly very high; human/cow and human/mouse show
varying levels of conservation. The human/chicken alignment also shows some conservation in the non-coding areas. The human/fugu alignmentshows conservation only within the first coding exon, and to a lesser degree within the regions upstream and downstream of that exon. (B) First
introns of cMet, comparison of CLUSTALW and MLAGAN alignments. We compared the alignment generated by LAGAN and CLUSTALW for the
first intron of the cMet gene in eight mammalian sequences (human, baboon, cat, dog, cow, pig, mouse, and rat). The alignments between allof the species except rodents were similar. VISTA plots of the projections to human and mouse are shown. CLUSTALW (top) misaligned the mousesequence around 4 Kb and 10 Kb, whereas MLAGAN (bottom) found significant conservation in these regions.

One of these challenges is the development of systems foraligning multiple long genomic sequences efficiently and reliably. This problem is becoming increasingly important as
several entire genomes from related species become available.
MLAGAN is suitable for high-throughput reliable multiple alignment of genomic sequences. MLAGAN is based on
LAGAN, which is designed for rapid and reliable alignment ofa pair of genomic sequences that may be very similar (e.g.,human and chimpanzee), or very distant (e.g., human andfugu). Both programs are available for use on a public web
server, http://lagan.stanford.edu, and are compatible with the
VISTA visualization tool. The source code for LAGAN and
MLAGAN is available from the authors.
LAGAN and MLAGAN assume that one has already identified apparent orthologous regions between two species, and726

Genome Researchwww.genome.org

that there are no genomic rearrangements. Identifying orthologous regions is challenging in the context of an automatic pipeline for whole-genome alignment. While genome
rearrangements have been researched at the level of genes(Sankoff 1999), aligning sequences in the presence of rearrangements is a direction for future research.
An interesting question in research on sequence comparisons has been whether local or global alignment algorithms are more appropriate for a given application. Whengenomic sequences come in several unordered fragments, orwhen genomic rearrangements are frequent, local alignments
are more appropriate, and less likely to miss true conservation. At the same time, genomes such as human, mouse, rat,and even fugu, share large blocks of essentially uninterrupted
synteny. As the number of sequences increases, pairwise local

LAGAN: Multiple Alignment of Genomic DNAalignments may be increasingly difficult to reconcile into an overall
picture of conservation. We believethat multiple global alignments oforthologous blocks provide a muchcleaner picture of the genome evolution, and will become increasingly important as the number ofavailable genomes increases.
Our results suggest that multiple alignments are better thanpairwise alignments at aligning
conserved exons between distantspecies. For that reason, a three-way
global alignment of the human,
mouse, and fugu genomes may be agood way to obtain accurate alignments of the human and fugu
genes.
Some of the most promising
methods for annotating biologicalfeatures such as genes and regula- Figure 3 Multiple alignment of a misannotated exon of the cMet gene. According to the MLAGAN
alignment, the RefSeq-annotated start position of the exon is not well conserved among the mammals.
tory elements are based on two- There is also a large gap (∼250 positions) in the alignment after the first 20 nucleotides of this exon,
species sequence comparisons.
generated by sequence in cat that has no clear homolog in any of the other species aligned. However,Such annotation methods use se- a consensus splice acceptor site and the start position of an Acembly gene prediction exon (exon 9 inMet.f and Met.e, exon 10 in Met.a and Met.h), beginning 54 nucleotides downstream of the RefSeq
quence conservation, plus the coexon 10, is neatly conserved across all the mammals studied. Additionally, no human ESTs correspondincidence of signature of the bioto the RefSeq annotation, whereas all ESTs (as annotated in the UCSC genome browser) align precisely
logical feature in the two species, aswith the Acembly gene prediction exon.
signals to detect the feature. Characteristic signatures such as open
reading frames and splice-site consensus sequences in genes,
k-long words (k-mers) that match with at most c differences
between the two sequences. Given a maximum distance d,nested nucleotide complementarities in non-coding RNA
and maximum shift s, two seeds that are x- and y-letters apartgenes, and motifs in regulatory elements, should be muchin the first and second sequences, respectively, can be chainedeasier to recognize whenever they are conserved across multogether if x ⱕ d, y ⱕ d, and | x ⳮ y | ⱕ s. A seed is chained totiple species. Multiple-comparison-based annotation methods
the single previous seed that creates the highest scoring chain
may prove powerful in annotating biological features withinamong all chains that end with this seed. CHAOS also supa multiple alignment.
ports a translation option, in which both nucleic sequences
Multiple sequence alignment will become an increasare translated in all six coding frames (three forward and three
ingly important tool for biological discovery as several gereverse), and all combinations of frames are compared in turn.
nomes suitable for cross-species comparison become availAmino acids are grouped (Stanfel 1996), and all amino acidsin the same group are considered equal.
able. Our systems should enable researchers to align long orthologous regions from several species, and design multiple
Scoring of Chainsalignment pipelines on a whole-genome scale.
After computing the maximal chains, CHAOS scores eachchain by using match and mismatch penalties for the letters
METHODSof each seed, and gap penalties proportional to | x ⳮ y | foreach pair of chained seeds. CHAOS throws away chains that
Global Pairwise Alignment, LAGAN
score below some threshold t. We augment this scoring
method, by adding a rapid rescoring step: Chains that score
LAGAN aligns a pair of genomic sequences in three main
below t are immediately thrown away. Chains that score
steps: (1) generation of local alignments, (2) construction of aabove t are rescored by performing ungapped extensions inrough global map, and (3) computation of the final global
both directions from each seed, and finding the optimal loalignment. These three steps are described in detail next.
cation to insert exactly one gap of size | x ⳮ y | . If the alignment is done on amino acid sequences, the rescoring is done
Generation of Local Alignments
using a BLOSUM (Henikoff and Henikoff 1992) matrix.
To generate local alignments between the two sequences,
LAGAN uses CHAOS, a method that finds local homologies
between two sequences, and chains them into a rough global
map (Brudno and Morgenstern 2002). Instead of CHAOS, any
efficient local alignment method can be used (details on how
to incorporate a different local aligner into LAGAN can befound in our website, http://lagan.stanford.edu). Here, wesummarize CHAOS and describe rapid rescoring, an improvement we have made to this algorithm.
The CHAOS algorithm works by chaining short words,
the seeds, which match between the two sequences. Given aword length k, and a degeneracy c, a (k, c)-seed is a pair ofConstruction of a Rough Global Map
LAGAN orders the local alignments produced by CHAOS intoa rough global map (Fig. 1B,C). A local alignment is a vector(b, e, b⬘, e⬘, s), representing begin and end positions of thealignment in each sequence, and the score of the alignment.
Given two local alignments A1 = (b1, e1, b1⬘, e1⬘, s1), A2 = (b2,e2, b2⬘, e2⬘, s2), we say A1 < A2 if and only if e1 < b2 and e1⬘ < b2⬘.
A chain of local alignments A1 < A2 < … < Ak, has score s1 + s2
+ … + sk. The optimal rough global map is the highest-scoring
chain, which can be computed using Sparse Dynamic Pro-

Genome Researchwww.genome.org

727

Brudno et al.

gramming in time O(n log n) where n is the total number of
local alignments (Eppstein et al. 1992).
Recursive Anchoring
When computing the local alignments with CHAOS, thechoice of parameters k (length of seeds), c (maximum degeneracy of seeds), and t (score threshold) presents a tradeoff
between speed (more restrictive parameters: higher k, lower c),
and sensitivity to detecting all significant alignments (morepermissive parameters: lower k, higher c, lower t). LAGAN uses
a recursive method similar to the one used in GLASS (Batzoglou et al. 2000), to achieve a combination of speed and sensitivity. During a recursive anchoring step, LAGAN calls
CHAOS with a restrictive set of parameters, computes a roughglobal map based on the resulting local alignments after thelatter are rescored, and then recursively calls CHAOS with amore permissive set of parameters in the regions between
each anchor of the global map.
Translated Anchoring
As an option, some of the recursive anchoring steps can betranslated. This option biases the aligner to perform best inprotein-coding regions and potentially worse in other conserved regions, and therefore as a default translated anchoringis disabled.

Computation of Global Alignment
To compute the final global alignment, LAGAN uses therough global map to limit the search area of dynamic programming. In standard Needleman-Wunsch, the search area
consists of a (M+1) ⳯ (N+1) matrix, where M and N are thelengths of sequences X and Y, respectively. The algorithm
produces the optimal global alignment, which is the highestscoring path from cell (0, 0) to cell (M, N); each cell in the pathis greater than the previous cell, by 1 in the x coordinate, ycoordinate, or both.
The Needleman-Wunsch algorithm evaluates the optimal intermediate alignment score for all (M+1) ⳯ (N+1) cells
in the alignment matrix, and therefore requires time O(MN).
If, however, it is known that the alignment passes through a
point (i, j), then there is no need to evaluate cells in the rectangles from (i+1, 0) to (M, jⳮ1), and from (0, j+1) to (iⳮ1, N).
LAGAN uses this idea: for every anchor in the rough global
map, starting at (i, j) and ending at (i⬘, j⬘), LAGAN limits thecomputation of Needleman-Wunsch within the following
three areas (where r is a parameter, typically 15): (1) the rectangle (0, 0) to (i + r, j + r), (2) the rectangle (i⬘ ⳮ r, j⬘ ⳮ r) to(M, N), (3) the band enclosed by the two diagonals from(i ⳮ r, j + r) to (i⬘ ⳮ r, j⬘ +r), and from (i + r, j ⳮ r) to (i⬘ + r,j⬘ ⳮ r; Fig. 4). In this sense the anchors in LAGAN are more
flexible than the anchors in MUMmer, AVID, and GLASS; inthese programs, the anchors are fixed, whereas in LAGAN
they provide only approximate locations by which the alignment should pass.
Memory-efficient Implementation
The search area of the limited-area dynamic programming
step consists of necks that enclose the anchors, and rectangles
that connect two consecutive anchors (Fig. 1C). In practice,
LAGAN performs the entire computation with memory proportional to the size of the largest rectangle. LAGAN achieves
this memory efficiency as follows: First, it allocates working
memory for one rectangle and the neck that follows it, andcomputes the Needleman-Wunsch matrix within this rectangle and neck. Second, LAGAN traces back all optimal alignments ending in the cells at the rightmost column of theneck; in practice, these alignments quickly converge upon asingle optimal alignment. Finally LAGAN deallocates all
working memory, except the memory necessary to keep the728

Genome Researchwww.genome.org

Figure 4 Limited area of dynamic programming around an anchor.
An anchor between (i, j) and (i⬘, j⬘) in the rough global alignment
limits the search area of Needleman-Wunsch. The alignment is required to pass through the diagonal from (i ⳮ r, j + r) to (i + r, j ⳮ r),stay within the diagonals from (i ⳮ r, j + r) to (i⬘ ⳮ r, j⬘ + r) and from(i + r, j ⳮ r) to (i⬘ + r, j⬘ ⳮ r), and exit through the diagonal from(i⬘ ⳮ r, j⬘ + r) to (i⬘ + r, j⬘ ⳮ r).

traced-back alignments. These three steps are then repeated
for the next rectangle and neck.

Running Time Analysis
The running time of LAGAN is dominated by the total number of cells in the rectangles between consecutive anchors. The
number of cells in the necks of constant width r is O[r · (M +N)] and is linear in the sequence lengths.
Suppose the anchoring step yielded n anchors, and let(x0, y0),…,(xn, yn) be the dimensions of each of the n + 1
rectangles. Letnn兺x and L = 兺 y

L1 =

2

i

i= 0

i

i= 0

denote the total length of the inter-anchor segments in eachsequence. At this stage we are effectively working with sequences of length L1 and L2, as we can assume the anchors
will be aligned in linear time and therefore ignore theirlength.
The total number of cells in these rectangles isn兺x y ,i i

i= 0

which can be rewritten asn兺 关共x − x兲共y − y兲 + xy + x y − xy 兴
i

i= 0

i

i

i

= n ⭈ cov共xi,yi兲 + xL1 + L2 y − nxy= n ⭈ cov共xi,yi兲 + L1L2 Ⲑ n + L1L2 Ⲑ n − L1L2 Ⲑ n= L1L2 Ⲑ n + n ⭈ cov共xi,yi兲

where x̄ and ȳ are the means of each distribution, andcov(xi,yi) is the covariance.
The first term of this expression depends only on theeffective lengths of the sequences and the total number ofanchors. If we assume a lower bound on acceptable anchor
density, then L1L2/n behaves linearly in sequence length because L1/n and L2/n are O(1). The second term is at most nLAGAN: Multiple Alignment of Genomic DNA

␴x␴y, where ␴ denotes the standard deviation of the corresponding distribution. The upper bound is achieved in case of
perfect correlation in the spacing of anchors in the two sequences. Assuming constant anchor density—a reasonable assumption for a fixed pair of organisms—this upper bound islinear in sequence length provided the standard deviation ␴of each distribution is constant. Thus, if the anchors arespaced approximately evenly, and with a constant density,
the running time will be linear in sequence length.

Global Multiple Alignment, MLAGANMLAGAN aligns multiple genomic sequences in two main
phases, namely (1) a progressive alignment phase based on
LAGAN pairwise alignments, and (2) an optional iterative improvement phase. Next we describe in detail the progressivealignment phase, the scoring method used in MLAGAN, andthe iterative improvement phase.

Progressive Multiple Alignment with Anchors
During progressive alignment, MLAGAN aligns the sequencesin the order of the given phylogenetic tree. For example,
MLAGAN aligns sequences from human, chimpanzee, mouse,rat, and chicken, in the following order: (1) (human, chimpanzee), (2) (mouse, rat), (3) (human/chimpanzee, mouse/rat), (4) (human/chimpanzee/mouse/rat, chicken). Eachalignment step merges two sequences or alignments into alarger alignment, effectively building a profile of all the sequences. The scoring method used is similar to CLUSTALW’s
method, with some significant differences in the treatment ofgaps. We explain our scoring method in detail in the next
subsection.
In Step 2.1, MLAGAN aligns two (multi-)sequences X andY, to generate multi-sequence X/Y, by performing limited-area
dynamic programming (LAGAN Step 3). This step assumes
that the rough global map between the two (multi-) sequences is known, as guaranteed by Step 2.2. In Step 2.2,
MLAGAN recomputes the rough global map between the
newly generated multi-sequence X/Y and each (multi)sequence Z, in two steps:
1. The anchors between X/Y and Z are computed as follows:First, all anchors in the rough global maps between X andZ, and between Y and Z, become anchors between X/Y andZ, with score equal to their original score. Second, eachanchor between X and Z that overlaps an anchor between
Y and Z, is reweighed with score equal to (s1 + s2) ⳯ I/U,
where s1, s2 are the scores of the (X, Z) and (Y, Z) anchors,respectively, I is the length of intersection, and U is thelength of union of the anchors (summed in X/Y and Z;
Fig. 5).
2. The rough global map between X/Y and Z is the highestweight chain of these anchors, computed using the Longest Increasing Subsequence algorithm.
In Step 2.3, MLAGAN selects the next two closest (multi)sequences in the phylogenetic tree, and goes to Step 2.1.

Multiple Alignment Scoring with Affine Gaps
Scoring schemes for multiple alignments is an open researcharea, with some innovative solutions published recently, suchas T-COFFEE (Notredame et al. 2000). A full probabilistic
framework for modeling multiple sequence alignments was
implemented in HANDEL (Holmes and Bruno 2001). Twoclassical models are sum-of-pairs, and consensus. The sum-ofpairs model scores the multiple alignment according to thesum of scores of all pairwise alignments. The consensus model
creates a consensus string by a (weighted) majority vote ateach position, and scores the multiple alignment according to
the sum of pairwise scores between the consensus sequenceand each individual sequence (Gusfield 1999). Each scoringmodel has advantages and disadvantages. The sum-of-pairs

Figure 5 Generation of anchors during progressive alignment.
Multi-sequence X/Y is aligned to sequence Z. Anchors between X andZ (top) and anchors between Y and Z (middle) are remapped to coordinates in the X/Y multi-sequence, and given a new score. Then, theLongest Increasing Subsequence algorithm is applied to select a subset of the remapped anchors, as the anchors between X/Y and Z.

model, for example, penalizes each additional gapped sequence less than previous ones. Given a gap penalty of g, in a
K-sequence alignment the first gap in a given position gets
penalized by (Kⳮ1) ⳯ g, the second gap by (Kⳮ3) ⳯ g, andthe ith gap by (Kⳮ2i+1) ⳯ g. As a result, once a few gaps are
already open it becomes advantageous to open more gaps
whenever they result in small improvements in the number ofmatches. An alternative is to penalize gaps per sequence, but
this is also problematic because a deletion in one sequence issignificantly favored over an insertion in one sequence, leading to artificial compression of the multiple alignment. The
consensus model does not properly model different lineages:
A column with seven ‘T’s and five ‘A’s aligned, scores equally
to one with seven ‘T’s, two ‘C’s, two ‘G’s, and one ‘A’. We use
a combination of the approaches: sum-of-pairs for scoring
substitutions, and consensus (scaled appropriately) for scoring gaps. Our approach is most similar to CLUSTALW, whichuses a sum-of-pairs approach for scoring substitutions, andheuristically weighted per-sequence penalties to score gaps.
A straightforward implementation of the consensus affine-gap model in multiple alignments leads to the stacking
effect: Because gap-open penalties are large compared tomatch and mismatch scores, often it is favorable to artificiallyopen additional gaps in order to stack the gap openings (Fig.
6). To overcome this problem we use a gap-end penalty equal
to the gap-open penalty (effectively, we divide equally thepenalty for opening a gap into gap open and gap end). As canbe seen from Figure 7, introducing such a gap-end penalty
eliminates the stacking effect.
We define the multiple alignment score precisely as follows. Given a multiple alignment of K sequences, let (Aij) be
the K ⳯ L alignment matrix: Ai1 … AiL is the ith aligned sequence, where Aij ∈ {A, C, G, T, –} (‘–’ denotes a gap). Define
the K ⳯ L matrix (Bij) over the four letters ‘N’ (nucleotide), ‘O’(gap open), ‘G’ (gap continue), and ‘C’ (gap close): Bij = ‘G’ in
all gapped positions (where Aij = ‘–’), except the ones opening
a gap, in which case Bij = ‘O’; Bij = ‘N’ in all nucleotide positions (where Aij ⫽ ‘–’), except the ones closing a gap, in whichcase Bij = ‘C’.
Let m, d, g, and c be the match, mismatch, gap-open, andgap-continue penalties, respectively. Define the function S(x,y), where x, y ∈ {A, C, G, T, –}, as follows: S(x, y) = 0, if x = ‘–’or y = ‘–’; otherwise, S(x, y) = m if x = y, and S(x, y) = d if x ⫽ y.
The function S(x, y) will be used for the sum-of-pairs matchand mismatch component of the score.

Genome Researchwww.genome.org

729

Brudno et al.

Figure 6 The stacking effect. (A) The correct multiple alignment of
seven sequences. (B) The multiple alignment of the same seven sequences, with the stacking effect. When only gap-open and gapcontinue penalties are used, the stacked alignment (B) incurs one
gap-open penalty less than the correct alignment (A), and therefore is
optimal. When gap-end penalties are used, the correct alignment (A)is optimal because it incurs two gap-end penalties less than thestacked alignment (B).

refinement: Each sequence is realigned under the constraintthat the realignment stays within radius r from the originalalignment. This procedure improves the alignment locally,
but does not allow large-scale changes. MLAGAN introduces alimited-area version of iterative refinement (Step 3 of Algorithm MLAGAN) that has the advantages of performing morework in the areas that need it, and of allowing larger-scale
adjustments: Each sequence is removed iteratively, and every
region of the removed sequence that improved significantly
the score of the multiple alignment becomes an anchor. For
each position i in the sequence, the cumulative contribution
si of positions 1,…,i to the existing multiple alignment scoreis calculated by summing the score of position i to si–1. Whensi falls below 0 at position i, si is set to 0. When si reaches athreshold T, si is set to 0, and an anchor is created at position
i. The removed sequence is then re-aligned to the remaining
ones using the LAGAN algorithm with these anchors. LAGAN
allows for small changes around the anchors, while doing afull alignment search in between the anchors. Therefore, theabove procedure makes small improvements in well-aligning
regions and allows for larger adjustments between them.

Training, Testing, and Default Parameters

We trained the parameters of CHAOS, LAGAN, and MLAGAN
on the stem cell leukemia (SCL) region from human, mouse,chicken, fugu, and zebrafish (Göttgens et al. 2002). During
training we did not use the ROSETTA dataset, or the CFTR
region, where we tested. One exception is a number of softth
Let Nj be the number of ‘N’s in the j column of (Bij).
ware bug fixes that were specific to the size of the dataset,
Similarly, define O j , G j , and C j . Define the function
such as memory leaks, which were corrected after applying
T(i) = min(Oi, KⳮOi) ⳯ (g + c) + min(Gi, KⳮGi) ⳯ c + min(Ci,LAGAN and MLAGAN to the CFTR sequences.
KⳮCi) ⳯ g. The function T(i) will be used for the consensusLAGAN and MLAGAN compute the rough global map
based gap component of the score. The multiple alignment
using recursive calls to CHAOS with parameters (k, c, t), wherescore is then:k is the k-mer size, c is the allowed degeneracy, and t is thethreshold (see Methods: Generation of Local Alignments).
SCOREMLAGAN关共Aij兲兴 =S共aij,aik兲
共K − 1兲T共i兲 +The default recursive steps are (12, 0, 30), (13, 1, 30), (8, 1, 30),
1ⱕiⱕL
1ⱕjⱕK j<kⱕKand (7, 1, 30), on sequences masked using RepeatMasker (Smitand Green). These steps are followed by a (7, 1, 30) step on the
Iterative Refinement With Anchors
unmasked sequences in order to anchor the conserved repeats. An optional (5, 0, 50) step on the translated sequences
One shortcoming of progressive alignment is that the initial
can be inserted between the (13, 1, 30) step and the (8, 1, 30)
pairwise alignments are fixed, and early errors cannot be corstep. The translated anchoring step is scored using a
rected later, when more information is available. A simple
BLOSUM62 matrix. This step is disabled by default, and was
example of this situation is shown in Figure 7, where sedisabled during the testing reported in Tables 1 and 2. The
quences 5 and 6 are initially misaligned. To fix such errors we
distance d parameter is 20 for genomic sequences and 8 foruse a technique known as iterative refinement.
comparison of peptides (translated). The shift s is 5 for all
In standard iterative refinement, as used by CLUSTALW,
runs.
each sequence is removed from the multiple alignment andDuring the limited-area computation of global alignre-aligned, yielding a better alignment. This process is rement (Step 3 of LAGAN, Fig. 4), the default radius is r = 15.
peated for a number of iterations, or until the score cannot be
The limited-area dynamic programming pass (Step 3 offurther improved (a local maximum is reached). Anson and
LAGAN) uses match m = +12 (or +18 for MLAGAN), mismatch
Myers (1997) introduced a limited-area version of iteratived = ⳮ8, gapopen g = ⳮ100, andgapcontinue c = ⳮ5. In MLAGAN,
the gapopen parameter gets split
equally, as explained in Section
3.2.2, into gapopen = ⳮ50 andgapend = ⳮ50.
All of the tests were run on a
2.3-GHz Pentium IV machine with
500MB of RAM, except for AVID on
the CFTR dataset, which due to thememory requirements was run on a
2-GHz Pentium IV machine with
1GB of RAM.
Most of the other programs wetested can take many different parameters. For all programs we used
the parameters that were either sugFigure 7 Part of a multiple alignment between eight mammalian sequences, before and after iteragested by the authors, or gave thetive refinement. Refinement improves the overall alignment by correctly placing the two T’s in the fifth
best results for that particular proand sixth sequence, and changing a 1-bp deletion in sequence 5 into an insertion in sequence 6.

兺冉

730

Genome Researchwww.genome.org

兺 兺冊LAGAN: Multiple Alignment of Genomic DNA

gram. Specifically: BLASTZ was run with a cutoff score of 2200
(K = 2200); DIALIGN used CHAOS anchoring, and was run
with the ‘ⳮnt’ (compare at the protein-coding level) option;
AVID version 0.91 was used, as it performed better than AVID
1.2 on the CFTR dataset. (AVID 1.2 performed slightly betteron the ROSETTA dataset.).

ACKNOWLEDGMENTS
M.B. was supported by an NSF Graduate Fellowship. G.M.C.
was supported by a Howard Hughes Medical Institute predoctoral fellowship. E.D.G. and the NISC Comparative Sequencing Program were supported by funds from the National Human Genome Research Institute. The following individuals
were key contributors within the NISC Comparative Sequencing Program: Jim Thomas (BAC isolation and mapping); Jeff
Touchman and Bob Blakesley (BAC sequencing); Gerry Bouffard, Steve Beckstrom-Sternberg, Pam Thomas, Jenny McDowell, and Baishali Maskeri (computational analyses). We thank
Berthold Göttgens and Michael Chapman for providing uswith the stem cell leukemia region, Burkhard Morgenstern forhelp with developing the CHAOS program, and Inna Dubchak and Alex Poliakov for testing LAGAN and providing uswith feedback.
The publication costs of this article were defrayed in partby payment of page charges. This article must therefore be
hereby marked “advertisement” in accordance with 18 USC
section 1734 solely to indicate this fact.

REFERENCES
Altschul, S.F., Gish, W., Miller, W., Myers, E.W., and Lipman, D.J.
1990. Basic local alignment search tool. J. Mol. Biol.
215: 403–410.
Altschul, S.F., Madden, T.L., Schoffer, A.A., Zhang, J., Zhang, Z.,
Miller, W., and Lipman, D.J. 1997. Gapped BLAST and
PSI-BLAST: A new generation of protein database search
programs. Nucleic Acids Res. 25: 3389–3402.
Anson, E.L. and Myers, E.W. 1997. Re-Aligner: A program forrefining DNA sequence multialignments. J. Comp. Biol.
4: 369–383.
Barton, G.J. and Sternberg, M.J.E. 1987. A strategy for the rapid
multiple alignment of protein sequences. J. Mol. Biol.
198: 327–337.
Batzoglou S., Pachter, L., Mesirov, J., Berger, B., and Lander, E.S.
2000. Human and mouse gene structure: Comparative analysis
and application to exon prediction. Genome Res. 10: 950–958.
Bonizzoni, P. and Vedova, G.D. 2001. The complexity of multiple
sequence alignment with SP-score that is a metric. TheoreticalComputer Science 259: 63–79.
Bray, N., Dubchak, I., and Pachter, L. 2003. AVID: A globalalignment program. Genome Res. 13: 97–102.
Brudno, M. and Morgenstern, B. 2002. Fast and sensitive alignmentof large genomic sequences. Proceeding of the IEEE Computer
Society Bioinformatics Conference (CSB).
Delcher, A.L., Kasif, S., Fleischman, R., Peterson, J., White, O., andSalzberg, S.L. 1999. Alignment of whole genomes. Nucleic Acids
Res. 27: 2369–2376.
Delcher, A.L., Phillippy, A., Carlton, J., and Salzberg, S.L. 2002. Fast
algorithms for large-scale genome alignment and comparison.
Nucleic Acids Res. 30: 2478–2483.
Dubchak, I., Brudno, M., Loots, G.G., Pachter, L., Mayor, C., Rubin,E.M., and Frazer, K.A. 2000. Active conservation of noncoding
sequences revealed by three-way species comparisons. Genome
Res. 10: 1304–1306.
Eppstein, D., Galil, Z., Giancarlo, R., and Italiano, G.F. 1992. Sparse
dynamic programming I: Linear cost functions. JACM
39: 546–567.
Gotoh, O. 1996. Significant improvement in accuracy of multiple
protein sequence alignments by iterative refinement as assessedby reference to structural alignments. J. Mol. Biol. 264: 823–838.
Göttgens, B., Barton, L.M., Chapman, M.A., Sinclair, A.M., Knudsen,
B., Grafham, D., Gilbert, J.G.R., Rogers, J., Bentley, D.R., andGreen, A.R. 2002. Transcriptional regulation of the stem cellleukemia gene (SCL)—Comparative analysis of five vertebrate
SCL loci. Genome Res. 12: 749–759.
Gusfield, D. 1999. Algorithms on strings, trees, and sequences: Computerscience and computational biology, pp. 351–353. Cambridge
University Press, Cambridge, UK.
Hardison, R., Chao, K.-M, Adamkiewicz, M., Price, D., Jackson, J.,
Zeigler, T., Stojanovic, N. and Miller, W. 1993. Positive andnegative regulatory elements of the rabbit ⑀-globin gene revealedby an improved multiple alignment program and functional
analysis. DNA Seq. 4: 163–176.
Hardison, R., Chao, K.-M., Schwartz, S., Stojanovic, N., Ganetsky,
M., and Miller, W. 1994. Globin Gene Server: A prototype E-maildatabase server featuring extensive multiple alignments and datacompilation for electronic genetic analysis. Genomics
21: 344–353.
Henikoff, S. and Henikoff, J.G. 1992. Amino acid substitution
matrices from protein blocks. Proc. Natl. Acad. Sci.
89: 10,915–10,919.
Höhl, M., Kurtz, S., and Ohlebusch, E. 2002. Efficient multiple
genome alignment. Bioinformatics 18: 312–320.
Holmes, I. and Bruno, W.J. 2001. Evolutionary HMMs: A Bayesian
approach to multiple alignment. Bioinformatics 17: 803–820.
Kent, W.J. 2002. BLAT—The BLAST-like alignment tool. Genome Res.
12: 656–664.
Kent, W.J. and Zahler, A.M. 2000. Conservation, regulation, synteny,and introns in a large-scale C. briggsae–C. elegans genomicalignment. Genome Res. 10: 1115–1125.
Mayor, C., Brudno, M., Schwartz, J.R., Poliakov, A., Rubin, E.M.,
Frazer, K.A., Pachter, L.S., and Dubchak, I. 2000. VISTA:
Visualizing global DNA sequence alignments of arbitrary length.
Bioinformatics 16: 1046.
Miller, W. 2001. Comparison of genomic DNA sequences: Solved
and unsolved problems. Bioinformatics 17: 391–397.
Morgenstern, B. 1999. DIALIGN 2: Improvement of thesegment-to-segment approach to multiple sequence alignment.
Bioinformatics 15: 211–218.
Morgenstern, B., French, K., Dress, A., and Werner, T. 1998.
DIALIGN: Finding local similarities by multiple sequencealignment. Bioinformatics 14: 290–294.
Needleman, S.B. and Wunsch, C.D. 1970. A general method
applicable to the search for similarities in the amino acidsequence of two proteins. J. Mol. Biol. 48: 443–453.
Ning, Z., Cox, A.J., and Mullikin, J.C. 2001. SSAHA: A fast searchmethod for large DNA databases. Genome Res. 11: 1725–1729.
Notredame, C., Higgins, D.G., and Heringa, J. 2000. T-coffee: A
novel method for fast and accurate multiple sequence alignment.
J. Mol. Biol. 302: 205–217.
Sankoff, D. 1999. Genome rearrangement with gene families.
Bioinformatics 15: 909–917.
Schwartz, S., Zhang, Z., Frazer, K.A., Smit, A., Riemer, C., Bouck, J.,
Gibbs, R., Hardison, R., and Miller, W. 2000. PipMaker—A web
server for aligning two genomic DNA sequences. Genome Res.
10: 577–586.
Simon A., Stone, E.A., and Sidow, A. 2002. Inference of functional
regions in proteins by quantification of evolutionary constraints.
Proc. Natl. Acad. Sci. 99: 2912–2917.
Smith, T.F. and Waterman, M.S. 1981. Identification of common
molecular subsequences. J. Mol. Biol. 147: 195–197.
Stanfel, L.E. 1996. A new approach to clustering the amino acids. J.
Theor. Biol. 183: 195–205.
Taylor, W. 1988. A flexible method to align large numbers of
biological sequences. J. Mol. Evol. 28: 161–169.
Thompson, J.D., Higgins, D.G., and Gibson, T.J. 1994. CLUSTALW:
Improving the sensitivity of progressive multiple sequencealignment through sequence weighting, position specific gap
penalties and weight matrix choice. Nucleic Acids Res.
22: 4673–4680.
Wang, L. and Jiang, T. 1994. On the complexity of multiple
sequence alignment. J. Comp. Biol. 1: 337–348.

WEB SITE REFERENCES
http://lagan.stanford.edu; LAGAN alignment server.
http://ftp.genome.washington.edu/RM/RepeatMasker.html;
RepeatMasker (Smit, A.F.A. and Green, P.).
Received October 23, 2002; accepted in revised form December 11, 2002.

Genome Researchwww.genome.org

731