New Flexible Approaches for Multiple Sequence Alignment
Tetsuo SHIBUYA and Hiroshi IMAI
Department of Information Science, Faculty of Science, University of Tokyo
7-3-1 Hongo, Bunkyo-ku, Tokyo 113, Japan
Tel: +81-3-3812-2111 ext.4102 Fax: +81-3-5800-6933

Keywords: multiple sequence alignment, algorithm, suboptimal, enumeration, parametricAbstract
The multiple sequence alignment problem is applicable and important in various elds in molecular biology such as the prediction of three dimensional structures of proteins and the inferenceof phylogenetic tree. However, the optimal alignment based on the scoring criterion is not always
the biologically most signi cant alignment. We here propose two exible and ecient approaches
to solve this problem.
One approach is to provide many suboptimal alignments as alternatives for the optimal one.
It has been considered almost impossible to investigate such suboptimal alignments of morethan two sequences because of the enormous size of the problem. We propose techniques forenumeration of suboptimal alignments using Eppstein algorithm. We also discuss what kindof suboptimal alignment is unnecessary to enumerate, and propose an ecient enumeration
algorithm to enumerate only necessary alignments.
The other approach is parametric analysis. The obtained optimal solution with xed parameters such as gap penalties is not always the biologically best alignment. Thus, it is required
to vary parameters and check how the optimal alignments change. The way to vary parameters has been studied well on the problem of two sequences, but not on the multiple alignment
problem because of the diculty of computing the optimal solution. We propose techniques forthis parametric multiple alignment problem, and examine the features of obtained alignments
by various parametric analyses.
For both approaches, this paper performs experiments on various groups of actual protein
sequences, and examines the eciency of these algorithms and property of sequence groups.

1 Introduction
The multiple alignment is a problem to obtain the alignment of multiple sequences with thehighest score based on some given scoring criterion between characters. This problem appears
in various elds of molecular biology such as the prediction of three dimensional structures ofproteins and the inference of phylogenetic tree.
The method using dynamic programming (DP) is well-known for the alignment problems(Dayho et al. 1978, Gotoh 1995, Spouge 1989, Waterman 1995). This method requires O(nd )time and space for d sequences of length at most n. This method can be applied when n is
not so large and d is 2 or 3, but for larger problems, it is impractical. The A3 algorithm isa well-known algorithm for the general optimization and search problems. This algorithm canreduce the search space dramatically if a powerful estimator is used. Thus the A3 algorithmwith upper bounding operation is proposed recently for computing the optimal alignment ofmultiple sequences (Ikeda & Imai 1994, Ikeda 1995).
The optimal alignment based on the scoring criterion is not always the biologically mostsigni cant alignment. We propose two methods to cope with this problem. One is to enumerate
suboptimal solutions and the other is to do parametric analysis.
A suboptimal alignment is an alignment whose score is close to the optimal one. In fact, incase of aligning two sequences, the suboptimal alignments problem is well-studied (Chao 1994,
Naor et al. 1993, Shibayama et al. 1993, Zuker 1991) and used for many applications suchas predicting protein structure and so on (Saqi et al. 1991, Saqi et al. 1992, Waterman etal. 1987). In multiple alignment problem, we can see suboptimal alignments of each pair ofsequences with these methods for only two sequences as in (Zuker 1991), but these are not theaccurate suboptimal alignments of all the sequences.
Enumeration of the suboptimal alignments had not been considered as very practical evenin the case of aligning two sequences (Naor et al. 1993, Zuker 1991). But such enumeration
becomes easier because a new ecient algorithm for the k shortest paths problem is proposedby Eppstein (Eppstein 1994, Shibuya et al. 1995). This algorithm enumerates the lengths of the
k shortest paths in O(k + n + m) time and space if we are given the shortest path tree from thesource or to the destination for any graph with non-negative m edges and n vertices. Even if
we have to output the paths themselves, this algorithm requires only time linear to the output
size, add to the time given above. Note that the shortest path tree can be constructed with DP
or A3 algorithm.
For this approach, we rst discuss the method to obtain E1 , which can be done with someextension of the A3 algorithm. E1 represents all aligned groups of residues in optimal andsuboptimal alignments which are at most 1 worse than the optimal. Furthermore, based on
this extended A3 algorithm and the Eppstein algorithm, we go on to discuss the methods forthe enumeration problem.
Number of suboptimal solutions is very large, and we should restrict the outputs to important solutions if we know what kind of solutions is important. Hence we discuss what kind of
suboptimal alignment is unnecessary to enumerate, and propose an ecient technique to enumerate only necessary alignments. This technique is so splendid that it remains output size
sensitive even though we restricted solutions only to the necessary alignments in our concept.
1For the other approach of parametric analysis, we rst review the basic techniques for parametric analysis (Gus eld et al. 1992, Huang et al. 1994, Vingron et al. 1994, Waterman 1994,
Waterman et al. 1992, Zimmer 1997). and propose new techniques for multiple alignments. As
for the techniques, we use Eppstein algorithm to examine all the optimal solutions for one xed
parameter, and upper bounding technique for the parametric alignment. In most of previous
works, they computed only one optimal solution for one xed parameter in parametric analysis.
We enumerate all the optimal solutions because the parametric analysis is the analysis of optimal
solutions and we consider we should examine all optimal solutions. Fortunately, it is reasonable
to obtain with Eppstein algorithm.
Then we do a parametric study on gap penalties, score tables and weight matrices. Relatedwith the weight matrices, we show the (enhanced) A3 algorithm is applicable for the weighted
multiple alignment problem. Weighted problem is considered only when the phylogenetic tree isgiven, but our approach enables more exible study of the weighted multiple alignment problem.
This technique may also be useful in constructing or tuning phylogenetic tree.
For both approaches, this paper performs experiments on various groups of actual protein
sequences, and examines the eciency and practicality of these algorithms and property ofsequence groups we use.

2

2 Preliminaries
In this section, we rst explain the de nition of the multiple alignment problem, and then survey
the exact algorithms for it brie y.
2.1 Multiple Sequence Alignment Problem

In this section, we describe what the multiple sequence alignment problem is. These are givenin this problem: A

xed set of alphabets 6 including a gap which is denoted by -.

 A score function between characters s : 6 2 6 ! R.

Each member in 6 except for a gap is called a character, and a nite string of characters iscalled a sequence. On the other hand, a nite string of members in 6 is called a padded sequence.
Then the set of ith elements of two or more padded sequences is called a ith column of the setof the sequences. Furthermore, a set of consecutive columns is called a region. A set of padded
sequences (S10 ; S20 ; . . . ; Sd0 ) is called a d-alignment or simply an alignment of (S1; S2 ; . . . ; Sd ) if
and only if all of the padded sequences have same length l, ith column contains at least one
character for any i(1  i  l), and (S10 ; S20 ; . . . ; Sd0 ) becomes (S1 ; S2; . . . ; Sd) if all gaps aredeleted. Furthermore, the (pairwise) projection Aij of an alignment A = (S10 ; S20 ; . . . ; Sd0 ) isde ned as the alignment obtained from (Si0; Sj0 ) by removing columns without any characters.
The score of a 2-alignment is de ned as the summation of scores of all columns obtained
directly with the score function. The score of a d-alignment is de ned as the summation of thescores of all the pairwise projections of the d-alignment. Finally, we can de ne the problem:the alignment problem is a problem of nding the alignment of given sequences with the largest
score.
In this de nition of the problem, the gap penalty is linear to sum of the length of gap
sequences. This kind of a gap penalty is called the linear gap penalty. If there is some starting
gap penalty which is costed to the starting of the gap sequence, it is called the ane gap penalty.
In this thesis, we mainly deal with the linear gap penalty.
The weighted sum-of-pairs multiple alignment problem (Altsuchul 1989, Gotoh 1995) is a
generalization of the simple sum-of-pairs multiple alignment problem described above. Thisversion of the problem is often used when the phylogenetic tree is given. In this problem, we
optimize sum of weighted score of each pairwise sequences alignment: we multiply the score of
the alignment of the ith and the j th sequence by wij . We call (wij ) a weight matrix. Note that(wii ) is always 0 for any i. We deal with this kind of problem only in the section 4.3.3.
The number of characters in 6 is 20 in the case of aligning protein sequences, and it is 4
in the case of aligning DNA sequences. The score function is usually given by a score table.
Table 1 shows the most famous and widely-used score table called PAM-250 (Altsuchul 1991,
Dayho et al. 1978). It is one of the variations of PAM matrices proposed by Dayho (Dayho
et al. 1978).
The multiple alignment problem can be easily transformed to the shortest path problem onsome grid-like directed acyclic graph with no negative edges. Let Sk be the kth sequence of d3

sequences to be aligned, and nk = O(n) be the length of Sk . Then suppose a directed acyclic
graph G = (V; E ) such that V = f(x1; . . . ; xd)jxi = 0; 1; . . . ; ni g and E = f(v; v + e)jv 2 V; e 2[0; 1]d; e 6= 0g. In this graph, a path from s = (0; . . . ; 0) to t = (n1 ; . . . ; nd ) corresponds to analignment of the sequences.
In the alignment problem of two sequences, the length of an edge is de ned from the score
table between characters, and the length of a path from s to t equals the score of the corresponding alignment. Figure 1 shows an example of it. In the graph, each diagonal edge corresponds tothe match of characters, and each horizontal or vertical edges corresponds to the inserted gaps.
In the multiple alignment problem, the sum of all the scores for alignments of pairwise sequences
is generally used as the score. This score of the alignment equals the length of the corresponding
path, if we de ne the length of each edge as the sum of the lengths of the corresponding edges
in the pairwise projections of the alignment. In this way, the longest path problem from s to t
in this graph is equivalent to the original alignment problem. This longest path problem can beeasily transformed to the shortest path problem by reversing the signs of the lengths (Gupta etal. 1995, Ikeda & Imai 1994, Ikeda 1995). Thus we can use shortest path algorithms for graphs
like above to compute the optimal alignment of given sequences.
2.2 Exact Algorithms for Alignment Problem

2.2.1 Dynamic Programming

The alignment graph presented in the section 2.1 is a layered graph and we can easily use
dynamic programming (DP) to obtain the shortest path (Dayho et al. 1978, Gotoh 1992,
Gotoh 1995, Spouge 1989, Waterman 1995) in time linear to the size of the graph size.

Algorithm 1 (Dynamic Programming) Let l(u; v) be the length of the edge (u; v) in thealignment graph, s be the source (0; 0; . . . ; 0) and t be the destination (n1 ; n2 ; . . . ; nd ).
1. Let p(s) be 0.

X n do the following:1
X x = i, compute the following value p(v), andFor all v = (x0 ; x1 ; . . . ; x ) such that

2. For i = 1 to i =

j d

j

d

k

1kd
let previous(v) be v 0 e which satis es this equation (1). Note that this p(v) equals to theshortest path length from s to v.

p(v) = e2[0min(p(v 0 e) + l(v 0 e; v));1]d ;e6=0(1)

3. We can obtain the shortest path by tracing back previous(v) from the destination t.

But the DP requires O(nd ) memory space where n is the length of the longest sequence andd is the number of sequences, because it needs to store all vertices of the graph in memory.
Thus the DP requires too much memory in the case of large d, and the DP can only deal withalignment of 2 or 3 sequences in ordinary case of aligning actual protein sequences.
42.2.2 Dijkstra Method
Dijkstra method is the most famous algorithm for the shortest path problem on graphs without
negative edges. The outline of this algorithm is as follows:

Algorithm 2 (Dijkstra Method) Let the graph in assumption be G = (V; E ), l(v; w) be thelength of the edge (v; w) which is always non-negative, s be the source and t be the destination.

1. Let p(v) be the potential of a vertex v, which is initialized as +1 except for the source s
whose potential is zero, and S be an empty set.
2. Add vertex v0 which has the minimum potential in V 0 S . Then, stop if v0 is t.
3. For all edges (v0 ; v ) in E, if p(v0 ) + l(v0 ; v) is smaller than p(v), replace p(v) with p(v0 ) +l(v0; v) and replace the path kept in v with the shortest path from s to v0 added with theedge (v0 ; v).
4. Go to step 2.

This algorithm requires O(m + n log n) time where m is the number of the edges and n is thenumber of the vertices in the graph, which is worse than DP, and this algorithm itself is not so
useful for alignment problem. But it will be remarkably more ecient if it is extended to theA3 algorithm as in the next section 2.2.3.

2.2.3 A3 algorithm
The A3 algorithm will not search the whole graph in nding the shortest path if a good estimate
for the shortest path length from each vertex to the destination t can be used.
The basic algorithm for the A3 algorithm is like following:

Algorithm 3 (A3 Algorithm) Let the graph in assumption be G = (V; E ), l(v; w) be thelength of the edge (v; w) which is always non-negative, s and t be the source and the destination,and h(v) be heuristic estimate for the length of the shortest path to t which is not longer than
the actual one.
1. Let p(v) be the potential of a vertex v, which is initialized as +1 except for the source s
whose potential is zero, and S be an empty set.
2. Add vertex v0 , whose value of p(v ) + h(v ) is smallest in V 0 S , to S . Then, stop if v0 is t.
3. For all edges (v0 ; v ) in E, if p(v0 ) + l(v0 ; v) is smaller than p(v), replace p(v) with p(v0 ) +l(v0; v) and replace the path kept in v with the shortest path from s to v0 added with theedge (v0 ; v). If v is in S , remove it from S .
4. Go to step 2.

In this algorithm, p(v) for vertex v in S is also the length of the shortest path from s to v.
p(v) + h(v) is the estimate for the shortest path from s to t via v. The searched vertices by theA3 algorithm is always within searched vertices by the Dijkstra method. In this way, the A35

algorithm can get the shortest path more e ectively. If the estimate h(v ) equals to the actual
shortest path length from v to t, this algorithm searches only on the shortest path.
The estimate h(v ) must not be longer than the shortest path length from v to t, because
the nal obtained path must not be longer than the other paths. If some of the estimate h(v ) is
longer than the actual shortest path, this algorithm is called A algorithm, which is often used
as an approximate algorithm for the shortest path problem.
In the A3 algorithm, the shortest path from s may not appear rst, and a shorter path may
be found in the future search, which is the reason of the removal of vertices from S in step 3.
It makes this algorithm rather inecient. This can be avoided if the estimator is dual feasible.
The de nition of `dual feasible' is as follows:

De nition 1 The estimator h for the shortest path to t is called dual feasible if and only if hsatis es the following constraint:8(u; v ) 2 El(u; v) + h(v)  h(u)(2)

If the estimator is dual feasible, the A3 algorithm can be easily translated to the Dijkstra
method by modifying the length of the edges (Hsu 1994, Ikeda et al. 1994, Ikeda 1995):

Theorem 1 (Ikeda Hsu et al.) Let h be a dual feasible estimator for s. The Dijkstra method

on a graph in which the length of the edge (u; v), or l(u; v) is replaced by l0 (u; v) as follows isequivalent to the A3 algorithm on the original graph.

l0(u; v) = l(u; v) + h(v) 0 h(u)(3)
Proof: l0(u; v) is non-negative because of dual feasibility of h. Thus, the shortest path can besearched with the Dijkstra method in the modi ed graph. Let p be the shortest path from s tov. Then the potential of v used in searching with the Dijkstra method on the modi ed graph is
described as follows:

X l (u; v)( )X l(u; v) + h(v) h(s)=p(v) =

0u;v 2p0(u;v)2pThis means the Dijkstra method on the new graph is equivalent to the A3 algorithm on theoriginal graph, because h(s) is constant.
2
Ikeda and Imai (Ikeda & Imai 1994) show the following estimator is very useful for alignment
problem in case d > 2. Let Gij be the corresponding graph to the alignment of Si and Sj , vij be
the corresponding vertex in Gij to v in G, and L3(u; v ) be the shortest path length from u to v.
Then h(v) = 1i<j d L3 (uij ; vij ) can be used as a powerful estimator for the multiple alignment
problem. This estimator is easily be shown to be dual feasible, i.e. l(u; v) + h(v)  h(u). Hence
the A3 algorithm can be applied as following.

P

Algorithm 4 (A3 algorithm for the Alignment Problem)

P

1. For each of i and j (1  i < j  d), apply DP to graph Gij from tij to calculate L3 (vij ; tij )for each vij in Vij . Then let h(v) be 1i<jd L3 (uij ; vij ).

6

2. Modify the length of edge (u; v) in G as follows, using h(v) above, and compute the shortest
path with Dijkstra method.

l0(u; v) = l(u; v) + h(v) 0 h(u)(4)

Note that the time and space used for the DP in the step 1 is negligible, if d is large. This A3
algorithm can deal with alignment problem of 5 to 6 normal sequences in reasonable time.
A vertex in the graph for the multiple alignment has 2d 0 1 edges going out from it, and theA3 algorithm examines all the descendant vertices and keeps in a heap the information about
all of them. If an upper bound L+ (s; t) for the s-t shortest path, which corresponds to the lower
bound of the score of the alignment, is given, the necessary space for the heap can be reduced(Ikeda 1995): we can ignore w such that L3 (s; v) + l(v; w) > L+ (s; t), when we examine thedescendant vertices of v. If the necessary space for the heap is reduced, the computing time of
the A3 algorithm will be also reduced. This is called the enhanced A3 algorithm.
Note that the branch-and-bound techniques implemented in MSA program (Gupta et al. 1995)
is equivalent to this enhanced A3 algorithm, and Araki et al. (Araki et al. 1993) showed thatthe estimated score computed directly by the score matrix is useful for the 2-alignment problem.

7

3 Enumeration of Suboptimal Alignments of Multiple Sequences
In this section, we deal with the rst approach to the exible alignment, that is enumerationof suboptimal solutions. At rst, we introduce the Eppstein algorithm (Eppstein 1994) which
is very ecient for enumeration of suboptimal paths in ordinary graphs. We then consider
how to obtain E1 eciently for multiple alignment problem. Furthermore, we also discuss how
to enumerate the suboptimal solutions, introducing the new notation of classes of suboptimal
solutions. Based on these algorithms, we also do experiments on actual protein sequences.
3.1 Eppstein Algorithm

Eppstein (Eppstein 1994) proposed an algorithm which nds implicitly the k shortest paths forthe graph G with non-negative m edges and n vertices regardless of cycles, in O(m + n + k) timeafter the shortest path tree is constructed. Eppstein (Eppstein 1994) also proposed an easieralgorithm of O(m + n log n + k) time. Note that before the proposition of Eppstein, the bestalgorithm for such a problem was O(k(m + n log n)).
In the algorithm, we use  (u; v) for the edge (u; v) as in equation (6). This  (u; v) denotes
how much longer the path will be using the edge (u; v ) than the optimal path by way of v, andtherefore this value is always non-negative.
If an edge (u; v) is on the shortest path tree, (u; v ) is zero, otherwise, it is called a sidetrackand (u; v) may not be zero. If we go along an s-t path p other than the shortest path, there
must be one or more sidetracks on p, and we de ne sidetrack(p) as the nearest sidetrack froms within them.
Let (tail(p); head(p)) be sidetrack(p). Then we can suppose a heap, in which the parent ofa path p is a path which is same as p from head(p) to t, but go along the shortest path froms to head(p) instead of using sidetrack(p). We de ne this parent of p as parent(p) and we callp a child of parent(p). The root of the heap is the shortest path, and all the paths from s tot appear in the heap once. In this heap, p is (sidetrack(p)) longer than parent(p). Figure 2
shows the path heap of the alignment graph of two sequences KN and QK.
We call a heap i-heap if the node of the heap has only i children at most (it is not required
to be balanced). The basic concept of the Eppstein algorithm is to modify this path heap to
4-heap, sharing as many nodes as possible. Figure 3 shows an example of this compact versionof path heap of the same alignment graph in Figure 2, in which some of the nodes are shared
and the number of nodes in it is reduced.
From this heap, we can obtain the k shortest paths in O(k) time (Frederickson 1993), orO(k log k) time in sorted form. The following is the outline of this algorithm:

Algorithm 5 (Eppstein)
1. Construct the shortest path tree from s to all the other vertices.
2. For each vertex v, construct HG (v ), that is, a 3-heap of sidetracks (u0 ; u), such that u ison the shortest path from s to v , ordered by (u0 ; u). Let the length from the root of HG (v)to a node n be (u; v) if n represents sidetrack (u; v ).

8

(a) For each vertex v, construct Hout (v), that is a 2-heap in which the root has only one
child, of sidetracks (v0 ; v) ordered by  (v0 ; v).
(b) For each vertex v, construct HT (v), that is, a 2-heap of vertices on the shortest pathfrom s to v ordered by the value  of the root of the heap made in step 2-(a).
(c) Merge Hout (v) and HT (v ) to make HG (v). Then let the length of the edge from node
n1 to node n2 in this heap be (n2 ) 0 (n1).
3. For each v in G, make an edge from each node in HG (v ) which represents a sidetrack(u0; u) to the root of HG (u0), and de ne the length of this new edge as the value of the root.
4. Make a new node for each v in G, and make an edge from this node to the root of HG (v).
Let the length of this edge be  of the root. Let this new graph be P (G).

Then we can nd a heap Hv (G) in P (G) for any v , considering the root as the node madein step 4 for v, and the value of a node as the length from the root to the node. There is a
one-to-one correspondence between the nodes in Hv (G) and the paths from v to t in G, and thek smallest nodes in this virtual heap Hv (G) correspond to the k shortest path. Moreover, we
can easily restore the path from the node of the heap, which can be done in O(n0 ) time wheren0 is the size of the output alignment.
In this algorithm, we can also compute the k shortest paths using the shortest path tree to t,which is the original method described in the Eppstein's paper (Eppstein 1994). Note that the
shortest path tree in step 1 is constructed generally by the Dijkstra method, but for problems
such as the alignment problem, we can also use DP. Eppstein showed the step 2 can be done inO(n + m) time, but this is a complicated algorithm and takes much time in practice, and weuse a far more easier algorithm that can be done in O(n log n + m) time, which is also proposedby Eppstein (Eppstein 1994): we make HG (v ) one by one from s to the other vertices along theshortest path tree, sharing as many nodes as possible.
3.2 Upper Bounding Technique for Computing E1

E1 is a set of vertices which are used by the s-t paths whose lengths are at most 1 longer than

the shortest path, and it corresponds to all aligned groups of residues in optimal and suboptimalalignments in original problem. The problem to compute it is well-studied (Ikeda 1995, Naoret al. 1993, Shibayama et al. 1993, Zuker 1991). Here we show how to compute this set E1
with A3 algorithm (Ikeda 1995). For any path p from s to t, the modi ed path length by theexpression (4) is only h(t) 0 h(s) longer than the original length. This value is not relevant top, thus E1 on the modi ed graph is same as the original one.
Theorem 2 Any paths from s to t on a graph in which the length of the edge (u; v), or l(u; v)
is replaced by l0 (u; v ) as in (4) are a constant shorter than those on the original graph.
Proof: Let p be a path from s to t, and h be a dual feasible estimator for the shortest path
length to t. Then the length of a path p or length0(p) in new graph is described by the lengthof p or length(p) in the original graph as follows:l0(u; v)length0(p) =X(u;v)2p9=

X l(u; v) + h(t) h(s)(u;v)2p0= length(p) + h(t) 0 h(s)(5)

According to this, all the paths on the new graph from s to t are h(s) 0 h(t), which is constant,shorter than those on the original graph.
2
Note that the following corollaries can be easily derived from this Theorem 2.

Corollary 1 E1 related to the paths from s to t on a graph in which the length of the edge(u; v), or l(u; v ) is replaced by l0(u; v) as in (4) are same as that on the original graph.

Corollary 2 The k shortest paths from s to t on a graph in which the length of the edge (u; v),or l(u; v ) is replaced by l0 (u; v) as in (4) are same as those on the original graph.

Hence, rst we modify the edge lengths with some dual feasible estimator, and then we canobtain E1 with the Dijkstra method as follows (Ikeda 1995) on this modi ed graph.

Algorithm 6 (E1)
1. Search from s by the Dijkstra method until the shortest path from s to t is discovered.
2. Search successively until a vertex v , to which the shortest path from s is more than 1
longer than the s-t shortest path, is discovered.
3. Modify the length of each edge (u; v ) to (u; v) as follows:(u; v) = l(u; v) + L3 (s; u) 0 L3 (s; v)(6)

4. Apply the Dijkstra method from t until a vertex from which the shortest path to t is longer
than 1 is discovered in this modi ed graph. E1 is the set of vertices searched in this step.

A vertex in the graph for the multiple alignment has 2d 0 1 edges going out from it, andthe Dijkstra algorithm examines all the descendant vertices and keeps the information about allof them. If an upper bound L+ (s; t) for the s-t shortest path is given, we can also reduce thenecessary space for heap as in the case of computing the optimal solution with enhanced A3
algorithm: we can ignore w such that L3(s; v) + l(v; w) > L+ (s; t) + 1, when we examine thedescendant vertices of v.
In general, such kind of an upper bound is dicult to obtain. However, we can use theactual shortest path length obtained in step 1 for the upper bound in step 2: we can ignore w
such that L3 (s; v) + l(v; w) > L3(s; t) + 1. Note that if we are given some upper bound of thesolution before computing the optimal, we can of course use it too.
3.3 Extending Eppstein Algorithm to Reduce Memory SpaceIn this section, we discuss how to enumerate eciently all the suboptimal alignments whosescores are at most 1 lower than the optimal one. The original Eppstein algorithm requires
searching all over the graph, and requires much memory. But it is evident that we only have to10

apply the Eppstein algorithm in the subset E1 after computing E1 as in the previous section,and then search the Eppstein's heap structure with the depth rst method.
However, if we use the easier O(n log n + m) algorithm in step 2 (b) of Eppstein algorithm(algorithm 5) in section 3.1, we do not have to compute E1 additionally.
First, we must take the step 1 and 2 in the section 3.2, using upper bounding technique. These
procedures cannot be skipped. After these procedures, we implement the Eppstein algorithm asfollows:

Algorithm 7 (Eppstein algorithm with A3)

1. Construct the Eppstein's heap structure only on the shortest path.
2. Search for suboptimal solutions which are at most 1 worse than the optimal (root) fromthe root of Hs (G) with depth rst search method. If we encounter HG (v ) which has notbeen constructed yet, we construct the heap structures of vertices on the shortest path froms to v for which we have not not constructed heaps yet.

With this method, when we nish enumerating all the suboptimal alignments, the set ofvertices for which we constructed the Eppstein heap is also E1 .

Theorem 3 The sets S of the vertices for which the Eppstein's heap structures are computedin the algorithm 7 is E1 .

Proof: We construct HG(v) for all the vertices in E1 in the algorithm 7. Thus we can easily

see that E1 2 S . A newly encountered vertex v for which HG (v) has not constructed is in E1 ,because we only search suboptimal paths which are 1 longer than the shortest path at most.
Add to that, to obtain HG (v) for some vertex v, we only have to compute HG (u) for each vertex
u on the shortest path from the source s to v in the easier method of the step 2(b) in algorithm 5,
which is also trivially in E1 . Thus, we do not construct Eppstein's heap structures for vertices
outside of E1 . Hence, we conclude that S = E1 .
2
Thus we do not have to compute E1 additionally. Notice that this technique can be also
used in general graphs other than the graphs for alignments.
3.4 Classi cation of Suboptimal Alignments

In this section, we classify suboptimal alignments. We introduce a notion of alignment class Di
as follows:

De nition 2 Di is a class of alignments which have i regions di erent from the optimal alignment, which is in D0 .
Figure 5 shows some examples of suboptimal multiple alignments of protein fragments. (a)is the optimal alignment, and (b), (c) and (d) are suboptimal alignments. The regions boundedby boxes are the regions which are di erent from the optimal alignment. (b) and (c) have only
one such region. On the other hand, (d) has two, both of which appear also in (b) or (c).
According to our notion of classi cation, the optimal alignment (a) is in the class D0 (andnone of the others are in this class), suboptimal alignments (b) and (c) are in the class D1, andsuboptimal alignment (d) is in the class D2 .
11Considering the fact that we can easily reconstruct (d) from (b) and (c), (d) is not soimportant as (b) nor (c) and is sometimes unnecessary to enumerate. Thus, it will be a goodnews if we can eciently enumerate only the alignments in the classes D0 and D1 .
3.5 Avoiding Unnecessary Alignments

The paths to which the alignments in the class Di corresponds branch o i times from theshortest path to which the alignment in the class D0 corresponds. Hence, we can consider avery easy branch-and-bound technique to avoid such alignments in Di (i  2) in enumeration of
suboptimal alignments: when we search the Eppstein's heap structure, if head(p) of s-t path p
is on the s-t shortest path and parent(p) is not the shortest path, ignore p and its all conceptual
children (de ned in section 3.1). Figure 6 shows the concept of this technique.
Recall that one of the good feature of the Eppstein algorithm is that we can obtain thesolutions in the time linear to the output size after having constructed the Eppstein heap. But
the above technique is not output sensitive, because it requires checking some of the nodes which
are in the class D2. If there are many such solutions, the computing time becomes rather large.
Thus we modi ed the Eppstein algorithm to overcome this problem. Here we show how toconstruct a heap structure whose nodes represent only alignments in classes D0 (the root or the
optimal solution) and D1:

Algorithm 8 (Modi ed Eppstein Algorithm)

1. Construct the shortest path tree from s to all the other vertices.
2. For each vertex v, construct HG (v), that is, a 3-heap of sidetracks (u0 ; u) ordered by  (u0 ; u)
as follows. Let the length from the root of HG (v) to a node n be (u; v) if n represents
sidetrack (u; v ).
(a) For each vertex v, construct Hout (v), that is a 2-heap in which the root has only one
child, of sidetracks (v0 ; v) ordered by  (v0 ; v).
(b) For each vertex v that is not on the s-t shortest path, construct HT (v), that is, a
2-heap of vertices on the s-v shortest path but not on the s-t shortest path ordered by
the value  of the root of the heap made in step 2-(a).
(c) For each vertex v on the s-t shortest path, construct HT (v ), that is, a 2-heap ofvertices on the s-v shortest path ordered by the value  of the root of the heap madein step 2-(a).
(d) Merge Hout (v) and HT (v ) to make HG (v). Then let the length of the edge from node
n1 to node n2 in this heap be (n2 ) 0 (n1).
3. For each v on the s-t shortest path, make a node Nv .
4. For each v in G, make an edge from each node in HG (v ) which represents a sidetrack(u0; u) to the root of HG (u0) if u0 is not on the s-t shortest path, and de ne the length of
this new edge as the value of the root. Otherwise, make an edge from the node to Nu .
0

5. Make a new node for each v in G, and make an edge from this node to the root of HG (v).
Let the length of this edge be  of the root. Let this new graph be P (G).

12

Once the heap was constructed, what we should do next is same as in the original Eppsteinalgorithm: we only have to search from the root of the virtual heap. Accordingly, we caneciently enumerate only the alignments in class D1 , and apparently it is output sensitivealgorithm once the heap is given. We can easily see that this algorithm can be also used with
A3 algorithm using the techniques in section 3.3.
We can also extend this algorithm for enumeration of alignments in class Di (i  c) for any
c in the same way, though the algorithm becomes much more complicated:

Algorithm 9
1. Construct the shortest path tree from s to all the other vertices.
2. For each vertex v, construct heaps H (i) (v) (0  i < c for the vertices on the s-t shortestGpath, and 0 < i  c for the others), a 3-heap of sidetracks (u0 ; u) as follows. Let the lengthfrom the root of HG(i) (v) to a node n be (u; v ) if n represents sidetrack (u; v).
(a) For each vertex v , construct H (i) (v) (0  i < c for the vertices on the s-t shortestout

path, and 0 < i  c for the others), that is a 2-heap in which the root has only one
child, of sidetracks (v0 ; v) ordered by  (v0 ; v).
(b) For each vertex v that is not on the s-t shortest path, construct c 0 2 heaps HT(i) (v)(0 < i < c), all of which are 2-heaps of vertices on the s-v shortest path ordered by
the value  of the root of the heap made in step 2-(a).
(c) For each vertex v that is not on the s-t shortest path, construct heap HT(c) (v), that is,
2-heap of vertices on the s-v shortest path but not on the s-t shortest path, orderedby the value  of the root of the heap made in step 2-(a).
(d) For each vertex v on the s-t shortest path, construct HT(i) (v) (0  i), that is, a 2-heapof vertices on the s-v shortest path ordered by the value  of the root of the heap madein step 2-(a).
(i) (v ) and H (i) (v ) to make H (i) (v). Then let the length of(e) For each i and v, merge HoutTGthe edge from node n1 to node n2 in these heaps be  (n2 ) 0 (n1 ).

3. For each v on the s-t shortest path, make a node Nv .

4. For each v in G, make an edge from each node in HG(i) (v) which represents a sidetrack(u0; u) to the root of HG(i) (u0) if u0 is not on the s-t shortest path, and de ne the length of
this new edge as the value of the root. If u0 is on the s-t shortest path, make an edge tothe root of HG(i+1) (u0 ) if i + 1 < c, or to Nu if i + 1 = c.
05. Make a new node for each v on the s-t shortest path, and make an edge from this node
to the root of HG(0) (v). Let the length of this edge be  of the root. Let this new graph be
P (G).

This algorithm requires O(c(n + m)) time to construct the heap structure. In contrast, thebranch-and-bound technique we mentioned at rst in this section can deal with this kind ofproblem easily: we only have to remember how many branches from the s-t shortest path the13parent path has, in searching suboptimal paths in the Eppstein heap structure. Thus, when cis large, which technique is more ecient may depends on cases. But, note that the problem
whose c is large is not so important as that whose c is small, i.e. 1 or 2.
3.6 Extracting Knowledge from Eppstein Heap

As mentioned by Eppstein (Eppstein 1994), the Eppstein heap has a good feature: some ofnumerical values for each suboptimal solution can be obtained in O(1) time with some simple
pre-process of O(jE1 j) time. In the case of multiple alignment problem, these can be obtained insuch an ecient way for example: the number of aligned groups in which all residues are same,the number of gaps, the score computed with another score table, the length of the alignment,
and so on. Our algorithm in the section 3.5 which enumerates only alignments in the class D1
have the same feature too.
Let us consider the case of computing the number of gaps in each obtained suboptimal
solution. Let the number of gaps p we want to know be gaps(p). Then the algorithm will belike this:Algorithm 10
1. For each vertex v in V , compute the number of gaps contained in the part of the s-v shortest
path. This can be done in O(n) time, by the depth rst search on the shortest path tree tot. Then let this value be g(v).
2. For each sidetrack (u; v ), compute following0(u; v):

0 (u; v) = g(u) 0 g(v)(7)

3. When we search in the path heap and obtained a path p (which represents some alignment),
we compute the gaps(p) as follows from the value of gaps(parent(p)).

gaps(p) = gaps(parent(p)) + 0(sidetrack(p))(8)

In the same way, we can get many kinds of values for each solution. But note that there issome values that cannot be obtained in O(1) time. For example, ane gap cost is one of them,though it is very important value.
3.7 Experimental Results

In this section, we examine the eciency of our approach and investigate the properties ofsuboptimal alignments through experiments. In the experiment, we used the PAM-250 matrix,
and linear gap penalty bx where x is the gap length and b is the minimum value in the PAM250 matrix, 08. All the experiments are done on Sun Ultra 1 workstation with 128 megabytememory.

14

3.7.1 Case with High Similarity
We rst did experiments on a group of 8 sequences with high similarity in Table 2. According
to it, the average scores per amino pair of these pairwise alignments are about 2:5 to 4. Add tothis, the optimal score of multiple alignment of all these 8 sequences is 33129 and its length is331298  2:59 (Table 3). These
456, thus the average score per amino pair of this alignment is 456
1(2)are higher than in the experiment in the next section 3.7.2. Figure 7 shows one of the optimalalignments of all the 8 sequences in Table 2. You can easily recognize the high similarity.
As for computing alignments of less than 8 sequences, we could apply the simple A3 algorithm. However, for alignment of the 8 sequences, we used the upper bounding technique(enhanced A3 ) because 128 megabyte memory is not enough for computing with the simple A3
algorithm: we used in the experiment the optimal solution as the upper bound to see the besteciency of this enhanced algorithm. In any case, we used the upper bounding technique after
the optimal solution is obtained.
According to Table 3, the DP takes a lot of time compared with the A3 algorithm when d issmall, but it is negligible when d is large. This table also shows that, the additional searchingtime required for computing suboptimal solutions is not so much as long as 1 is not much larger
than in these experiments: it requires at most twice the time in total as in the case of computing
only the optimal alignment in these experiments if 1  40.
Figure 8, Table 4 and Table 5 show the results of enumerating the suboptimal alignments.
Figure 8(a) shows that there are enormous number of suboptimal alignments, and the number
increases exponentially as 1 increases. However, in Figure 8(b), we can see the number ofsuboptimal solutions is dramatically reduced by ignoring alignments in class Di (i  2). The
number of the alignments enumerated in this way is only 0.0003% (d = 4) to 0.4% (d = 8) of
all the alignments in case 1 = 30 (see Table 5): it seems dicult to check signi cance of all
the suboptimal alignments at most 10 worse than the optimal, but in our method, we can do it.
Accordingly, the enumeration time is also reduced drastically(see Table 5).
Figure 9 shows an example of suboptimal alignments. The score of it is 33139, which is only
10 worse than the optimal. In this gure, * in the rst line indicates the regions di erent from theoptimal alignment in Figure 7. This alignment has 5 such regions, which means this alignment isin the class D5 . The existence of such alignments makes enumeration of suboptimal alignments
more dicult. It is the reason that our enumeration approach which avoids enumerating suchalignments is very ecient.
Observing Figure 8(a), the number of the suboptimal alignments seems to be similar andirrelevant to d. It is an interesting fact, but this comparison is unfair. The number must becompared between the cases which have same value of (1d) : we must consider 1 per amino pair.
2
For example, it is all right to compare the case 1 = 28 (d = 8) and the case 1 = 10 (d = 5).
In this case, the number of suboptimal alignments in the former case is 7168718
16112  444:9 times asthat of the latter case.
According to Table 4, jE1 j and the size of Eppstein heap for this size of 1 is not so large.
Thus, the enumeration time in Table 5(b) is small, though it includes times for constructing theheap. In Figure 8(b), the number of suboptimal alignments in class D1 increases much whend = 2; 7; 8 compared with other cases. The reason of this is seen in Table 4. The jE1 j in cases
15d = 2; 7; 8 is larger than the others: there may be many alignments which have a large region

di erent from the optimal. On the other hand, in Figure 8(a), the number of the alignments
in case d = 4 is large compared with others, but much of these must be combinations of small
number of `necessary' alignments.

3.7.2 Case with Low Similarity
We next did experiments on 5 globin sequences as in Table 6. According to Table 6, the averagescores per amino acid pair of pairwise alignments of them are about 0:2 to 1:3. The score of the
optimal multiple alignment of these 5 sequences is 543 and its length is 165, thus the averagescore per amid acid pair of this alignment is 1655431(5)  0:33, which is lower than the previous case.
2
Figure 10 shows the optimal alignment of the sequences in Table 6. It shows the dissimilarity
compared with the case of the previous experiments.
Figure 11 and Table 7 show the result of the experiments. According to Table 7, the searchingtime by simple A3 algorithm is far longer than in the previous experiments for same d, though
the length is short. It is because the estimator of the A3 algorithm is not so powerful in casewith low similarity.
According to Figure 11, the number of alignments in this experiment is also drasticallyreduced as in the experiments in the section 3.7.1 by ignoring alignments in class Di (i  2):
the number of alignments in D1 is only 0:004% of that of all the suboptimal alignments in cased = 5 and 1 = 20.
As mentioned by Zuker (Zuker 1991), the alignments with low score are not always insignificant. In general, if the lengths of sequences to be aligned are short, the size of E1 will besmall. However, the size of E1 is larger than in the previous experiments for same d and 1.
Hence we can conclude that sequences we use in this experiment is not so signi cant as in theprevious experiment. In this way, we can use the size of E1 as a factor of the signi cance of thealignment.

16

4 Parametric Analysis of Multiple Sequence Alignment
In this section, we describe the techniques for parametric analysis of multiple sequence alignment
problem. Furthermore, we also do experiments on actual protein sequences.
4.1 Basic Techniques

In this section, we describe basic methods to check how the optimal solution varies as theparameters such as gap penalties change. The easiest approach for this kind of problem is
to change the parameter little by little and check the optimal solution, but we cannot know
how little we should change the parameter. Recently the techniques for parametric analysis
are developed (Gus eld et al. 1992, Huang et al. 1994, Vingron et al. 1994, Waterman 1994,
Waterman 1995,waterman et al. 1992, Zimmer 1997). In those previous works, they also did
parametric analysis which deal with more than one parameters, but algorithms for them are notso ecient as the one-parameter case and it will often be nonsense if the parameters are notrelated each other. Thus we deal with only one parameter at one time in this thesis.
We consider the case in which the score of some alignment Ai is expressed with parameterp as follows:si(p) = a(Ai) + b(Ai) 1 p(9)
Gap penalty satis es this expression for example.
From here, we explain how to divide 1-parameter (1-dimensional) space to regions in which
the optimal alignments are always same. Let ai be a(Ai ) and bi be b(Ai ). Let pi and pj be
the values of the parameter which satis es pi < pj and has di erent optimal solutions. Let thealignment Ai be the alignment with smallest value of b among the optimal alignments at p = pi
and Aj be the alignment with largest value of b among the optimal alignments at p = pj . Then
0 aj
this two alignments Ai and Aj has the same score at p = pij = 0 abi 0
bj . If the optimal scorei
at p = pij equals to si (pij ) = sj (pij ), there are only two regions between pi and pj . Otherwise,
we can apply the same technique recursively (i.e. apply between pi and pij and between pij andpj ) to obtain such division. Figure 12 shows an example of this procedure.
Letting n be the number of regions which we want to obtain, we only have to compute theoptimal solutions 2n 0 1 times. Thus we can eciently do parametric analysis in the case of one
parameter.
In the algorithm, the alignments with the largest or smallest value of b among the optimalalignments at some xed parameter are required. These can be easily obtained by some lexicographical extension of DP (Vingron et al. 1994, Waterman 1994, Waterman 1995, Waterman etal. 1992). This technique can be applied also to the (enhanced) A3 algorithm, but the aim of
the parametric analysis is to examine all the optimal solutions. Accordingly, it is not preferable
to ignore most optimal solutions if there are many. Thus, for this purpose, we enumerate theoptimal solutions by the Eppstein algorithm in the section 3.1 and pick up the solutions withthe largest or smallest value of b.

17

4.2 Upper Bounding Technique for Parametric Alignment

As we stated in the section 2.2, the A3 algorithm will be more ecient if some upper bounding
value for the optimal solution is given (it is called the enhanced A3 algorithm). In the parametricalignment problem, si (pij ) = sj (pij ) in the section 4.1 is a lower bound of the score of the optimalalignment at p = pij . Thus it can be used as this upper bounding value in computing the optimalalignments at p = pij with A3 algorithm. Note that the enhanced A3 algorithm will show best
performance if the si (pij ) = sj (pij ) is the optimal score at p = pij , which always happens at thenal stage of the parametric analysis.
4.3 Experimental Results

In this section, we do parametric analysis on groups of actual protein sequences. We use thefamous PAM-250 matrix for score matrix in the experiments in sections 4.3.1 and 4.3.3. Here,all the experiments are also done on Sun Ultra 1 workstation with 128 megabyte memory insections 4.3.1 and 4.3.3 and on Sun Ultra 1 workstation with 256 megabyte memory in section
4.3.2.

4.3.1 Parametric Gap Penalty
Previous works on parametric analysis mainly dealt with gap penalty, because it is a very
important factor of the alignment problem (Gus eld et al. 1992, Huang et al. 1994, Vingron
et al. 1994, Waterman 1994, Waterman et al. 1992). But these works dealt only with the2-alignment problem.
We did parametric analysis of gap penalty using the top d sequences (d < 7) in Table 2. In
general, the most popular gap penalty is the minimum value in the score matrix, which is 08
in this PAM-250 case. We did parametric analysis for d-sequence alignment (2  d  6) with
gap penalty between 02 and 016.
Table 8 shows the result of the experiment. In Table 8, the rst row of each entry of d shows
the boundaries of the regions, but several of the ends are not the boundaries: the ends with
0 in #Max and #Min entry are not boundaries. The second row shows the numbers of the
optimal solutions at the value. The last two rows show the numbers of optimal solutions withlargest/smallest value of b in the section 4.1. Thus, these values equal to the number of the
optimal solutions between the boundaries.
According to the table, it is observed that the intervals of the regions become smaller (i.e. theoptimal solution is not stable), as the penalty increases regardless of d. It also shows that there
are not so much di erence between di erent d's, which means we can do parametric analysis
as easily as in the 2-alignment case. The table also shows that there are more than 1 optimal
solution in all cases in the experiments.
Figure 13 shows the number of visited nodes by A3 algorithm in computing the all the
optimal alignments under various gap penalties. According to this gure, the number of thevisited nodes increases drastically as the gap penalty increases especially when gap penalty islarger than 04. This is the reason why we analyzed gap penalty only up to 02:5 or 03:5 whend  5: the required space was too large to compute when the gap penalty is around 02.
18

We think this increase of the search space is caused by the instability of the optimal solution,
based on many various experiments. It will be more clear in the next section 4.3.2. In general,
it is known that the number of required space for A3 algorithm is large if the similarity among
the group is low, probably for same origin.
We also did experiments on 10 TNF- sequences in Table 9, which are very similar to eachother. Table 11 shows the result of it. In the experiments, we examined between gap penalties
016 and 02:5, for d = 2; 4; 6; 8; 10 sequences. It is between 016 and 04:5 in case of d = 10,
because of computational diculty of obtaining the optimal solution with gap penalties near to
0. There are fewer regions than the EF-1 case. The reason of this may be the high similarity
among this group.

4.3.2 Parametric Score Matrix
In this section, we deal with the parametric analysis of score matrix. Score matrix is very
important parameter and it deeply a ects the quality of the output alignment. Thus parametricanalysis of it is not only very important but also requires much care to deal with.
There are (n + 2)(2 n 0 1) parameters, where n is the number of characters, to change in thescore matrix, thus what we can do is very limited simple analysis. Some of the previous works
dealt with this topic. For example, Vingron et al. (Vingron et al. 1994) analyzed the varying
solutions as they added some constant (which is the parameter to change) to each element in
the score matrix. But it is not so interesting because those constants does not have any meaningat all. We should analyze in more practical and interesting way.
We implemented a program to analyze how the optimal solutions will change as the score
(2)
(1)
(2)
matrix changes linearly between two score matrices (s(1)
ij ) and (sij ): we used (p 1 sij +(1 0 p) 1 sij )as the score matrix and considered p as the parameter to be changed. Note that our program
can of course deal with the analysis what Vingron et al. did (Vingron et al. 1994).
At rst, we did experiments on 8 sequences (fragments) of rhodopsin superfamily shown intable 12. Figure 14 shows the optimal alignment of 5 sequences taken from table 12. But, in the
biologists' view, it is not the best alignment (Blanck et al. 1987, Hargrave 1991). According tothem, the K's (lycine) which are marked with * at the rst line in gure 14 must be aligned.
There must be various methods to cope with this kind of problems, but one of the simplest
method is to check the optimal alignment obtained with a score table whose score for K-K
is large. But, too large score for K-K makes irrelevant K's aligned and the alignment will be
unnatural. Thus, we should choose as small score as possible in the constraint that we canobtain a biologically good solution. In this way, we can align reasonably those K's.
In this point of view, we did parametric analysis on these sequences, letting the score forK-K be the parameter to change. In the experiment, we let the scores other than it are same as
those in the PAM-250 score matrix in Table 1. In the PAM-250 score matrix, the score for K-K
is 5, hence we only examined in the region where its score is larger than 5.
Table 14 shows the result of this experiment. The rst row of each entry shows the score
for K-K which is boundary of regions except for several ends with - in #Max and #Min entries.
The second row shows the number of the optimal solutions and the other two rows shows thenumber of the optimal solutions with largest/smallest b in the section 4.1. In the table, * denotes19the smallest score for K-K which induces the optimal alignment whose K's are aligned desirably.
Figure 15 shows the alignment of 5 sequences obtained with smallest score for K-K such that K's
are aligned desirably. This alignment is biologically much better than that obtained with theordinary PAM-250 score matrix.
Figure 16 shows how the number of visited nodes by the A3 algorithm change as the score
for K-K changes. This gure shows the case of 4 and 5 sequences. In this gure, we can easily
see that the number is large where the boundary in the parametric analysis exists. Thus, as we
mentioned in the section 4.3.1, we conclude that computing the optimal solution will be moredicult if the solution is instable.
It is also interesting that in the case of the alignment of 4 rhodopsin sequences, the numberof visited nodes is xed when the score for K-K is larger than 84. It is because the all possible
K's are aligned and thus letting the score larger makes no di erence.
PAM-t score table is computed based on the evolutionary permutation probability of amidacids per period proportional to t (Dayho et al. 1978), and we can consider various PAM score
tables such as PAM-120, PAM-250 and so on. Altschul (Altsuchul 1991) suggested that it isbetter to align sequences with two or three di erent PAM scores based on statistical analysis.
Thus parametric analysis on this t is useful. Because PAM-t converges into some matrix as t
becomes larger, parametric analysis on t can be approximated by the linear parametric analysis
between PAM matrices if t is large enough. But note that, because of the same reason, theoptimal solutions obtained with di erent large t's are often same. Note that the exact parametricanalysis on this t is dicult and remains as one of future works.
Table 15 shows the result of this kind of experiment using 6 TNF- sequences in Table 9.
(1)
The rst row of each entry shows the value p of (p 1 s(1)ij + (1 0 p) 1 sij ) which is boundary ofregions except for several ends with - in #Max and #Min entries, where (s(1)) is PAM-250 scorematrix and (s(2)) is PAM-320 matrix. The second row shows the number of the optimal solutions
and the other two rows shows the number of the optimal solutions with largest/smallest b in
the section 4.1. Notice that there are no other optimal solutions other than those at both ends
between the PAM-250 score matrix and the PAM-320 score matrix.

4.3.3 Parametric Weight Matrix
In this section, we deal with weighted version of the multiple alignment problem as we mentioned
in the section 2.1. Computing the optimal solution of this problem by the (enhanced) A3
algorithm is rather easy: all we have to do is using h(v) = 1i<j d wij 1 L3(uij ; vij ) as theestimator. Thus we can use same techniques as in previous experiments.
A weight matrix for aligning sequences whose phylogenetic tree is known can be made if divergence between sequences are given (Altsuchul 1989). But what should we do if the divergence
are ambiguous? In such case, parametric analysis between reasonable two weight matrices helps.
Thus, parametric analysis of weight matrix may helps tuning parameters of a phylogenetic tree.
There are (d 0 2)(2 d + 1) parameters to change in the weight matrix, thus what we can do isvery limited simple analysis. We implemented a program to analyze how the optimal solutions
change as weight matrix changes linearly between two weight matrices, like in the case of thesection 4.3.2.

P

20

We did experiments between following two weight matrices of (wij ) and (wij(16;n) ):(0i=ji 6= jwij = 1w(p;n) =
ij

(p w
1

wij

i = n or j = nijotherwise

(10)
(11)

In this equation, (wij ) corresponds to the simple sum-of-pairs multiple alignment, and wij(p;n) increases the importance of nth sequence to p times as the simple sum-of-pairs multiple alignment.
If biologically good alignment is discovered in the experiment, we can estimate the importanceof the sequence which was increased.
Table 16 shows experiment results using the top 6 EF-1 sequences in Table 2. The rst
column is the name of the sequence whose importance was increased. The rst row of eachentry shows the value of p of wij(p;n) which is boundary of regions except for several ends with- in #Max and #Min entries. The second row shows the number of the optimal solutions andthe other two rows shows the number of the optimal solutions with largest/smallest b in thesection 4.1.
In this experiment, we notice that the optimal solutions will change even when only p = 1:33
in some of the cases (cases of Tha and Ent). It means we should take more care of the weightmatrix. This experiment also show that there are more than 1 optimal solution in all the casesin this experiment. In the experiment, the number of the regions are not too large to deal with
(6 to 10 in this experiment). This means this approach is very reasonable to take.
We also did experiments on 5 EF-2 sequences in Table 17, which are very similar to each otherand very long. The experiments are done also between weight matrices of (wij ) and (wij(16;n) ). In
the experiments, we increased importance of one of the sequences. Table 18 shows the results.

21

5 Conclusions
In multiple alignment problem, it is often said that the optimal solution obtained with only thescoring criteria is not always the biologically best one. Thus, in this thesis, we took two exible
approaches to exible multiple alignment problem to overcome this problem. One approach isby suboptimal analysis and the other is by parametric analysis.
We rst investigated method for enumerating suboptimal alignments of multiple sequences.
Previous works related to this subject could not deal with enumeration in practice because of
the enormous amount of the suboptimal alignment.
At rst, we introduced the Eppstein algorithm (Eppstein 1994). We then proposed ecient
upper bounding technique to obtain E1 eciently for multiple alignment problem using A3
algorithm. Using this technique, we extended Eppstein algorithm to reduce the memory space.
Furthermore, we classi ed the suboptimal alignments into classes Di based on the numberof the di erent regions from the optimal solution. We showed that the suboptimal alignments
in classes Di (i  2) can be easily constructed from the alignments in the class D0 (the optimal
solution) and the class D1, and suggested that the alignments in the classes Di (i  2) is notnecessary to enumerate. In this point of view, we proposed an algorithm to enumerate only the
suboptimal alignments in the class D1 . We also discussed how to obtain information from the
Eppstein's heap structure. We showed that some of the important values of the alignments canbe obtained in constant time after some preprocess.
Based on these algorithms and observations, we did experiments on actual protein sequences
to see the eciency of our algorithms and to examine the property of the sequences. We did
experiments on both groups with high similarity and with low similarity and we saw our method
was very dramatical.
Next, we investigated parametric optimization of the multiple alignment problem. The
previous works studied this topic only on 2-alignment case, but multiple alignment case is also
very important.
We rst described the basic technique of parametric analysis, which is a technique to divide
the parameter space into regions where the optimal solution is always same. Then we discussed
how to use upper bounding technique of the enhanced A3 algorithm in this parametric analysis.
In multiple alignment problem, there are three kind of parameters, i.e. gap penalty, scorematrix, and weight matrix. We did parametric experiments on each of them, and examined
the property of the multiple alignment problem in the view of parametric optimization. In
the experiments, we showed that the parametric analysis of multiple alignment problem is aspractical as that of 2-alignment.
As for future works, there are several problems to do. At rst, we should develop more exible
method, such as parametric analysis of suboptimal solutions. Using ane gap cost which makes
the problem much dicult also remains as one of future works. Algorithms we proposed are not
only for the alignment problem and can also be applied to many other optimization problems,
which we should do as a future work too.

22

Acknowledgement
This work was supported in part by the Grant-in-Aid for Scienti c Research on Priority Areas,\Genome Science" from the Ministry of Education, Science, Sports and Culture of Japan.

23

References
Altsuchul S. F. 1991. Amid Acid Substitution Matrices from an Information Theoretic
Perspective, J. Mol. Biol., 219, 555-565.
Altsuchul S. F., Carroll R. J., and Lipman D. J. 1989. Weights for Data Related by a Tree, J.
Mol. Biol., 207, 647-653.
Altschul S. F. and Erickson B. W. 1986. Optimal Sequence Alignment Using Ane Gap Costs,
Bull. Math. Biol., 48, 603-616.
Araki S., Goshima M., Mori S., Nakashima H., Tomita S., Akiyama Y., and Kanehisa M. 1993.
Application of Parallelized DP and A3 Algorithm to Multiple Sequence Alignment, Proc.
Genome Informatics Workshop IV, 94-102.
Blanck A. and Oesterhelt D. 1987. The Halo-opsin Gene. II. Sequence, Primary Structure ofHalorhodopsin and Comparison with Bacteriorhodopsin, EMBO J., 6, 265-273.
Chao K. M. 1994. Computing All Suboptimal Alignments in Linear Space, Proc. 5th
Symposium on Combinatorial Pattern Matching, 31-42. LNCS 807, Springer-Verlag, New York.
Dayho M. O., Schwartz R. M., and Orcutt B. C. 1978. Atlas of Protein Sequence andStructure (M. O. Dayho ed.), 5, suppl. 3, 345-352, National Biomedical Research Foundation,
Washington D. C.
Eppstein D. 1994. Finding the k Shortest Paths, Proc. 25th IEEE Annual Symposium on
Foundation of Computer Science, 154-165.
Frederickson G. N. 1993. An Optimal Algorithm for Selection in a Min-Heap, Information andComputation, 104, 197-214.
Gotoh O. 1982. An Improved Algorithm for Matching Biological Sequences, J. Mol. Biol. 162,705-708.
Gotoh O. 1995. A Weighting System and Algorithm for Aligning Many PhylogeneticallyRelated Sequences, Comput. Applic. Biosci., 11, 543-551.
Gotoh O. 1993. Optimal Alignment between Groups of Sequences and its Application toMultiple Sequence Alignment, Comput. Applic. Biosci., 9, 361-370.
Gotoh O. 1995. Parity, 10, No. 12, 13-19, in Japanese.

24Gracy J., Chiche L. and Sallantin J. 1993. Improved alignment of weakly homologous protein
sequences using structural information, Protein Engineering 6, 821-829.
Gupta S. K., Kececioglu J. D., and Scha er A. A. 1995. Improving the Practical Space andTime Eciency of the Shortest-paths Approach to Sum-of-pairs Multiple Sequence Alignment,
J. Comput. Biol., 2, No. 3, 459-472.
Gus eld D., Bakasubramanian K., and Naor D. 1992. Parametric Optimization of Sequence
Alignment, Proc. 3rd ACM-SIAM Annual Symposium on Discrete Algorithms, 432-439.
Hargrave P. A. 1991. Current Opinion on Structural Biology, 1, 575-581.
Hsu M. 1994. A Study on the Shortest-Path Algorithm for Route Navigation, A Master's
Thesis, Department of Information Science, University of Tokyo.
Huang X., Pevzner P. A., and Miller W. 1994. Parametric Recomputing in Alignment Graphs,
Proc. 5th Annual Symposium on Combinatorial Pattern Matching, 87-101, LNCS 807,
Springer-Verlag, New York.
T. Ikeda, 1995. Applications of the A3 Algorithm to Better Routes Finding and Multiple
Sequence Alignment, A Master's Thesis, Department of Information Science, University ofTokyo.
Ikeda T., Hsu M., Imai H., Nishimura S., Shimoura H., Hashimoto T., Temmoku K., andMitoh K. 1994. A Fast Algorithm for Finding Better Routes by AI Search Techniques, Proc.
IEEE International Confrence on Vehicle Navigation and Information System, 90-99.
Ikeda T. and Imai H. 1994. Fast A3 Algorithms for Multiple Sequence Alignment, Proc.
Genome Informatics Workshop V, 90-99.
Imai H. and Ikeda T., 1995. k-group Multiple Alignment Based on A* Search, Proc. Genome
Informatics Workshop VI, 9-18.
Naor D. and Brutlag D., 1993. On suboptimal alignments of biological sequences, Proc. 4th
Symposium on Combinatorial Pattern Matching, 179-196. LNCS 684, Springer-Verlag, NewYork.
Saqi M. A. and Sternberg M. J. 1991. A simple method to generate non-trivial alternatealignments of protein sequences, J. Mol. Biol., 219, 727-732.
Saqi M. A., Bates P. A., and Sternberg M. J. E. 1992. Towards an automatic method ofpredicting protein structure by homology: an evaluation of suboptimal sequence alignments,
Protein Engineering, 5, 305-311.
25

Shibayama G. and Imai H. 1993. Finding K -best Alignment of Multiple Sequences, Proc.
Genome Informatics Workshop IV, 120-129.
Shibuya T., Ikeda T., Imai H., Nishimura S., Shimoura H., and Tenmoku K. 1995. Finding a
Realistic Detour by AI Search Techniques, Proc. 2nd Intelligent Tranportation Systems, 4,2037-2044.
Shibuya T. and Imai H. 1997. Enumerating Suboptimal Alignments of Multiple BiologicalSequences Eciently, Proc. Paci c Symposium on Biocomputing, 409-420.
Shibuya T. and Imai H. 1997. New Flexible Approaches for Multiple Sequence Alignment,
Proc. 1st Annual International Conference on Computational Molecular Biology, 267-276.
Shibuya T. and Imai H. 1996. Parametric Alignment of Multiple Biological Sequences, Proc.
Genome Informatics Workshop 1996, 41-50.
Shibuya T., Imai H., Nishimura S., Shimoura H., and Tenmoku K. 1996. Detour Queries inGeographical Databases for Navigation and Related Algorithm Animations, Proc.
International Symposium on Cooperative Database Systems for Advanced Applications.
Shirai Y. and Tsuji J. 1982. Arti cial Intelligence, Iwanami Course: Information Science, 22,
Iwanami, Tokyo, in Japanese.
Spouge J. L. 1989. Speeding Up Dynamic Programming Algorithms for Finding Optimal
Lattice Paths, SIAM J. Appl. Math., 49, 1552-1566.
Vingron M. and Waterman M. S. 1994. Sequence Alignment and Penalty Choices: Review ofConcepts, Case Studies and Implications, J. Mol. Biol., 235, 1-12.
Waterman M. S. 1995. Introduction to Computational Biology: Maps, Sequences andGenomes, Chapman & Hall.
Waterman M. S. 1994. Parametric and Ensemble Sequence Alignment Algorithms, Bull. Math.
Biol., 56, 743-767.
Waterman M. S. and Eggert M. 1987. A New Algorithm for Best Subsequence Alignmentswith Application to tRNA-rRNA Comparisons, J. Mol. Biol., 197, 723-728.
Waterman M. S., Eggert M., and Lander E. 1992. Parametric Sequence Comparisons, Proc.
Natural Academy of Science, USA, 89, 6090-6093.

26

Zimmer R. and Lengauer T. 1997. Fast and Numerically Stable Parametric Alignment of
Biosequences, Proc. 1st Annual International Conference on Computational Molecular Biology,344-353.
Zuker M. 1991. Suboptimal sequence alignment in molecular biology. Alignment with error
analysis, J. Mol. Biol., 221, 403-420.

27

28

C
S
T
P
A
G
N
D
E
Q
H
R
K
M
I
L
V
F
Y
W

Table 1: PAM-250 score matrix

12
0 2
02 1 303 1 0 6
02 1 1 1 2
03 1 0 01 1 5
04 1 0 01 0 0 2
05 0 0 01 0 1 2 405 0 0 01 0 0 1 3 405 01 01 0 0 01 1 2 2 403 01 01 0 01 02 2 1 1 3 6
04 0 01 0 02 03 0 01 01 1 2 605 0 0 01 01 02 1 0 0 1 0 3 505 02 01 02 01 03 02 03 02 01 02 0 0 6
02 01 0 02 01 03 02 02 02 02 02 02 02 2 506 03 02 03 02 04 03 04 03 02 02 03 03 4 2 602 01 0 01 0 01 02 02 02 02 02 02 02 2 4 2 404 03 03 05 04 05 04 06 05 05 02 04 05 0 1 2 01
0 03 03 05 03 05 02 04 04 04 0 04 04 02 01 01 02
08 02 05 06 06 07 04 07 07 05 03 2 03 04 05 02 06
C S T P A G N D E Q H R K M I L V

97 10
0 0 17F Y W

Table 2: Sequences of EF-TU and EF-1 to be aligned and their scores of pairwise sequencealignments. We use the top d sequences in this table in the experiments.
Sequences
SpeciesProtein
Halobacterium marismortui (Hal) EF-TU
Methanococcus vannielii (Met) EF-TU
Thermoplasma acidophilum(Tha) EF-1
Thermococcus celer(Thc) EF-1
Sulfolobus acidocaldarius (Sul) EF-1
Entamoeba histolytica (Ent) EF-1
Plasmodium falciparum (Pla) EF-1
Stylonychia lemnae(Sty) EF-1

Length
421
428
424
428
435
430
443
446

29

Pairwise Scores
Met Tha Thc Sul Ent Pla Sty
1329 1314 1221 1109 1099 1000 971
1336 1247 1150 1176 1087 1045
1311 1261 1233 1063 1072
1132 1130 1049 991
1192 1131 1099
1584 1551
1636

Table 3: Searching time (sec) by the A3 algorithm in the experiment on the d sequences of
EF-TU and EF-1 . In case d = 8, the enhanced A3 utilizing the optimal score is used. Notethat only DP is used in case d = 2.
best score
Pre-process DP
Search (optimal)
Search (1 = 10)
Search (1 = 20)
Search (1 = 30)
Search (1 = 40)

d=2 d=3 d=4 d=5 d=6 d=7 d=8
1329
0.32
-

3970
1.00
0.18
0.18
0.22
0.27
0.33

7709 12314 18101 24912 33129
4.30 7.23 11.1 16.0 20.5
0.52 3.35 19.6 426 5427
0.60 3.63 20.9 439 5686
0.73 4.17 23.1 462 6735
0.93 5.00 26.9 498 8027
1.23 6.22 31.5 552 9623

Table 4: Size of E1 and Eppstein's heap structure in the experiments on d sequences of EF-TU
and EF-1 . The heap size in the table does not include the number of nodes made in step 4 of
Eppstein algorithm, which is same as jE1 j.
1
10
20
30
40

jE10 jheap size
jE20 jheap size
jE30 jheap size
jE40 jheap size

d=2 d=3 d=4 d=5 d=6 d=7 d=8503 485437 277
1101 595
7184 1010
1447 946
12983 4949
2011 2528
17934 25861

513 553
411 503609 689
1266 1701
817 901
3417 3997
1170 1249
8648 10036

534 579 540
454 674 404
691 799 784
1750 2604 2672
864 1246 1316
3594 7552 8539
1156 1973 2254
7785 18407 23973

Table 5: Enumerating time (sec) when 1 = 30 in the experiment on the d sequences EF-TU
and EF-1 . (a) is the case enumerating all the suboptimal alignments, and (b) is the caseenumerating alignments in class D0 (optimal) and D1. The time of constructing the Eppstein's
heap structure is included in the time below.
(a) #alignments
time (sec)(b) #alignments
time (sec)d=2 d=3d=4

d=5 d=6 d=7

d=838047513 8804702 327522816 85923864 20689104 49633652 13857237

152.87 54.98 1830.18 510.63 124.55 292.50 109.00
6968 1695
1117 2176 1659 41791 60589
0.23 0.13
0.32 0.95 2.10 9.83 29.62
30

Table 6: Globin sequences to be aligned and their scores of pairwise sequence alignments.
globinLumbricus terrestris - AIII
Aplysia limacina
Busycon canaliculatum
Chironomus thummi thummi - VIIA
Chironomus thummi thummi - IIIa

(Lum)
(Apl)
(Bus)
(Ct7)
(Ct3)

Length
157
146
147
145
151

Apl
29

Bus
15
126

Ct7
35
177
111

Ct3
41
140
64
191

Table 7: The best score, searching time (sec) by simple A3 algorithm and the size of E1 in theexperiment on the d globin sequences.
best score
Pre-Process DP
Search (optimal)
Search (1 = 10)
Search (1 = 20)
Search (1 = 30)
Search (1 = 40)
jE10 j
jE20 j
jE30 j
jE40 jd=2 d=3 d=4 d=5
29
0.05
357
776
1149
1544

31

103 354 543
0.27 0.52 0.83
0.73 7.40 837
0.85 8.13 865
0.93 9.43 909
1.22 11.22 964
1.63 14.13 1080
508 380 554
1184 866 1159
2403 1575 2448
3839 2669 4569

Table 8: The result of the experiment on parametric gap penalty using EF-1 sequences.
d=2

d=3

d=4

d=5

d=6

Gap penalty#Solutions
#Max
#Min
Gap penalty#Solutions
#Max
#Min
Gap penalty#Solutions
#Max
#Min
Gap penalty#Solutions
#Max
#Min
Gap penalty#Solutions
#Max
#Min

016
4
016
8
016
16
016
2
016
4
-

05

12
8
4
03:5
16
8
8
08
32
16
16
07:5
4
2
2
06:5
16
4
4

03

24
16
8
03
24
16
8
03:83
32
16
16
04
4
2
2
04:5
8
4
4

02:5

192
8
16
02:75
32
16
16
03:5
32
16
16
03:38
4
2
2
04
8
4
4

02
576
32
8
02:5
72
16
16
02:5
32
16
16
03:17
4
2
2
03:5
4
-

02:2

48
32
16
02:33
48
32
16
03
4
2
2

02
256
96
32
02:25
160
128
32
02:88
12
4
2

02
4608
384
128
02:75
8
4
4

02:5
24
4
4

Table 9: TNF- sequences used in the experiment
Species
Homo sapiens (HSp)
Mus musculus (MM)Sus scrofa (SS)Ovis orientalis aries (OOAp)
Bos primigenius taurus (BPT)
Equus caballus (EC)
Oryctolagus cuniculus (OC)Rattus norvegicus (RN)
Homo sapiens (HSi)Ovis orientalis aries (OOAi)

Proteintumor necrosis factor (TNF- ) precursor
tumor necrosis factor (TNF- ) precursor
tumor necrosis factor (TNF- ) precursor
tumor necrosis factor (TNF- ) precursor
tumor necrosis factor (TNF- ) inhibitor
tumor necrosis factor (TNF- ) precursor
tumor necrosis factor (TNF- ) precursor
tumor necrosis factor (TNF- ) precursor
tumor necrosis factor (TNF- ) inhibitor
tumor necrosis factor (TNF- ) inhibitor

32

Length
233
235
232
234
233
234
234
235
233
233

Table 10: Scores of pairwise alignments of TNF- sequences
HSp
MM
SS
OOAp
BPT
EC
OC
RN
HSi
OOAi

MM
955

SS
990
919

OOAp
924
872
980

BPT
935
872
990
1059

EC
1013
945
982
916
933

OC
949
952
906
840
865
962

RN
943
1124
912
856
854
931
946

HSi
1041
952
974
915
926
1013
961
940

OOAi
932
866
988
1145
1067
906
845
850
922

Table 11: The result of the experiment on parametric gap penalty using TNF- .
d=2

d=4

d=6

d=8

d = 10Gap penalty#Solutions
# Max
# Min
Gap penalty#Solutions
# Max
# Min
Gap penalty#Solutions
# Max
# Min
Gap penalty#Solutions
# Max
# Min
Gap penalty#Solutions
# Max
# Min

-16
6
-16
12
-16
6
-16
1
-16
1
-

-2.5
18
12
6
-6.67
24
12
12
-4.95
12
6
6
-5.57
2
1
1
-5.28
2
1
1

33

-3.75
24
12
12
-3.3
18
12
6
-4.5
2
1
1
-4.5
1
-

-2.88
24
12
12
-2.75
24
12
12
-3.36
3
2
1

-2.67
24
12
12
-2.5
36
24
12
-2.92
4
2
2

-2.5
36
24
12

-2.75
4
2
2

-2.5
2
-

Table 12: Sequences of Rhodopsin Superfamily used in the experiments
SpeciesHalobacterium sp.
Homo Sapiens
Homo Sapiens
Gallus gallus
Homo Sapiens
Homo Sapiens
Drosophila melanogaster
Drosophila melanogasterProteinLengthHalorhodopsin precursor 39
Green-sensitive opsin
48
Red-sensitive opsin
48
Rhodopsin
55
Blue-sensitive opsin
54
Rhodopsin
55
Opsin RH3
53
Opsin RH4
53

(Hal)
(HSg)
(HSr)
(GGd)
(HSb)
(HSd)
(DM3)
(DM4)

Table 13: Scores of pairwise sequences of Rhodopsin Superfamily
HSg HSr GGd HSb
Hal -35 -32 -79 -61
HSg
256 108 100
HSr
107 100
GGd
147
HSb
HSd
DM3

34

HSd DM3 DM4-83 -31 -26
111 2
-10
110 7
-5
275 3018
150 34
35
33
21
277

Table 14: The result of the experiment on parametric score matrix using rhodopsin sequences.
* denotes the smallest score for K-K which induces the optimal alignment whose K's are aligned
desirably.

d=2
d=3
d=4
d=5
d=6
d=7
d=8

Score of K-K 5 53* 128#Solutions 1 6
5
#Max
- 1
#Min
- 5
Score of K-K 5 55.5 128#Solutions 2 7
5
#Max
- 2
#Min
- 5
Score of K-K 5 38 58.5* 128#Solutions 4 36
48
12
#Max
- 32
12
#Min
- 4
32
Score of K-K 5 24 89* 124 128#Solutions 24 32
20
28 16
#Max
- 8
12
16 #Min- 24
8
12 Score of K-K 5 20 109.5* 128#Solutions 24 28
6
2
#Max
- 4
2
#Min
- 24
4
Score of K-K 5 16 33.33 64#Solutions 4 8
6
2
#Max
- 4
2
#Min
- 4
4
Score of K-K 5 16.5 32#Solutions 8 16
8
#Max
- 8
#Min
- 8
-

Table 15: The result of the linear parametric analysis between PAM-250 and PAM-320 using 6
TNF- sequencesp

0 0.5 1
d = 6 #Solutions 6 12 9
#Max
- 6 3
#Min
- 6 6

35

Table 16: The result of the experiment on parametric weight matrix using EF-1 sequences.
Hal

Met

Tha

Thc

Sul

Ent

Weight#Solutions
#Max
#Min
Weight#Solutions
#Max
#Min
Weight#Solutions
#Max
#Min
Weight#Solutions
#Max
#Min
Weight#Solutions
#Max
#Min
Weight#Solutions
#Max
#Min

1
4
1
4
1
4
1
4
1
4
1
4
-

1:33
16
4
4
2:25
8
4
4
2:33
8
4
4
1:8
8
4
4
2
12
4
4
1:33
8
4
4

3:17
12
4
8
4: 2
16
8
8
5: 5
8
4
4
4
12
2
4
3:75
8
4
4
3
8
4
4

3
16
8
4
3
12
8
4
3
8
4
4
3
8
4
4
3
8
4
4
1:5
8
4
4

4:5
8
4
4
4:4
16
8
8
7:67
8
4
4
5
12
4
2
5
12
8
4
3:66
8
4
4

4:75
8
4
4
4:5
16
8
8
8
16
4
4
5:33
8
4
4
6
16
8
8
5
8
4
4

5
48
8
4
12
16
8
8
16
4
6
8
4
4
6:43
16
8
8
6
8
4
4

6:57
16
8
8
14:5
16
8
8

16
8
16
24
16
8

10:33
8
4
4
8
16
8
8
10:33
12
8
4

12
8
4
4
10:25
16
8
8
13
24
16
8

14:44
8
4
4
14
16
8
8
16
16
-

16
12
8
4
16
8
-

Table 17: EF-2 (translation elongation factor eEF-2) sequences to be aligned and their pairwisescoresSequences
Species
Homo sapiens (HS)Cricetinae gen. sp. (CGS)
Drosophila melanogaster (DM)
Dictyostelium discoideum (DD)
Saccharomyces cerevisiae (SC)Protein
eEF-2
eEF-2
eEF-2
eEF-2
eEF-2

Length
858
858
844
830
842

Pairwise Scores
CGS DM
DD
SC
4216 3409 2463 2941
3392 2456 2931
2416 2959
2446

Table 18: The result of the experiment on parametric weight matrix using EF-2.
HS

Weight#Solutions
#Max
#Min

1
16
-

1.33
32
16
16

2
32
8
16

36

3
16
8
8

3.14
16
8
8

6
16
8
8

9
24
16
8

16
24
-

37

s

A

T

G

C

A

C

T

t
Figure 1: The graph for the alignment of two sequences ATGC and ACT. The s-t path in the baldline represents the alignment of ATGC- and A--CTK
Q

NThe optimal alignment

K

K

N

Q
K
Shortest path tree

Figure 2: The path heap of the alignment graph of two sequences KN and QK
38

K

N

Q

The optimal alignment

K

Figure 3: The compact path heap of the alignment graph of two sequences KN and QK. Some of
the nodes are shared, and the number of nodes in it is reduced.

sEt
Figure 4: E1 is a set of vertices which are used by the s-t paths whose lengths are at most 1
longer than the shortest path
39

REAFSQAIWRATFAQVPESRSLFKR==
ADFLV-ALF-EKFPDSANFFADFKGKSKNG-S-LLFGLLFKTYPDTKKHFKHFDLAAVF-TAYPDIQARFPQFAGK-DVAS
GSGVE-ILY-FFLNKFPGNFPMFKKLG(a) The optimal alignment
REA
ADF
KNG
LAA
GSG

FSQ AIWRATFAQVPESRSLFKR==
LV- ALF-EKFPDSANFFADFKGKS-S- LLFGLLFKTYPDTKKHFKHFD-VF TAYPDIQARFPQFAGK-DVAS-VE ILY-FFLNKFPGNFPMFKKLGREAFSQAIWRATFAQVPESRS LF KR==ADFLV-ALF-EKFPDSANFFA DF KGKSKNG-S-LLFGLLFKTYPDTKK HF KHFD
LAAVF-TAYPDIQARFPQFAG -K DVASGSGVE-ILY-FFLNKFPGNFP MF KKLG(b) A suboptimal alignment(c) Another suboptimal alignment

REA FSQ AIWRATFAQVPESRS LF KR==ADF LV- ALF-EKFPDSANFFA DF KGKSKNG -S- LLFGLLFKTYPDTKK HF KHFD
LAA -VF TAYPDIQARFPQFAG -K DVAS
GSG -VE ILY-FFLNKFPGNFP MF KKLG

(d) Unnecessary alignment to check
Figure 5: Examples of suboptimal alignments of multiple protein sequences40

root = the optimal alignmentaConceptual path heap

d

b

c
e

all of these alignments have
two or more regions different
from the optimal alignment

a

b

d

c

e

Figure 6: An example of a conceptual path heap. We can easily construct a heap which does
not contain unnecessary alignments such as c and its all descendants.

41

1
80
Hal
MS-DEQHQNLAIIGHVDHGKSTLVGRLLYETGSVPEHVIEQHKEEAEEKGKGGFEFAYVMDNLAEERERGVTIDIAHQEFMet
MAKTKPILNVAFIGHVDAGKSTTVGRLLLDGGAIDPQLIVRLRKEAEEKGKAGFEFAYVMDGLKEERERGVTIDVAHKKFTha
MASQKPHLNLITIGHVDHGKSTLVGRLLYEHGEIPAHIIEEYRKEAEQKGKATFEFAWVMDRFKEERERGVTIDLAHRKF
Thc
MAKEKPHINIVFIGHVDHGKSTTIGRLLFDTANIPENIIKKFE-EMGEKGK-SFKFAWVMDRLKEERERGITIDVAHTKF
Sul
MS-QKPHLNLIVIGHVDHGKSTLIGRLLMDRGFIDEKTVKEAEEAAKKLGKDSEKYAFLMDRLKEERERGVTINLSFMRF
Ent
MPKEKTHINIVVIGHVDSGKSTTTGHLIYKCGGIDQRTIEKFEKESAEMGKGSFKYAWVLDNLKAERERGITIDISLWKF
Pla
Sty MGKEKTHINLVVIGHVDSGKSTTTGHIIYKLGGIDRRTIEKFEKESAEMGKGSFKYAWVLDKLKAERERGITIDIALWKF
MPKEKNHLNLVVIGHVDSGKSTSTGHLIYKCGGIDKRTIEKFEKEAAEMGKGSFKYAWVLDKLKAERERGITIDIALWNF
81
160
Hal
STDTYDFTIVDCPGHRDFVKNMITGASQADNAVLVVAA-D---D-GV-QP-QTQEHVFLARTLGIGELIVAVNKMD-L-V
Met
PTAKYEVTIVDCPGHRDFIKNMITGASQADAAVLVVNVDDA--KSGI-QP-QTREHVFLIRTLGVRQLAVAVNKMD-T-V
Tha
ETDKYYFTLIDAPGHRDFVKNMITGTSQADAAILVISARDG--E-GV-ME-QTREHAFLARTLGVPQMVVAINKMDATSP
Thc
ETPHRYITIIDAPGHRDFVKNMITGASQADAAVLVVAV-T---D-GV-MP-QTKEHAFLARTLGINNILVAVNKMD-M-V
Sul
ETRKYFFTVIDAPGHRDFVKNMITGASQADAAILVVSAKKGEYEAGMSAEGQTREHIILSKTMGINQVIVAINKMDLADTEnt
ETSKYYFTIIDAPGHRDFIKNMITGTSQADVAILIVAAGTGEFEAGISKNGQTREHILLSYTLGVKQMIVGVNKMD-A-I
Pla
ETPRYFFTVIDAPGHKDFIKNMITGTSQADVALLVVPADVGGFDGAFSKEGQTKEHVLLAFTLGVKQIVVGVNKMD-T-V
Sty ETAKSVFTIIDAPGHRDFIKNMITGTSQADAAILIIASGQGEFEAGISKEGQTREHALLAFTMGVKQMIVAVNKMDDKSV
161
240
Hal
DYGESEYKQVVEEV-KDLLTQVRFDSENAKFIPVSAFEGDNIAEESEHTGWYDGEILLEALNELPAPEPPTDAPLRLPIQMetNFSEADYNELKKMIGDQLLKMIGFNPEQINFVPVASLHGDNVFKKSERNPWYKGPTIAEVIDGFQPPEKPTNLPLRLPIQThaPYSEKRYNEVKADA-EKLLRSIGFK-D-ISFVPISGYKGDNVTKPSPNMPWYKGPTLLQALDAFKVPEKPINKPLRIPVE
Thc
NYDEKKFKAVAEQV-KKLLMMLGYK-N-FPIIPISAWEGDNVVKKSDKMPWYNGPTLIEALDQMPEPPKPTDKPLRIPIQ
Sul
PYDEKRFKEIVDTV-SKFMKSFGFDMNKVKFVPVVAPDGDNVTHKSTKMPWYNGPTLEELLDQLEIPPKPVDKPLRIPIQEnt
QYKQERYEEIKKEI-SAFLKKTGYNPDKIPFVPISGFQGDNMIEPSTNMPWYKGPTLIGALDSVTPPERPVDKPLRLPLQPlaSty KYSEDRYEEIKKEV-KDYLKKVGYQADKVDFIPISGFEGDNLIEKSDKTPWYKGRTLIEALDTMQPPKRPYDKPLRIPLQ
NWDQGRFIEIKKEL-SDYLKKIWLQPRQDPFIPISGWHGDNMLEKSPNMPWFTGSTLIDALDALDQPKRPKDKPLRLPLQ
241
320
Hal
DVYTISGIGTVPVGRVETGILNTGDNVSFQPSD-V----S-GEVKTVEMHHEEVPKAEPGDNVGFNVRGVGKDDIRRGDVMetDVYTITGVGTVPVGRVETGIIKPGDKVVFEPAG-A----I-GEIKTVEMHHEQLPSAEPGDNIGFNVRGVGKKDIKRGDVTha
DVYSITGIGTVPVGRVETGVLKPGDKVIFLPAD-K----Q-GDVKSIEMHHEPLQQAEPGDNIGFNVRGIAKNDIKRGDVThc
DVYSIKGVGTVPVGRVETGVLRVGDVVIFEPASTIFHKPIQGEVKSIEMHHEPMQEALPGDNIGFNVRGVGKNDIKRGDVSulEVYSISGVGVVPVGRIESGVLKVGDKIVFMPVG-K----I-GEVRSIETHHTKIDKAEPGDNIGFNVRGVEKKDVKRGDVEnt
DVYKISGIGTVPVGRVETGILKPGTIVQFAPSG-V----S-SECKSIEMHHTALAQAIPGDNVGFNVRNLTVKDIKRGNVPlaGVYKIGGIGTVPVGRVETGILKAGMVLNFAPSA-V----V-SECKSVEMHKEVLEEARPGDNIGFNVKNVSVKEIKRGYVSty DVYKIGGIGTVPVGRVETGLLKPGMVLTFAPMN-I----T-TECKSVEMHHESLTEAEPGDNVGFTVKNLSVKDLRRGYV
321
400
Hal
CGPADDPPSVA--ET-FQAQIVVMQHPSVITEGYTPVFHAHTAQVACTVESIDKKIDPSSGEVAE-ENPDFIQNGDAAVVMet
LGHTTNPPTVA--TD-FTAQIVVLQHPSVLTDGYTPVFHTHTAQIACTFAEIQKKLNPATGEVLE-ENPDFLKAGDAAIVTha
CGHLDTPPTVV--KA-FTAQIIVLNHPSVIAPGYKPVFHVHTAQVACRIDEIVKTLNPKDGTTLK-EKPDFIKNGDVAIV
Thc
AGHTNNPPTVVRPKDTFKAQIIVLNHPTAITVGYTPVLHAHTLQVAVRFEQLLAKLDPRTGNIVE-ENPQFIKTGDSAIVSul
AGSVQNPPTVA--DE-FTAQVIVIWHPTAVGVGYTPVLHVHTASIACRVSEITSRIDPKTGKEAE-KNPQFIKAGDSAIVEnt
ASDAKNQPAVG-CED-FTAQVIVMNHPGQIRKGYTPVLDCHTSHIACKFEELLSKIDRRTGKSMEGGEPEYIKNGDSALVPla
Sty ASDTKNEPAKG-CSK-FTAQVIILNHPGEIKNGYTPLLDCHTSHISCKFLNIDSKIDKRSGKVVE-ENPKAIKSGDSALVASDSKNDPAKD-TTN-FLAQVIVLNHPGQIQKGYAPVLDCHTAHIACKFDEIESKVDRRSGKVLE-EEPKFIKSGEAALV
456
401
Hal
TVRPQKPLSIEPSSEIPELGSFAIRDMGQTIAAGKV--L---G---V-NE-----R
Met
KLIPTKPMVIESVKEIPQLGRFAIRDMGMTVAAGMA--I---Q---VTAKN----KTha
KVIPDKPLVIEKVSEIPQLGRFAVLDMGQTVAAGQC--I---D---L-EK-----R
Thc
VLRPTKPMVIEPVKEIPQMGRFAIRDMGQTVAAGMV--I---S---I-QKA----ESul
KFKPIKELVAEKFREFPALGRFAMRDMGKTVGVGVI--I---D---VKPRKVE-VK
Ent
KIVPTKPLCVEEFAKFPPLGRFAVRDMKQTVAVGVV--K---A---V-T------P
Pla
Sty SLEPKKPMVVETFTEYPPLGRFAIRDMRQTIAVGIINQLKRKNLGAVTAKAPA-KK
RMVPQKPMCVEAFNQYPPLGRFAVRDMKQTVAVGVIKEVVKKEQKGMVTKAAQKKK

Figure 7: The optimal alignment of the 8 EF-TU and EF-1 sequences42#solutions
d=2
d=3
d=4
d=5

100000000

10000000

d=6
d=7
d=8

1000000

100000

10000

1000

100

10

1

Delta
0

10

20

30

(a) enumerating all the alignments#solutions
d=2
d=3
d=4
d=5
d=6
d=7
d=8

100000000

10000000

1000000

100000

10000

1000

100

10

1

Delta
0

10

20

30

40

(b) enumerating alignments in classes D0 and D1
Figure 8: Number of the suboptimal alignments of d sequences of EF-TU and EF-1 whosescores are at most 1 worse than the optimal. (a) is the case enumerating all the alignments.
(0  1  30) (b) is the case enumerating alignments in classes D0 and D1 . (0  1  40)43

1 ****80
Hal
MSDEQ-HQNLAIIGHVDHGKSTLVGRLLYETGSVPEHVIEQHKEEAEEKGKGGFEFAYVMDNLAEERERGVTIDIAHQEFMet
MAKTKPILNVAFIGHVDAGKSTTVGRLLLDGGAIDPQLIVRLRKEAEEKGKAGFEFAYVMDGLKEERERGVTIDVAHKKFTha
MASQKPHLNLITIGHVDHGKSTLVGRLLYEHGEIPAHIIEEYRKEAEQKGKATFEFAWVMDRFKEERERGVTIDLAHRKF
Thc
MAKEKPHINIVFIGHVDHGKSTTIGRLLFDTANIPENIIKKFE-EMGEKGK-SFKFAWVMDRLKEERERGITIDVAHTKF
Sul
MS-QKPHLNLIVIGHVDHGKSTLIGRLLMDRGFIDEKTVKEAEEAAKKLGKDSEKYAFLMDRLKEERERGVTINLSFMRF
Ent
MPKEKTHINIVVIGHVDSGKSTTTGHLIYKCGGIDQRTIEKFEKESAEMGKGSFKYAWVLDNLKAERERGITIDISLWKF
Pla
Sty MGKEKTHINLVVIGHVDSGKSTTTGHIIYKLGGIDRRTIEKFEKESAEMGKGSFKYAWVLDKLKAERERGITIDIALWKF
MPKEKNHLNLVVIGHVDSGKSTSTGHLIYKCGGIDKRTIEKFEKEAAEMGKGSFKYAWVLDKLKAERERGITIDIALWNF

81
**
**
Hal
STDTYDFTIVDCPGHRDFVKNMITGASQADNAVLVVAAD----D-GV-QP-QTQEHVFLARTLGIGELIVAVNKMD--LVMet
PTAKYEVTIVDCPGHRDFIKNMITGASQADAAVLVVNVDDA--KSGI-QP-QTREHVFLIRTLGVRQLAVAVNKMD--TV
Tha
ETDKYYFTLIDAPGHRDFVKNMITGTSQADAAILVISARDG--E-GV-ME-QTREHAFLARTLGVPQMVVAINKMDATSP
Thc
ETPHRYITIIDAPGHRDFVKNMITGASQADAAVLVVAVT----D-GV-MP-QTKEHAFLARTLGINNILVAVNKMD--MV
Sul
ETRKYFFTVIDAPGHRDFVKNMITGASQADAAILVVSAKKGEYEAGMSAEGQTREHIILSKTMGINQVIVAINKMDLADTEnt
ETSKYYFTIIDAPGHRDFIKNMITGTSQADVAILIVAAGTGEFEAGISKNGQTREHILLSYTLGVKQMIVGVNKMD--AIPlaETPRYFFTVIDAPGHKDFIKNMITGTSQADVALLVVPADVGGFDGAFSKEGQTKEHVLLAFTLGVKQIVVGVNKMD--TV
Sty ETAKSVFTIIDAPGHRDFIKNMITGTSQADAAILIIASGQGEFEAGISKEGQTREHALLAFTMGVKQMIVAVNKMDDKSV

161
240
Hal
DYGESEYKQVVEEV-KDLLTQVRFDSENAKFIPVSAFEGDNIAEESEHTGWYDGEILLEALNELPAPEPPTDAPLRLPIQMetNFSEADYNELKKMIGDQLLKMIGFNPEQINFVPVASLHGDNVFKKSERNPWYKGPTIAEVIDGFQPPEKPTNLPLRLPIQThaPYSEKRYNEVKADA-EKLLRSIGFK-D-ISFVPISGYKGDNVTKPSPNMPWYKGPTLLQALDAFKVPEKPINKPLRIPVE
Thc
NYDEKKFKAVAEQV-KKLLMMLGYK-N-FPIIPISAWEGDNVVKKSDKMPWYNGPTLIEALDQMPEPPKPTDKPLRIPIQ
Sul
PYDEKRFKEIVDTV-SKFMKSFGFDMNKVKFVPVVAPDGDNVTHKSTKMPWYNGPTLEELLDQLEIPPKPVDKPLRIPIQEnt
QYKQERYEEIKKEI-SAFLKKTGYNPDKIPFVPISGFQGDNMIEPSTNMPWYKGPTLIGALDSVTPPERPVDKPLRLPLQPlaSty KYSEDRYEEIKKEV-KDYLKKVGYQADKVDFIPISGFEGDNLIEKSDKTPWYKGRTLIEALDTMQPPKRPYDKPLRIPLQ
NWDQGRFIEIKKEL-SDYLKKIWLQPRQDPFIPISGWHGDNMLEKSPNMPWFTGSTLIDALDALDQPKRPKDKPLRLPLQ

241
320
Hal
DVYTISGIGTVPVGRVETGILNTGDNVSFQPSD-V----S-GEVKTVEMHHEEVPKAEPGDNVGFNVRGVGKDDIRRGDVMetDVYTITGVGTVPVGRVETGIIKPGDKVVFEPAG-A----I-GEIKTVEMHHEQLPSAEPGDNIGFNVRGVGKKDIKRGDVTha
DVYSITGIGTVPVGRVETGVLKPGDKVIFLPAD-K----Q-GDVKSIEMHHEPLQQAEPGDNIGFNVRGIAKNDIKRGDVThc
DVYSIKGVGTVPVGRVETGVLRVGDVVIFEPASTIFHKPIQGEVKSIEMHHEPMQEALPGDNIGFNVRGVGKNDIKRGDVSulEVYSISGVGVVPVGRIESGVLKVGDKIVFMPVG-K----I-GEVRSIETHHTKIDKAEPGDNIGFNVRGVEKKDVKRGDVEnt
DVYKISGIGTVPVGRVETGILKPGTIVQFAPSG-V----S-SECKSIEMHHTALAQAIPGDNVGFNVRNLTVKDIKRGNVPla
Sty GVYKIGGIGTVPVGRVETGILKAGMVLNFAPSA-V----V-SECKSVEMHKEVLEEARPGDNIGFNVKNVSVKEIKRGYVDVYKIGGIGTVPVGRVETGLLKPGMVLTFAPMN-I----T-TECKSVEMHHESLTEAEPGDNVGFTVKNLSVKDLRRGYV

321
**
400
Hal
CGPADDPPSVA--ET-FQAQIVVMQHPSVITEGYTPVFHAHTAQVACTVESIDKKIDPSSGEVAEE-NPDFIQNGDAAVVMet
LGHTTNPPTVA--TD-FTAQIVVLQHPSVLTDGYTPVFHTHTAQIACTFAEIQKKLNPATGEVLEE-NPDFLKAGDAAIVTha
CGHLDTPPTVV--KA-FTAQIIVLNHPSVIAPGYKPVFHVHTAQVACRIDEIVKTLNPKDGTTLKE-KPDFIKNGDVAIV
Thc
AGHTNNPPTVVRPKDTFKAQIIVLNHPTAITVGYTPVLHAHTLQVAVRFEQLLAKLDPRTGNIVEE-NPQFIKTGDSAIVSul
AGSVQNPPTVA--DE-FTAQVIVIWHPTAVGVGYTPVLHVHTASIACRVSEITSRIDPKTGKEAEK-NPQFIKAGDSAIVEnt
ASDAKNQPAVG-CED-FTAQVIVMNHPGQIRKGYTPVLDCHTSHIACKFEELLSKIDRRTGKSMEGGEPEYIKNGDSALVPlaASDTKNEPAKG-CSK-FTAQVIILNHPGEIKNGYTPLLDCHTSHISCKFLNIDSKIDKRSGKVVEE-NPKAIKSGDSALV
Sty ASDSKNDPAKD-TTN-FLAQVIVLNHPGQIQKGYAPVLDCHTAHIACKFDEIESKVDRRSGKVLEE-EPKFIKSGEAALV401*******Hal
TVRPQKPLSIEPSSEIPELGSFAIRDMGQTIAAGKV--L------GV-NE-----R
Met
KLIPTKPMVIESVKEIPQLGRFAIRDMGMTVAAGMA--I------QVTAKN----KTha
KVIPDKPLVIEKVSEIPQLGRFAVLDMGQTVAAGQC--I------DL-EK-----R
Thc
VLRPTKPMVIEPVKEIPQMGRFAIRDMGQTVAAGMV--I------SI-QKA----ESul
KFKPIKELVAEKFREFPALGRFAMRDMGKTVGVGVI--I------DVKPRKVE-VK
Ent
KIVPTKPLCVEEFAKFPPLGRFAVRDMKQTVAVGVV--K------AV-TP======Pla
Sty SLEPKKPMVVETFTEYPPLGRFAIRDMRQTIAVGIINQLKRKNLGAVTAKAPA-KK
RMVPQKPMCVEAFNQYPPLGRFAVRDMKQTVAVGVIKEVVKKEQKGMVTKAAQKKK

Figure 9: An example of suboptimal alignments of the 8 EF-TU and EF-1 sequences

1KKQCGVLEGLKVKSEWGRAYGSGHDREAFSQAIWRATFAQVPESRSLFKRVHGDH-TS--DPA-FIAHAERVLGGLDIAISTL
83
Lum
Apl
S-LSAAEADL-AGKSWAPVFAN-KN--ANGADFLVALFEKFPDSANFFADFKGKSVADIKASPKLRDVSSRIFTRLNEFVNNABus
G-LDGAQKTA-LKESWKVLGADGPTMMKNGSLLFGLLFKTYPDTKKHFKHFDDATFAAMDTTGVGKAHGVAVFSGLGSMICSI
Ct7
Ct3 APLSADQASL-VKSTWAQV----RN--S-EVEILAAVFTAYPDIQARFPQFAGKDVASIKDTGAFATHAGRIVGFVSEIIALI
V-ATPAMPSM-TDAQVAAVKGDWEKIKGSGVEILYFFLNKFPGNFPMFKKL-GNDLAAAKGTAEFKDQADKIIAFLQGVIEKL
84
165
Lum
DQP--A-TLKEELDHLQVQHEGRKIPDNYFDAFKTAILHVVAAQLGERCYSNNEEI-HDAIACDGFARVLPQVLERGIKGHH
Apl
ANA--G-KMSAMLSQFAKEHVGFGVGSAQFENVRSMFPGFVAS-VA----APPAGA-DAAWT-KLFGLII-DAL---KAAGABusDDD--D-CVBGLAKKLSRNHLARGVSAADFKLLEAVFKZFLDE--A----TQRKAT-DAQKD-AD-GALL-TML---IKAHV
Ct7
Ct3 GNESNAPAVQTLVGQLAASHKARGISQAQFNEFRAGLVSYVSS-NV----AWNAAA-ESAWT-AGLDNIF-GLL---FAA-L
GSDM-G-GAKALLNQLGTSHKAMGITKDQFDQFRQALTELLGN-LG---FGGNIGAWNATVD-LMFHVIF-NAL---DGTPV

Figure 10: One of the optimal alignments of the 5 globin sequences44#solutions

d=2
d=3
d=4

100000000

d=510000000
1000000
100000
10000
1000
100
10
1
0

20 Delta

10(a) enumerating all the alignments#solutions

d=2
d=3
d=4
d=5

100000000
10000000
1000000

100000
10000
1000
100
10
1
0

10

20

30

40 Delta(b) enumerating alignments in classes D0 and D1
Figure 11: Number of the suboptimal alignments of d globin sequences whose scores are at
most 1 worse than the optimal alignment. (a) is the case enumerating all the alignments.
(0  1  20) (b) is the case enumerating alignments in classes D0 and D1 . (0  1  40)
45

Score

1
4
3
5

6

7
2

p2

p7

p5 p6 p3 p4 p1

Parameter

Figure 12: An example of division of 1-parameter space. In this case, there are 4 regions between
p1 and p2.

46

#Visited Nodes
d=3
d=4800000d=5
d=6600000400000200000

0
-15.0

-10.0

-5.0

0.0

Gap Penalty

Figure 13: Number of visited nodes by A3 algorithm for various gap penalties.

Figure 14: The optimal alignment of 5 rhodopsin sequences based on PAM-250
1**55
Hal ==GLALVQSVGVTSWAYS-VLDVF-AK-YVF---AF-I-LLR-WV-ANN---ER=
HSg NPGYPFHPLMAALPAFFAKSATIYNPVIYVFMNRQFRNCILQ-LF-GKK----V=
HSr NPGYAFHPLMAALPAYFAKSATIYNPVIYVFMNRQFRNCILQ-LF-GKK----V=
GGd NQGSDFGPIFMTIPAFFAKSSAIYNPVIYIVMNKQFRNCMITTLCCGKNPLGDED
HSb NRNHGLDLRLVTIPSFFSKSACIYNPIIYCFMNKQFQACIMK-MVCGKAMTDESD

47

Figure 15: The optimal alignment of 5 rhodopsin sequences when the score for K-K is between
89 and 124
1*
60
Hal GLALVQSVGVTSWAYSVLDVFAK--YVF-A--FILL-R--W---VAN---NER=======
HSg N----PGYPFHPLMAALPAFFAKSATIYNPVIYVFMNRQFRNC-ILQLF-GKK-V=====
HSr N----PGYAFHPLMAALPAYFAKSATIYNPVIYVFMNRQFRNC-ILQLF-GKK-V=====
GGd N----QGSDFGPIFMTIPAFFAKSSAIYNPVIYIVMNKQFRNCMITTLCCGKNPLGDE-D
HSb N----RNHGLDLRLVTIPSFFSKSACIYNPIIYCFMNKQFQAC-IMKMVCGKA-MTDESD

48

#Visited Nodes
5000

4500

4000

3500

3000

2500

2000

1500

1000

500

0

Score for K-K
0.00

20.00

40.00

60.00

80.00

100.00120.00(a) Case of the alignment of 4 rhodopsin sequences#Visited Nodes
100000

90000

80000

70000

60000

50000

40000

30000

20000

10000

0

Score for K-K
0.00

20.00

40.00

60.00

80.00

100.00120.00(a) Case of the alignment of 5 rhodopsin sequences
Figure 16: Number of visited nodes by A3 algorithm for various scores for K-K. The dotted lines
denote the boundaries of the regions where the optimal solution is same.
49