Tuesday, November 5, 2013

Sequence Alignment using Dynamic Programming

Problem :
Input :
2 strings or sequence of characters i.e. X = x1x2...xm, Y = y1y2...yn

Goal - To define the similarity matrix between the 2 sequences with best alignment. This is based on penalty score, called Needleman/Wunsch score.

Eg. consider the sequence - AGGGCT and AGGCA , best alignment here is:
AGGGCT
 |  |  |      |
AGG - CA

Penalty
Penalty = penalty of match+penalty of difference + penalty of space

In the problem all the penalties are given.


For this example, the two sequences to be globally aligned are G A A T T C A G T T A (sequence #1)
G G A T C G A (sequence #2)
So M = 11 and N = 7 (the length of sequence #1 and sequence #2, respectively)
A simple scoring scheme is assumed where
  • Si,j = 1 if the residue at position i of sequence #1 is the same as the residue at position j of sequence #2 (match score); otherwise
  • Si,j = 0 (mismatch score)
  • w = 0 (gap penalty) 
Applying the DP
 Now lets apply the DP here. Suppose X is x1x2...xm, and Y=y1y2...yn. In this case, the last array position can have 3 possibilities (in top and bottom sequence respectively):
  1. xm and yn are aligned
  2. xm is aligned with space / gap
  3. space / gap  is aligned with yn
(Note that space cannot align with space, as it will be useless)

Optimal substructure
Let X' = X- {xm} , Y' = Y - {yn}
If case 1  holds, then induced alignment of X' and Y' as X' and Y' are of same size.
Similarly for case 2 and case 3.

Relevant sub problems 
Here sub problems get smaller, by plucking the characters from either 1st string or 2nd string or both. Hence we need to track how many characters we have plugged, and hence we need 2D array to keep track tof 2 things.


Optimal solution for penalty:
Let P denotes the Penalty here.
Lets consider the case 1 - xm and yn are aligned
Mij = Sij  + M(i-1)(j-1)

Case2 - xm is aligned with space / gap
Mij = w + M(i-1,j)

Case 3 - space / gap  is aligned with yn
Mij = w + M(i,j-1)

Mi,j = MAXIMUM[
     Mi-1, j-1 + Si,j (match/mismatch in the diagonal),
     Mi,j-1 + w (gap in sequence #1),
     Mi-1,j + w (gap in sequence #2)]
 
In case the penalties are such that we are given penalty Sij of mismatch is higher than penalty of match, then we have to use minimum function.
 
Three steps in dynamic programming
  1. Initialization
  2. Matrix fill (scoring)
  3. Traceback (alignment)
Initialization Step-edge cases
The first step in the global alignment dynamic programming approach is to create a matrix with M + 1 columns and N + 1 rows where M and N correspond to the size of the sequences to be aligned.
Since this example assumes there is no gap opening or gap extension penalty, the first row and first column of the matrix can be initially filled with 0.(Note in normal problem if Sgap is not 0, then use it some as i.gap)


Matrix Fill Step

One possible (inefficient) solution of the matrix fill step finds the maximum global alignment score by starting in the upper left hand corner in the matrix and finding the maximal score Mi,j for each position in the matrix. In order to find Mi,j for any i,j it is minimal to know the score for the matrix positions to the left, above and diagonal to i, j. In terms of matrix positions, it is necessary to know Mi-1,j, Mi,j-1 and Mi-1, j-1.
For each position, Mi,j is defined to be the maximum score at position i,j; i.e.
Mi,j = MAXIMUM[
     Mi-1, j-1 + Si,j (match/mismatch in the diagonal),
     Mi,j-1 + w (gap in sequence #1),
     Mi-1,j + w (gap in sequence #2)]
Note that in the example, Mi-1,j-1 will be red, Mi,j-1 will be green and Mi-1,j will be blue.
Using this information, the score at position 1,1 in the matrix can be calculated. Since the first residue in both sequences is a G, S1,1 = 1, and by the assumptions stated at the beginning, w = 0. Thus, M1,1 = MAX[M0,0 + 1, M1, 0 + 0, M0,1 + 0] = MAX [1, 0, 0] = 1.
A value of 1 is then placed in position 1,1 of the scoring matrix.
Since the gap penalty (w) is 0, the rest of row 1 and column 1 can be filled in with the value 1. Take the example of row 1. At column 2, the value is the max of 0 (for a mismatch), 0 (for a vertical gap) or 1 (horizontal gap). The rest of row 1 can be filled out similarly until we get to column 8. At this point, there is a G in both sequences (light blue). Thus, the value for the cell at row 1 column 8 is the maximum of 1 (for a match), 0 (for a vertical gap) or 1 (horizontal gap). The value will again be 1. The rest of row 1 and column 1 can be filled with 1 using the above reasoning.
Now let's look at column 2. The location at row 2 will be assigned the value of the maximum of 1(mismatch), 1(horizontal gap) or 1 (vertical gap). So its value is 1.
At the position column 2 row 3, there is an A in both sequences. Thus, its value will be the maximum of 2(match), 1 (horizontal gap), 1 (vertical gap) so its value is 2.
Moving along to position colum 2 row 4, its value will be the maximum of 1 (mismatch), 1 (horizontal gap), 2 (vertical gap) so its value is 2. Note that for all of the remaining positions except the last one in column 2, the choices for the value will be the exact same as in row 4 since there are no matches. The final row will contain the value 2 since it is the maximum of 2 (match), 1 (horizontal gap) and 2(vertical gap).
Using the same techniques as described for column 2, we can fill in column 3.
After filling in all of the values the score matrix is as follows:

Traceback Step

After the matrix fill step, the maximum alignment score for the two test sequences is 6. The traceback step determines the actual alignment(s) that result in the maximum score. Note that with a simple scoring algorithm such as one that is used here, there are likely to be multiple maximal alignments.
The traceback step begins in the M,J position in the matrix, i.e. the position that leads to the maximal score. In this case, there is a 6 in that location.


Traceback takes the current cell and looks to the neighbor cells that could be direct predacessors. This means it looks to the neighbor to the left (gap in sequence #2), the diagonal neighbor (match/mismatch), and the neighbor above it (gap in sequence #1). The algorithm for traceback chooses as the next cell in the sequence one of the possible predacessors. In this case, the neighbors are marked in red. They are all also equal to 5.


Since the current cell has a value of 6 and the scores are 1 for a match and 0 for anything else, the only possible predacessor is the diagonal match/mismatch neighbor. If more than one possible predacessor exists, any can be chosen. This gives us a current alignment of

    (Seq #1)      A 
                  |
    (Seq #2)      A
So now we look at the current cell and determine which cell is its direct predacessor. In this case, it is the cell with the red 5.

 


The alignment as described in the above step adds a gap to sequence #2, so the current alignment is

    (Seq #1)     T A
                   |
    (Seq #2)     _ A
Once again, the direct predacessor produces a gap in sequence #2.
After this step, the current alignment is

      (Seq #1)     T T A
                       |
                   _ _ A 
Continuing on with the traceback step, we eventually get to a position in column 0 row 0 which tells us that traceback is completed. One possible maximum alignment is :
Giving an alignment of :
          G A A T T C A G T T A
          |   |   | |   |     | 
          G G A _ T C _ G _ _ A 
 
An alternate solution is:
Giving an alignment of :
          G _ A A T T C A G T T A
          |     |   | |   |     | 
          G G _ A _ T C _ G _ _ A

 
There are more alternative solutions each resulting in a maximal global alignment score of 6. Since this is an exponential problem, most dynamic programming algorithms will only print out a single solution.

Time complexity

O(m n)
m = size of sequence X, and n = size of Y

 

0 comments:

Post a Comment