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
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):
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)
Three steps in dynamic programming
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)
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.
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:
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
The alignment as described in the above step adds a gap to sequence #2, so the current alignment is
After this step, the current alignment is
Giving an alignment of :
Giving an alignment of :
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
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)
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):
- xm and yn are aligned
- xm is aligned with space / gap
- space / gap is aligned with yn
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
- Initialization
- Matrix fill (scoring)
- Traceback (alignment)
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.
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).
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.
(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.
(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 :
G A A T T C A G T T A
| | | | | |
G G A _ T C _ G _ _ A
An alternate solution is:
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
