Advertisement
CasualGamer

Smith Waterman

Jan 5th, 2020
227
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 4.70 KB | None | 0 0
  1. import numpy as np
  2.  
  3. from cell import cell
  4.  
  5. def alignment_nw(a,b):
  6.     if a==b:
  7.         return 3
  8.     else:
  9.         return -3
  10.  
  11. def recursion(paths,D):
  12.     new_paths = []
  13.     returnvalue = paths
  14.     keep_going = False
  15.     for path in paths:
  16.         if path[-1].get_value() == 0:
  17.             continue
  18.         pointers = path[-1].get_pointers()
  19.         for pointer in pointers:
  20.             if pointers.size==0:
  21.                 continue
  22.             indices = np.where(D == pointer)
  23.             path = np.append(path,D[indices])
  24.             keep_going = True
  25.             if len(new_paths) == 0:
  26.                 new_paths = [path]
  27.             else:
  28.                 new_paths.append(path)
  29.     if keep_going:
  30.         returnvalue=recursion(new_paths,D)
  31.     return returnvalue
  32.  
  33. def main(A,B):
  34.     A = np.array([ a for a in A ])
  35.     B = np.array([ b for b in B ])
  36.     A = np.insert(A,0,'-')
  37.     B = np.insert(B,0,'-')
  38.     D = np.array([[cell() for i in range(B.size)] for j in range(A.size)])
  39.     for i in range(1,A.size):
  40.         D[i][0].set_value(0)
  41.     for i in range(1,B.size):
  42.         D[0][i].set_value(0)
  43.     for i in range(1,A.size):
  44.         for j in range(1,B.size):
  45.             D[i][j].set_value(0)
  46.             candidates = [D[i-1][j-1].get_value()+alignment_nw(A[i],B[j]), D[i-1][j].get_value()-2,D[i][j-1].get_value()-2]
  47.             index_array = np.argwhere(candidates == np.amax(candidates)).flatten().tolist()
  48.             if candidates[index_array[-1]]>0:
  49.                 D[i][j].set_value(candidates[index_array[-1]])
  50.             for k in index_array:
  51.                 if k == 0:
  52.                     D[i][j].add_pointer(D[i-1][j-1])
  53.                 elif k == 1:
  54.                     D[i][j].add_pointer(D[i-1][j])
  55.                 elif k == 2:
  56.                     D[i][j].add_pointer(D[i][j-1])
  57.     print(D)
  58.     #Tracing arrows back
  59.     path = [[np.amax(D)]]
  60.     paths = recursion(path,D)
  61.     print(paths)
  62.  
  63.     previous_step_indices = np.array([-1,-1])
  64.     alignment = ['','']
  65.     score = 0
  66.     for id,path in enumerate(paths):
  67.         for step in path:
  68.             indices = np.where(D == step)
  69.             if(previous_step_indices[0]==-1):
  70.                 previous_step_indices=indices
  71.                 continue
  72.             elif indices[0]<previous_step_indices[0] and indices[1]<previous_step_indices[1]:
  73.                 alignment = [A[previous_step_indices[0][-1]]+alignment[0],B[previous_step_indices[1][-1]]+alignment[1]]
  74.                 if A[previous_step_indices[0][-1]]==B[previous_step_indices[1][-1]]:
  75.                     score=score+1
  76.                 else:
  77.                     score=score-1
  78.             elif indices[0]<previous_step_indices[0]:
  79.                 alignment = [A[previous_step_indices[0][-1]]+alignment[0],'-'+alignment[1]]
  80.                 score=score-1
  81.             elif indices[1]<previous_step_indices[1]:
  82.                 alignment = ['-'+alignment[0],B[previous_step_indices[1][-1]]+alignment[1]]# np.char.add(['-',B[indices[1]]],alignment)#('-'+alignment[0],B[indices[1]]+alignment[1])
  83.                 score=score-1
  84.             previous_step_indices=indices
  85.         print('\npath number: {} score: {}\n{}\n{}\n'.format(id+1,score,alignment[0],alignment[1]))
  86.         alignment = ['','']
  87.         score=0
  88.     return
  89.  
  90. main('ACGTC','AGTCA')
  91. #main('GCATGCU','GATTACA')
  92. #main('GGTTGACTA','TGTTACGG')
  93.  
  94.  
  95.  
  96. ###################BACKUP#######################
  97.  
  98. #def main(A,B):
  99. #    A = np.array([ a for a in A ])
  100. #    B = np.array([ b for b in B ])
  101. #    A = np.insert(A,0,'-')
  102. #    B = np.insert(B,0,'-')
  103. #    cell = np.dtype([('symbol', np.unicode_, 3), ('value', np.int32)])
  104. #    D = np.zeros((A.size,B.size),dtype=cell)#([('H', 1), ('V', 0)], dtype=cell)
  105. #    for i in range(1,A.size):
  106. #        D[i][0]=('',-i)
  107. #    for i in range(1,B.size):
  108. #        D[0][i]=('',-i)
  109. #    for i in range(1,A.size):
  110. #        for j in range(1,B.size):
  111. #            candidates = [D[i-1][j-1]['value']+alignment_nw(A[i],B[j]), D[i-1][j]['value']-1,D[i][j-1]['value']-1]
  112. #            #print('intex i: {}, j: {}, candidates: {}'.format(i,j,candidates))
  113. #            index_array = np.argwhere(candidates == np.amax(candidates)).flatten().tolist()
  114. #            print(index_array)
  115. #            D[i][j]['value']=candidates[index_array[-1]]
  116. #            for k in index_array:
  117. #                if k == 0:
  118. #                    D[i][j]['symbol']+="\\"
  119. #                elif k == 1:
  120. #                    D[i][j]['symbol']+='^'
  121. #                elif k == 2:
  122. #                    D[i][j]['symbol']+='<'
  123. #            #print(index_array)
  124. #            #print('A[i-1]: {},B[j-1]: {},B[j]: {})'.format(A[i-1],B[j-1],B[j]))
  125. #    print(D)
  126. #    #Tracing arrows back
  127. #    print(D[3][4]['symbol'])
  128. #    return
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement