MATH-4961 / MATH-6961 / CSCI-4964 / CSCI-6964

Alignment: Definitions and algorithms

FROM DOT PLOTS TO ALIGNMENT

Let ${\bf A} = a_1a_2a_3 \dots a_m$ and ${\bf B} = b_1b_2b_3 \dots b_n$ be 2 biomolecular sequences of the same type; DNA, RNA or protein. In a dot plot, ${\bf A}$ is placed from top to bottom in the 0^th column of an m+1 x n+1 matrix. ${\bf B}$ is placed from left to right in the 0^th row of the same matrix. ``Dots'' plotted in row i and column j indicate at least that a<sub>i</sub> = b<sub>j</sub>. If the dot plot uses a ``window'' and ``stringency'' method with w > 1, or a k-word method with k > 1, then not all matches are plotted in general.

Another way to write down the information in a dot plot is to write the two sequences from left to right, one above the other, and to use line segments to join a_i to b_j whenever there is a dot in row i and column j. The dot plot in Figure 5 would yield the image in figure 1. This can be regarded as a (bipartite) graph, with the two sequence elements acting as the vertices. We can call this a ``dot plot graph''.

  

Figure 1: The ``Alignment'' corresponding to the dot plot in Figure 5 (dot plot notes).

The graph corresponding to the dot plot in Figure 6 would be very complicated. The idea of ``alignment'' is to take a subset of this graph, yielding an ``alignment graph'', or what is normally called an ``alignment of two sequences''. There are 2 requirements that must be satisfied:

1.

Each vertex has at most one edge.

2.

No 2 edges cross each other.

The graph in Figure 1 violates both of these conditions. Note that both solid and dotted lines have been used in this Figure. Taking only the solid lines yields a valid alignment. The same is true for the dotted lines. Note that these 2 valid alignments are not the only ones that can be achieved by deleting edges. The whole theory of sequence alignment is to decide what edges to use.

This leads us to the first definition of sequence alignment. This definition will evolve as the course progresses.

DEFINITION: An alignment between sequences ${\bf A}$ and ${\bf B}$ is a set of integer pairs, $\{(i_h,j_h)\}_{h=1}^{k}$ that satisfy:

1.

a_ih = b_jh for $1 \leq h \leq k$, (this condition is temporary and will be dropped later)

2.

i_h < i_h+1 for $1 \leq h < k$ and

3.

j_h < j_h+1 for $1 \leq h < k$.

Residues a_ih and b_ih are said to be ``aligned'' or ``matched'' with each other. Later on, we will drop the requirement that matched pairs must be equal.

We now change the look of the alignment graph by agreeing to write the 2 sequences so that the edges are always vertical lines. This is equivalent to saying that matched residues are always ``directly over one another'', or ``in the same column''. In order to accomplish this, gaps must be left in 1 or both sequences, sometimes at the ends. These gaps are usually written explicitly, using the symbol ``-''. For example, an alignment of
${\bf A} =$ AATAACATTCACGCGTCTGC and
${\bf B} =$ AATATTGACGCATTCTGCAT appears below.

AATAACATTCACGC-GTCTGC--
|||   ||| ||||  |||||  
AAT---ATTGACGCATTCTGCAT

Written like this, we can describe the alignment in evolutionary terms. We might imagine that ${\bf A}$ and ${\bf B}$ are the descendants of a common ancestor, ${\bf U}$. The sequence ${\bf U}$ changed into ${\bf A}$ and ${\bf B}$ through a series of ``primitive evolutionary events''. We define these to be:

1.

Substitution: One residue changes to another.

2.

Insertion: A residue is added to a sequence.

3.

Deletion: A residue is deleted from a sequence.

Sequence elements that do not change are called ``conserved''. In an alignment, matched pairs of residues are inferred to be conserved. This may not be true, since the same change may occur in both.

When a_i appears above b_j and $a_i \neq b_j$, we know that a change has occurred. It is possible that more than 1 change has occurred.

When a_i appears above a gap, we know that either a_i has been inserted into the original sequence and appears in ${\bf A}$, or that a_i was in ${\bf U}$ and was deleted from ${\bf B}$. Similarly when a gap appears above b_j, we know that either b_j has been inserted into the original sequence and appears in ${\bf B}$, or that b_j was in ${\bf U}$ and was deleted from ${\bf A}$.

Thus a gap in an alignment indicates that either an insertion or a deletion has occurred. We cannot deduce when a residue is either inserted or deleted from both sequences. The inference of an ancestral sequence is always error prone, and is approached using statistical methods. A series of consecutive ``-''s is usually called a (single) gap. The number of dashes is referred to as the size, or length of the gap. Dashes at the end of a sequence are called ``end gaps''. Because of the insertion/deletion ambiguity, gaps are often called ``indels''.

SEQUENCE EDITING

Consider the process of turning ${\bf A}$ into ${\bf B}$ using a series of substitutions, insertions and deletions. We call this ``editing'' ${\bf A}$ into ${\bf B}$.We start by writing ${\bf A}$ from left to right without any gaps. We place a copy of ${\bf A}$ below it. We make a single change in this copy. If a substitution is made to a<sub>i</sub>, then a<sub>i</sub> is changed in the copy. If a<sub>i</sub> is deleted, then this character is replaced by a ``-'' in the copy. If an insertion occurs after a_i, then this element is added to the copy and a ``-'' is added at this position in the original sequence. We continue in this manner, getting k+1 rows for k edit operations. Every time an element is inserted, a gap character is added in the same position all the way up to the top. Eliminating all middle rows gives us an alignment of ${\bf A}$ with ${\bf B}$.

Conversely, an alignment is a record of editing. Matched residues are not edited. a_i over b_j means that b_j is substituted for a_i. a_i over a gap means that a_i is deleted. A gap over b_j indicates that b_j has been inserted.

MINIMUM EDIT DISTANCE

It is clear that we can edit ${\bf A}$ to ${\bf B}$ using at most m+n edit operations. That is, we can delete all the residues of ${\bf A}$, one by one, and then add all the residues of ${\bf B}$, one by one. In fact, we could do this in at most $\max\{m,n\}$ operations, using $\min\{m,n\}$ substitutions followed by $\max\{m,n\}-\min\{m,n\}$ insertions or deletions. This is the best that could be done to turn ACACACACACA into GTGGTGTTGTGTGT.

ACACACACACA---

GTGGTGTTGTGTGT

In general, we can do better.

Let ${\cal D}({\bf A},{\bf B})$ be the minimum number of edits required to turn ${\bf A}$ into ${\bf B}$.

Proposition: ${\cal D}$ is a metric

Proof:

• If ${\bf A} = {\bf B}$, then ${\cal D}({\bf A},{\bf B}) = 0$, since no editing is required.
• ${\cal D}({\bf A},{\bf B}) = {\cal D}({\bf B},{\bf A})$, since any series of edits turning ${\bf A}$ into ${\bf B}$ can be performed in reverse to turn ${\bf B}$ into ${\bf A}$.
• If ${\bf A}, {\bf B}$ and ${\bf C}$ are 3 sequences, then one can edit ${\bf A}$ into ${\bf B}$ in ${\cal D}({\bf A},{\bf C}) + {\cal D}({\bf B},{\bf C})$ steps by first editing ${\bf A}$ into ${\bf C}$ and then editing ${\bf C}$ into ${\bf B}$. Thus ${\cal D}({\bf A},{\bf B}) \leq {\cal D}({\bf A},{\bf C}) + {\cal D}({\bf B},{\bf C})$

How can we compute ${\cal D}({\bf A},{\bf B})$? For $1 \leq i \leq m$, let ${\bf A}_i$ be the i^th prefix of ${\bf A}$. That is ${\bf A}_i = a_1a_2\dots a_i$. Similarly, ${\bf B}_j = b_1b_2\dots b_j$ for $1 \leq j \leq n$. Let $D_{i,j} = {\cal D}({\bf A}_i,{\bf B}_j)$. We define ${\bf A}_0$ and ${\bf B}_0$ to be the empty string. It is clear that $D_{i,j} = \max\{i,j\}$ if i or j are 0. That is, D_i,0 = i and D_0,j = j. Suppose that we know D_h,k for all $0 \leq h \leq i$ and $0 \leq k \leq j$, but not for (i,j). Can we compute D_i,j?

In fact, if i > 0 and j > 0, then

$$D_{i,j} = \min\{D_{i-1,j}+1 , D_{i,j-1}+1, D_{i-1,j-1}+1-\delta(a_i,b_j)\}$$

Proof:

• We can edit ${\bf A}_i$ to ${\bf A}_{i-1}$ in 1 step by deleting a_i. We can then edit ${\bf A}_{i-1}$ to ${\bf B}_{j}$ in D_i-1,j steps. This gives us an edit from ${\bf A}_i$ to ${\bf B}_{j}$ in D_i-1,j+1 steps.
• We can edit ${\bf B}_{j}$ to ${\bf B}_{j-1}$ in 1 step by deleting b_j. We can then edit ${\bf B}_{j-1}$ to ${\bf A}_i$ in D_i,j-1 steps. This gives us an edit from ${\bf B}_{j}$ to ${\bf A}_i$ in D_i,j-1+1 steps.
• We and edit ${\bf A}_{i-1}$ to ${\bf B}_{j-1}$ in D_i-1,j-1 steps. We can apply the same edit operations to ${\bf A}_i$ by leaving a_i untouched. We end up with ${\bf B}_{j-1}$ with a_i added to the end. If a_i = b_j, there is nothing more to do. The ``cost'' is $0 = 1 - \delta(a_i,b_j)$ Otherwise, we can change a_i to b_j in $1 = 1 - \delta(a_i,b_j)$ step. Thus we have transformed ${\bf A}_{i-1}$ to ${\bf B}_{j-1}$ in $D_{i-1,j-1}+1-\delta(a_i,b_j)$ steps.

Thus

$$D_{i,j} \leq \min\{D_{i-1,j}+1 , D_{i,j-1}+1, D_{i-1,j-1}+1-\delta(a_i,b_j)\}$$

The equality follows from the fact that at least 1 of the following must occur in transforming ${\bf A}_i$ to ${\bf B}_{j}$

• a_i is deleted. This step costs 1. Transforming ${\bf A}_{i-1}$ to ${\bf B}_{j}$ costs another D_i-1,j.
• b_j is inserted. This step costs 1. Transforming ${\bf A}_i$ to ${\bf B}_{j-1}$ costs another D_i,j-1.
• Neither a_i nor b_j is deleted. Changing a_i to b_j costs $1-\delta(a_i,b_j)$. Editing the rest costs D_i-1,j-1.

Working in term of ``edit distance'' is cumbersome. It is easier to work with alignments. If ${\cal A}$ is the sequence alphabet, let $d: {\cal A}' \times {\cal A}' \rightarrow [0,\infty)$, where ${\bf A}'$ is ${\cal A} \cup \{-\}$. That is, we add the dash character to the alphabet. We call d the ``(edit) cost function''. Assume that $a, b \in {\cal A}$. We assume that d(a,a) = 0. d will usually be symmetric, but it does not have to be. If it is not symmetric, then ${\cal D}$ will not be a metric.

• If $a \neq b$, d(a,b) > 0 and d(a,b) is called a substitution cost.
• d(a,-)>0 is call a deletion cost (the cost of deleting a) or an insertion cost (the cost of inserting a).
• d(-,b)>0 is call an insertion cost (the cost of inserting b) or a deletion cost (the cost of deleting b).

At this point, we will change the original definition of alignment by dropping the first condition that a_ih=b_jh for all matched pairs. The d function will make that distinction for us.

Let $(i_1,j_1), (i_2,j_2), (i_3,j_3), \dots, (i_k,j_k)$ be an alignment between ${\bf A}$ and ${\bf B}$. Put this alignment into ``alignment format'':

$$\begin{array}{lccccc} {\bf A}' = & a'_1 & a'_2 & a'_3 & \dot... ...\ {\bf B}' = & b'_1 & b'_2 & b'_3 & \dots b'_{k'} \end{array}$$

The a'_is and b'_js are either residues or gap characters. Deleting the gaps from ${\bf A}'$ and ${\bf B}'$ gives us back ${\bf A}$ and ${\bf B}$. Then the cost of the alignment is given by $$\sum_{h=1}^{k'} d(a'_h,b'_h)$$. This also gives us a ``generalized edit distance'', and a generalized minimum edit distance. If d(a,-) = d(-,b) = 1 and d(a,b) = 1 for $a \neq b$, then we have the old definition. ${\cal D}({\bf A},{\bf B})$ is the minimum edit distance between ${\bf A}$ and ${\bf B}$, which is the same as the minimum alignment cost. What is D_i,j now? We'll give a formula and proof thinking about alignments, not edits.
Clearly, $D_{i,0} = \sum_{h=1}^{i} d(a_h,-)$ and $D_{0,j} = \sum_{h=1}^{j} d(-,b_h)$. Suppose that we know D_h,k for all $0 \leq h \leq i$ and $0 \leq k \leq j$, but not for (i,j). Then:

$$\begin{array}{cl} D_{i,j} = \min\{ & D_{i-1,j}+d(a_i,-) , D_{i,j-1}+d(-,b_j), \\ & D_{i-1,j-1}+d(a_i,b_j)\} \end{array}$$

This recursion is called the ``fill algorithm''.

Proof: An alignment between ${\bf A}_i$ and ${\bf B}_{j}$ can be created by

• taking a minimum cost alignment between ${\bf A}_{i-1}$ and ${\bf B}_{j}$ and adding a_i/- to the end,
• taking a minimum cost alignment between ${\bf A}_i$ and ${\bf B}_{j-1}$ and adding -/b_j to the end,
• taking a minimum cost alignment between ${\bf A}_{i-1}$ and ${\bf B}_{j-1}$ and adding a_i/b_j to the end.

This establishes the $\leq$ part. To establish equality, consider an optimal alignment of ${\bf A}_i$ and ${\bf B}_{j}$. The last 2 characters are either a_i/-, -/b_j or a_i/b_j. This yields equality for cases 1, 2 or 3, respectively.

EXAMPLE TTATGGACTT versus CTTGGCTAGG. d = 2 and 1 for all substitutions and indels, respectively.

		C	T	T	G	G	C	T	A	G	G
	0	1	2	3	4	5	6	7	8	9	10
T	1	2	$\cdot$1	$\cdot$2	3	4	5	$\cdot$6	7	8	9
T	2	3	$\cdot$2	$\cdot$1	2	3	4	$\cdot$5	6	7	8
A	3	4	3	2	3	4	5	6	$\cdot$5	6	7
T	4	5	$\cdot$4	$\cdot$3	4	5	6	$\cdot$5	6	7	8
G	5	6	5	4	$\cdot$3	$\cdot$4	5	6	7	$\cdot$6	$\cdot$7
G	6	7	6	5	$\cdot$4	$\cdot$3	4	5	6	$\cdot$7	$\cdot$6
A	7	8	7	6	5	4	5	6	$\cdot$5	6	7
C	8	$\cdot$7	8	7	6	5	$\cdot$4	5	6	7	8
T	9	8	$\cdot$7	$\cdot$8	7	6	5	$\cdot$4	5	6	7
T	10	9	$\cdot$8	$\cdot$7	8	7	6	$\cdot$5	6	7	8

The D matrix after the fill algorithm is executed. The dots, $\cdot$, indicate positions where a_i = b_j.

		C	T	T	G	G	C	T	A	G	G
	0	1	2	3	4	5	6	7	8	9	10
T	1
T	2
A	3
T	4
G	5
G	6
A	7
C	8
T	9
T	10

The D matrix before filling. The 0^th row and column have been defined.

The value of D_i,j can be determined if D_h,k is known for $0 \leq h \leq i$, $0 \leq k \leq j$, $(h,k) \neq (i,j)$. There are different ``fill orders'' that can be used.

Column fill :

  for (j=1; j<=m; j++){
     for (i=1; i <=n; i++){
        Define D[i,j];
     }
  }

Row fill :

  for (i=1; i<=n; i++){
     for (j=1; j <=m; j++){
        Define D[i,j];
     }
  }

The set

$${\rm codia}_d = \{(i,j) \: \vert \: 1 \leq i \leq m, \: 1 \leq j \leq n, \: j+i=d\}$$

is called the ``d^th co-diagonal''. d varies from 2 to m+n. If $(i,j) \in {\rm dia}_d$, then D_i,j can be computed using values of D_h,k for (h,k) in ${\rm codia}_{d'}$, d' < d.

Diagonal fill :

   for (d=2; d<=m+n; d++){
      for (i=min{1,d-n}; i<=min{m,d-1}; i++){
         j = d - i;
         Define D[i,j];
      }
   }

Note that with the current algorithm, we only need to keep the present and previous column (row) to execute the fill algorithm in column (row) fill order. For the diagonal fill, we need the 2 previous co-diagonals.

TRACEBACK : After the ``fill algorithm'' is executed, an alignment can be constructed using the ``traceback algorithm''. The alignment is constructed from right to left; that is, backwards. This is because we start with D_m,n, and determine which of the 3 cases resulted in the value of D_m,n. This process continues. A traceback algorithm is explained below and a slight variant is illustrated in Figure 2.

1.

Set $i=m, \: j=n$. The alignment is empty. Proceed to 2.

2.

Test for:

(a)

D_i,j = D_i-1,j + d(a_i,-). Proceed to 3.

(b)

D_i,j = D_i,j-1 + d(-,b_j). Proceed to 4.

(c)

D_i,j = D_i-1,j-1 + d(a_i,b_j). Proceed to 5.

3.

Set i=i-1. Add a_i/- to the alignment. Proceed to 6.

4.

Set j=j-1. Add -/b_j to the alignment. Proceed to 6.

5.

Set $i=i-1, \: j=j-1$. Add a_i/b_j to the alignment. Proceed to 6.

6.

Test for i=0 or j=0:

(a)

$i>0, \: j>0$. Proceed to 2

(b)

$i=0, \: j>0$. Add $--{\dots}-/b_1b_2{\dots}b_j$ to the alignment. Stop.

(c)

$i>0, \: j=0$. Add $a_1a_2{\dots}a_i/--{\dots}-$ to the alignment. Stop.

(d)

$i=0, \: j=0$. Stop.

  

Figure 2: In this variant of the traceback algorithm, only the matched pairs, a_i/b_j are listed in the alignment. The gaps are inferred.

		C	T	T	G	G	C	T	A	G	G
	0	1	2	3	4	5	6	7	8	9	10
T	1	$\uparrow$2	1	2	3	4	5	6	7	8	9
T	2	3	$\nwarrow$2	1	2	3	4	5	6	7	8
A	3	4	$\uparrow$3	2	3	4	5	6	5	6	7
T	4	5	4	$\nwarrow$3	4	5	6	5	6	7	8
G	5	6	5	4	$\nwarrow$3	4	5	6	7	6	7
G	6	7	6	5	4	$\nwarrow$3	4	5	6	7	6
A	7	8	7	6	5	$\uparrow$4	5	6	5	6	7
C	8	7	8	7	6	5	$\nwarrow$4	5	6	7	8
T	9	8	7	8	7	6	5	$\nwarrow$4	$\leftarrow$5	$\leftarrow$6	$\leftarrow$7
T	10	9	8	7	8	7	6	5	6	7	$\uparrow$8

The execution of the traceback is illustrated above. The symbols: ``$\uparrow$'', ``$\leftarrow$'' and ``$\nwarrow$'' mean:

1.

Move up (Set i=i-1),

2.

Move left (Set j=j-1) and

3.

Move up and left. (Set i=i-1 and j=j-1), respectively.

The resulting alignment is:

-TTATGGACT---T
  | ||| ||
C-T-TGG-CTAGG-

GAP PENALTIES :

We have already introduced the notion of a gap as a series of insertions or deletions in one sequence with respect to another in the context of an alignment. The ``size'' of a gap is the number of indels. Our definition of a gap allows for consecutive indels in one or both sequences. In the alignment above, -/C and T/- are consecutive entries that contain indels. We call this a single gap. Similarly, the four last columns: -/A, -/G, -/G and T/- constitute 1 gap, not 2. (This is not standard usage.)

The cost of a gap often depends only on its size, s. A function $g:\{0,1,2,\dots\} \rightarrow [0,\infty)$, which gives the cost g(s) to a gap of size s, is called a ``gap cost'' or ``gap penalty'' function. g(0) is always assumed to be 0. In our general definition, d(a,-) = d(-,a) can take on different values for each letter of the (DNA or amino acid) alphabet. If the value of d(a,-) is the same $\forall a$, the a gap of size s has cost $\beta s$, where $\beta = d(a,-)$. A gap penalty function of the form $g(s)=\alpha+\beta s$ is called an ``affine gap penalty function'', where $\alpha$ and $\beta$ are non-negative constants.

Alignment cost with gap cost functions.

In this case, d(a_i,b_j) still gives the cost of substituting a_i with b_j, but d(a_i,-) and d(-,b_j) are no longer used. The cost of an alignment cannot be determined by adding up one summand for each column. If there are $\Delta_s$ gaps of size s in an alignment, then the cost of an alignment between ${\bf A}$ and ${\bf B}$ is

$$\sum_{h=1}^{k}d(a_{i_h},b_{j_h}) + \sum_{s=1}^{\infty} \Delta_s g(s)$$

At this time, we add 2 ``fake matches'' to every alignment to take care of the end gaps. This means that we do not need to list the end gaps as special cases. We will assume that every alignment (not the ``local'' alignments that will be considered later) begins with (i₀,j₀) = (0,0) and ends with (i_k+1,j_k+1) = (m+1,n+1). We assume that these ``fake matches'' have 0 weight. Note that a gap occurs whenever i_h-i_h-1+j_h-j_h-1>2, and the gap size is i_h-i_h-1+j_h-j_h-1-2, for $1 \leq h \leq k+1$. Thus:

$$\sum_{h=1}^{k+1}g(i_h-i_{h-1}+j_h-j_{h-1}-2) = \sum_{s=1}^{\infty} \Delta_s g(s).$$

There are algorithms to compute minimum cost alignments for such gap functions, but we will not consider them. Instead, we will deal with these gap functions when we look at ``maximum similarity'' alignment.

PATH GRAPHS

Alignments can be visualized in so-called ``path graphs'' that are related to dot plots. Alignments are constructed by eye using dot plots. Path graphs allow one to define the set of all possible alignments in a formal, yet intuitive way.

One type of path graph is what I call the ``truck route'' path graph. I use this term because the path is connected. It cannot hop from one location to another. Truck route path graphs formalize the representation that shows matches and gaps explicitly.

An m+1 x n+1 square grid of points is constructed. The grid can described as $\{(i,j) \; \vert \; 0 \leq i \leq m, 0 \leq j \leq m\}$, where i refers to row, and j to column. Two distinct points (i,j) and (i',j') have an edge between them if and only if $\max \{ \vert i-i'\vert, \vert j-j'\vert \} = 1$.

Any collection of connected edges from (0,0) to (m+1,n+1) that moves from each node (i-1,j-1) to either (i,j-1), (i,j) or to (i-1,j), describes an alignment.

When a diagonal edge, $(i-1,j-1) \mapsto (i,j)$ is used, we say that a(i) is matched with b(j).

When a vertical edge, $(i-1,j-1) \mapsto (i,j-1)$ is used, we say that a(i) is inserted with respect to sequence ${\bf b}$, or that a(i) is in a gap.

Similarly, when a horizontal edge, $(i-1,j-1) \mapsto (i-1,j)$ is used, we say that b(j) is inserted with respect to sequence ${\bf a}$, or that b(j) is in a gap.

Figure 3 illustrates an alignment together with its path graph. These kinds of path graphs are the most useful ones for looking at multiple sequence alignments. They are not good for describing alignments with affine gap penalties.

  

Figure 3: An alignment of 2 DNA sequences together with its ``truck route'' path graph. Every column in the alignment corresponds to an edge in the graph.

  

Figure 4: An alignment of 2 DNA sequences together with its ``airline route'' path graph. Matched pairs of residues correspond to vertices in path graph. Gaps are not shown explicitly, but are inferred when the path jumps more than 1 unit in either direction. The yellow edges are part of the ``traceback path'' that will be discussed later.

In contrast, the ``airline route'' path graph uses an m x n grid, and matched pairs correspond to vertices, not edges. A path consists of a series of vertices, as defined above. Gaps are inferred when either i<sub>l+1</sub> - i<sub>l</sub> > 1 or j<sub>l+1</sub> - j<sub>l</sub> > 1. The alignments that have been introduced in the course up to now are best visualized using ``airline route'' path graphs.

Consider the m x n array of integer pairs: ${\cal G}_{m,n} = \{(i,j) \: \vert \: 0 \leq i \leq m, \: 0 \leq j \leq n\}$. Consider all line segments joining (i,j) to (i-1,j), (i,j) to (i,j-1) and (i,j) to $(i-1,j-1), \: \forall \: (i,j) \in {\cal G}_{m,n}$. (We exclude negative indices.) A ``path'' from (0,0) to (m,n) is a sequence of these connected line segments, joining distinct pairs $(0,0)=(i_0,j_0), (i_1,j_1), (i_2,j_2), (i_3,j_3), \dots (i_k,j_k)= (m,n)$, where $0 \leq i_h - i_{h-1} \leq 1$ and $0 \leq j_h - j_{h-1} \leq 1, \: \forall \: 0<h\leq k$. There is a 1:1 correspondence between paths and alignments, with the segment:

• (i,j) to (i-1,j-1) corresponding to a_i/b_j in an alignment,
• (i,j) to (i-1,j) corresponding to a_i/- in an alignment and
• (i,j) to (i,j-1) corresponding to -/b_j in an alignment.

SIMILAR (DNA) SEQUENCES

If 2 (DNA) sequences are known to be fairly similar, say 90% or better and n and m are large (say 10,000 or larger) and ``roughly the same size'', then much time (and space) can be saved in alignment.

Diagonal band fill

The set

$${\rm dia}_d = \{(i,j) \: \vert \: 1 \leq i \leq m, \: 1 \leq j \leq n, \: j-i=d\}$$

is called the ``d<sup>th</sup> diagonal''. d varies from 1-m to n-1. If ${\bf A}$ and ${\bf B}$ are roughly of the same length and are fairly similar, then a minimum cost alignment will lie ``close to'' the 0^th diagonal, ${\rm dia}_0$. This means that for column fill, it suffices to consider only values of i that satisfy $\vert j-i\vert \leq bw$ when filling the j^th column, where bw (band width) is a constant. In practice, bw might be 100 when m and n are many thousands. The values of D_i,j just outside the boundaries are assumed to be $+\infty$ (ie., a large number). There is a risk of computing an alignment that is not minimum cost. In practice, bw can be increased and new alignments can be computed until the minimum cost no longer changes. When m and n are not roughly equal, the bounds for i in the j^th column can be computed by demanding that $\vert i-[\frac{jm}{n}]\vert \leq bw$.

MAXIMUM SIMILARITY ALIGNMENT

As before, ${\cal A}$ is the sequence alphabet, and $w: {\cal A} \times {\cal A} \rightarrow {\bf R}$. We call w the ``weight function'' that gives a measure of similarity. For example, $w(a,b)=\delta(a,b)$ would be a natural and simple weight function. So would $w(a,b)=2\delta(a,b)-1$. We assume that w is symmetric.

We can think of w as the inverse of d. However, gap costs are treated the same as before. The definition of alignment remains the same. This time, we try to maximize the similarity of an alignment, so that the gap penalties (gap costs) are subtracted, not added. If there are $\Delta_s$ gaps of size s in an alignment, then the cost of an alignment between ${\bf A}$ and ${\bf B}$ is

$$\sum_{h=1}^{k}w(a_{i_h},b_{j_h}) - \sum_{s=1}^{\infty} \Delta_s g(s)$$

If $g(s) = {\beta}s$, then we can define $w(a,-) = w(-,b) = -\beta$ for all letters a and b. In this case, the alignment similarity is $$\sum_{h=1}^{k} w(a'_h,b'_h)$$, where a' and b' are characters in ${\cal A}' = {\cal A} \cup \{-\}$. In this case, we define $${\cal S}({\bf A},{\bf B})$$ to be the maximum similarity of all alignments between ${\bf A}$ and ${\bf B}$.

Sufficient condition for ${\cal D}({\bf A},{\bf B})$ to be equivalent to ${\cal S}({\bf A},{\bf B})$.

Suppose that w and d are score functions for sequence alignment by similarity maximization and distance minimization, respectively. Let the gap penalty for distance minimization be given by some function g. We wish to construct a maximum similarity alignment.

Suppose, as well, that w(a,b) + d(a,b) = C for some constant C, and that $g(s) \geq \frac{Cs}{2}$. Then distance minimization alignment is equivalent to similarity maximizing alignment.

Proof: For any alignment, there are k matched pairs and a total of m+n-2k elements in gaps, including possible end gaps.

Let $\Delta_l$ be the number of gaps of size l. Then

$$m+n-2k = \sum_{l=1}^{\infty} l \Delta_{l}.$$

Thus

$$\begin{array}{cl} & \displaystyle \sum_{h=1}^{k} d(a_{i_h},b... ... \sum_{l=1}^{\infty} (g(l) -Cl/2)\Delta_l \right ]. \end{array}$$

Hence, ${\cal D}({\bf A},{\bf B}) = C(m+n)/2 - {\cal S}({\bf A},{\bf B})$ if we take g(s) - Cs/2 to be the similarity gap penalty (non-negative).

Example The familiar case where $d(a',b') = 1 - \delta(a',b')$ gives what is called the ``Levenstein distance'' between 2 sequences. In this case, g(s) = s. If we let $w(a,b)=\delta(a,b)$, which seems ``natural'', then C=1, so g'(s) = s/2 as a gap penalty for maximum similarity alignment would give identical results.

For a general gap function, g, how can we compute ${\cal S}({\bf A},{\bf B})$ and a maximum similarity alignment?

Let S_i,j be the maximum similarity of an alignment between ${\bf A}_i$ and ${\bf B}_{j}$. Note that S_0,j = -g(j) and that S_i,0 = -g(i). What sort of recursion can be used in this case? The answer is given below.

$$\begin{array}{ll} S_{i,j} = \max \left\{ \right. & S_{i-1,j-... ...0 \end{array} } \left\{S_{i-l,j-k}-g(k+l) \right\}. \end{array}$$

Proof : An optimal alignment either ends with a match between a_i and b_j, or with an end gap. The S_i-1,j-1+w(a_i,b_j) term adds the weight of an a_i/b_j match to an optimal alignment between ${\bf A}_{i-1}$ and ${\bf B}_{j-1}$. This either extends a match or closes an end gap. The second term considers all possible ``last matched pair'' cases.

While this formula is straightforward, computation time is ${\cal O}(m^{2}n^{2})$, since one must backtrack for each (i,j) pair. We can do better when g has a special form.

Case 1. $g(s) = {\beta}s$, where $\beta$ is a positive constant. This is the ``linear gap penalty'' case.

$$S_{i,j} = \max \left \{ S_{i-1,j-1}+w(a_i,b_j), S_{i-1,j}- \beta, S_{i,j-1} - \beta \right \} .$$

Case 2. $g(s)=\alpha+\beta s$, where $\alpha$ and $\beta$ are non-negative constants. (Both are not 0.) This is the ``affine gap penalty'' case. The linear gap penalty recursion breaks down because when we add an unpaired a_i to an optimal alignment between ${\bf A}_{i-1}$ and ${\bf B}_{j}$ or b_j to an optimal alignment between ${\bf A}_i$ and ${\bf B}_{j-1}$ we don't know if we are creating a new gap, with a cost of $\alpha + \beta$, or just extending an existing one, with a cost of only $\beta$. This can be solved in a number of ways. They all involve the use of at least one extra array that must be filled.

``Zuker's method'': Let S¹_i,j be the maximum similarity between ${\bf A}_i$ and ${\bf B}_{j}$ that ends with a gap. That is, a_i and b_j are not matched with each other. Let S²_i,j be the maximum similarity between ${\bf A}_i$ and ${\bf B}_{j}$ that ends with an a_i/b_j match. Clearly

$${\cal S}({\bf A},{\bf B}) = \max \{ S^{1}_{m,n} , S^{2}_{m,n} \}$$

Initial conditions.

• $S^{1}_{0,0} = -\infty$. S¹_i,0 = g(i) for i>0 and S¹_0,j = g(j) for j > 0.
• S²_0,0 = 0, $S^{2}_{i,0} = -\infty$ for i>0 and $S^{2}_{0,j} = -\infty$ for j>0.

Recursions:

$$\begin{array}{lcll} S^{1}_{i,j} & = & \max \{ & S^{1}_{i-1,j... ...& S^{1}_{i-1,j-1}, S^{2}_{i-1,j-1} \} + w(a_i,b_j). \end{array}$$

Proof : We consider cases. For S¹_i,j:

• a_i is unpaired and b_j is paired.
  • b_j is paired to a_i-1. This gives us an end gap of size 1, and so the score is S²_i-1,j - g(1).
  • b_j is paired to a_i' for i' < i-1. This gives us an end gap of size >1, and so the score is $S^{1}_{i-1,j} - \beta$ since a_i-1 cannot be paired.
• a_i is paired and b_j is unpaired.
  • a_i is paired to b_j-1. This gives us an end gap of size 1, and so the score is S²_i,j-1 - g(1).
  • a_i is paired to b_j' for j' < j-1. This gives us an end gap of size >1, and so the score is $S^{1}_{i,j-1} - \beta$ since b_j-1 cannot be paired.
• Both a_i and b_j are unpaired. In this case, the score is $S^{1}_{i-1,j} - \beta = S^{1}_{i,j-1} - \beta$.

For S²_i,j, a_i is paired with b_j. a_i-1 and b_j-1 are either paired with each other, or they are not.

Alignment examples.

 AACAGGCAGGCCATT--AGGGAG--CTACGG-GTTGTGCCGACCGTCCAAT---AATGGTCGGCCTCTC
 |: |:|  :|:|||:  || :||  || ||: :||:||  ||   |||:|    |||  |:::| |      
 AGAAAG--AGTCATCTGAGCAAGGACT-CGATATTATG--GA---TCCGA-AGAAAT--TTAACGT---

In this example, identity is indicated by |, transitions ( $A{\leftrightarrow}G$, $C{\leftrightarrow}T$) by :, and transversions ( $A{\leftrightarrow}C$, $A{\leftrightarrow}T$, $G{\leftrightarrow}C$, $G{\leftrightarrow}T$) by blanks, respectively.

Example: ${\bf A}$ = AGGCTACGG and ${\bf B}$ = AGGGACTCGAT. w = 10, 2 and -5 for identity, transitions and transvertions, respectively.

Optimal alignment:

 AGG--CTACGG-
 |||  || ||: 
 AGGGACT-CGAT

S¹:

		A	G	G	G	A	C	T	C	G	A	T
	$-\infty$	-11	-12	-13	-14	-15	-16	-17	-18	-19	-20	-21
A	-11	-12	-1	-2	-3	-4	-5	-6	-7	-8	-9	-10
G	-12	-1	-2	9	8	7	6	5	4	3	2	1
G	-13	-2	9	8	19	18	17	16	15	14	13	12
C	-14	-3	8	19	18	17	16	17	16	15	14	13
T	-15	-4	7	18	17	16	17	16	27	26	25	24
A	-16	-5	6	17	16	15	16	27	26	25	24	25
C	-17	-6	5	16	15	16	15	26	25	26	25	24
G	-18	-7	4	15	14	15	26	25	26	25	36	35
G	-19	-8	3	14	15	14	25	24	25	36	35	38

S²:

		A	G	G	G	A	C	T	C	G	A	T
	0	$-\infty$	$-\infty$	$-\infty$	$-\infty$	$-\infty$	$-\infty$	$-\infty$	$-\infty$	$-\infty$	$-\infty$	$-\infty$
A	$-\infty$	10	-9	-10	-11	-4	-20	-21	-22	-16	-9	-25
G	$-\infty$	-9	20	9	8	-1	-9	-10	-11	3	-6	-14
G	$-\infty$	-10	9	30	19	10	2	1	0	14	5	-3
C	$-\infty$	-18	-7	4	25	14	28	19	26	10	9	15
T	$-\infty$	-19	-8	3	14	20	19	38	21	21	10	24
A	$-\infty$	-5	-2	9	20	27	15	14	33	29	36	20
C	$-\infty$	-21	-10	1	12	15	37	18	37	28	24	38
G	$-\infty$	-15	4	15	26	17	11	32	21	47	30	20
G	$-\infty$	-16	3	14	25	28	12	21	27	36	49	31

Up to now, we have been dealing with what is called ``global alignment''. When dealing with edit distance, we must transform one sequence into another. As we dealt with minimum cost alignment and maximum similariy alignment, we always kept the notion that we are dealing with all of ${\bf A}$ and ${\bf B}$. In particular, we penalize ``end gaps''. We will now consider another class of alignment problem.

Local sequence alignment :

If $1 \leq i_1 \leq i_2 \leq m$ and $1 \leq j_1 \leq j_2 \leq n$, let ${\bf A}_{i_1,i_2} = a_{i_1}a_{i_{1}+1}{\dots}a_{i_2}$ and let ${\bf B}_{j_1,j_2} = b_{j_1}b_{j_{1}+1}{\dots}b_{j_2}$. Let

$${\cal H}({\bf A},{\bf B}) = \max \left \{ 0, \right. \begin{... ...cal S}({\bf A}_{i_1,i_2},{\bf B}_{j_1,j_2}) \left. \right \} .$$

For a general gap function, g, how can we compute ${\cal H}({\bf A},{\bf B})$ and a maximum similarity alignment? Such an alignment is called a ``local alignment''.

Let H_i,j be the maximum similarity of an alignment between ${\bf A}_i$ and ${\bf B}_{j}$ that includes either the (a_i,b_j) match or else counts an end gap, or is empty.

Then

$${\cal H}({\bf A},{\bf B}) = \begin{array}{c} \max \\ i,j \end{array} \{ H_{i,j} \}.$$

Note that H_i,0 = H_0,j = 0, since these values can only be $\max \{0,-g(i)\}$ or $\max \{0,-g(j)\}$. The recursion for i>0 and j>0 is:

$$\begin{array}{ll} H_{i,j} = \max \left\{0, \right. & H_{i-1,... ...0 \end{array} } \left\{H_{i-l,j-k}-g(k+l) \right\}. \end{array}$$

Proof : An optimal alignment either ends with a match between a_i and b_j, or with an end gap. The H_i-1,j-1+w(a_i,b_j) term adds the weight of an a_i/b_j match to an optimal local alignment between ${\bf A}_{i-1}$ and ${\bf B}_{j-1}$. This either extends a match or closes an end gap. The second term considers all possible ``last matched pair'' cases. If the last 2 terms are negative, then the 0 option (no alignment) can be used.

While this formula is straightforward, computation time is ${\cal O}(m^{2}n^{2})$, since one must backtrack for each (i,j) pair. As with S, we can do better when g has a special form. In these forms, the algorithm is named ``Smith-Waterman'' after Temple Smith and Michael Waterman who first published it.

Case 1. $g(s) = {\beta}s$, where $\beta$ is a positive constant. This is the ``linear gap penalty'' case.

$$\begin{array}{ll} H_{i,j} = \max \left \{0, \right. & H_{i-1... ...ft. H_{i-1,j}- \beta, H_{i,j-1} - \beta \right \} . \end{array}$$

Case 2. Affine gap penalty. $g(s)=\alpha+\beta s$.

``Zuker's method'': Let H¹_i,j be the maximum similarity between ${\bf A}_i$ and ${\bf B}_{j}$ that ends with a gap. That is, a_i and b_j are not matched with each other. Let H²_i,j be the maximum similarity between ${\bf A}_i$ and ${\bf B}_{j}$ that ends with an a_i/b_j match. It is understood that either H¹_i,j or H²_i,j will be 0 (no alignment) if there is no alignment with a non-negative score.

Clearly

$${\cal H}({\bf A},{\bf B}) = \begin{array}{c} \max \\ i,j \end{array}\{ H^{2}_{i,j} \}$$

The initial conditions are simple. H¹_i,0 = H¹_0,j = H²_i,0 = H²_0,j = 0.

Recursions:

$$\begin{array}{lcll} H^{1}_{i,j} & = & \max \{0, & H^{1}_{i-1... ...i-1,j-1}+w(a_i,b_j), H^{2}_{i-1,j-1}+w(a_i,b_j) \}. \end{array}$$

Proof : We consider cases. For H¹_i,j:

• There is no alignment with positive score. Otherwise, there is such an alignment.
• a_i is unpaired and b_j is paired.
  • b_j is paired to a_i-1. This gives us an end gap of size 1, and so the score is H²_i-1,j - g(1).
  • b_j is paired to a_i' for i' < i-1. This gives us an end gap of size >1, and so the score is $H^{1}_{i-1,j} - \beta$ since a_i-1 cannot be paired.
• a_i is paired and b_j is unpaired.
  • a_i is paired to b_j-1. This gives us an end gap of size 1, and so the score is H²_i,j-1 - g(1).
  • a_i is paired to b_j' for j' < j-1. This gives us an end gap of size >1, and so the score is $H^{1}_{i,j-1} - \beta$ since b_j-1 cannot be paired.
• Both a_i and b_j are unpaired. In this case, the score is $H^{1}_{i-1,j} - \beta = H^{1}_{i,j-1} - \beta$.

For H²_i,j, if there is an alignment with positive score, then a_i is paired with b_j. a_i-1 and b_j-1 are either paired with each other, or they are not.

Alignment examples. Identity, transition and transversion weights are 1, 0 and -1, respectively. $\alpha = \beta = 1$.
${\bf A}$ = TCTTCTCCAAGGCGTTAACT and
${\bf B}$ = AACTTCGTTTGAGGCTTCTT

        10   
CTTC-TCCAAGGC
|||| |:::||||
CTTCGTTTGAGGC
      10

H¹

		A	A	C	T	T	C	G	T	T	T	G	A	G	G	C	T	T	C	T	T
	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
T	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
C	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
T	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1
T	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	1	1	0
C	0	0	0	0	0	1	0	2	1	0	0	0	0	0	0	0	0	1	0	2	2
T	0	0	0	0	0	0	2	1	1	1	0	1	0	0	0	0	0	0	2	1	3
C	0	0	0	0	0	0	1	1	1	1	1	0	0	0	0	0	0	0	1	3	2
C	0	0	0	0	0	0	0	1	1	1	1	1	0	0	0	0	0	0	0	2	3
A	0	0	0	0	0	0	1	0	1	1	1	0	1	0	0	0	0	0	1	1	2
A	0	0	0	0	0	0	0	1	0	0	0	1	0	2	1	0	0	0	0	0	1
G	0	0	0	0	0	0	0	0	0	0	0	0	2	1	3	2	1	0	0	0	0
G	0	0	0	0	0	0	0	0	0	0	0	0	1	3	2	4	3	2	1	0	0
C	0	0	0	0	0	0	0	0	0	0	0	0	0	2	4	3	5	4	3	2	1
G	0	0	0	0	0	0	0	0	0	0	0	0	0	1	3	5	4	4	3	2	1
T	0	0	0	0	0	0	0	0	0	1	0	0	0	0	2	4	4	4	5	4	3
T	0	0	0	0	0	0	0	0	1	0	2	1	0	0	1	3	4	5	5	5	4
A	0	0	0	0	0	0	0	0	0	2	1	1	0	0	0	2	3	5	5	4	4
A	0	0	0	0	0	0	0	0	0	1	1	0	1	1	0	1	2	4	4	4	3
C	0	0	0	0	1	0	0	0	0	0	0	1	1	0	0	0	1	3	3	3	3
T	0	0	0	1	0	2	1	0	0	0	0	0	0	0	0	0	0	2	3	2	4

H²

		A	A	C	T	T	C	G	T	T	T	G	A	G	G	C	T	T	C	T	T
	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
T	0	0	0	0	1	1	0	0	1	1	1	0	0	0	0	0	1	1	0	1	1
C	0	0	0	1	0	1	2	0	0	1	1	0	0	0	0	1	0	1	2	0	1
T	0	0	0	0	2	1	1	1	1	1	2	0	0	0	0	0	2	1	1	3	1
T	0	0	0	0	1	3	1	0	2	2	2	1	0	0	0	0	1	3	1	2	4
C	0	0	0	1	0	1	4	0	0	2	2	1	0	0	0	1	0	1	4	1	2
T	0	0	0	0	2	1	1	3	3	2	3	1	0	0	0	0	2	1	1	5	3
C	0	0	0	1	0	2	2	1	3	3	2	2	0	0	0	1	0	2	2	2	5
C	0	0	0	1	1	0	3	1	1	3	3	1	1	0	0	1	1	0	3	2	3
A	0	1	1	0	0	0	0	3	0	0	2	3	2	1	0	0	0	0	0	2	1
A	0	1	2	0	0	0	0	1	2	0	0	2	4	2	1	0	0	0	0	0	1
G	0	0	1	1	0	0	0	1	0	1	0	1	2	5	3	0	0	0	0	0	0
G	0	0	0	0	0	0	0	1	0	0	0	1	1	3	6	2	1	0	0	0	0
C	0	0	0	1	0	0	1	0	1	0	0	0	0	0	2	7	4	3	3	1	0
G	0	0	0	0	0	0	0	2	0	0	0	1	0	1	3	3	6	4	3	2	1
T	0	0	0	0	1	1	0	0	3	1	1	0	0	0	0	3	6	7	4	4	3
T	0	0	0	0	1	2	1	0	1	4	2	0	0	0	0	2	5	7	7	6	5
A	0	1	1	0	0	0	1	1	0	0	3	2	2	0	0	0	2	4	6	6	5
A	0	1	2	0	0	0	0	1	0	0	1	3	3	2	0	0	1	2	4	5	5
C	0	0	0	3	0	0	1	0	1	0	1	0	2	2	1	1	1	2	5	4	5
T