-
Notifications
You must be signed in to change notification settings - Fork 0
/
Modified_Needleman_Wunch.py
70 lines (61 loc) · 1.81 KB
/
Modified_Needleman_Wunch.py
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
import pandas as pd
import numpy as np
gap_penalty=-1
match_score=1
mismatch_score=0
def initialize(r,c):
mat=np.zeros((r+1,c+1))
return(mat)
def match(seq1,seq2):
mat=np.ones((len(seq1),len(seq2)))
for i in range(len(seq1)):
for j in range(len(seq2)):
if(seq1[i]==seq2[j]):
mat[i,j]=match_score
else:
mat[i,j]=mismatch_score
return(mat)
def backtrack(mat,seq1,seq2,match_table):
a_a=""
a_b=""
i=len(seq1)
j=len(seq2)
while(i>0 or j>0):
if(i>0 and j>0 and mat[i,j]==mat[i-1,j-1]+match_table[i-1,j-1]):
a_a = seq1[i-1] + a_a
a_b = seq2[j-1] + a_b
i = i-1
j = j-1
elif(i>0 and mat[i,j]==mat[i-1,j]+gap_penalty):
a_a = seq1[i-1] + a_a
a_b = '_' + a_b
i =i-1
else:
a_a = '_' + a_a
a_b = seq2[j-1] + a_b
j =j-1
return(a_a,a_b)
def matrix(mat,match_table):
for i in range(1,mat.shape[0]):
for j in range(1,mat.shape[1]):
if(i==mat.shape[0]-1 or j==mat.shape[1]-1):
a=mat[i-1,j]
b=mat[i,j-1]
else:
a=mat[i-1,j]+gap_penalty
b=mat[i,j-1]+gap_penalty
c=mat[i-1,j-1]+match_table[i-1,j-1]
mat[i,j]=max(a,b,c)
return(mat)
def improve_n_w(seq1,seq2):
seq1=seq1.lower()
seq2=seq2.lower()
mat=initialize(len(seq1),len(seq2))
match_table=match(seq1,seq2)
mat=matrix(mat,match_table)
print(mat)
alignment=backtrack(mat,seq1,seq2,match_table)
return(alignment)
print("Result is : ")
print("Input is ('acactg','acactgatcg')")
print(improve_n_w('acactg','acactgatcg'))