Needleman Wunsch Algorithm in C++
The Needleman-Wunsch algorithm is used in the field of bioinformatics to align proteins or sequences of nucleotides. The Needleman-Wunsch algorithm was the first application of dynamic programming to compare biological sequences. Saul B. Needleman and Christian D. Wunsch gave the Needleman-Wunsch algorithm, published in 1970. This algorithm divides the large problem into smaller subproblems, and the solution of subproblems is used to find a solution that is optimal for a more significant problem. Sometimes, it is also called the optimal matching algorithm and the global alignment technique. This algorithm is still used for optimal global alignment, especially when the quality of the worldwide alignment is essential.
Consider an example of two DNA sequences: GCATGCG and GATTACA.
Construct a table as shown in the figure. Start the first string at the top of the third column and the other at the start of the third row.
G | C | A | T | G | C | G | ||
G | ||||||||
A | ||||||||
T | ||||||||
T | ||||||||
A | ||||||||
C | ||||||||
A |
The next step is to decide the score of each pair of letters. For the above example, one possible alignment can:
12345678 GCATG-CG G-ATTACA
These letters may be matched, mismatched, or matched to a gap.
- Match: The two letters are said to match if the current index is the same.
- Mismatch: The two letters are said to be mismatched if the current index is different.
- Indel (Insertion and Deletion): The best alignment consists of one letter aligning to a gap in the other string.
Each situation is assigned a score, and the total score of the pairings is the score of the whole alignment candidate.
- Match: +1
- Mismatch or Indel: -1
The score would be zero for the above alignment:
GCATG-CG G-ATTACA +-++--+- à 1 * 4 + (- 1) * 4 = 0
Now start with a zero in the second row and second columns and move with cells row by row, calculate the score of each cell. The calculation of the score is by comparing the scores of the cells neighbouring the left cell, top or top-left.
- The path from the top or left cell represents an indel pairing, so take the left and top cells' scores and add the indel score to each.
- The diagonal path represents a match/mismatch, so take the score of the top-left diagonal cell and add the score for a match if the corresponding bases (letters) in the row and column match or the score for mismatch if they do not.
G | C | A | T | G | C | G | ||
0 | −1 | −2 | −3 | −4 | −5 | −6 | −7 | |
G | −1 | |||||||
A | −2 | |||||||
T | −3 | |||||||
T | −4 | |||||||
A | −5 | |||||||
C | −6 | |||||||
A | −7 |
The first case of the Needleman Wunsch Algorithm with existing that are in all 3 directions. It is the intersection of the first letters, in this case, G and G.
G | ||
0 | −1 | |
G | −1 | X |
There are three possible sums:
- The score of the top-left neighbour is score 0. G and G have matching pairs, so 1 is added to the score, and then the score is 0 + 1.
- The score of the top neighbour is – 1, and moving from there, that shows as indel (insertion and deletion), so add a score for insertion and deletion (-1) + (-1) = (-2).
- The score of the left neighbour has a score of – 1, which shows an insertion and deletion that also gives a result of (-2).
The highest candidate is 1 for 1,-2, and -2, so x is replaced by 1 in the cell.
G | ||
0 | −1 | |
G | −1 | 1 |
The cell that gives the result of the highest candidate score must be corrected. In the further example, the diagonal moves for X and Y show a mismatch.
G | C | ||
0 | −1 | −2 | |
G | −1 | 1 | X |
A | −2 | Y |
X:
- Top: (−2) + (−1) = (−3)
- Left: (+1) + (−1) = (0)
- Top-Left: (−1) + (−1) = (−2)
Y:
- Top: (1) + (−1) = (0)
- Left: (−2) + (−1) = (−3)
- Top-Left: (−1) + (−1) = (−2)
For X and Y, the highest score is 0:
G | C | ||
0 | −1 | −2 | |
G | −1 | 1 | 0 |
A | −2 | 0 |
T | G | |
T | 1 | 1 |
A | 0 | X |
- Top: (1) + (−1) = (0)
- Top-Left: (1) + (−1) = (0)
- Left: (0) + (−1) = (−1)
For this case, all directions reaching the highest score must be noted as possible origin cells in the finished diagram in Figure 1, for example, the cell in row and column 7.
In this way, the table is filled with scores of all the possible candidates that are aligned; the score on the bottom suitable cell shows the alignment score for the alignment that is best.
Tracing Arrows Back to Origin:
Select a path form the bottom suitable back cell to the top left cell following the direction of the arrows. From that path, a sequence is made with the following rules.
A diagonal arrow shows a match or mismatch, so the letter of the columns and the row letter of the original cell align.
A horizontal or vertical arrow represents an insertion and deletion (indel). Vertical arrows will align a gap ("- ") to the row letter (the "side" sequence), and horizontal arrows align the gap to the letter of the column that is the top sequence.
By following these steps,
G → CG → GCG → -GCG → T-GCG → AT-GCG → CAT-GCG → GCAT-GCG
A → CA → ACA → TACA → TTACA → ATTACA → -ATTACA → G-ATTACA
↓
(branch) → TGCG → -TGCG → ...
→ TACA → TTACA → ...
System of Scoring:
The scoring system that is simplest is that it gives a value for every match, mismatch and indel (insertion or deletion). The score for the match is 1; mismatch is -1, and indel (insertion or deletion) is – 1. Thus, the lower the alignment, the larger the score, and the larger the distance is edited; for this system of scoring, one expects a high score. Another scoring system could be:
- Match = 1
- Indel (insertion or deletion) = - 2
- Mismatch = - 1
For the system, the alignment score will show the distance that is edited between the two strings. A different systems of scoring can be devised for various situations; for example, if gaps are assumed to be very bad for the alignment, one may use a different system of scoring that penalizes the gaps heavily, such as
- Match = 1
- Mismatch = - 1
- Indel (insertion or deletion) = - 2
Similarity Matrix:
A more complex system of scoring feature values not only for the type of alternation but also for the letters that are included. For example, a match between A and A may be given 1, but if T and T are a match, it can be a score of 4. To show all the possible combinations of letter and their score, we used a similarity matrix that is represented as.
A | G | C | T | |
A | 1 | −1 | −1 | −1 |
G | −1 | 1 | −1 | −1 |
C | −1 | −1 | 1 | −1 |
T | −1 | −1 | −1 | 1 |
Each row in the similarity matrix shows a switch from one of the letters that the cell matches. Hence, this shows all the matches and mismatches that are possible. It should be noted that all the matches will go along the diagonal. Also, not all the tables need to be filled. If the T – T = 4 rule is to be implemented for the above similarity matrix:
A | G | C | T | |
A | 1 | −1 | −1 | −1 |
G | −1 | 1 | −1 | −1 |
C | −1 | −1 | 1 | −1 |
T | −1 | −1 | −1 | 4 |
Other matrices used for scoring have been constructed statistically, which give weights to actions that are different and appropriate to a particular situation. It is important to have weighted scoring matrices in alignment with the protein sequence due to the changing frequency of the various amino acids. Two large families of matrices used for scoring are PAM and BLOSUM.
Gap Penalty:
During the alignment of sequences, there are gaps, i.e., insertion or deletions (indels); sometimes, they can be large. In terms of biology, a large gap is more likely due to a large deletion as to oppose multiple single deletions. The score of two indels (insertion or deletion) can be worse than one large indel (insertion or deletion). The common way to do this is via a large gap-start score for each letter, which enhances the insertion or deletion. For instance, new-indel can cost – 5, and extend-indel may cost – 1. Doing this, an alignment such as
GAAAAAAT G - - A – A – T
As many equal alignments, some with small alignments will now align as:
GAAAAAAT GAA - - - - T
Or any other that may have 4 long gaps.
Advanced Presentation of the Algorithm:
The aligned characters' scores are given by a similarity matrix. Here, S (a, b) is the character similarity of a and b. Linear gap penalty is used by it and called as d.
For instance, if the matrix of similarity was
A | G | C | T | |
A | 10 | −1 | −3 | −4 |
G | −1 | 7 | −5 | −3 |
C | −3 | −5 | 9 | 0 |
T | −4 | −3 | 0 | 8 |
Then, the alignment will be:
AGACTAGTTAC CGA - - - GACGT
Having a gap penalty of – -5, that would have the score given as
S (A, C) + S (G, G) + S (A, A) + (3 × d) + S (G, G) + S (T, A) + S (T, C) + S (A, G) + S (C, T) = −3 + 7 + 10 − (3 × 5) + 7 + (−4) + 0 + (−1) + 0 = 1
For finding the alignment that has the highest score, a two-dimensional array is used. The ith row and jth column are denoted as Fij. In sequence A, there is only one row for each character, and in sequence B, there is only one column for each character. To align a sequence of size n and m, the memory used is O (nm). Hirschberg's algorithm only stores the subset of the array in memory and uses O (min {n, m}) space.
As the algorithm moves further, the Fij will be updated to an optimal score for the alignment of the first i = 0, …, n characters in A and B, and the first j = 0, …, m characters. The optimality principle is given as follows:
- Base condition
F0j = d * j
Fi0 = d * i
- The recursion equation is given based on the optimality principle:
Fij = max (Fi – 1, j – 1 + S (Ai, Bj), Fi, j-1 + d, Fi-1,j + d)
The algorithm for computing the F matrix is given as follows:
d ← Gap penalty score for i = 0 to length(A) F(i,0) ← d * i for j = 0 to length(B) F(0,j) ← d * j for i = 1 to length(A) for j = 1 to length(B) { Match ← F(i−1, j−1) + S(Ai, Bj) Delete ← F(i−1, j) + d Insert ← F(i, j−1) + d F(i,j) ← max(Match, Insert, Delete) }
Once the F matrix is calculated, the entry Fnm gives the output of the maximum score among all the alignments that are possible. To evaluate the alignment that actually gives the score, start from the bottom right cell and make a comparison of the value with the three possible sources (Match, Insert, and Delete) to see which it came from. If there is a match, then Ai and Bj are aligned. If there is a delete, then the gap is aligned with Ai, and if insert, then Bj is aligned with a gap.
Implementation of Needleman Wunsch Algorithm in C++:
#includeusing namespace std; void getMinimumPenalty(string x, string y, int pxy, int pgap) { int i, j; int m = x.length(); int n = y.length(); int dp[n+m+1][n+m+1] = {0}; for (i = 0; i <= (n+m); i++) { dp[i][0] = i * pgap; dp[0][i] = i * pgap; } for (i = 1; i <= m; i++) { for (j = 1; j <= n; j++) { if (x[i - 1] == y[j - 1]) { dp[i][j] = dp[i - 1][j - 1]; } else { dp[i][j] = min({dp[i - 1][j - 1] + pxy , dp[i - 1][j] + pgap , dp[i][j - 1] + pgap }); } } } int l = n + m; i = m; j = n; int xpos = l; int ypos = l; int xans[l+1], yans[l+1]; while ( !(i == 0 || j == 0)) { if (x[i - 1] == y[j - 1]) { xans[xpos--] = (int)x[i - 1]; yans[ypos--] = (int)y[j - 1]; i--; j--; } else if (dp[i - 1][j - 1] + pxy == dp[i][j]) { xans[xpos--] = (int)x[i - 1]; yans[ypos--] = (int)y[j - 1]; i--; j--; } else if (dp[i - 1][j] + pgap == dp[i][j]) { xans[xpos--] = (int)x[i - 1]; yans[ypos--] = (int)'_'; i--; } else if (dp[i][j - 1] + pgap == dp[i][j]) { xans[xpos--] = (int)'_'; yans[ypos--] = (int)y[j - 1]; j--; } } while (xpos > 0) { if (i > 0) xans[xpos--] = (int)x[--i]; else xans[xpos--] = (int)'_'; } while (ypos > 0) { if (j > 0) yans[ypos--] = (int)y[--j]; else yans[ypos--] = (int)'_'; } int id = 1; for (i = l; i >= 1; i--) { if ((char)yans[i] == '_' && (char)xans[i] == '_') { id = i + 1; break; } } cout << "Minimum Penalty in aligning the genes = "; cout << dp[m][n] << " "; cout << "The aligned genes are : "; for (i = id; i <= l; i++) { cout<<(char)xans[i]; } cout << " "; for (i = id; i <= l; i++) { cout << (char)yans[i]; } return; } int main(){ string gene1 = "AGGGCT"; string gene2 = "AGGCA"; int misMatchPenalty = 3; int gapPenalty = 2; getMinimumPenalty(gene1, gene2, misMatchPenalty, gapPenalty); return 0; }
Output:
Minimum Penalty in aligning the genes = 5 The aligned genes are : AGGGCT A_GGCA
Complexity of Needleman Wunsch Algorithm:
Space Complexity: O (n * m)
Time Complexity: O (n * m)