两个矩阵
一个是惩罚矩阵：
image.png
#生成惩罚矩阵
A=np.ones([5,5],dtype=np.int16)*-2+np.eye(5,dtype=np.int16)*5
A[0,0]=0
name_list=["_","A","C","G","T"]
punish_matrix=pd.DataFrame(A,columns=name_list,index=name_list)
这个的作用是为打分提供个参考，将每个碱基的对应情况进行打分
打分矩阵：
我们看下原理图：
image.png
s来自惩罚矩阵，w为惩罚分，H来自打分矩阵
假设我有两条序列：
str_one="ACTTTCA"
str_two="AGGCTTA"
str_one="ACTTTCA"
str_two="AGGCTTA"
index=["_"]+[i for i in str_one]
columns=["_"]+[i for i in str_two]
score_matrix=pd.DataFrame(np.zeros([len(str_one)+1,len(str_two)+1]),index=index,columns=columns)
# 设置行列名，以-，和你的序列为行列名，一条为行名，另一条为列名
punish=-3 #惩罚分
for i in range(len(str_one)+1):
for j in range(len(str_two)+1):
if i==0 or j ==0:
score_matrix.iloc[i,j]=0+punish*i+punish*j
else:
insert=score_matrix.iloc[i,j-1]+punish
delect=score_matrix.iloc[i-1,j]+punish
match=score_matrix.iloc[i-1,j-1]+punish_matrix.loc[str_one[i-1],str_two[j-1]]
score_matrix.iloc[i, j]=max(insert,delect,match) #选择最大的值作为填充
image.png
对于这个矩阵来说，我们要找到最底角的值开始回溯
对于这个矩阵来说是3.0
摘自知乎
那么我们先找到右下角的那个值，通过右下角这个值进行路径的回溯，即获得行列的全长即可：
i=len(str_one)
j=len(str_two)
while (i > 0 or j > 0) :
if (i > 0 and j > 0 and score_matrix.iloc[i, j]==score_matrix.iloc[i-1,j-1]+punish_matrix.loc[str_one[i-1],str_two[j-1]]):
s1 += str_one[i - 1]
s2 += str_two[j - 1]
i -= 1
j -= 1
elif (j > 0 and score_matrix.iloc[i, j]==score_matrix.iloc[i,j-1]+punish):
s1 += '-'
s2 += str_two[j - 1]
j -= 1
elif (i > 0 and score_matrix.iloc[i, j]==score_matrix.iloc[i-1,j]+punish):
s1 += str_one[i - 1]
s2 += '-'
i -= 1
全部代码
import numpy as np
import pandas as pd
A=np.ones([5,5],dtype=np.int16)*-2+np.eye(5,dtype=np.int16)*5
A[0,0]=0
name_list=["_","A","C","G","T"]
punish_matrix=pd.DataFrame(A,columns=name_list,index=name_list)
str_one="ACTTTCA"
str_two="AGGCTTA"
s1 = ''
s2 = ''
index=["_"]+[i for i in str_one]
columns=["_"]+[i for i in str_two]
score_matrix=pd.DataFrame(np.zeros([len(str_one)+1,len(str_two)+1]),index=index,columns=columns)
punish=-3
for i in range(len(str_one)+1):
for j in range(len(str_two)+1):
if i==0 or j ==0:
score_matrix.iloc[i,j]=0+punish*i+punish*j
else:
insert=score_matrix.iloc[i,j-1]+punish
delect=score_matrix.iloc[i-1,j]+punish
match=score_matrix.iloc[i-1,j-1]+punish_matrix.loc[str_one[i-1],str_two[j-1]]
score_matrix.iloc[i, j]=max(insert,delect,match)
i=len(str_one)
j=len(str_two)
# i，j 均为最大值，即满足矩阵右下角的元素的角标
while (i > 0 or j > 0) :
if (i > 0 and j > 0 and score_matrix.iloc[i, j]==score_matrix.iloc[i-1,j-1]+punish_matrix.loc[str_one[i-1],str_two[j-1]]):
s1 += str_one[i - 1]
s2 += str_two[j - 1]
i -= 1
j -= 1
elif (j > 0 and score_matrix.iloc[i, j]==score_matrix.iloc[i,j-1]+punish):
s1 += '-'
s2 += str_two[j - 1]
j -= 1
elif (i > 0 and score_matrix.iloc[i, j]==score_matrix.iloc[i-1,j]+punish):
s1 += str_one[i - 1]
s2 += '-'
i -= 1
print(score_matrix)
print(s1[::-1] + '\n' + s2[::-1])
结果是：
image.png
考虑空间：
如果序列太长，那么空间复杂度会大大加大，我们可以采用中点法，将长序列分割成两段，然后进行比对
image.png
如上图，它的打分矩阵就不需要完全填充完，分成两段后只需要将蓝色部分的打分举证填充完毕即可
image.png
即绿色部分
import numpy as np
import pandas as pd
import argparse
A=np.ones([5,5],dtype=np.int16)*-2+np.eye(5,dtype=np.int16)*5
A[0,0]=0
name_list=["_","A","C","G","T"]
punish_matrix=pd.DataFrame(A,columns=name_list,index=name_list)
def mid_align(str_one,str_two):
s1 = ''
s2 = ''
index = ["_"] + [i for i in str_one]
columns = ["_"] + [i for i in str_two]
score_matrix = pd.DataFrame(np.zeros([len(str_one) + 1, len(str_two) + 1]), index=index, columns=columns)
punish = -3
for i in range(len(str_one) + 1):
for j in range(len(str_two) + 1):
if i == 0 or j == 0:
score_matrix.iloc[i, j] = 0 + punish * i + punish * j
else:
insert = score_matrix.iloc[i, j - 1] + punish
delect = score_matrix.iloc[i - 1, j] + punish
match = score_matrix.iloc[i - 1, j - 1] + punish_matrix.loc[str_one[i - 1], str_two[j - 1]]
score_matrix.iloc[i, j] = max(insert, delect, match)
i = len(str_one)
j = len(str_two)
# i，j 为矩阵最大的，即满足矩阵“右下角"元素
while (i > 0 or j > 0):
if (i > 0 and j > 0 and score_matrix.iloc[i, j] == score_matrix.iloc[i - 1, j - 1] + punish_matrix.loc[
str_one[i - 1], str_two[j - 1]]):
s1 += str_one[i - 1]
s2 += str_two[j - 1]
i -= 1
j -= 1
elif (j > 0 and score_matrix.iloc[i, j] == score_matrix.iloc[i, j - 1] + punish):
s1 += '-'
s2 += str_two[j - 1]
j -= 1
elif (i > 0 and score_matrix.iloc[i, j] == score_matrix.iloc[i - 1, j] + punish):
s1 += str_one[i - 1]
s2 += '-'
i -= 1
list = [s1[::-1],s2[::-1]]
return list
if __name__ == "__main__":
parser = argparse.ArgumentParser(description="Multi-value dictionary")
parser.add_argument("-s1", "--s1", required=True, type=str, help="First sequence")
parser.add_argument("-s2", "--s2", required=True, type=str, help="Second sequence")
Args = parser.parse_args()
A = np.ones([5, 5], dtype=np.int16) * -2 + np.eye(5, dtype=np.int16) * 5
A[0, 0] = 0
name_list = ["_", "A", "C", "G", "T"]
punish_matrix = pd.DataFrame(A, columns=name_list, index=name_list)
str_one = Args.s1
str_two = Args.s2
n = len(str_one)
m = len(str_two)
i_mid = n // 2 #取中点，然后将序列分为两段
string_one = str_one[0:i_mid] #分割序列
string_two = str_two[0:i_mid]
string_three = str_one[i_mid:len(str_one)]
string_four = str_two[i_mid:len(str_two)]
print(mid_align(string_one,string_two)[0] + mid_align(string_three,string_four)[0])
print(mid_align(string_one,string_two)[1] + mid_align(string_three,string_four)[1]) #拼接序列
结果为：
image.png
