Skip to content
Permalink
Branch: master
Find file Copy path
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
65 lines (48 sloc) 1.42 KB
from blosum62 import *
import numpy as np
from numpy import argmin
from fasta import readFASTA
# this is an implementation of the needleman-wunsch algorithm
# and it should handle the affine gap situation
def mismatch_penalty(A,B):
if A!=B:
return -blosum62[A+B]
else:
return -blosum62[A+B]
s=readFASTA("rosalind_gcon.txt")
A=s.values()[0]
B=s.values()[1]
#A='MMM'
#B='YYYYYYY'
nA=len(A)+1
nB=len(B)+1
def gap(x,start=10,step=1):
return start + step*x
dpV=np.zeros((nB,nA),dtype=float)
dpG=np.zeros((nB,nA),dtype=float)
dpE=np.zeros((nB,nA),dtype=float)
dpF=np.zeros((nB,nA),dtype=float)
e=gap(1)
f=gap(2)-gap(1)
i=1
j=1
dpG[j,i]=mismatch_penalty(B[0],A[0])
dpE[j,i]=gap(1)+gap(1)
dpF[j,i]=gap(1)+gap(1)
i=1
for j in xrange(2,nB):
dpG[j,i]=mismatch_penalty(B[j-1],A[0])+gap(j-1)
dpE[j,i]=min(dpE[j-1,i]+f,dpV[j-1,i]+e)
dpF[j,i]=gap(j)+gap(1)
j=1
for i in xrange(2,nA):
dpG[j,i]=mismatch_penalty(B[0],A[i-1])+gap(i-1)
dpE[j,i]=gap(i)+gap(1)
dpF[j,i]=min(dpF[j,i-1]+f,dpV[j,i-1]+e)
for i in xrange(2,nA):
for j in xrange(2,nB):
dpG[j,i]=dpV[j-1,i-1]+mismatch_penalty(B[j-1],A[i-1])
dpE[j,i]=min(dpE[j-1,i]+f,dpV[j-1,i]+e)
dpF[j,i]=min(dpF[j,i-1]+f,dpV[j,i-1]+e)
dpV[j,i]=min(dpE[j,i],dpF[j,i],dpG[j,i])
print dpV[nB-1,nA-1]
You can’t perform that action at this time.