Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
#!/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"))