Vol. 19 Suppl. 1 2003, pages i54–i62
DOI: 10.1093/bioinformatics/btg1005

BIOINFORMATICS

Glocal alignment: finding rearrangements during
alignment
Michael Brudno 1,∗,†, Sanket Malde 1,†, Alexander Poliakov 2,
Chuong B. Do 1, Olivier Couronne 2, Inna Dubchak 2 andSerafim Batzoglou 1
1 Department of Computer Science, Stanford University, Stanford, CA 94305 and2 Life Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA

94720, USA
Received on January 6, 2003; accepted on February 20, 2000

ABSTRACT
Motivation: 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
biological features between distant species. The two mainclasses of pairwise alignments are global alignment,
where one string is transformed into the other, and localalignment, where all locations of similarity between thetwo strings are returned. Global alignments are less prone
to demonstrating false homology as each letter of one
sequence is constrained to being aligned to only one letterof the other. Local alignments, on the other hand, can copewith rearrangements between non-syntenic, orthologous
sequences by identifying similar regions in sequences;this, however, comes at the expense of a higher false
positive rate due to the inability of local aligners to take
into account overall conservation maps.
Results: In this paper we introduce the notion of glocalalignment, a combination of global and local methods,
where one creates a map that transforms one sequence
into the other while allowing for rearrangement events. We
present Shuffle-LAGAN, a glocal alignment algorithm thatis based on the CHAOS local alignment algorithm and the
LAGAN global aligner, and is able to align long genomic
sequences. To test Shuffle-LAGAN we split the mouse
genome into BAC-sized pieces, and aligned these pieces
to the human genome. We demonstrate that ShuffleLAGAN compares favorably in terms of sensitivity andspecificity with standard local and global aligners. From thealignments we conclude that about 9% of human/mouse
homology may be attributed to small rearrangements, 63%of which are duplications.
Availability: Our systems, supplemental information, andthe alignment of the human and mouse genomes using∗ To whom correspondence should be addressed.
†

These authors contributed equally to the work.

i54

Shuffle-LAGAN are available at
http://lagan.stanford.edu/glocal.
Contact: glocal@cs.stanford.edu

1 INTRODUCTION
The availability of complete vertebrate genomes, such as
the human, mouse, and fugu (Lander et al., 2001; Venter etal., 2001; Waterson et al., 2002; Aparicio et al., 2002) has
pushed comparative genomics into a new era. Comparing
genomic sequences from related species is a fruitful sourceof biological insight, as functional elements such as genes
and regulatory sites tend to exhibit significant sequencesimilarity across related genomes, whereas regions thatare not functional tend to be non-conserved (Dubchak etal., 2000; Göttgens et al., 2002). Performing alignment
on a whole genome scale introduces new challenges, as
it is no longer feasible to do manual post-processing of thealignments to ensure their correctness.
The availability of multiple genomic sequences also
makes it possible to learn more about the mutations
that DNA undergoes during evolution, and thus create
better alignment programs. In particular, most current
programs penalize insertions, deletions and substitutionsof individual base pairs between two sequences (we callthese events simple edits). It is known, however, that DNA
also mutates by various rearrangement events, such astranslocations (a subsegment is removed and inserted ina different location but the same orientation), inversions(a subsegment of DNA is removed from the sequenceand then inserted back in the same location but opposite
orientation), duplications (a copy of a subsegment isinserted into the sequence, the original subsegment isunchanged), or a combination of the above (See Fig. 1). In
this paper we explore local rearrangement events, which
we define to be longer than 100 base pairs (bp) and shorter
than 100 thousand bp (Kbp).
Of the two main methods for pairwise alignment—c Oxford University Press 2003; all rights reserved.
Bioinformatics 19(1) Glocal alignment: finding rearrangements during alignment

Fig. 1. Common rearrangement events. Top left corner corresponds start of both sequences, diagonal lines to areas of homology.
(A) No rearrangements (consistent homology). (B) Translocation. (C) Inversion. (D) Duplication.

global alignment, which shows how one sequence canbe transformed into another using a combination of
the simple edits, and local alignment, which identifies
local similarities between regions of sequences—neither
handles rearrangement events satisfactorily. Globalalignment algorithms such as Needleman and Wunsch(1970), Dialign (Morgenstern, 1999), MUMmer (Delcher
et al., 1999, 2002), Avid (Bray et al., 2003), and LAGAN(Brudno et al., 2003) do not handle these events at all: themap between the two sequences that a global alignment
algorithm creates must be monotonically increasing inboth sequences. While local alignment methods such as
Smith and Waterman (1981), BLAST (Altschul et al.,1990, 1997), CHAOS (Brudno and Morgenstern, 2002)
and BLASTZ (Schwartz et al., 2003) are able to identify
homology in the presence of rearrangements between two
sequences, they do not suggest how the two sequences
could have evolved from their common ancestor. Localalignment algorithms also dependent on a cutoff score:
only those alignments that score above a threshold are
accepted. Deciding on this cutoff can be difficult: if it is
set too high, some significant hits will be missed, and if it
is set too low there will be too many false positives. Also,in the case where both sequences have n paralogs (copies)
of a particular gene or feature, local aligners return n 2
local alignments between all of the pairs, whereas a
simple global alignment may more clearly reflect theevolutionary process.
Most work on detecting large scale rearrangements forwhole genomes has been done at the level of genes andchromosomes, not genomic sequences. Here each gene isrepresented by a unique letter, and a string of these letters
is a ‘chromosome.’ Hannenhalli and Pevzner (1995) were
the first to present an algorithm to find the minimal numberof inversions required to transform one chromosome
into another. This method has been extended to otheroperations, such as translocations, block interchanges andtranspositions (Hannenhalli, 1996; Christie, 1996; Bafna
and Pevzner, 1998). More recently Pevzner and Tesler(2003) have investigated the rearrangements that haveoccurred between the human and mouse genome, andfound that significant stretches of human/mouse homology
are due to micro-rearrangements, of size less than 1
million base pairs (Mbp).
Varré et al. (1999) proposed a distance metric between
two DNA sequences that models various rearrangement
events. Their algorithm, called Tnt1, builds a second
sequence from an initially empty string using insertions
and copying of blocks from the first sequence. The
distance between the two strings is defined to be theKolmogorov complexity of the program that builds thesecond sequence. This algorithm has several shortcomings, the most notable being its inability to handle simple
edit operations.
An alignment tool for genomic sequences that models
rearrangements must meet several criteria, the most important of which is speed: an aligner for genomic sequences
should be able to rapidly align sequences of length 1 Mbp
or more. Likewise, it is important to use an objective function that correctly penalizes various rearrangement eventsand simple edits. In this work we introduce the problem ofglocal alignment, a hybrid of global and local alignments,and specify what features a glocal aligner should have.
We also present Shuffle-LAGAN (SLAGAN), a glocalaligner capable of quickly aligning long genomic sequences. SLAGAN is based on the CHAOS and LAGAN
aligners (Brudno and Morgenstern, 2002; Brudno et al.,2003), and has been used to align the mouse genome,
split into BAC-sized pieces, to the human genome. We
show that SLAGAN, when compared with the regular
LAGAN global aligner, has better sensitivity and similar
specificity. When compared to the BLASTZ local aligner,
SLAGAN shows better specificity at a modest cost insensitivity. Using the results produced by SLAGAN, we
discuss the characteristics of the local rearrangement
events that occured since the divergence of the human andmouse genomes.
i55

M.Brudno et al.

2 ALGORITHMS
In this section we present a definition of glocal alignment,
and present SLAGAN, a glocal alignment algorithm that
we have implemented. SLAGAN is more sensitive thanregular global aligners, yet it has better specificity thanlocal aligners. Finally SLAGAN is able to rapidly alignlong genomic sequences of length greater than 1 Mbp. The
total time to align the human and mouse genomes using
SLAGAN and the technique of Couronne et al. (2003) ona modern cluster was about 12 hours of wall clock time
(25 CPU-days).
2.1 Glocal alignment
A glocal alignment between two sequences is a series ofoperations that transform one sequence into the other. We
believe the necessary set of operations includes insertions,deletions, point mutations, inversions, translocations andduplications. Each operation incurs a penalty, and the total
edit distance between the two sequences is the sum of
the penalties. It may also be advantageous to penalize
not only individual operations, but also combinations; forinstance an inverted translocation should possibly havea smaller penalty than the sum of an inversion penalty
and a translocation penalty. Finally, a glocal alignment
algorithm should be symmetric in the sequence order:
the resulting alignment should not depend on the orderin which the sequences are given. One of the main
shortcomings of the SLAGAN algorithm presented below
and the aforementioned Tnt1 algorithm is that neither
is symmetric in the sequence order. This is especially
noticeable when one aligns duplicated regions, as both
SLAGAN and Tnt1 will currently report duplications in
only one of the sequences.
This list of allowed operations is not meant to becomplete: as we learn more about the evolution of DNA
sequences, additional operations may become necessary.
It is also difficult to find good parameters for the various
penalties, as currently there is no sound mathematicalbasis for glocal alignments. Both these areas are promising
topics for future work.
2.2 The Shuffle-LAGAN algorithm
The SLAGAN algorithm consists of three distinct stages.
During the first stage the local alignments between thetwo sequences are found using the CHAOS tool (Brudno
and Morgenstern, 2002; Brudno et al., 2003). Second,
the maximal scoring subset of the local alignments undercertain gap penalties is picked to form a 1-monotonic
conservation map. It is the structure of this map thatmakes SLAGAN different from standard anchored globalaligners. Finally, the local alignments in the conservation
map that can be part of a common global alignment
are joined into maximal consistent subsegments, which
are aligned using the LAGAN global aligner (Brudno eti56al., 2003). See Figure 2 for a graphical overview of thealgorithm. The exact parameters used in the algorithm are
specified in Appendix 1.
2.2.1 Generation of local alignments. To generate local alignments between the two sequences, SLAGAN uses
CHAOS, a method that finds small matching words withdegeneracy, and chains them into local alignments. Here,
we only summarize the CHAOS algorithm. (See Brudno
and Morgenstern, 2002; Brudno et al., 2003).
CHAOS works by chaining short words, the seeds,
which match between the two sequences. Given a word
length k, and a degeneracy c, a (k, c)-seed is a pair of klong words (k−mers) that match with at most c differences
between the two sequences. Given a maximum distance d,and maximum shift s, two seeds that are x- and y-lettersapart in the first and second sequences, respectively, canbe chained together if x  d, y  d, and |x − y|  s.
A seed is chained to the single previous seed that creates
the highest scoring chain among all chains that end withthis seed. CHAOS also supports a translation option, in
which both nucleic sequences are translated in all 6 coding
frames (3 forward and 3 reverse), and all combinationsof frames are compared in turn. Amino acids are grouped(Stanfel, 1996), and all amino acids in the same group are
considered equal. Finally, the chains are extended using
the standard ungapped BLAST (Altschul et al., 1990)
extension, until the score drops below a certain threshold.
After computing the chains, CHAOS scores each chain
by a rapid scoring mechanism: ungapped extensions are
performed in both directions from each seed, and the
optimal location to insert a gap of length exactly |x − y|is found. The resulting alignment is scored using astandard Needleman–Wunsch edit metric with penalties
described in Appendix 1. If the alignment is done on
amino acid sequences, the rescoring is done using a
BLOSUM62 matrix (Henikoff and Henikoff, 1992).
2.2.2 Building the 1-monotonic conservation map.
Most tools for rapid global alignment start with a setof local alignments, which they resolve into a ‘roughglobal map’ (Delcher et al., 1999; Batzoglou et al., 2000;
Bray et al., 2003; Brudno et al., 2003). The rough global
map must be non-decreasing in both sequences. In order
to allow our alignment algorithm to catch rearrangements, we relax this assumption to allow the map to be
non-decreasing in only one sequence, without putting
any restrictions on the second sequence. We call this a
1-monotonic conservation map. Here we show how tocreate a 1-monotonic conservation map under the affine
gap model.
For this portion of the algorithm we represent each
local alignment L generated in the previous section as
L = (star t1 , end 1 , star t 2 , end 2 , score, strand), a vector
with six fields: the start and end positions of the local

Glocal alignment: finding rearrangements during alignment

Fig. 2. Overview of the SLAGAN Algorithm. (A) The local alignment between the two sequences are generated using CHAOS. (B) The
highest scoring 1-monotonic map (indicated in bold) is found. (C) The maximal consistent subsegments of the 1-monotonic map (dashedboxes) are aligned using LAGAN.

alignment in the two sequences, the score of the alignment,
and the strand on which the alignment occurs in the second
sequence. We define local alignments on the positive
strand to have both start positions less than their respective
end positions; on the negative strand, the end position insequence 2 is smaller than the start position.
Consider two local alignments, L 1 and L 2 . We call
L 1 and L 2 1-monotonic if L 2 .star t1 > L 1 .end1 .
We call these alignments consistent if (1) they are 1monotonic, (2) they are both on the same strand, and(3) L 2 .star t 2 > L 1 .end 2 for alignments on the positive
strand or L 2 .star t 2 < L 1 .end 2 for alignments on thenegative strand. An ordered list of local alignments[L 1 . . . L k ] is 1-monotonic or consistent if for any pairof local alignments L i , L j if i < j then L i and L j
are 1-monotonic or consistent, respectively. Intuitively,a list of local alignments is 1-monotonic if it is strictlyincreasing in the first sequence. The list is consistent if thealignments in it are strictly increasing in both sequences
and all on the positive strand, or strictly increasing in thefirst sequence and strictly decreasing in the second andall are on the negative strand. Equivalently, a list of localalignments is consistent if all of them can be chained
into a single global alignment. A consistent, ordered listof local alignments will consequently be referred to as aconsistent chain, or a consistent subchain.
The problem of finding the highest scoring consistentsubset of local alignments can be solved using standard
dynamic programming in time O(n 2 ), by chaining eachalignment L j to the highest scoring consistent subchain
ending in some alignment L i such that L i and L j areconsistent. Eppstein and colleagues (Eppstein et al., 1992)
have shown how to use the sparse dynamic programming technique to speedup the algorithm to run in timeO(n log n). The key observation is that at any point in the
dynamic programming, the set of optimal consistent subchains considered so far partitions the two-dimensional
space of (star t 1 , star t2 ) coordinates into ‘influence areas’ such that any local alignment in a particular influence
area must be chained to a particular previous optimal
consistent subchain. By considering local alignments
in order of increasing star t 1 , we need only perform a
one-dimensional search to find the appropriate influence
area for chaining each local alignment; more specifically,
this may be accomplished in O(log n) time by storing
the active influence areas sorted by diagonal number(star t 1 –star t 2 ) in a balanced binary tree (see Eppstein etal., 1992 for details). We generalize this result to finding
the highest scoring 1-monotonic subset of alignments, also
under the affine gap model. Here we present an overview
of the algorithm (a more detailed description is in our
online supplement, http://lagan.stanford.edu/glocal/).
Within the new model, the assumptions that all localalignments in a chain lie on the same strand and must be
consistent are relaxed. In particular, when chaining localalignment L 2 to a subchain ending in L 1 , L 1 may be on
either the positive or negative strand (+/−), L 1 may comeeither before (L 1 .end 2 < L 2 .star t 2 ) or after (L 1 .end 2 >L 2 .star t 2 ) L 2 in sequence 2 (+/−), and L 2 may be on
either the positive or negative strand (+/−). We represent
the eight possible combinations of L 1 .strand, direction,and L 2 .strand by 3 bits, in that order.
The original Eppstein O(n log n) algorithm dealt withthe special case where an alignment on the positive strand
is always joined to a preceding subchain on the same
strand (+ + +). In the generalized model, however, analignment need not have the same strand or direction as
the subchain to which it is joined; more precisely, fori57

M.Brudno et al.

each local alignment L 2 on a given strand L 2 .strand, wemust consider joining it to consistent subchains that areeither on the positive or negative strand (L 1 .strand) andthat either precede or follow it in sequence 2 (direction).
Finding the optimal consistent subchain with which to join
L 2 in each of the four possible combinations of L 1 .strand
and direction may be done in O(log n) time using thebalanced binary trees as in the original Eppstein algorithm.
Taking the maximum over the four cases then gives the
optimal chaining for L 2 . Since the gap penalties usedand hence the influence regions depend on L 2 .strand
as well, a total of eight balanced binary trees (for eachcombination of L 1 .strand, direction, and L 2 .strand)
must be maintained, of which only the four consistent
with each L 2 .strand are used in determining the proper
chaining for L 2 . In a sense, we may regard this as runningeight parallel versions of Eppstein’s algorithm at once,building a single common chain.
In practice we use only four different gap penalties, as
we eliminate the symmetric cases, e.g. if both strands are
positives and chained in the positive direction (+ + +),they should be penalized the same way as negative strands
chained in the negative direction (− − −). One can thinkof the four gap penalties as corresponding to the regular
gap penalties (+ + + and − − −), inversions (+ + −and − − +), translocated inversions (+ − − and − + +)and translocations (+ − + and − + −). The gap penalty
charged for all transitions consists of three parts: (1)the gap open penalty is charged for the consistent cases(+ + + and − − −) if the two segments are on different diagonals (L 1 .end 1 –L 1 .end 2  =L 2 .star t 1 –L 2 .star t 2 ),and is always charged for the other cases (in these cases
it is the inversion, translocation or the inverted translocation penalty); (2) the gap continue penalty is equalto|(L 1 .end 1 –L 1 .end 2 )–(L 2 .star t1 –L 2 .star t 2 )| ×
constant; and (3) the distance between two alignments
is defined to be min(|L 1 .end 1 –L 2 .star t 1 |, |L 1 .end 2 –L 2 .star t 2 |), and it is also penalized as (distance ×
constant). Finally, in order to not pay both an inversion
and a translocated inversion penalty for a simple inversion,
we add a ‘momentum heuristic’: if a majority of the last
five elements in the chain have a different strand than thecurrent local alignment, and the case we are considering isa translocated inversion (+ − − or − + +), the regular gap
penalty is charged instead. This is done in order to charge
inversions a smaller penalty than translocated inversions(see Appendix 1 for the exact penalties used).
2.2.3 Aligning consistent subsegments. Recall that two
local alignments are considered to be consistent if they
can both be a part of a global alignment. Once we havea 1-monotonic conservation map it is straight-forward togenerate the maximal consistent subsegments of the map
by simply sorting all of the local alignments in the 1i58

monotonic map by their start1 coordinates, taking thefirst alignment to be the start of a consistent subsegment,
and adding additional local alignments while they are
all consistent. As soon as an alignment is found to beinconsistent with the current subsegment, we start a newsubsegment.
Because the edges of the local alignment may not be
reliable borders of the end of homology (there may be
additional conservation which was not strong enoughto meet the local alignment criterion) the start and end
positions of consistent subsegments are expanded in thefirst sequence to the borders of the previous and next
consistent subsegments, respectively. The borders are also
expanded in the second subsequence, in proportion to
how much they were expanded in the first sequence.
The expanded consistent subsegments can now be aligned
using any global alignment algorithm.
Note that because of the expansion, the consistentsubsegments will now overlap in the first sequence. In
order to get a true glocal alignment between the twosequences, it is now necessary to clip the overlappingalignments so that they do not overlap in the first sequence.
This can be done using two linear passes over thealignments, one forward and one backward, to find the
optimal point at which to end one alignment and start theother.

3 RESULTS
We tested SLAGAN by aligning the human and mouse
genomes using the whole genome alignment techniqueused in the Berkeley Genome Pipeline (Couronne et al.,2003). In this technique the mouse genome is split up
into contigs of 250 Kbp. The potential human orthologs
for each contig are found using the BLAT aligner (Kent,
2002). The human sequence is then extended around the
BLAT anchor, and aligned to the mouse contig using
the tested aligner. If the aligner being used is a globalaligner and the BLAT hits fall on both strands, then thealigner is called both with the original mouse contig anda reverse-complemented copy, making it possible to catchinversions. However, a global aligner in this context would
not be able to deal with small scale translocations or
duplications on the same strand. When using SLAGAN
about 40% of all local alignments created by CHAOS are
eliminated while creating the 1-monotonic conservation
map, indicating that the homology map is ‘cleaned up’
compared to just having all local alignments. This allows
for more sensitive settings when generating the localalignments.
3.1 Quality of SLAGAN alignments
In order to evaluate the quality of SLAGAN alignments
we have used the metrics developed to evaluate thevarious alignment programs for the mouse genome paper

Glocal alignment: finding rearrangements during alignment

Table 1. Sensitivity and specificity of Shuffle-LAGAN compared to a localaligner (BLASTZ) and a global aligner (LAGAN)
Dataset

BLASTZLAGAN

all human v. all mousehuman chr 20 v. all mousehuman chr 20 v. mouse chr 2% non-orthologous

39.6
40.5
37.2
8.1%

36.5
41.6
41.4
0.5%

Shuffle-LAGAN
38.2
42.5
42.4
0.2%

The first three rows are the total coverage of the human dataset byalignments with the mouse dataset. The last line is the percentage of
non-onhologous alignments on human chromosome 20(Waterson et al., 2002). In particular, we consider thesensitivity of an alignment program to be the percentage
of base pairs in the alignment it produces that meet aparticular Smith–Waterman scoring threshold. The exact
matrices we use (‘total’ and ‘tight’) are identical to thoseused to evaluate the mouse genome aligners (see Appendix
1). To test the specificity of the aligners we applied thetechnique used to evaluate the specificity of BLASTZ: thecoverage of human chromosome 20.
Because human chromosome 20 is considered to becompletely ortholgous to mouse chromosome 2, very little
difference should be seen between the coverage of humanchromosome 20 by mouse chromosome 2 and by thewhole mouse genome. The results are summarized inTable 1. The overall conclusion is that while SLAGAN
may not be quite as sensitive as BLASTZ on the wholegenome scale, it is able to align a larger percentage ofthe orthologous regions of human chromosome 20 andhas higher specificity as shown by the lower percentage
of non-orthologous alignments. As compared to regular
LAGAN, SLAGAN is more sensitive, especially foraligning genes and nearby areas (see Table 2), which are
the areas of the genome that are more likely to undergo
duplications in order to evolve new function. SLAGAN
appears to be slightly more specific than LAGAN, though
the difference may not be statistically significant.

3.2 Analysis of rearrangements
After aligning the whole human and mouse genomes witha glocal aligner, it becomes possible to characterize theextent to which BAC-sized chunks in the human genome
have maintained a linear order, and by this whether
global alignments are an appropriate method to analyze
BAC-sized chunks of the genome. Our results (Table 2)
indicate that overall, the homologies found on the wholegenome level by global and glocal alignments are very
similar. There is, however, a noticeable difference when
one considers the alignment of genes and gene-related
elements, such as coding regions, UTRs and elements
close to genes, where one sees a significant improvementof about 2% in total coverage. This result suggeststhat as much as 2% of the gene coding regions in thehuman genome may have evolved by local translocationor duplication since the human/mouse divergence. The
evolution of genes by local duplication is supported by
anecdotal evidence of co-location of homologous genes
such as the globin cluster (Flint et al., 2001).
Once the 1-monotonic map is found, it is possible toclassify the various types of rearrangements on it. We
generate the highest scoring consistent map of all the localalignments above a threshold CHAOS score in order todetermine whether the main conservation is on the positiveor negative strand. Using this map, we classify any localalignment that lies on the opposite strand as an inversion.
We label as a translocation any local alignment that is on
the main strand, but not part of the consistent chain, as
well as any alignment on the weak strand that could not
be part of the main consistent chain had it been on themain strand. Finally, if a particular alignment covers atleast 70% of another alignment in the mouse sequence,
the lower scoring of the two local alignments is labeled
a duplication. Because SLAGAN is not symmetric only
the duplications in mouse are found. Thus, the number of
duplications may be underestimated by as much as a factorof two. We would also not be able to locate translocations
and duplications that have been split between two contigs.
Table 3 summarizes the different rearrangement events,
their proportion in the human genome, and the level towhich they are conserved, both overall and per base pair.
It is notable that duplications as a whole score lower
per base pair (42.7) than sequences that have undergone
other rearrangements (54.4). This could be the result of
duplicated sequences being free to mutate to evolve newfunction, while sequences rearranged but not duplicatedbeing more constrained to their previous function. Of thenon-duplicated sequences, simple inversions tend to be theshortest (average length of 122 bp), and have the highest
score per base pair (55.3), followed by translocations(188 bp, 46.8).

4 DISCUSSION
Sequence alignment is one of the oldest and most successful applications of computer science to biology. Despite
the considerable advances achieved after several decadesof research in this area, many important challenges
remain. One of these challenges is the developmentof systems for aligning sequences in the presence ofrearrangements. This is becoming increasingly importantwhen automatically comparing whole genomes of related
species.
We developed SLAGAN, a system suitable for highthroughput reliable alignment of genomic sequences in
the presence of rearrangements. SLAGAN is based on
LAGAN, a pairwise aligner that is designed for rapid andreliable alignment of a pair of genomic sequences that
may be very similar (e.g. human and chimpanzee), or veryi59

M.Brudno et al.

Table 2. Coverage of various features of the human genome by BLASTZ, LAGAN, and Shuffle-LAGAN. The tight matrix is the one from Schwartz et al.
(2003). Coding, UTR and Upstream500 refer to the corresponding areas of Refseq genes

BLASTZ
Coverage
Overall
Coding (CDS)
UTRs
Upstream200

LAGAN
Tight1
5.6
92.5
39.6/26.02
28.3

Total
39.5
98.2
86.1/85.92
85.2

Shuffle-LAGAN

Total
36.49
93.18
82.46
77.71

Tight
5.24
88.31
28.78
26.56

Total
38.17
95.46
84.35
80.48

Tight
5.45
90.34
29.16
27.1

1 These results use chromosome 20 rather than the whole genome as the human sequence, as numbers for the latter were not published. Whole genome
numbers tend to be slightly lower. 2 BLASTZ authors did not publish the overall coverage of UTRs, the 5 UTR/3 UTR coverage is given instead.

Table 3. Classification of local rearrangements found between the human and mouse genomes by Shuffle-LAGAN and LAGAN.

Main Strand

ConsistentTranslocated only
Duplicated only
Transloc. & Dupl.
Total

Inverted

Total

% of
LengthAvg
LengthScore
per bp% ofLengthAvg
LengthScore
per bp% ofLengthAvg
LengthScore
per bp

90.8
1.0
2.3
1.5
95.4

143
150
274
183
145

59.9
53.6
44.8
44.8
59.3

1.1
1.5
0.2
1.9
4.7

133
159
249
267
185

55.9
54.0
39.4
39.2
47.6

91.8
2.5
2.3
3.5
263Mbp

143
155
272
222
133

59.9
53.8
44.2
41.7
58.1

% of Length is the percentage of the total length (in the human sequence) of all local alignments (263 Mbp). Avg Length is the length of the average localalignment. Score per bp is the CHAOS score divided by length

distant (e.g. human and zebrafish). SLAGAN source codeis available at http://lagan.stanford.edu/glocal. The sourcecode for LAGAN is also available from the authors.
The introduction of the concept of glocal alignment suggests several new research directions, of which multiple
glocal alignment is perhaps the most natural. The scoring
scheme that we use in the SLAGAN algorithm is simplistic, and it is likely that a more sophisticated scheme would
result in better alignments. One of the major drawbacksof the SLAGAN algorithm is that it is not symmetric in
the sequence order, and thus misses duplications in the sequence that is constrained to be 1-monotonic. We suspectthat chaining with rearrangements may be an NP-complete
problem if one does not require the chain to be monotonic
in one of the two sequences. It should also be possible tocreate glocal alignment algorithms using completely different edit models. This work represents a ‘first stab’ at anopen problem, and that more work in this area is necessary.

ACKNOWLEDGEMENTS
M.B. was supported by the NSF Graduate Fellowship. The
authors would like to thank Arend Sidow, Gregory Cooper,and Eugene Davydov for useful conversations, and Eric
Green for providing the CFTR sequence data that we used
to train our parameters.
i60

REFERENCES
Aparicio,S. et al. (2002) Whole genome shotgun assembly andanalysis of the genome of Fugu rubripes. Science, 297, 1301–
1310.
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.
Bafna,V. and Pevzner,P. (1998) Sorting by transpositions. SIAM J.
Discrete Math., 11, 224–240.
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.
Bray,N., Dubchak,I. and Pachter,L. (2003) AVID: a global alignment program. Genome Res., 13, 97–102.
Brudno,M., Do,C.B., Cooper,G.M., Kim,M.F. and Davydov,E. (2003) NISC comparative Sequencing Program,
Green,E.D., Sidow,A. and Batzoglou,S. LAGAN and MultiLAGAN: efficient tools for large-scale multiple alignmentof genomic DNA. Genome Res., 13, 721–731.
Brudno,M. and Morgenstern,B. (2002) Fast and sensitive alignmentof large genomic sequences. Proceeding of the IEEE Computer
Society Bioinformatics Conference (CSB).
Chiaramonte,F., Yap,V.B. and Miller,W. (2002) Scoring pairwise

Glocal alignment: finding rearrangements during alignment

genomic sequence alignments. In Proceedings of the PacificSymposium on Biocomputing (PSB).
Christie,D.A. (1996) Sorting permutations by block interchanges.
Inf. Process. Lett., 60, 165–169.
Couronne,O., Poliakov,A., Bray,N., Ishkhanov,T., Ryaboy,D., Rubin,E., Pachter,L. and Dubchak,I. (2003) Strategies and tools forwhole-genome alignments. Genome Res., 13, 73–80.
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.
Flint,J., Tufarelli,C., Peden,J., Clark,K., Daniels,R.J., Hardison,R.,
Miller,W., Philipsen,S., Tan-Un,K.C., McMorrow,T. et al. (2001)
Comparative genome analysis delimits a chromosomal domain
and identifies key regulatory elements in the alpha globin cluster.
Hum. Mol. Genet., 10, 371–382.
Gö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
cell leukemia gene (SCL)—comparative analysis of five vertebrate SCL loci. Genome Res., 12, 749–759.
Hannenhalli,S. (1996) Polynomial Algorithm for Computing
Translocation Distance between Genomes. DAMATH: Discrete
Applied Mathematics and Combinatorial Operations Researchand Computer Science.
Hannenhalli,S. and Pevzner,P.A. (1995) Transforming cabbage into
turnip (polynomial algorithm for sorting signed permutations
by reversals). Proceedings of the 27th Annual Symposium on
the Theory of Computing (STOC 95). Las Vegas, Nevada, pp.
178–189.
Henikoff,S. and Henikoff,J.G. (1992) Amino acid substitution
matrices from protein blocks. Proc. Natl Acad. Sci. USA, 89,
10 915–10 919.
Kent,W.J. (2002) BLAT–the BLAST-Like Alignment Tool. GenomeRes., 12, 656–664.
Lander,E.S. et al. (2001) Initial sequencing and analysis of thehuman genome. Nature, 409, 860–921.
Morgenstern,B. (1999) DIALIGN 2: improvement of the segment–
to–segment approach to multiple sequence alignment. Bioinformatics, 15, 211–218.
Needleman,S.B. and Wunsch,C.D. (1970) A general method applicable to the search for similarities in the amino acid sequence oftwo proteins. J. Mol. Biol., 48, 443–453.
Pevzner,P. and Tesler,G. (2003) Genome rearrangements in mammalian evolution: lessons from human and mouse genomes.
Genome Res., 13, 37–45.
Schwartz,S., Kent,W.J., Smit,A., Zhang,Z., Baertsch,R., Hardison,R.C., Haussler,D. and Miller,W. (2003) Human-mouse alignments with BLASTZ. Genome Res., 13, 103–107.
Simon,A., Stone,E.A. and Sidow,A. (2002) Inference of functional

regions in proteins by quantification of evolutionary constraints.
PNAS, 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.
Varré,J.S., Delahaye,J.P. and Rivals,E (1999) Transformation distances: a family of dissimilarity measures based on movementsof segments. Bioinformatics, 15, 194–202.
Venter,J.C. et al. (2001) The sequence of the human genome.
Science, 291, 1304–1351.
Waterson,G.A. et al. (2002) Initial sequencing and comparative
analysis of the mouse genome. Nature, 420, 520–562.

APPENDIX 1: DEFAULT PARAMETERS
The initial local alignments are generated using CHAOS
using seeds of size 11 with one degeneracy. The chains arethen rescored with matches and mismatches scored using
the match matrices from Chiaramonte et al. (2002), a gap
open penalty of –50 and a gap continuation penalty of –25.
All local alignments scoring above 2000 are returned, andare used to create the 1-monotonic conservation map.
When creating the 1-monotonic conservation map,
we penalize all inversions with a cost of—1000— 1 ×
distance, and translocations with—2000—2.5 × distance,
where the distance is the offset between the two alignments, as defined in section 2.2.2. In order to prevent
false chaining of two segments that are consistent we
use a small gap continuation penalty of −0.5, with nogap opening penalties, as they are not meaningful when
one is dealing with gaps between alignments, rather than
individual base pairs.
Finally for reasons of efficency we have arbitrarily
limited the expansion of consistent subsegments to only
25 kilobase pairs at each end. A whole genome alignment
pipeline often tries to align non-homologous regions, andbecause LAGAN runs slowly in such cases it is necessary
to limit the input sequence lengths. LAGAN is run withits default parameters. We do not currently perform thefinal merging of the overlapping alignments. In orderto build the 1-monotonic map we use repeat masked(Repeatmasker, Smit and Green, unpublished) sequences,
but we use the unmasked sequences for the final LAGAN
alignment, allowing us to align conserved repeats without
being misled by non-homologous repeats when building
the conservation map.
All of the parameters were hand trained on the Cystic
Fibrosis Transmembrane Conductance Regulator region(CFTR), representing ∼1.8 megabases of human chromosome 7 and its ortholog in mouse (Thomas et al.
manuscript in preparation, 2003).
The whole genome alignments were run on a 24 nodefarm of two processor 2.3 GHz Pentium 4 machines.
For evaluation we used human genome June 2002 freezei61

M.Brudno et al.

Total:
Gap open 500, Gap extend 25, threshold 2500
A
C
G
T
A
90
−100
−50
−100
C
−100
110
−100
−50
G
−50
−100
110
−100
T
−100
−50
−100
90

Fig. 3. Matrices and cutoffs.

(http://genome.ucsc.edu/) and mouse genome MGSC
v3 freeze (http://www.ncbi.nlm.nih.gov/genome/guide/mouse/index.html). The matrices and cutoffs used tocompute coverage are given in Figure 3.

i62

Tight:
Gap open 2000, Gap extend 50, threshold 3400
A
C
G
T
A
100
−200
−100
−200
C
−200
100
−200
−100
G
−100
−200
100
−200
T
−200
−100
−200
100