两个不同长度序列的 Needleman-wunsch 算法

Question

我写了一段代码，用 Needleman-Wunsch 算法比对两个序列。对于两个相同长度的序列，它进展顺利。当我使用两个不同长度的序列时，出现错误：string index out of range。谁能帮帮我？

那是我的代码：

import numpy as np 
s1='ACGGAAGGACAGGGAGATAGCGCGGGGAACGAGGAAGAGGGGGGATTCCAGGGAAGGGTAAT' 
s2='TCCGTATAATAGTGCTGTACTAAGCAAATTTATAGTTCTCTAGAAAGTGCCCGC'


M=len(s1)
N=len(s2)
match=+2
gap=-1
mismatch=-2
matrix=np.zeros((N+1,M+1))

#builduing the score matrix 

# 1. putting horizontal and vertical gaps in columns and rows 
for i in range(1,N+1):
    matrix[i,0]=matrix[i-1,0]+gap
for j  in range(1,M+1):
    matrix[0,j]=matrix[j-1,0]+gap
    
# 2. putting in the score/values in the score matrix 
    
for i in range(1,N+1):
       for j in range(1,M+1):
           if s1[i-1] == s2[j-1]:
               score1 = matrix[i-1,j-1] + match
           else:
               score1 = matrix[i-1,j-1] + mismatch
           score2 = matrix[i,j-1] + gap
           score3 = matrix[i-1,j] + gap
           matrix[i,j] = max(score1,score2,score3)
#create a traceback matrix with D=diagonal, V=vertikal, H=horizontal 
trace_mat=np.zeros((N+1,M+1),dtype=str)
for i in range(1,N+1):
    trace_mat[i,0]= 'V'
for j  in range(1,M+1):
    trace_mat[0,j]='H'

非常感谢你的帮助!!!

Answer 1

给出更新后的代码

import numpy as np 
s1='ACGGAAGGACAGGGAGATAGCGCGGGGAACGAGGAAGAGGGGGGATTCCAGGGAAGGGTAAT' 
s2='TCCGTATAATAGTGCTGTACTAAGCAAATTTATAGTTCTCTAGAAAGTGCCCGC'

N=len(s1)
M=len(s2)
match=+2
gap=-1
mismatch=-2
matrix=np.zeros((N+1,M+1))

#builduing the score matrix 

# 1. putting horizontal and vertical gaps in columns and rows 
for i in range(1,N+1):
    matrix[i,0]=matrix[i-1,0]+gap
for j  in range(1,M+1):
    matrix[0,j]=matrix[0,j-1]+gap
    
# 2. putting in the score/values in the score matrix 
    
for i in range(1,N+1):
       for j in range(1,M+1):
           if s1[i-1] == s2[j-1]:
               score1 = matrix[i-1,j-1] + match
           else:
               score1 = matrix[i-1,j-1] + mismatch
           score2 = matrix[i,j-1] + gap
           score3 = matrix[i-1,j] + gap
           matrix[i,j] = max(score1,score2,score3)
#create a traceback matrix with D=diagonal, V=vertikal, H=horizontal 
trace_mat=np.zeros((N+1,M+1),dtype=str)
for i in range(1,N+1):
    trace_mat[i,0]= 'V'
for j  in range(1,M+1):
    trace_mat[0,j]='H'

如果您所做的其他一切都是正确的，这应该可以正常工作。但是，我想指出你做错了什么：

行 'M=len(s1)'、'N=len(s2)' 是我指出的错误之一。这表示 N 指向第二个字符串，M 指向第一个字符串。因此，矩阵定义为matrix[len(s2)][len(s2)]。但是，在嵌套的for循环中，您使用'i'作为矩阵的行来定义s1字符串，而j作为列作为矩阵来定义字符串s2。但是，这是不正确的，因为您将 s1 定义为矩阵的行。
如果我声明s1为矩阵的行，s2为矩阵的列，那么下面的循环要重新考虑，

    for j in range(1,M+1):
        matrix[0,j]=matrix[j-1,0]+gap

我认为不应该这样做，只有修复第 1 点中所述的错误才能解决 string index out of range。但是，我认为这个循环段应该写成如下

    for j  in range(1,M+1):
        matrix[0,j]=matrix[0,j-1]+gap

循环的修改初始化了矩阵的起始垂直值。假设采用动态规划方法，我相信这也应该得到解决。

两个不同长度序列的 Needleman-wunsch 算法

Needleman-wunsch algorithm for two sequences of different length

python

alignment