Attempt to fix the second algorithm.

.gitignore

+*.swp
+balibase/

alignementseq.py

@@ -113,73 +113,102 @@ def align2steps(seq1, seq2, d=8, e=4):
     mm = [[0 for _ in range(m+1)] for _ in range(n+1)]   # scores matrix
     mx = [[0 for _ in range(m+1)] for _ in range(n+1)]
     my = [[0 for _ in range(m+1)] for _ in range(n+1)]
-    path = [[0 for _ in range(m+1)] for _ in range(n+1)]   # last step for optimal score
+    pathm = [[0 for _ in range(m+1)] for _ in range(n+1)]   # last step for optimal score
+    pathx = [[0 for _ in range(m+1)] for _ in range(n+1)]
+    pathy = [[0 for _ in range(m+1)] for _ in range(n+1)]
+    
+    lower_bound = - (n+m) * max(e,d)
     
     # first line and col
     for i in range(1, n+1):
-        mx[i][1] = my[i][1] = -d -(i-1)*e
-        path[i][0] = (i-1, 0)
+        mx[i][0] = -d -(i-1)*e
+        pathx[i][0] = pathx
+        # Best alignment ending with no gap but with 0 letters of one seq does
+        # not exist... but is still used in the formulas below
+        mm[i][0] = lower_bound
+        # Same for best alignment with 0 letter of Y ending on a gap on the X
+        # side with a letter of Y
+        my[i][0] = lower_bound
+        
     for j in range(1, m+1):
-        mx[1][j] = my[1][i] = -d - (j-1)*e
-        path[0][j] = (0, j-1)
+        my[0][j] = -d - (j-1)*e
+        pathy[0][j] = pathy
+        mm[0][j] = lower_bound
+        mx[0][j] = lower_bound
     
     # fill table
-    for i in range(2, n+1):
-        for j in range(2, m+1):
+    for i in range(1, n+1):
+        for j in range(1, m+1):
             a = x[i-1]
             b = y[j-1]
             
             # find max for M
             s1 = mm[i-1][j-1] + score(a,b)
             s2 = mx[i-1][j-1] + score(a,b)
-            s3 = my[i][j-1] + score(a,b)
+            s3 = my[i-1][j-1] + score(a,b)
             
             if s1 >= s2:
                 if s1 >= s3:
                     mm[i][j] = s1
-                    path[i][j] = (i-1, j-1)
+                    pathm[i][j] = pathm
                 else:
                     mm[i][j] = s3
-                    path[i][j] = (i, j-1)
+                    pathm[i][j] = pathy
             else:   # s2 > s1
                 if s2 >= s3:
                     mm[i][j] = s2
-                    path[i][j] = (i-1, j)
+                    pathm[i][j] = pathx
                 else:
                     mm[i][j] = s3
-                    path[i][j] = (i, j-1)
+                    pathm[i][j] = pathy
             
             # find max for I_x
             s4 = mm[i-1][j] - d
             s5 = mx[i-1][j] - e
             if s4 >= s5:
                 mx[i][j] = s4
+                pathx[i][j] = pathm
             else:
                 mx[i][j] = s5
+                pathx[i][j] = pathx
             
             # find max for I_y
             s6 = mm[i][j-1] - d
             s7 = my[i][j-1] - e
             if s6 >= s7:
                 my[i][j] = s6
+                pathy[i][j] = pathm
             else:
                 my[i][j] = s7
-    
-    bestScore = mm[n][m]
+                pathy[i][j] = pathy
+
+
+    bestScore, bestPath = max(
+            [(mm,pathm),(mx,pathx),(my,pathy)],
+            key=(lambda t:t[0][n][m]))
+    bestScore = bestScore[n][m]
     alx = ""
     aly = ""
     i, j = n, m
+    path = bestPath
     while i > 0 and j > 0:
-        i2, j2 = path[i][j]
-        if i==i2:
-            alx+='-'
-        else:
-            alx+=x[i-1]
-        if j==j2:
-            aly+='-'
+        prev_path = path[i][j]
+        if path == pathm:
+            i -= 1
+            j -= 1
+            alx += x[i]
+            aly += y[j]
+        elif path == pathy:
+            j -= 1
+            alx += '-'
+            aly += y[j]
         else:
-            aly+=y[j-1]
-        i,j = i2, j2
+            i -= 1
+            alx += x[i]
+            aly += '-'
+        path = prev_path
+
+
     while i>0:
         aly+='-'
         alx+=x[i-1]
@@ -198,8 +227,10 @@ def test():
 #    print(align("chat", "cat"))
 #    print(align("chat", "cgat"))
 #    print(align("chat", "at"))
-    s, x, y = align(seqs[0], seqs[1])
-    print(s,x,y)
+    s, x, y = align2steps(seqs[0], seqs[1])
+    print(s)
+    print(x.seq)
+    print(y.seq)
     with open("balibase/us/1.fasta", "w") as fd:
         Bio.SeqIO.write((x, y), fd, "fasta")
 test()