cmp_nw/needleman-wunsch.py

   1 #!/usr/bin/python
   2
   3 from numpy import *
   4 import sys
   5
   6 def NeedlemanWunsch(seqA, seqB, match, mismatch, gap):
   7         a = len(seqA)
   8         b = len(seqB)
   9         F = zeros((a+1,b+1))
  10         for i in range(0, a):
  11                 F[i,0] = gap * i
  12         for j in range(0, b):
  13                 F[0,j] = gap * j
  14
  15         for i in range(1, a+1):
  16                 for j in range(1, b+1):
  17                         F[i,j] = -10000
  18
  19                         if (seqA[i-1] == seqB[j-1]):
  20                                 chscore = match
  21                         else:
  22                                 chscore = mismatch
  23                         if (F[i-1,j-1] + chscore > F[i,j]):
  24                                 F[i,j] = F[i-1,j-1] + chscore
  25
  26                         if (F[i-1,j] + gap > F[i,j]):
  27                                 F[i,j] = F[i-1,j] + gap
  28                         if (F[i,j-1] + gap > F[i,j]):
  29                                 F[i,j] = F[i,j-1] + gap
  30                         # print i, j, F[i,j]
  31         # print "** ", a, b, F[a,b]
  32         return F[a,b]
  33
  34 print NeedlemanWunsch(sys.argv[1],sys.argv[2],1,0,-1)