Bioinformatics

Experiments with Python and Bioinformatics

global.py

@@ -0,0 +1,139 @@
+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
+
+
+