程式能實現什麼

a.完成gap值的自定義輸入以及兩條需比對序列的輸入
b.完成得分矩陣的計算及輸出
c.輸出序列比對結果
d.使用matplotlib對得分矩陣路徑的繪製

一、實現步驟

1.使用者輸入步驟

a.輸入自定義的gap值
b.輸入需要比對的鹼基序列1（A,T,C,G）換行表示輸入完成
b.輸入需要比對的鹼基序列2（A,G）換行表示輸入完成

輸入（示例）：

2.程式碼實現步驟

1.獲取到使用者輸入的gap，s以及t
2.呼叫構建得分矩陣函式，得到得分矩陣以及方向矩陣
3.將得到的得分矩陣及方向矩陣作為引數傳到回溯函式中開始回溯得到路徑，路徑儲存使用的是全域性變數，存的仍然是方向而不是座標位置減少儲存開銷，根據全域性變數中儲存的方向將比對結果輸出。

全域性比對程式碼如下：

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,9,0],[-4,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,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,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,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,j-1,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,i-1,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,t)
      finalList.pop()

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

#畫圖函式
def drawArrow(mark,a,b,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,0)
        lX = lX - 1
      elif m == 1:
        arrow(ax,1)
        lX = lX - 1
        lY = lY + 1
      elif m == 2:
        arrow(ax,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,gap)
  print("The scoring matrix is as follows：")     #輸出得分矩陣
  print(mark)
  remount(mark,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,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,3)) #direction矩陣用來儲存得分矩陣中得分來自的方向  第一個表示左方 第二個表示左上 第三個表示上方 1表示能往哪個方向去
                    #由於得分可能會來自多個方向，所以使用三維矩陣存
  for i in range(1,threeGrade)           #finalGrade表示最終來自三個方向上的得分
      if max(finalGrade) >= 0:                  #如果該最大值是大於0的則 選取三個方向上的最大得分存入得分矩陣 否則不對矩陣進行修改
        mark[i][j] = max(finalGrade)
        for k in range(0,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,t):
  if direction[i][j][0] == 1 :
    if all(direction[i][j-1] == [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,t)
      finalList.pop()            #該方向走完後將這個方向彈出 繼續下一輪判斷 下面兩個大的判斷同理
  if direction[i][j][1] == 1 :
    if all(direction[i-1][j-1] == [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,t)
      finalList.pop()
  if direction[i][j][2] == 1 :
    if all(direction[i-1][j] == [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,t)
      finalList.pop()
#畫箭頭函式
def arrow(ax,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,int(mark[i-1][j-1]))
  lX = my + 1
  lY = a - mx + 1
  for n in range(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,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,t) #呼叫回溯函式
  print(finalOrder1)
  print(finalOrder2)
  drawArrow(mark,j)

二、實驗結果截圖

1.全域性比對

2.區域性比對

到此這篇關於詳解基於python的全域性與區域性序列比對的實現(DNA)的文章就介紹到這了,更多相關python全域性與區域性序列比對內容請搜尋我們以前的文章或繼續瀏覽下面的相關文章希望大家以後多多支援我們！

總結

本次實驗使用動態規劃對全域性序列比對進行了實現，自己卡的最久的地方是回溯以及畫圖的時候。剛開始在實現回溯的過程中，老是找不準回溯的條件以及將所有的路徑都記錄下來的方法，最後是使用的方向矩陣，也就是重新定義一個與得分矩陣等大的矩陣（但是這個矩陣是三維），存放的是每個位置能夠回溯的方向，第一個數值表示左邊，第二個表示左上，第三個表示上方，為0時表示當前方向不能回溯，沒有路徑，為1時表示能回溯，當該位置的所有能走的方向都走完時即可返回。將所有路徑記錄下來的方法是定義全域性變數，當有路徑能夠走到終點時便將這條路徑存放入該全域性變數中。
繪圖的時候使用的是matplotlib中的散點圖，然後將每個點的得分以註釋的形式標記在該點的右上角，並用箭頭將路徑繪出。不得不說的是，這個圖確實太醜了，我學識淺薄，也沒想到能畫出這個圖的更好的方法，還希望老師指點。
總的來說這次實驗經歷的時間還比較長，主要是因為python也沒有很熟悉，很多函式也是查了才知道，然後視覺化更是瞭解的少，所以畫出來的圖出奇的醜，還有回溯的時候也是腦子轉不過彎來，所以要學習的東西還有很多，需要更加努力。
本次實驗還能夠有所改進的地方是：
1.把兩個比對演算法結合，讓使用者能夠選擇使用哪種比對方式。
2.作出一個更好看的介面，增加使用者體驗感。
3.把圖畫的更美觀。
（老丁已閱，USC的同學們謹慎借鑑）