Sequence Alignment

Key ideas in solution:

• First, the alignment problem can be expressed as a graph problem:

  

  • The diagonal edges correspond to character matches (whether identical or not).
  • Each down edge corresponds to a gap in string 2.
  • Each across edge corresponds to a gap in string 1.
  • Goal: find the longest (max value) path from top-left to bottom-right.
    => longest-path problem in Directed Acyclic Graph.
  • Can be solved in O(mn) time (using topological sort)
    (m = length of string1, n = length of string2).
  • Suppose we label the vertices (i, j) where i is the row, and j is the column.
  • The best path to an intermediate node is via one of its neighbors: NORTH, WEST or NORTHWEST.

    

• Dynamic programming approach:
  • Building an actual graph (adjacency matrix etc) is unnecessary.
  • Let D_i,j
    = value of best path to node (i,j) from (0, 0).
    = score of optimal match of first i chars in string 1 to first j chars in string 2.
  • Define (for clarity)

    Ci-1,j-1	=	Di-1,j-1 + exact-match-score
    Ci-1,j	=	Di-1,j + string1-gap-score
    Ci,j-1	=	Di,j-1 + string2-gap-score

  • Thus,

    Di,j	=	max (Ci-1,j-1, Ci-1,j, Ci,j-1)
    	=	max (Di-1,j-1 + exact-match-score, Di-1,j + string1-gap-score, Di,j-1 + string2-gap-score)

  • Initial conditions:
    • D_0,0 = 0
    • D_i,0 = D_i-1,0 + gap-score for i-th char in string 1.
    • D_0,j = D_0,j-1 + gap-score for j-th char in string 2.

Implementation:

• Pseudocode:

  Algorithm: align (s1, s2)
  Input: string s1 of length m, string s2 of length n
  
  1.    Initialize matrix D properly;
  
        // Build the matrix D row by row.
  2.    for i=1 to m
  3.      for j=1 to n
  
            // Initialize max to the first of the three terms (NORTH).
  4.        max = D[i-1][j] + gapScore2 (s2[j-1])
  
            // See if the second term is larger (WEST).
  5.        if max < D[i][j-1] + gapScore1 (s1[i-1]) 
  6.          max = D[i][j-1] + gapScore1 (s2[i-1]) 
  7.        endif
  
            // See if the third term is the largest (NORTHWEST).
  8.        if max < D[i-1][j-1] + matchScore (s1[i-1], s2[j-1])
  9.          max = D[i-1][j-1] + matchScore (s1[i-1], s2[j-1])
  10.       endif
  
  11.        D[i][j] = max
  
  12.     endfor
  13.   endfor
  
        // Return the optimal value in bottom-right corner.
  14.   return D[m][n]
  
  Output: value of optimal alignment
     

• Explanation:
  • The 2D array D[i][j] stores the D_i,j values.
  • The method matchScore returns the value in matching two characters.
  • The method gapScore1 returns the value of matching a particular character in string 1 with a gap.
    (gapScore2 is similarly defined).
  • Note: D[i][0] and D[0][j] represent initial conditions.
  • Note: there is one more row than the length of string 1.
    (likewise for string 2).

• Computing the actual alignment:
  • A separate 2D array tracks which term contributed to the maximum for each (i, j).
  • After computing all of D, traceback and obtain path (alignment).

In-Class Exercise: (Solution) Download this template and implement the above algorithm. Gap and match costs are given.

Analysis:

• The double for-loop tells the story: O(mn) to compute D_m,n.

• Note: computing the actual alignment requires tracing back: O(m + n) time (length of a path).

• Space required: O(mn).

An improvement (in space):

• Genes are often longer than 10,000 bases (letters) long
  => space required is prohibitive.

• First, note that we only need successive rows in computing D
  => space requirement can be reduced to O(n).

• However, O(mn) space is still required to construct the actual alignment.

• A different approach:
  • Consider the "middle" character in string 1
    => the character at position m / 2.
  • In the optimal alignment, this character will either align with some j-th character in string 1, or a gap.
    => we can try all possible j values.
  • To try these alignments, compute the m/2-th row (in O(mn) time).
  • Output the optimal alignment found.
  • Recurse on either side of the optimal alignment.

• It is possible to show: O(mn) time and O(n) space.

Local Alignment - The Smith-Waterman Algorithm

Consider a global alignment example:

• An alignment is forced across the whole table.
  => each position is scored.

• However, the score in one particular section may be higher
  => local alignments may have higher scores.

• Goal in local alignment: find best possible local alignment.