global.py

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=0,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)
print nA, nB
print '*'
for i in xrange(1,nA):
    for j in xrange(1,nB):
         if dpV[j,i]<dpE[j,i]:
             continue
    
        #x=min(dpV[j,i],dpE[j,i],dpF[j,i])
print '*'

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)
print '*'           
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)

print '*' 

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 '*'
#traceback 

i=nA-1
j=nB-1
R=[]
S=[]
s=argmin([dpE[j,i],dpF[j,i],dpG[j,i]])

while i>=1 and j>=1:
            
    if s==2:
        R.append(A[i-1])
        S.append(B[j-1])
        #print '**',i,j,s,dpE[j,i],dpF[j,i],dpG[j,i],''.join(reversed(R)),''.join(reversed(S))
        i=i-1
        j=j-1
    
    elif s==0:
        R.append('-')
        S.append(B[j-1])
        #print '**',i,j,s,dpE[j,i],dpF[j,i],dpG[j,i],''.join(reversed(R)),''.join(reversed(S))
        j=j-1
        while  j>=1  and dpE[j,i]+gap(j)-gap(j-1)<dpV[j,i]+gap(1):
            S.append(B[j-1])
            R.append('-') 
            #print '*',i,j,s,dpE[j,i],dpF[j,i],dpG[j,i],''.join(reversed(R)),''.join(reversed(S))
            j=j-1
       
       
    elif s==1:
        S.append('-')
        R.append(A[i-1])
        print '**',i,j,s,dpE[j,i],dpF[j,i],dpG[j,i],''.join(reversed(R)),''.join(reversed(S))
        i=i-1
        while i>=1 and dpF[j,i]+gap(i)-gap(i-1)<dpV[j,i]+gap(1):
            S.append('-')
            R.append(A[i-1]) 
            #print '*',i,j,s,dpE[j,i],dpF[j,i],dpG[j,i],''.join(reversed(R)),''.join(reversed(S))
            i=i-1
        

   
    s=argmin([dpE[j,i],dpF[j,i],dpG[j,i]])


if i==0:
    while j>0:
        R.append('-')
        S.append(B[j-1])
        j=j-1

elif j==0:
    while i>0:
        S.append('-')
        R.append(A[i-1])
        i=i-1


Rx= ''.join(reversed(R))
Sx=''.join(reversed(S))

print -dpV[nB-1][nA-1]
print Rx
print Sx