蛋白质一级结构全局比对 Needleman-Wunsch 算法的 Python 实现

有关全局对比算法的基本思想可以看这篇文章，这里仅仅用 Python 浅浅实现一下。程序在最后。禁止转载。以下面这两个短肽为例说明一下程序的思路：

SIVK
SSIIVP

在这里我们设计了两个矩阵：matrix1为得分矩阵，用来计算和记录两序列各个配对的得分；matrix2为溯回矩阵，用于记录各配对得分的来源，以及最终对比结果的输出。

Initialize_Matrix部分初始化上面两个矩阵，初始化结果如下：

[[  0.   0.   0.   0.   0.   0.]
 [  0.   0.  -3.  -6.  -9. -12.]
 [  0.  -3.   0.   0.   0.   0.]
 [  0.  -6.   0.   0.   0.   0.]
 [  0.  -9.   0.   0.   0.   0.]
 [  0. -12.   0.   0.   0.   0.]
 [  0. -15.   0.   0.   0.   0.]
 [  0. -18.   0.   0.   0.   0.]] 
 
 [[0. 0. 0. 0. 0. 0.]
 [0. 0. 2. 2. 2. 2.]
 [0. 1. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0. 0.]]

Needleman_Wunsch部分中s代表得分(score)，d代表来源(direction)。这一部分需要接收对比序列seq1第j个残基和seq2第i个残基在矩阵中的位置(i,j)，其匹配得分w以及其上/左/左上方的最终得分(su,sl,sul)，根据这些输入及得分规则算出最终得分s和来源d。

blosum为配对分数的一个简易库（因为嫌麻烦所以只写了需要用到的，如果是核酸序列那就更简单了），并且可以查询得到配对分数。

check_d用来将matrix2来进行溯回以输出对比序列。

而历遍seq1和seq2的部分整合到了main()里面，由于matrix1索引为0的行(列)是留出来写序列的（后来发现好像numpy的matrix似乎写不进字符串），索引为1的行(列)是用来空对比的（或者说是gap对比），所以我们seq1第j个残基和seq2第i个残基 (j, i) 对应到matrix里面就成了 [i+2, j+2]，(后面的说明先咕一会

import numpy as np

def Initialize_Matrix(seq1,seq2,gap):
    matrix1 = np.zeros((len(seq2)+2,len(seq1)+2))
    matrix2 = np.zeros((len(seq2)+2,len(seq1)+2))
    matrix1[1, 1:] = np.linspace(0, (len(matrix1[1, 1:]) - 1) * gap, num=len(matrix1[1, 1:]), endpoint=True)
    matrix1[1:, 1] = np.linspace(0, (len(matrix1[1:, 1]) - 1) * gap, num=len(matrix1[1:, 1]), endpoint=True)
    matrix2[1, 2:] = 2
    matrix2[2:, 1] = 1
    return matrix1,matrix2

def Needleman_Wunsch(gap,w,su,sl,sul):
    s = max(sul+w,su+gap,sl+gap)
    if sul+w == s:
        d = 3
    elif sl+gap == s:
        d = 2
    else:
        d = 1
    return [s,d]

def blosum(res1,res2):
    blosum_bank = {'AS':1,'AI':-1,'AV':0,'AP':-1,'IS':-2,'II':4,
                   'IV':3,'IP':-3,'VS':-2,'VV':4,'VP':-2,
                   'KS':0,'KI':-3,'KV':-2,'KP':-1,'SP':-1,'SS':4}
    w = blosum_bank.get(res1+res2,blosum_bank.get(res2+res1))
    return w

def check_d(matrix2,seq1,seq2):
    i = -1
    j = -1
    alain_seq1 = ''
    alain_seq2 = ''
    while i >= -(len(seq2)+1) and j >= -(len(seq1)+1):
        #print(alain_seq1)
        if matrix2[i][j] == 3:
            alain_seq1 = seq1[j] + alain_seq1
            alain_seq2 = seq2[i] + alain_seq2
            i -= 1
            j -= 1
        elif matrix2[i][j] == 1:
            alain_seq1 = '-' + alain_seq1
            alain_seq2 = seq2[i] + alain_seq2
            i -= 1
        elif matrix2[i][j] == 2:
            alain_seq1 = seq1[j] + alain_seq1
            alain_seq2 = '-' + alain_seq2
            j -= 1
        elif matrix2[i][j] == 0:
            print('done')
            break
        else:
            print('wrong')
            break
    return alain_seq1,alain_seq2


def main():
    seq1 = 'SIVK'
    seq2 = 'SSIIVP'
    gap = -3
    matrix1,matrix2 = Initialize_Matrix(seq1,seq2,gap)
    #print(matrix1,'\n',matrix2)
    for j in range(len(seq1)+len(seq2)):
        for i in range(j+1):
            if j-i <= len(seq1)-1 and i <= len(seq2)-1:
                result = Needleman_Wunsch(gap,blosum(seq1[j-i],seq2[i]),matrix1[i+1][j-i+2],
                                          matrix1[i+2][j-i+1],matrix1[i+1][j-i+1])
                matrix1[i+2][j-i+2] = result[0]
                matrix2[i+2][j-i+2] = result[1]
                #print(j-i-1,i-1,seq1[j-i-1],seq2[i-1])
                #print(matrix1)
            else:
                print('wrong')
                continue
    alain_seq1, alain_seq2 = check_d(matrix2,seq1,seq2)
    print(matrix1,'\n\n',matrix2,'\n\n',alain_seq1,'\n',alain_seq2)

main()

最终输出结果：

done

[[  0.   0.   0.   0.   0.   0.]
 [  0.   0.  -3.  -6.  -9. -12.]
 [  0.  -3.   4.   1.  -2.  -5.]
 [  0.  -6.   1.   2.  -1.  -2.]
 [  0.  -9.  -2.   5.   5.   2.]
 [  0. -12.  -5.   2.   8.   5.]
 [  0. -15.  -8.  -1.   6.   6.]
 [  0. -18. -11.  -4.   3.   5.]] 

 [[0. 0. 0. 0. 0. 0.]
 [0. 0. 2. 2. 2. 2.]
 [0. 1. 3. 2. 2. 2.]
 [0. 1. 3. 3. 3. 3.]
 [0. 1. 1. 3. 3. 2.]
 [0. 1. 1. 3. 3. 2.]
 [0. 1. 1. 1. 3. 3.]
 [0. 1. 1. 1. 1. 3.]] 

 -S-IVK 
 SSIIVP

但是这个方法还存在一些问题，就比如：如果在su,sl,sul计算得到的最总分数里存在两个相同的最大值，也就是存在多条回溯路径，那么该怎么实现多个输出呢，或者说设计一个优先级的判断……

等有兴趣再更