生物信息-009-双序列局部比对及算法

在看局部比对的算法之前，先看看序列 a 和序列 b 这两条序列，谁和序列 c 更相似一些。AC 做全局比对的得分是-46 分。BC 做全局比对的得分是 8 分。单从分数上看，b跟 c 更相似，因为得分高。但是如果只看红框里这部分，a 和 c 显然比 b 和 c 对得要好。如果只比较红框里这部分的话，a 和 c 的得分是 24 分，比 b 和 c 的得分要高。这个例子告诉我们，对于像这样一长一短的两条序列，比较局部比比较全长更有意义。这就是为什么除了全局比对，还有局部比对



局部比对的算法和全局比对很相似，只是在选最大值时通过增加了第四个元素“0”，来达到比对局部的效果。序列 p 和序列 q，一长一短，其他输入值跟全局比对的一样



局部比对的计算公式在全局比对的基础上增加了第四个元素“0”。得分矩阵初始值仍是0，但第一行和第一列与全局比对不同，全是 0



s(1,4)的计算：上面格 s(0,4)+gap=0+-5=-5，左边格 s(1,3)+gap=4+-5=-1，

斜上格 s(0,3)+w(1,4)=0+0=0，还有一个 0。max(-5,-1,0,0)=0。这个 0 和 s(1,1)

的 0 是不一样的。这个 0 应该画上斜上的箭头，因为它可以来自公式中的 0，也可以来自斜上格。而 s(1,1)的 0 没有箭头因为它只来自公式中的 0。两种情况虽然都是 0，但来源不同，一定要通过箭头标识清楚。



按照公式，填充满整个得分矩阵。与全局比对不同，局部比对的得分不是在右下角，而是在整个矩阵中找最大值。这个最大值才是局部比对的最终得分，他可能出现在任何一个位置。这次箭头追溯也不是从右下角到左上角，而是从刚刚找到的最大值开始追溯到没有箭头为止。追溯箭头终止的位置也可以是得分矩阵中的任何一个位置



python代码实现

def introduce():
print("*********************************************")
print("Welcome to use short sequence alignment tool!")
print("Author : sunchengquan")
print("input1: long sequence!!!!!")
print("input2: short sequence!!!!!")
print("*********************************************")
print('\n')

def matrix(seq1, seq2):
input_1 = []
input_2 = []
for i in range(len(seq1)):
input_1.append(i * 0)
for i in range(len(seq2)):
input_2.append(i * 0)
return input_1, input_2

def matchBase(base1, base2):
base = base1+base2
if base1 == base2:
match_score = {"AA": 10, "GG": 7, "CC": 9, "TT": 8}
if base1 == 'A':
return match_score['AA']
elif base1 == 'G':
return match_score['GG']
elif base1 == 'C':
return match_score['CC']
else:
return match_score['TT']
else:
mismatch = {"AG":-1,"GA":-1,"AC":-3,"CA":-3,"AT":-4,"TA":-4,"GC":-5,"CG":-5,"GT":-3,"TG":-3,"CT":0,"TC":0}
if base  == "AG" or base  == "GA":
return mismatch['AG']
elif base  == "AC" or base  == "CA":
return mismatch['AC']
elif base  == "AT" or base  == "TA":
return mismatch['AT']
elif base  == "GC" or base  == "CG":
return mismatch['GC']
elif base  == "GT" or base  == "TG":
return mismatch['GT']
elif base  == "CT" or base  == "TC":
return mismatch['CT']

def getScore(i, j, result_match,score):

num_s1 = "(" + str(j - 1) +"," +str(i - 1) + ")"
num_si = "(" + str(j - 1) +"," + str(i) + ")"
num_sj = "(" + str(j) + "," + str(i - 1) + ")"

score_1 = score[num_s1] + result_match
score_i = score[num_si]
score_j = score[num_sj]
score_max = max(score_1, score_i, score_j,0)

a = "(" + str(j) + "," + str(i) + ")"

if score_max == score_1:
con[a] = score_max
elif score_max == score_i:
a_i[a] = score_max
else:
b_j[a] = score_max

score[a] = score_max

def getPath(j, i,seq1, seq2,flag):
a = "(" + str(j) + "," + str(i) + ")"
score_res1 = con.get(a)

if score_res1 != None:
res1.append(seq1[i])
res2.append(seq2[j])
if j == 0:
return res2
res_j = getPath(j - 1, i - 1, seq1, seq2, flag)
else:
if j == 0:
return res2
res_j = getPath(j, i - 1, seq1, seq2, flag)
return res_j

def run(seq1, seq2):
score = {}
input_1, input_2 = matrix(seq1, seq2)
for i in range(len(input_1)):
s = "(0," + str(i) + ")"
score[s] = input_1[i]

for i in range(len(input_2)):
s = "(" + str(i) + ",0)"
score[s] = input_2[i]

for j in range(len(seq2)-1):
j += 1
for i in range(len(seq1)-1):
i += 1
result_match = matchBase(seq1[i], seq2[j])
getScore(i, j, result_match, score)

flag = True
res_j = getPath(len(seq2) - 1, len(seq1) - 1, seq1, seq2, flag)
return res_j,score

if __name__ == "__main__":

flag = True
while flag:
introduce()
seq1 = input("Please input long sequence:")
seq2 = input("Please input short sequence:")
seq1 = ">" + seq1.upper()
seq2 = ">" + seq2.upper()
a_i = {}
b_j = {}
con = {}
res1 = []
res2 = []

res_j ,score= run(seq1, seq2)
res_j.reverse()
res1.reverse()
print(" ".join(res1))
print(" ".join(res_j))
print("局部比对得分：%s" %(max(score.values())))
tmp = input("press n to exit：[y/n]")
if tmp.strip() == "n":
flag = False
else:
flag = True

运行结果

*********************************************
Welcome to use short sequence alignment tool!
Author : sunchengquan
input1: long sequence!!!!!
input2: short sequence!!!!!
*********************************************

Please input long sequence:ACGTC
Please input short sequence:CG
C G
C G
局部比对得分：16
press n to exit：[y/n]

部分来源于 ：山东大学 基础医学院 生物信息学课件http://www.crc.sdu.edu.cn/bioinfo