pythonに基づくグローバル・シーケンス対の実装(DNA)


プログラムは何を実現できますか?
a.gap値のカスタム入力とシーケンスに対する入力を完了する
b.得点マトリックスの計算と出力を完了する
c.出力シーケンス比結果
d.ポイントマトリックスの経路をmatplotlibで描画する
一、実現ステップ
1.ユーザ入力ステップ
a.カスタムのgap値を入力する
b.入力は対の塩基配列1(A,T,C,G)を必要とし、改行は入力完了を表します。
b.入力は対の塩基配列2(A,T,C,G)を必要とし、改行は入力完了を表します。
入力(例):
在这里插入图片描述
2.コード実現手順
1.ユーザ入力のgap、s及びtを取得する。
2.得点行列関数の構築を呼び出し、得点行列と方向行列を得る
3.得られた得点行列と方向行列をパラメータとしてトレース関数に戻していく経路です。経路格納は大域変数を使用しています。保存されているのは座標位置ではなく、方向です。格納オーバーヘッドは大域変数に格納されている方向によって比較結果を出力します。
4.グローバル変数に格納されている方向によって、matplotlibを使って経路を描き出します。
グローバルマッチングコードは以下の通りです。

import matplotlib.pyplot as plt
import numpy as np

#        finalList          finalOrder1,finalOrder2        finalRoad              
def createList():
  global finalList
  global finalOrder1
  global finalOrder2
  global finalRoad
  finalList = []
  finalOrder1 = []
  finalOrder2 = []
  finalRoad = []


#  A G C T       ,              
def createDic():
  dic = {'A':0,'G':1,'C':2,'T':3}
  return dic
#      
# A G C T
def createGrade():
  grade = np.matrix([[10,-1,-3,-4],
            [-1,7,-5,-3],
            [-3,-5,9,0],
            [-4,-3,0,8]])
  return grade

#              
def getGrade(a,b):
  dic = createDic() #              
  grade = createGrade() #     grade
  return grade[dic[a],dic[b]]

#                 、    gap 
def createMark(s,t,gap):
  a = len(s)             #      a,b
  b = len(t)
  mark = np.zeros((a+1,b+1))     #         
  direction = np.zeros((a+1,b+1,3)) #direction                                            1         
                    #             ,          

  direction[0][0] = -1        #                       -1
  mark[0,:] = np.fromfunction(lambda x, y: gap * (x + y), (1, b + 1), dtype=int)   #  gap            
  mark[:,0] = np.fromfunction(lambda x, y: gap * (x + y), (1, a + 1), dtype=int)   #  gap            
  for i in range(1,b+1):
    direction[0,i,0] = 1
  for i in range(1, a + 1):
    direction[i, 0, 2] = 1

  for i in range(1,a+1):
    for j in range(1,b+1):
      threeMark = [mark[i][j-1],mark[i-1][j-1],mark[i-1][j]]     #threeMark                         
      threeGrade = [gap,getGrade(s[i-1],t[j-1]),gap]         #threeGrade                           
      finalGrade = np.add(threeMark,threeGrade)           #finalGrade              
      mark[i][j] = max(finalGrade)                  #                  
      #                 ,           
      for k in range(0,len([y for y,x in enumerate(finalGrade) if x == max(finalGrade)])):
        directionList = [y for y,x in enumerate(finalGrade) if x == max(finalGrade)]
        direction[i][j][directionList[k]] = 1
  return mark,direction

#                                       
def remount(mark,direction,i,j,s,t):
  if direction[i][j][0] == 1 :
    if direction[i][j-1][0] == -1:      #                     
      finalList.append(0)          #              
      finalList.reverse()          #                
      index1 = 0              #         s                   
      index2 = 0              #         t   
      for k in finalList:
        if k == 0 :
          finalOrder1.append("-")
          finalOrder2.append(t[index2])
          index2 += 1
        if k == 1 :
          finalOrder1.append(s[index1])
          finalOrder2.append(t[index2])
          index1 += 1
          index2 += 1
        if k == 2 :
          finalOrder1.append(s[index1])
          finalOrder2.append("-")
          index1 += 1
      finalList.reverse() #             
      finalRoad.append(np.array(finalList)) #                      
      finalList.pop()            #              
      return
    else :
      finalList.append(0)          #                  ,     
      remount(mark,direction,i,j-1,s,t)
      finalList.pop()            #                                
  if direction[i][j][1] == 1 :
    if direction[i-1][j-1][0] == -1:

      finalList.append(1)
      finalList.reverse()          #                 
      index1 = 0              #          s                   
      index2 = 0              #          t   
      for k in finalList:
        if k == 0 :
          finalOrder1.append("-")
          finalOrder2.append(t[index2])
          index2 += 1
        if k == 1 :
          finalOrder1.append(s[index1])
          finalOrder2.append(t[index2])
          index1 += 1
          index2 += 1
        if k == 2 :
          finalOrder1.append(s[index1])
          finalOrder2.append("-")
          index1 += 1
      finalList.reverse() #             
      finalRoad.append(np.array(finalList)) #                      
      finalList.pop()
      return
    else :
      finalList.append(1)
      remount(mark,direction,i-1,j-1,s,t)
      finalList.pop()
  if direction[i][j][2] == 1 :
    if direction[i-1][j][0] == -1:
      finalList.append(2)
      finalList.reverse()           #                 
      index1 = 0                #          s                   
      index2 = 0                #          t   
      for k in finalList:
        if k == 0 :
          finalOrder1.append("-")
          finalOrder2.append(t[index2])
          index2 += 1
        if k == 1 :
          finalOrder1.append(s[index1])
          finalOrder2.append(t[index2])
          index1 += 1
          index2 += 1
        if k == 2 :
          finalOrder1.append(s[index1])
          finalOrder2.append("-")
          index1 += 1
      finalList.reverse()          #             
      finalRoad.append(np.array(finalList)) #                      
      finalList.pop()
      return
    else :
      finalList.append(2)
      remount(mark,direction,i-1,j,s,t)
      finalList.pop()

#     
def arrow(ax,sX,sY,aX,aY):
  ax.arrow(sX,sY,aX,aY,length_includes_head=True, head_width=0.15, head_length=0.25, fc='w', ec='b')

#    
def drawArrow(mark, direction, a, b, s, t):
  #a s    4  b t    6
  fig = plt.figure()
  ax = fig.add_subplot(111)
  val_ls = range(a+2)
  scale_ls = range(b+2)
  index_ls = []
  index_lsy = []
  for i in range(a):
    if i == 0:
      index_lsy.append('#')
    index_lsy.append(s[a-i-1])
  index_lsy.append('0')
  for i in range(b):
    if i == 0:
      index_ls.append('#')
      index_ls.append('0')
    index_ls.append(t[i])
  plt.xticks(scale_ls, index_ls)      #     
  plt.yticks(val_ls, index_lsy)
  for k in range(1,a+2):
    y = [k for i in range(0,b+1)]
    x = [x for x in range(1,b+2)]
    ax.scatter(x, y, c='y')
  for i in range(1,a+2):
    for j in range(1,b+2):
      ax.text(j,a+2-i,int(mark[i-1][j-1]))
  lX = b+1
  lY = 1
  for n in range(0,len(finalRoad)):
    for m in (finalRoad[n]):
      if m == 0:
        arrow(ax,lX,lY,-1,0)
        lX = lX - 1
      elif m == 1:
        arrow(ax,lX,lY,-1,1)
        lX = lX - 1
        lY = lY + 1
      elif m == 2:
        arrow(ax, lX, lY, 0, 1)
        lY = lY + 1
    lX = b + 1
    lY = 1
  ax.set_xlim(0, b + 2) #        ,   [0,1]
  ax.set_ylim(0, a + 2) #        ,   [0,1]
  ax.set_aspect('equal') # x  y    
  plt.show()
  plt.tight_layout()
if __name__ == '__main__':
  createList()
  print("Please enter gap:")
  gap = int(input())       #  gap         tip:                          
  print("Please enter sequence 1:")
  s = input()           #            
  print("Please enter sequence 2:")
  t = input()           #            
  a = len(s)           #  s   
  b = len(t)           #  t   
  mark,direction = createMark(s,t,gap)
  print("The scoring matrix is as follows:")     #      
  print(mark)
  remount(mark,direction,a,b,s,t) #      
  c = a if a > b else b     #                     
  total = int(len(finalOrder1)/c)
  for i in range(1,total+1):   #        
    k = str(i)
    print("Sequence alignment results "+k+" is:")
    print(finalOrder1[(i-1)*c:i*c])
    print(finalOrder2[(i-1)*c:i*c])
  drawArrow(mark, direction, a, b, s, t)

ローカルマッチングコードは以下の通りです。

import matplotlib.pyplot as plt
import numpy as np
import operator
#                           0
#        finalList          finalOrder1,finalOrder2       
def createList():
  global finalList
  global finalOrder1
  global finalOrder2
  global finalRoad
  finalList = []
  finalOrder1 = []
  finalOrder2 = []
  finalRoad = []


#  A G C T       ,              
def createDic():
  dic = {'A':0,'G':1,'C':2,'T':3}
  return dic
#      
# A G C T
def createGrade():
  grade = np.matrix([[10,-1,-3,-4],
            [-1,7,-5,-3],
            [-3,-5,9,0],
            [-4,-3,0,8]])
  return grade

#              
def getGrade(a,b):
  dic = createDic() #              
  grade = createGrade() #     grade
  return grade[dic[a],dic[b]]

#                 、    gap 
def createMark(s,t,gap):
  a = len(s)             #      a,b
  b = len(t)
  mark = np.zeros((a+1,b+1))     #         
  direction = np.zeros((a+1,b+1,3)) #direction                                            1         
                    #             ,         
  for i in range(1,a+1):
    for j in range(1,b+1):
      threeMark = [mark[i][j-1],mark[i-1][j-1],mark[i-1][j]]     #threeMark                         
      threeGrade = [gap,getGrade(s[i-1],t[j-1]),gap]         #threeGrade                           
      finalGrade = np.add(threeMark,threeGrade)           #finalGrade              
      if max(finalGrade) >= 0:                  #         0                                
        mark[i][j] = max(finalGrade)
        for k in range(0,len([y for y,x in enumerate(finalGrade) if x == max(finalGrade)])): #                 ,           
          directionList = [y for y,x in enumerate(finalGrade) if x == max(finalGrade)]
          direction[i][j][directionList[k]] = 1
  return mark,direction

#                                       
def remount(mark,direction,i,j,s,t):
  if direction[i][j][0] == 1 :
    if all(direction[i][j-1] == [0,0,0]):      #                     
      finalList.append(0)          #              
      finalList.reverse()          #                
      index1 = i              #         s                   
      index2 = j-1              #         t   
      for k in finalList:
        if k == 0 :
          finalOrder1.append("-")
          finalOrder2.append(t[index2])
          index2 += 1
        if k == 1 :
          finalOrder1.append(s[index1])
          finalOrder2.append(t[index2])
          index1 += 1
          index2 += 1
        if k == 2 :
          finalOrder1.append(s[index1])
          finalOrder2.append("-")
          index1 += 1
      finalList.reverse()
      finalRoad.append(np.array(finalList)) #                      
      finalList.pop()            #              
      return
    else :
      finalList.append(0)          #                  ,     
      remount(mark,direction,i,j-1,s,t)
      finalList.pop()            #                                
  if direction[i][j][1] == 1 :
    if all(direction[i-1][j-1] == [0,0,0]):
      finalList.append(1)
      finalList.reverse()           #                 
      index1 = i-1              #          s                   
      index2 = j-1              #          t   
      for k in finalList:
        if k == 0 :
          finalOrder1.append("-")
          finalOrder2.append(t[index2])
          index2 += 1
        if k == 1 :
          finalOrder1.append(s[index1])
          finalOrder2.append(t[index2])
          index1 += 1
          index2 += 1
        if k == 2 :
          finalOrder1.append(s[index1])
          finalOrder2.append("-")
          index1 += 1
      finalList.reverse()
      finalRoad.append(np.array(finalList)) #                      
      finalList.pop()
      return
    else :
      finalList.append(1)
      remount(mark,direction,i-1,j-1,s,t)
      finalList.pop()
  if direction[i][j][2] == 1 :
    if all(direction[i-1][j] == [0,0,0]):
      finalList.append(2)
      finalList.reverse()           #                 
      index1 = i-1                #          s                   
      index2 = j                #          t   
      for k in finalList:
        if k == 0 :
          finalOrder1.append("-")
          finalOrder2.append(t[index2])
          index2 += 1
        if k == 1 :
          finalOrder1.append(s[index1])
          finalOrder2.append(t[index2])
          index1 += 1
          index2 += 1
        if k == 2 :
          finalOrder1.append(s[index1])
          finalOrder2.append("-")
          index1 += 1
      finalList.reverse()
      finalRoad.append(np.array(finalList)) #                      
      finalList.pop()
      return
    else :
      finalList.append(2)
      remount(mark,direction,i-1,j,s,t)
      finalList.pop()
#     
def arrow(ax,sX,sY,aX,aY):
  ax.arrow(sX,sY,aX,aY,length_includes_head=True, head_width=0.15, head_length=0.25, fc='w', ec='b')

#    
def drawArrow(mark, direction, a, b, s, t,mx,my):
  #a s    4  b t    6
  fig = plt.figure()
  ax = fig.add_subplot(111)
  val_ls = range(a+2)
  scale_ls = range(b+2)
  index_ls = []
  index_lsy = []
  for i in range(a):
    if i == 0:
      index_lsy.append('#')
    index_lsy.append(s[a-i-1])
  index_lsy.append('0')
  for i in range(b):
    if i == 0:
      index_ls.append('#')
      index_ls.append('0')
    index_ls.append(t[i])
  plt.xticks(scale_ls, index_ls)      #     
  plt.yticks(val_ls, index_lsy)
  for k in range(1,a+2):
    y = [k for i in range(0,b+1)]
    x = [x for x in range(1,b+2)]
    ax.scatter(x, y, c='y')
  for i in range(1,a+2):
    for j in range(1,b+2):
      ax.text(j,a+2-i,int(mark[i-1][j-1]))
  lX = my + 1
  lY = a - mx + 1
  for n in range(0,len(finalRoad)):
    for m in (finalRoad[n]):
      if m == 0:
        arrow(ax,lX,lY,-1,0)
        lX = lX - 1
      elif m == 1:
        arrow(ax,lX,lY,-1,1)
        lX = lX - 1
        lY = lY + 1
      elif m == 2:
        arrow(ax, lX, lY, 0, 1)
        lY = lY + 1
    lX = b + 1
    lY = 1
  ax.set_xlim(0, b + 2) #        ,   [0,1]
  ax.set_ylim(0, a + 2) #        ,   [0,1]
  ax.set_aspect('equal') # x  y    
  plt.show()
  plt.tight_layout()

if __name__ == '__main__':
  createList()
  print("Please enter gap:")
  gap = int(input())       #  gap         tip:                          
  print("Please enter sequence 1:")
  s = input()           #            
  print("Please enter sequence 2:")
  t = input()           #            
  a = len(s)           #  s   
  b = len(t)           #  t   
  mark,direction = createMark(s,t,gap)
  print("The scoring matrix is as follows:")     #      
  print(mark)
  maxDirection = np.argmax(mark) #        
  i = int(maxDirection/(b+1))
  j = int(maxDirection - i*(b+1))
  remount(mark,direction,i,j,s,t) #      
  print(finalOrder1)
  print(finalOrder2)
  drawArrow(mark, direction, a, b, s, t, i, j)
二、実験結果のスクリーンショット
1.グローバルマッチング
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
2.ローカルマッチング
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
ここではpythonに基づく大局と局部シーケンス比の実現(DNA)に関する記事を紹介します。より多くの関連python大局と局部シーケンス比の内容を紹介します。以前の文章を検索してください。または下記の関連記事を引き続きご覧ください。これからもよろしくお願いします。
締め括りをつける
今回の実験は動態計画を使って全体シーケンスを比較して実現しました。自分のカードの一番長いところは遡及及び絵を描く時です。遡ってはいけない条件を常に探して、すべての経路を記録する方法です。最後に使用する方向行列です。つまり、得点行列などの大きな行列を再定義します。位置ごとに遡る方向です。一番目の数値は左、二つ目は左上、第三は上を表しています。0の場合は現在の方向はさかのぼらないこと、パスがないこと、1の場合は遡ることができることを示しています。この位置のすべての歩く方向が終わったら戻れます。すべてのパスを記録する方法はグローバル変数を定義し、パスがあればそのパスをグローバル変数に保存することです。
絵を描くときは、matplotlibの散点図を使用して、各点の得点をその点の右上隅に注釈で表示し、矢印で経路を描画します。言わなければならないのは、この図は確かに醜すぎて、学識が浅くて、この絵を描くことができるとは思いませんでした。
全体としては、今回の実験の経験はまだ長いです。pythonもよく知られていないので、多くの関数が調べてから分かります。そして、可視化がもっとよく分かりません。だから、絵は奇抜な醜さがあります。
今回の実験で改善されたところは:
1.2つの対アルゴリズムを結合して、ユーザにどのような対を使うかを選択させます。
2.よりよく見るインターフェースを作って、ユーザーの体験感を高める。
3.絵をより美しくする。
(丁さんはもう読みました。USCの学生たちは参考にしてください。)