reeps

Sequence alignment

May 18th, 2018
166
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 1.46 KB | None | 0 0
  1. import numpy
  2. import scipy
  3. from string import *
  4.  
  5. #return value of each variation
  6. def Diagonal(n1, n2, pt):
  7.     if (n1 == n2):
  8.         return pt['MATCH']
  9.     else:
  10.         return pt['MISMATCH']
  11. #the function gets the optional elements of the alignment matrix and returns pointer's elements
  12. def Pointers( di, ho, ve):
  13.     pointer = max(di, ho, ve)
  14.     if (di == pointer):
  15.         return 'D'
  16.     elif (ho == pointer):
  17.         return 'H'
  18.     else:
  19.         return 'V'
  20.  
  21. def NW(s1, s2):
  22.     match =1   
  23.     mismatch = -1
  24.     gap = -2
  25.     penalty = { 'MATCH': match, 'MISMATCH' : mismatch, 'GAP': gap}
  26.     n = len(s1) + 1 #dimension of the matrix columns
  27.     m = len(s2) + 1  #dimension of the matrix rows
  28.     al_mat = numpy.zeros((m, n), dtype = int) # alignments matrix with zeros
  29.     p_mat = numpy.zeros((m, n), dtype = str) # alignments matrix with zeros
  30. # checking on gaps
  31.     for i in range(m):
  32.         al_mat[i][0] = penalty['GAP'] * i
  33.         p_mat[i][0] = 'V'
  34.     for j in range(n):
  35.         al_mat[0][j] = penalty['GAP'] * j
  36.         p_mat[0][j] = 'H'
  37.     p_mat[0][0] = 0 #return the first element of pointer matrix back to 0
  38.     for i in range(1, m):
  39.         for j in range(1, n):
  40.             di = al_mat[i - 1][j - 1] + Diagonal(s1[j - 1], s2[i - 1], penalty)
  41.             #value of matching on diagonal
  42.             ho = al_mat[i][j - 1] + penalty['GAP']
  43.             #horizontal
  44.             ve = al_mat[i - 1][j] + penalty['GAP']
  45.             #vertical
  46.             al_mat[i][j] = max(di, ho, ve)
  47.             p_mat[i][j] = Pointers(di, ho, ve)
  48.     print numpy.matrix(al_mat)
  49.     print numpy.matrix(p_mat)
  50.  
  51. NW('TCGCA', 'TCCA')
Advertisement
Add Comment
Please, Sign In to add comment