在序列比对的时候，有全局比对和局部比对两种方法，其中，Needleman-Wunsch比对算法是其中的一个很经典的全局比对算法。下面将用python从头实现，将考虑match，mismatch，gap和gap是否连续的因素。
先确定打分策略，先考虑match，mismatch和gap的分数。
现将match定为1分，mismatch -1分，gap是-2分。
现有两个序列，一个是ATCG，一个是TCG，以下是初始分值矩阵：
first_score
第一列和第一行表示一开始比对到gap上，因此计入gap的分数。
接着计算每行每列的分数，比如第三行第三列，可以得到三个分数，同一行列的为一个gap：
next
从三种路径中选择分数最高的作为这个位置的分数。接下来每个位置都要进行这样的计算，直到填满这个表格。
form
同时可以得到最佳路径：
path
以下用python实现：
- 比较两个碱基是否一致
# 比较两个碱基分数， gaps的分数考虑在cal_score
def diff(first, second):
if first == second:
return match
else:
return mismatch
- 计算分数，并考虑多段gap的情况。一般来说，结果应该更倾向于仅出现一段长gap而不是多个短gap。所以将第一次出现的gap定位-2分，第二次及之后连续出现的gap定为一个更小的分数。 即变量extend。
def cal_score(seq1, seq2):
nrow = len(seq1)
ncol = len(seq2)
scores = np.zeros((nrow, ncol))
path = np.zeros((nrow, ncol))
# 初始化第一行第一列
for index1 in range(nrow):
scores[index1,0] = gap * index1
path[index1,0] = 0 #上面
for index2 in range(ncol):
scores[0,index2] = gap * index2
path[0, index2] = 2 #左侧
for num1 in range(1, nrow):
for num2 in range(1, ncol):
# 得到上面，斜上和左侧的结果
last_score = [scores[num1 - 1, num2], scores[num1 - 1, num2 - 1], scores[num1, num2 - 1]]
change_score = diff(seq1[num1], seq2[num2])
current_score = []
if path[num1 - 1 , num2] == 0:
current_score.append(scores[num1-1, num2] + extend) # 上面
else:
current_score.append(scores[num1-1, num2] + gap)
current_score.append(scores[num1-1, num2 - 1] + change_score) # 斜上
if path[num1 , num2 -1] == 2:
current_score.append(scores[num1, num2-1] + extend) # 左侧
else:
current_score.append(scores[num1, num2-1] + gap)
current_index = current_score.index(max(current_score)) # 当前索引，不是0就gap
scores[num1, num2] = max(current_score)
path[num1, num2] = current_index
return scores, path
- 回溯路径并输出结果
def cal_seq(scores, path):
index1 = len(seq1) - 1
index2 = len(seq2) -1
top = ''
middle = ''
bottom = ''
while True:
if path[index1, index2] == 1:
top += seq1[index1]
bottom += seq2[index2]
if seq1[index1] == seq2[index2]:
middle += '|'
else:
middle += ' '
index1 -= 1
index2 -= 1
elif path[index1, index2] == 0:
top += seq1[index1]
bottom += '-'
middle += ' '
index1 -= 1
else:
top += '-'
bottom += seq2[index2]
middle += ' '
index2 -= 1
top_num = len(top) - top.count('-')
bottom_num = len(bottom) - bottom.count('-')
if top_num == max(len(seq1), len(seq2))-1 or bottom_num == max(len(seq1), len(seq2))-1:
break
return top, middle, bottom
4.命令行解析与主函数
import argparse
import numpy as np
def getArgs(argv=None):
parser = argparse.ArgumentParser()
parser.add_argument("-seq1", help = "", required = True)
parser.add_argument("-seq2", help = "", required = True)
args = parser.parse_args()
return args
if __name__ == '__main__':
args = getArgs()
seq1 = "_" + args.seq1
seq2 = "_" + args.seq2
match=5 # 或者别的分数
mismatch=-5
gap=-10
extend = -0.5
scores,path = cal_score(seq1, seq2)
top, middle, bottom = cal_seq(scores, path)
print(top[::-1])
print(middle[::-1])
print(bottom[::-1])
调用：python extend_alignment.py -seq1 ATCGATGGTATATATCGATC -seq2 ATCGATGAGTATAT
输出：
output
如有错误，欢迎指出！
