import numpy as np from random import randint DIAG UP LEFT def matrix

 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
import numpy as np
from random import randint
DIAG, UP, LEFT = 0, 1, 2
def matrix(seq1, seq2):
match, mismatch, indel = 1, -1, -1
n, m = len(seq1), len(seq2)
s, ptr = np.zeros((n+1, m+1)), np.zeros((n+1, m+1), dtype=int)
for i in range(1, n+1): s[i,0] = indel * i
for j in range(1, m+1): s[0,j] = indel * j
ptr[0,1:], ptr[1:,0] = LEFT, UP
for i in range(1,n+1):
for j in range(1,m+1):
if seq1[i-1] == seq2[j-1]: s[i,j], ptr[i,j] = s[i-1,j-1] + match, DIAG
else: s[i,j], ptr[i,j] = s[i-1,j-1] + mismatch, DIAG
if s[i-1, j] + indel > s[i,j]: s[i,j], ptr[i,j] = s[i-1,j] + indel, UP
if s[i, j-1] + indel > s[i,j]: s[i,j], ptr[i,j] = s[i, j-1] + indel, LEFT
return s, ptr
def trace(seq1, seq2, s, ptr) :
align1, align2, i, j, cur = '', '', len(seq1), len(seq2), lambda: ptr[i, j]
while i > 0 or j > 0:
if cur() == DIAG: align1, align2, i, j = seq1[i-1] + align1, seq2[j-1] + align2, i-1, j-1
elif cur() == LEFT: align1, align2, j = "-" + align1, seq2[j-1] + align2, j-1
elif cur() == UP: align1, align2, i = seq1[i-1] + align1, "-" + align2, i-1
return align1, align2
def needleman_alg(seq1, seq2) :
s, ptr = matrix(seq1, seq2)
alignment = trace(seq1, seq2, s, ptr)
return alignment, s[len(seq1), len(seq2)]
def random_DNA_sequence(length):
nucleotides = ['A', 'T', 'G', 'C']
return ''.join([nucleotides[randint(0,3)] for i in range(length)])
seq1 = random_DNA_sequence(20)
seq2 = random_DNA_sequence(20)
print needleman_alg(seq1, seq2)