Initial

alignementseq.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Jan 15 15:02:16 2019
+
+@author: ariane.delrocq
+"""
+
+import Bio.SubsMat.MatrixInfo
+
+mat = Bio.SubsMat.MatrixInfo.blosum62
+
+
+def score(a, b, gap=-8, id=1, mut=-1):
+    if not a or not b:
+        return gap
+    else:
+        a = a.upper()
+        b = b.upper()
+        try:
+            return mat[a,b]
+        except KeyError:
+            return mat[b,a]
+    
+#    elif a==b:
+#        return id
+#    else:
+#        return mut
+
+def align(x,y):
+    n=len(x)
+    m=len(y)
+    mat = [[0 for _ in range(m+1)] for _ in range(n+1)]   # scores matrix
+    path = [[0 for _ in range(m+1)] for _ in range(n+1)]   # last step for optimal score
+    
+    # first line and col
+    for i in range(1, n+1):
+        mat[i][0] = mat[i-1][0] + score('', x[i-1])
+        path[i][0] = (i-1, 0)
+    for j in range(1, m+1):
+        mat[0][j] = mat[j-1][0] + score('', y[j-1])
+        path[0][j] = (0, j-1)
+    
+    # fill table
+    for i in range(1, n+1):
+        for j in range(1, m+1):
+            a = x[i-1]
+            b = y[j-1]
+            s1 = mat[i-1][j-1] + score(a,b)
+            s2 = mat[i-1][j] + score('', a)
+            s3 = mat[i][j-1] + score('', b)
+            if s1 >= s2:
+                if s1 >= s3:
+                    mat[i][j] = s1
+                    path[i][j] = (i-1, j-1)
+                else:
+                    mat[i][j] = s3
+                    path[i][j] = (i, j-1)
+            else:   # s2 > s1
+                if s2 >= s3:
+                    mat[i][j] = s2
+                    path[i][j] = (i-1, j)
+                else:
+                    mat[i][j] = s3
+                    path[i][j] = (i, j-1)
+    
+    bestScore = mat[n][m]
+    alx = ""
+    aly = ""
+    i, j = n, m
+    while i > 0 and j > 0:
+        i2, j2 = path[i][j]
+        if i==i2:
+            alx+='-'
+        else:
+            alx+=x[i-1]
+        if j==j2:
+            aly+='-'
+        else:
+            aly+=y[j-1]
+        i,j = i2, j2
+    while i>0:
+        aly+='-'
+        alx+=x[i-1]
+        i -=1
+    while j>0:
+        alx+='-'
+        aly+=y[j-1]
+        j -=1
+    alx = alx[::-1]
+    aly = aly[::-1]
+    return bestScore, alx, aly
+
+print(align("chat", "cat"))
+print(align("chat", "cgat"))
+print(align("chat", "at"))