1.序列比較演算法（全域序列比對及區域序列比對的python實作）
- 前言
- 演算法思想介紹
- 實作功能及實作方法
- 運行結果演示
- 源代碼
- 遇到的問題及總結
前言
階段性地完成了DNA序列比較演算法，還有很多不足和需要完善的地方有待日后改進，
演算法思想介紹
一個很詳細完整的演算法介紹
雙序列全域比對及演算法
Needleman-Wunsch 演算法：動態規劃法
輸入值：兩條序列、替換記分矩陣以確定不同字母間的相似度得分，以及空位罰分
雙序列區域比對演算法
區域比對的計算公式在全域比對的基礎上增加了第四個元素“0”，得分矩陣初始值仍是0，但第一行和第一列與全域比對不同，全是0，
實作功能及實作方法
- 使用已定義的記分矩陣
- 設定需比對的序列、gap大小及要進行的比對選擇
- 列印最高得分的序列比對結果
方法：倒序查找最終路進行序列比對 - 列印得分矩陣及比對路徑
方法：使用遞回和堆疊記錄最終路徑
運行結果演示
-
雙序列全域比對
輸入序列：atcggtac；aatcgta
輸入序列：atcggt；aacg
-
雙序列區域比對
輸入序列：atccga；tcga
輸入序列：acgtc；cg
源代碼
#!/usr/bin/env python
# -*- coding:utf-8 -*-
from numpy import *
import copy
from matplotlib import pyplot as plt
from pandas import *
#創建全域比對得分矩陣
def GlobalScoreMatrix(m,n,w,replace,s,path,senquence1,senquence2,gap):
for i in range(m):
for j in range(n):
#判斷s（0,0）這一特殊情況
if i == 0 and j == 0:
s[i][j] = 0
elif i-1 < 0:#判斷第一行的特殊情況
s[i][j] = s[i][j - 1] + gap
path[i,j,0] = 1
elif j-1 < 0:#判斷第一列的特殊情況
s[i][j] = s[i - 1][j] + gap
path[i,j,1] = 1
else:
#獲取最大值
s[i][j] = max(s[i - 1][j - 1] + w[replace[senquence2[i - 1]]][replace[senquence1[j - 1]]],
s[i - 1][j] + gap, s[i][j - 1] + gap)
#記錄最大值來的方向
if s[i - 1][j - 1] + w[replace[senquence2[i - 1]]][replace[senquence1[j - 1]]] == s[i][j]:
path[i,j,2] = 1
if s[i - 1][j] + gap == s[i][j]:
path[i,j,1] = 1
if s[i][j - 1] + gap == s[i][j]:
path[i,j,0] = 1
#創建區域比對得分矩陣
def LocalScoreMatrix(m,n,w,replace,s,path,senquence1,senquence2,gap):
for i in range(1,m):
#區域矩陣第一行及第一列均為0，不需要再初始化
for j in range(1,n):
#獲取最大值,與全域比對不同之處在于選取最大值時將0加入選擇
s[i][j] = max(0,s[i - 1][j - 1] + w[replace[senquence2[i - 1]]][replace[senquence1[j - 1]]],
s[i - 1][j] + gap, s[i][j - 1] + gap)
#記錄最大值來的方向，若最大值為0則不必記錄
if s[i - 1][j - 1] + w[replace[senquence2[i - 1]]][replace[senquence1[j - 1]]] == s[i][j]:
path[i,j,2] = 1
if s[i - 1][j] + gap == s[i][j]:
path[i,j,1] = 1
if s[i][j - 1] + gap == s[i][j]:
path[i,j,0] = 1
#尋找全域序列比對路徑
def FindGlobalPath(i,j,path,OnePath,LastGlobalPath):
#遞回終止條件：回到原點（0，0）
if i == 0 and j == 0:
OnePath.append((i, j))
#將OnePath進行深拷貝再加入至最終路徑串列LastGlobalPath中
LastGlobalPath.append(copy.deepcopy(OnePath))
#將該點出堆疊
OnePath.pop()
else:
for k in range(3):
#判斷每個點來的方向
if path[i,j,k] == 1:
#下標0處記錄從左來的方向
if k == 0:
#將該點入堆疊
OnePath.append((i,j))
#進行遞回記錄下一個點
FindGlobalPath(i,j - 1,path,OnePath,LastGlobalPath)
#遞回回傳后將該點出堆疊，記錄另一方向
OnePath.pop()
#下標1處記錄從上來的方向
elif k == 1:
OnePath.append((i, j))
FindGlobalPath(i - 1, j, path,OnePath,LastGlobalPath)
OnePath.pop()
#下標2處記錄從左上來的方向
else:
OnePath.append((i, j))
FindGlobalPath(i - 1, j - 1, path,OnePath,LastGlobalPath)
OnePath.pop()
# 尋找區域序列比對路徑
def FindLocalPath(i, j, path, OnePath, LastLocalPath):
#遞回終止條件：某個沒有記錄方向的點
if all(path[i][j] == [0, 0, 0]):
OnePath.append((i, j))
# 將OnePath進行深拷貝再加入至最終路徑串列LastLocalPath中
LastLocalPath.append(copy.deepcopy(OnePath))
# 將該點出堆疊
OnePath.pop()
else:
for k in range(3):
# 判斷每個點來的方向
if path[i, j, k] == 1:
# 下標0處記錄從左來的方向
if k == 0:
# 將該點入堆疊
OnePath.append((i, j))
# 進行遞回記錄下一個點
FindLocalPath(i, j - 1, path, OnePath, LastLocalPath)
# 遞回回傳后將該點出堆疊，記錄另一方向
OnePath.pop()
# 下標1處記錄從上來的方向
elif k == 1:
OnePath.append((i, j))
FindLocalPath(i - 1, j, path, OnePath, LastLocalPath)
OnePath.pop()
# 下標2處記錄從左上來的方向
else:
OnePath.append((i, j))
FindLocalPath(i - 1, j - 1, path, OnePath, LastLocalPath)
OnePath.pop()
#輸出比對后的序列
def ShowContrastResult(LastPath,senquence1,senquence2):
#依次輸出每條路徑
for k,aPath in enumerate(LastPath):
rowS = ''
colS = ''
#每條路徑倒序遍歷
for i in range(len(aPath) -1,0,-1):
#方向從左邊來
if aPath[i][0] == aPath[i - 1][0]:
rowS += senquence1[aPath[i - 1][1] - 1]
colS += '-'
#方向從上面來
elif aPath[i][1] == aPath[i - 1][1]:
colS += senquence2[aPath[i - 1][0] -1]
rowS += '-'
#方向從左上來
else:
rowS += senquence1[aPath[i - 1][1] - 1]
colS += senquence2[aPath[i - 1][0] - 1]
#依次輸出每條路的序列比對結果
print("======比對結果",k+1,"======")
print("序列1:",rowS)
print("序列2:",colS)
# 判斷是否為最終路徑中的點
def judgePath(point, LastPath):
for aPath in LastPath:
if point in aPath:
return True
return False
# 畫出路徑圖
def ShowPaths(senquence1, senquence2, LastPath):
s1 = "0" + senquence1
s2 = "0" + senquence2
# 列索引
col = list(s1)
# 行索引
row = list(s2)
# 設定畫布大小
fig = plt.figure(figsize=(5, 5))
ax = fig.add_subplot(111, frameon=True, xticks=[], yticks=[], )
the_table = plt.table(cellText=s, rowLabels=row, colLabels=col, rowLoc='right',
loc='center', cellLoc='bottom right', bbox=[0, 0, 1, 1])
# 設定表格文本字體大小
the_table.set_fontsize(8)
# 畫出每個點的路徑圖
for i in range(m):
for j in range(n):
for k in range(3):
if path[i, j, k] == 1: # 畫出記錄的方向
# 下標0處記錄從左來的方向
if k == 0:
if judgePath((i, j), LastPath): # 若某點在在最終路徑中
# 畫出紅色箭頭
plt.annotate('', xy=(j / n, (2 * m - 2 * i - 1) / (2 * (m + 1))),
xytext=((2 * j + 1) / (2 * n), (2 * m - 2 * i - 1) / (2 * (m + 1))),
arrowprops=dict(arrowstyle="->", color='r', connectionstyle="arc3"))
else:
# 未在最終路徑中則畫出黑色箭頭
plt.annotate('', xy=(j / n, (2 * m - 2 * i - 1) / (2 * (m + 1))),
xytext=((2 * j + 1) / (2 * n), (2 * m - 2 * i - 1) / (2 * (m + 1))),
arrowprops=dict(arrowstyle="->", connectionstyle="arc3"))
# 下標1處記錄從上來的方向
elif k == 1:
if judgePath((i, j), LastPath):
plt.annotate('', xy=((2 * j + 1) / (2 * n), (2 * m - 2 * i - 1) / (2 * (m + 1))),
xytext=((2 * j + 1) / (2 * n), (m - i) / (m + 1)),
arrowprops=dict(arrowstyle="<-", color='r', connectionstyle="arc3"))
else:
plt.annotate('', xy=((2 * j + 1) / (2 * n), (2 * m - 2 * i - 1) / (2 * (m + 1))),
xytext=((2 * j + 1) / (2 * n), (m - i) / (m + 1)),
arrowprops=dict(arrowstyle="<-", connectionstyle="arc3"))
# 下標1處記錄從上來的方向
elif k == 2:
if judgePath((i, j), LastPath):
plt.annotate('', xy=((2 * j + 1) / (2 * n), (2 * m - 2 * i - 1) / (2 * (m + 1))),
xytext=(j / n, (m - i) / (m + 1)),
arrowprops=dict(arrowstyle="<-", color='r', connectionstyle="arc3"))
else:
plt.annotate('', xy=((2 * j + 1) / (2 * n), (2 * m - 2 * i - 1) / (2 * (m + 1))),
xytext=(j / n, (m - i) / (m + 1)),
arrowprops=dict(arrowstyle="<-", connectionstyle="arc3"))
plt.show()
#將堿基轉換為集合下標
replace = {'A':0,'G':1,'C':2,'T':3}
#構造替換計分矩陣
w = [[10,-1,-3,-4],[-1,7,-5,-3],[-3,-5,9,0],[-4,-3,0,8]]
#定義需要比對的序列
senquence1 = input("請輸入序列1：").upper()
senquence2 = input("請輸入序列2：").upper()
#定義輸入的gap
gap = int(input("請輸入gap："))
choise = int(input("請選擇要進行的序列比對(1-全域序列比對 2-區域序列比對) : "))
# 獲取序列的長度
m = len(senquence2) + 1
n = len(senquence1) + 1
#構建m*n全0矩陣
s = zeros((m,n))
#記錄每個點的方向，下標0處存盤從左來的方向，下標1處存盤從上來的方向，下標2處存盤從左上來的方向
#初始值均為0，若存在從某方向上來則將其對應下標的值置為1
path = zeros((m,n,3))
#記錄每條路徑
OnePath = []
#記錄所有全域序列比對路徑
LastGlobalPath = []
#記錄所有區域序列比對路徑
LastLocalPath = []
if choise == 1:#進行全域序列比對
#構建得分矩陣
GlobalScoreMatrix(m,n,w,replace,s,path,senquence1,senquence2,gap)
FindGlobalPath(m-1,n-1,path,OnePath,LastGlobalPath)
ShowContrastResult(LastGlobalPath,senquence1,senquence2)
ShowPaths(senquence1, senquence2, LastGlobalPath)
elif choise == 2:#進行區域序列比對
#構建得分矩陣
LocalScoreMatrix(m, n, w, replace, s, path, senquence1, senquence2, gap)
row = where(s == np.max(s))[0]
# 獲取得分矩陣中最大值的列索引
col = where(s == np.max(s))[1]
for i in range(len(row)):#依次尋找不同區域比對的比對路徑并輸出比對結果
FindLocalPath(row[i], col[i], path, OnePath, LastLocalPath)
ShowContrastResult(LastLocalPath, senquence1, senquence2)
ShowPaths(senquence1, senquence2, LastLocalPath)
else:
print("無效選擇！")
遇到的問題及總結
-
思維誤區： 最初在存盤最終路徑的問題上，認為在每次路徑遞回至初始結點時才將其入堆疊，導致遇到分岔路則無法記錄前一條路的完整路徑
經過高人指點后解決：每次到達一個結點就將其入堆疊，遞回結束后將其出堆疊 -
思維誤區2：最初以采用存盤每個點的前一點坐標為存盤方式，導致所有得分矩陣只能存盤一條路徑
再次經過高人指點后解決：改存盤方式為存盤每個點計算得來的方向 -
遞回終止條件存盤最終路徑 -涉及問題：list的賦值、淺拷貝、深拷貝
Python List的賦值方法 -
畫出路徑矩陣表格 -涉及問題：表格與畫布大小不一致且導致無法確定箭頭坐標
解決方式：使用bbox引數 -
獲取得分矩陣中最大值的索引 - 涉及問題：獲取where結果的索引值
-
區域比對遞回終止條件 - 涉及問題：串列比較 any() all()
總結
這次的作業因為拖延很晚才開始，遇到的問題也絕不僅僅只有以上幾個，最深刻的還是最開始的思維誤區，直接卡到崩潰，但其他問題也都耗費了或多或少的時間才解決，更有無數小問題除錯了無數遍才得以做出這個半成品，其實還有很多待完善的地方，比如輸入的字符判斷，比如一致度和相似度的計算，比如圖形化界面的撰寫，寫出來的代碼也還有很多不成熟的地方，但是因為時間問題，暫時就到這了，有時間再改進咯，
轉載請註明出處，本文鏈接：https://www.uj5u.com/qianduan/159729.html
標籤：其他