序列比对：Needleman-Wunsch 算法

序列比对：Needleman-Wunsch 算法

在比对两个序列时，不仅要考虑完全匹配的字符，还要考虑一个序列中的空格或间隙（或者，相反地，要考虑另一个序列中的插入部分）和不匹配，在序列比对中，需要找到最优的比对（最优比对大致是指要将匹配的数量最大化，将空格和不匹配的数量最小化）。如果要更正式些，您可以确定一个分数，为匹配的字符添加分数、为空格和不匹配的字符减去分数。

全局和局部序列比对:

全局序列比对尝试找到两个完整的序列 S1 和 S2 之间的最佳比对。如S1=GCCCTAGCG S2=GCGCAATG 如果设定每个匹配字符为1分，每个空格为-2分，每个不匹配为-1分，则下面的比对就是全局最优比对：S1'=GCCCTAGCG S2'=GCGC_AATG,连字符“_”代表空格。在 S2' 中有五个匹配字符，一个空格（或者反过来说，在 S1' 中有一个插入项），有三个不匹配字符。这样得到的分数是 (5×1) + (1×-2) + (3×-1) = 0，这是能够实现的最佳结果。

Needleman-Wunsch 算法：

该算法用来计算全局比对。它的思路与 LCS 算法相似。这个算法也使用二维表格，一个序列沿顶部展开，一个序列沿左侧展开。而且也能通过以下三个途径到达每个单元格：1.来自上面的单元格，代表将左侧的字符与空格比对。2.来自左侧的单元格，代表将上面的字符与空格比对。3.来自左上侧的单元格，代表与左侧和上面的字符比对（可能匹配也可能不匹配）。该单元格的值来自于一下3个中的最大值：（1）上方的值-2 （2）左边的值-2 （3）如果该单元格所在的行于所在的列对应的字符相等，则为左上值加1，否则为左上值-1。

#coding=utf_8_decode
def Needleman_Wunsch(str1,str2):
#字符串为空则返回
if str1=='' or str2=='':
return  ''
#字符串长度
m=len(str1)
n=len(str2)
#初始化
lcs=[[i*(-2)] for i in range(0,m+1)]
lcs[0]=[j*(-2) for j in range(0,n+1)]
#
for i in range(m):
for j in range(n):
lcs[i+1].append(
max(
lcs[i][j]+(1 if str1[i] == str2[j] else -1),
lcs[i][j+1]-2,
lcs[i+1][j]-2,
)
)
#for i in range(m+1):
#print lcs[i]
i=m-1
j=n-1
common_substr1 = u''
common_substr2 = u''
common_substr1 = u"%s%s" % (str1[i], common_substr1)
common_substr2 = u"%s%s" % (str2[j], common_substr2)
#回溯
while True:
if i == 0 and j == 0:
break
if str1[i] == str2[j]:
if lcs[i-1][j-1]+1>lcs[i-1][j]-2 and lcs[i-1][j-1]+1>lcs[i][j-1]-2:
i = i - 1
j = j -1
common_substr1 = u"%s%s" % (str1[i], common_substr1)
common_substr2 = u"%s%s" % (str2[j], common_substr2)

else:
if lcs[i][j+1] > lcs[i+1][j]:
i = i-1
common_substr1 = u"%s%s" % (str1[i], common_substr1)
common_substr2 = u"%s%s" % ('-', common_substr2)

else:
j = j-1
common_substr1 = u"%s%s" % ('-', common_substr1)
common_substr2 = u"%s%s" % (str2[j], common_substr2)

else:
if lcs[i-1][j-1]+1>lcs[i-1][j]-2 and lcs[i-1][j-1]+1>lcs[i][j-1]-2:
i = i - 1
j = j -1
common_substr1 = u"%s%s" % (str1[i], common_substr1)
common_substr2 = u"%s%s" % (str2[j], common_substr2)

else:
if lcs[i][j+1] > lcs[i+1][j]:
i = i-1
common_substr1 = u"%s%s" % (str1[i], common_substr1)
common_substr2 = u"%s%s" % ('-', common_substr2)

else:
j = j-1
common_substr1 = u"%s%s" % ('-', common_substr1)
common_substr2 = u"%s%s" % (str2[j], common_substr2)
print common_substr1
print common_substr2

Needleman_Wunsch("AGCACACA","ACACTA")