Implemented the Needleman-Wunsch algorithm

Author: alevchuk

pairwise_alignment.py renamed to alignment.py

@@ -26,57 +26,7 @@ def match_score(alpha, beta):
     else:
         return mismatch_penalty
 
-def water(seq1, seq2):
-    m, n = len(seq1), len(seq2)  # length of two sequences
-    
-    # Generate DP table and traceback path pointer matrix
-    score = zeros((m+1, n+1))      # the DP table
-    pointer = zeros((m+1, n+1))    # to store the traceback path
-    
-    P = 0
-    max_score = P        # initial maximum score in DP table
-    
-    # Calculate DP table and mark pointers
-    for i in range(1, m + 1):
-        for j in range(1, n + 1):
-            score_up = score[i-1][j] + gap_penalty
-            score_down = score[i][j-1] + gap_penalty
-            score_diagonal = score[i-1][j-1] + match_score(seq1[i-1], seq2[j-1])
-            score[i][j] = max(0,score_up, score_down, score_diagonal)
-            if score[i][j] == 0:
-                pointer[i][j] = 0 # 0 means end of the path
-            if score[i][j] == score_up:
-                pointer[i][j] = 1 # 1 means trace up
-            if score[i][j] == score_down:
-                pointer[i][j] = 2 # 2 means trace left
-            if score[i][j] == score_diagonal:
-                pointer[i][j] = 3 # 3 means trace diagonal
-            if score[i][j] >= max_score:
-                max_i = i
-                max_j = j
-                max_score = score[i][j];
-    
-    
-    align1, align2 = '', ''    # initial sequences
-    
-    i,j = max_i,max_j    # indices of path starting point
-    
-    #traceback, follow pointers
-    while pointer[i][j] != 0:
-        if pointer[i][j] == 3:
-            align1 = align1 + seq1[i-1]
-            align2 = align2 + seq2[j-1]
-            i = i - 1
-            j = j - 1
-        elif pointer[i][j] == 2:
-            align1 = align1 + '-'
-            align2 = align2 + seq2[j-1]
-            j = j - 1
-        elif pointer[i][j] == 1:
-            align1 = align1 + seq1[i-1]
-            align2 = align2 + '-'
-            i = i - 1
-    
+def finalize(align1, align2):
     align1 = align1[::-1]    #reverse sequence 1
     align2 = align2[::-1]    #reverse sequence 2
 
@@ -112,3 +62,107 @@ def water(seq1, seq2):
     print align1
     print symbol
     print align2
+
+
+def needle(seq1, seq2):
+    m, n = len(seq1), len(seq2)  # length of two sequences
+    
+    # Generate DP table and traceback path pointer matrix
+    score = zeros((m+1, n+1))      # the DP table
+   
+    # Calculate DP table
+    for i in range(0, m + 1):
+        score[i][0] = gap_penalty * i
+    for j in range(0, n + 1):
+        score[0][j] = gap_penalty * j
+    for i in range(1, m + 1):
+        for j in range(1, n + 1):
+            match = score[i - 1][j - 1] + match_score(seq1[i-1], seq2[j-1])
+            delete = score[i - 1][j] + gap_penalty
+            insert = score[i][j - 1] + gap_penalty
+            score[i][j] = max(match, delete, insert)
+
+    # Traceback and compute the alignment 
+    align1, align2 = '', ''
+    i,j = m,n # start from the bottom right cell
+    while i > 0 and j > 0: # end toching the top or the left edge
+        score_current = score[i][j]
+        score_diagonal = score[i-1][j-1]
+        score_up = score[i][j-1]
+        score_left = score[i-1][j]
+
+        if score_current == score_diagonal + match_score(seq1[i-1], seq2[j-1]):
+            align1 += seq1[i-1]
+            align2 += seq2[j-1]
+            i -= 1
+            j -= 1
+        elif score_current == score_left + gap_penalty:
+            align1 += seq1[i-1]
+            align2 += '-'
+            i -= 1
+        elif score_current == score_up + gap_penalty:
+            align1 += '-'
+            align2 += seq2[j-1]
+            j -= 1
+
+    # Finish tracing up to the top left cell
+    while i > 0:
+        align1 += seq1[i-1]
+        align2 += '-'
+        i -= 1
+    while j > 0:
+        align1 += '-'
+        align2 += seq2[j-1]
+        j -= 1
+
+    finalize(align1, align2)
+
+def water(seq1, seq2):
+    m, n = len(seq1), len(seq2)  # length of two sequences
+    
+    # Generate DP table and traceback path pointer matrix
+    score = zeros((m+1, n+1))      # the DP table
+    pointer = zeros((m+1, n+1))    # to store the traceback path
+    
+    max_score = 0        # initial maximum score in DP table
+    # Calculate DP table and mark pointers
+    for i in range(1, m + 1):
+        for j in range(1, n + 1):
+            score_diagonal = score[i-1][j-1] + match_score(seq1[i-1], seq2[j-1])
+            score_up = score[i][j-1] + gap_penalty
+            score_left = score[i-1][j] + gap_penalty
+            score[i][j] = max(0,score_left, score_up, score_diagonal)
+            if score[i][j] == 0:
+                pointer[i][j] = 0 # 0 means end of the path
+            if score[i][j] == score_left:
+                pointer[i][j] = 1 # 1 means trace up
+            if score[i][j] == score_up:
+                pointer[i][j] = 2 # 2 means trace left
+            if score[i][j] == score_diagonal:
+                pointer[i][j] = 3 # 3 means trace diagonal
+            if score[i][j] >= max_score:
+                max_i = i
+                max_j = j
+                max_score = score[i][j];
+    
+    align1, align2 = '', ''    # initial sequences
+    
+    i,j = max_i,max_j    # indices of path starting point
+    
+    #traceback, follow pointers
+    while pointer[i][j] != 0:
+        if pointer[i][j] == 3:
+            align1 += seq1[i-1]
+            align2 += seq2[j-1]
+            i -= 1
+            j -= 1
+        elif pointer[i][j] == 2:
+            align1 += '-'
+            align2 += seq2[j-1]
+            j -= 1
+        elif pointer[i][j] == 1:
+            align1 += seq1[i-1]
+            align2 += '-'
+            i -= 1
+
+    finalize(align1, align2)

needle

@@ -0,0 +1,15 @@
+#!/usr/bin/env python
+
+import sys, string
+import alignment
+
+# Command-line arguments
+f1 = open(sys.argv[1], 'r')
+seq1 = f1.readline()
+seq1 = string.strip(seq1)
+
+f2 = open(sys.argv[2], 'r')
+seq2 = f2.readline()
+seq2 = string.strip(seq2)
+
+alignment.needle(seq1, seq2)

needle-test

@@ -0,0 +1,6 @@
+#!/bin/bash
+
+./needle <(echo ALEXNILDCQGAAA%LIN@%) <(echo PYTHONNILDCQGHHHH%LINJJBB)
+echo
+
+./needle <(echo ONNALEXNILDCQGAAA%LIN@%BBSDF) <(echo PYTHONNILDCQGHHHH%LINJJBB)

water

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 
 import sys, string
-import pairwise_alignment
+import alignment
 
 # Command-line arguments
 f1 = open(sys.argv[1], 'r')
@@ -12,4 +12,4 @@ f2 = open(sys.argv[2], 'r')
 seq2 = f2.readline()
 seq2 = string.strip(seq2)
 
-pairwise_alignment.water(seq1, seq2)
+alignment.water(seq1, seq2)