【序列匹配算法】Smith Waterman算法 v2.0

Smith Waterman算法是用于两条序列局部匹配的经典算法。与Needleman Wunsch算法不同的是，Smith Waterman算法侧重于寻找局部最优匹配而非全局匹配。

阅读本文前强烈建议先阅读我的上一篇文章，因为重复的地方笔者没有再写一遍了：

【序列匹配算法】Needleman Wunsch算法 v2.0_Domainmain的博客-CSDN博客

那么在本文中，我只介绍与Needleman Wunsch不同的地方，其余相同的地方可以见上述文章。（想直接看完整代码请跳到最后）

首先再点明我们要匹配的序列：

序列1 AACGTACTCAAGTCT
序列2 TCGTACTCTAACGAT

第一个要修改的地方就是打分系统。在Needleman中，笔者举例的计分方式如下：（可见上一篇文章）

从左边来的格子算出黄色格子的得分为-2-2=-4；
从上边来的格子算出黄色格子的得分为-2-2=-4；
从对角线格子计算，由于黄色格子对应的A与T不匹配，所以分数为0-6=-6

这三个分数中取最大值-4作为黄色格子的值。

那么在Smith Waterman中，求最大值时还需要与0作比较。即第(i,j)格子的分数为：

而在Needleman中没有0这一项。所以在初始化矩阵的时候，就应该只有全0的矩阵而不应该有负数：

def ini_matrix(l1: list, l2: list, gap):"""初始化罚分矩阵l1: 序列一列表l2: 序列二列表gap: 空位得分"""# 获取序列长度构建初始矩阵n1 = len(l1)n2 = len(l2)score_matrix = np.zeros((n1+1, n2+1))# 返回矩阵return score_matrix

在打分时，也要注意和0比较：

# %% 计分
def score(matrix: np.array, l1: list, l2: list, match, mismatch, gap):"""计算矩阵得分matrix: 初始化的矩阵l1: 序列一列表l2: 序列二列表match: 匹配得分mismatch: 不匹配时得分gap: 空位得分"""# 循环计分for i in range(1, len(l1)+1):for j in range(1, len(l2)+1):# 计算三类分值from_left = matrix[i][j - 1] + gap  # 从左到右空位from_above = matrix[i - 1][j] + gap  # 从上到下空位if l1[i-1] == l2[j-1]:  # 对角线from_diag = matrix[i - 1][j - 1] + match  # 匹配else:from_diag = matrix[i - 1][j - 1] + mismatch  # 不匹配# 比较并赋分matrix[i][j] = max(from_left, from_above, from_diag, 0)return matrix

第二个就是回溯的起点和终点。在Needleman中，回溯的起点在矩阵右下角，终点在矩阵左上角。而Smith Waterman的回溯起点在矩阵中最大数值的位置，回溯过程中遇到左、左上、上这三个格子均为0时停止：

如图所示，在这个例子中，右侧标记的93为矩阵最大值，就是回溯起点；左侧标记的9就是一种可能的回溯终点，可以发现它的回溯邻居（左、上、左上）都为0。

回溯的方法同笔者上一篇文章，只用修改起点和条件即可：

# %% 回溯
def trace_back(res: np.array, l1: list, l2: list, match, mismatch, gap):"""回溯矩阵获得匹配结果索引res: 结果矩阵l1: 序列一列表l2: 序列二列表match: 匹配得分mismatch: 不匹配时得分gap: 空位得分"""path = []    # 最终所有路径# 找到矩阵中的最大值并入栈i = np.where(res == np.max(res))[0][0]j = np.where(res == np.max(res))[1][0]m_stack = [(i, j)]  # 主栈a_stack = []  # 辅助栈while m_stack:  # 当主栈非空时# 检查是否到终点row = m_stack[-1][0]col = m_stack[-1][1]if res[row-1][col-1] == res[row-1][col] == res[row][col-1] == 0:# 所有邻居都是0，到终点，存储索引路径，依次弹出栈顶path.append(m_stack.copy())m_stack.pop()else:if len(m_stack) > len(a_stack):# 检查主辅栈长度是否一致，不一致则添加新邻居a_stack.append([])if l1[row - 1] == l2[col - 1] and res[row][col] == res[row-1][col-1]+match:a_stack[-1].append((row-1, col-1))elif res[row][col] == res[row-1][col-1]+mismatch:a_stack[-1].append((row-1, col-1))if res[row][col] == res[row-1][col]+gap:a_stack[-1].append((row-1, col))if res[row][col] == res[row][col-1]+gap:a_stack[-1].append((row, col-1))# 检测辅栈栈顶列表是否为空，不空则可以访问邻居elif a_stack[-1] != []:m_stack.append(a_stack[-1].pop())# 辅助栈为空，则同时退栈一个elif a_stack[-1] == []:a_stack.pop()m_stack.pop()return path

其余的代码块就与Needleman几乎一致了。浅浅展示一下结果吧：

那么结合本次修改的地方，调整主函数的输出，最终完整的代码如下：

# %% 导入包
import numpy as np
import matplotlib.pyplot as plt
import time# %% 导入序列
def file(path: str):"""导入文件中的序列path: 文件路径"""with open(path) as file_obj:seq = file_obj.read()# 返回序列return seq# %% 初始化
def str_to_list(a: str, b: str):"""将序列字符串转换为单个字符列表a: 序列一b: 序列二"""return list(a), list(b)def ini_matrix(l1: list, l2: list, gap):"""初始化罚分矩阵l1: 序列一列表l2: 序列二列表gap: 空位得分"""# 获取序列长度构建初始矩阵n1 = len(l1)n2 = len(l2)score_matrix = np.zeros((n1+1, n2+1))# 返回矩阵return score_matrix# %% 计分
def score(matrix: np.array, l1: list, l2: list, match, mismatch, gap):"""计算矩阵得分matrix: 初始化的矩阵l1: 序列一列表l2: 序列二列表match: 匹配得分mismatch: 不匹配时得分gap: 空位得分"""# 循环计分for i in range(1, len(l1)+1):for j in range(1, len(l2)+1):# 计算三类分值from_left = matrix[i][j - 1] + gap  # 从左到右空位from_above = matrix[i - 1][j] + gap  # 从上到下空位if l1[i-1] == l2[j-1]:  # 对角线from_diag = matrix[i - 1][j - 1] + match  # 匹配else:from_diag = matrix[i - 1][j - 1] + mismatch  # 不匹配# 比较并赋分matrix[i][j] = max(from_left, from_above, from_diag, 0)return matrix# %% 回溯
def trace_back(res: np.array, l1: list, l2: list, match, mismatch, gap):"""回溯矩阵获得匹配结果索引res: 结果矩阵l1: 序列一列表l2: 序列二列表match: 匹配得分mismatch: 不匹配时得分gap: 空位得分"""path = []    # 最终所有路径# 找到矩阵中的最大值并入栈i = np.where(res == np.max(res))[0][0]j = np.where(res == np.max(res))[1][0]m_stack = [(i, j)]  # 主栈a_stack = []  # 辅助栈while m_stack:  # 当主栈非空时# 检查是否到终点row = m_stack[-1][0]col = m_stack[-1][1]if res[row-1][col-1] == res[row-1][col] == res[row][col-1] == 0:# 所有邻居都是0，到终点，存储索引路径，依次弹出栈顶path.append(m_stack.copy())m_stack.pop()else:if len(m_stack) > len(a_stack):# 检查主辅栈长度是否一致，不一致则添加新邻居a_stack.append([])if l1[row - 1] == l2[col - 1] and res[row][col] == res[row-1][col-1]+match:a_stack[-1].append((row-1, col-1))elif res[row][col] == res[row-1][col-1]+mismatch:a_stack[-1].append((row-1, col-1))if res[row][col] == res[row-1][col]+gap:a_stack[-1].append((row-1, col))if res[row][col] == res[row][col-1]+gap:a_stack[-1].append((row, col-1))# 检测辅栈栈顶列表是否为空，不空则可以访问邻居elif a_stack[-1] != []:m_stack.append(a_stack[-1].pop())# 辅助栈为空，则同时退栈一个elif a_stack[-1] == []:a_stack.pop()m_stack.pop()return path# %% 将回溯结果转化为碱基
def trace_to_base(match: list, l1: list, l2: list): """将回溯获得的索引转化为对应的碱基并打印match: 回溯索引元组列表l1: 序列一列表l2: 序列二列表"""# 初始化seq_match1 = [l1[match[0][0] - 1]]seq_match2 = [l2[match[0][1] - 1]]connect = []  # 连接符if seq_match1 == seq_match2:connect.append("|")else:connect.append("*")# 依次转换for index in range(1, len(match)):if match[index][0] != match[index - 1][0] and match[index][1] != match[index-1][1]:seq_match1.append(l1[match[index][0] - 1])seq_match2.append(l2[match[index][1] - 1])if l1[match[index][0] - 1] == l2[match[index][1] - 1]:connect.append("|")else:connect.append("*")elif match[index][0] == match[index - 1][0]:seq_match1.append("-")seq_match2.append(l2[match[index][1] - 1])connect.append(" ")elif match[index][1] == match[index - 1][1]:seq_match1.append(l1[match[index][0] - 1])seq_match2.append("-")connect.append(" ")# 转换为字符串seq_match1 = "".join(seq_match1)seq_match2 = "".join(seq_match2)connect = "".join(connect)return seq_match1, seq_match2, connect# %% 可视化
def matrix_visual(res_matrix: np.array, l1: list, l2: list, index: list, plot_val=False):"""可视化得分矩阵res_matrix: 得分矩阵l1: 序列一列表l2: 序列二列表index: 回溯的索引列表plot_val: 逻辑值,是否将矩阵数值绘制在图中,默认为False"""# 处理index列表path = set()for i in index:path = set(i).union(path)# 可视化得分矩阵fig, ax = plt.subplots()visual_matrix = plt.imshow(res_matrix)plt.colorbar(visual_matrix)# 修改轴l1.insert(0, '0')l2.insert(0, '0')ax.set_xticks(np.arange(0, len(l1)), minor=False)ax.set_yticks(np.arange(0, len(l2)), minor=False)ax.xaxis.set_ticks_position('top')ax.xaxis.set_label_position('top')ax.set_xticklabels(l1)ax.set_yticklabels(l2)plt.xlabel("Sequence 1")plt.ylabel("Sequence 2")# 是否画出标记if plot_val:for i in range(len(l2)):for j in range(len(l1)):if (j, i) in path:ax.text(i-0.35, j+0.1,res_matrix[j][i], fontsize=6, color="#FF49F5")else:ax.text(i-0.35, j+0.1, res_matrix[j][i], fontsize=6)# 紧密排布plt.tight_layout()plt.show()# %% 主函数
def sw(seq1: str, seq2: str, match, mismatch, gap, plot=False, plot_val=False):"""主函数seq1: 第一条序列seq2: 第二条序列match: 匹配得分mismatch: 不匹配时得分gap: 空位得分plot: 逻辑值,是否画出得分矩阵,默认不画出plot_val: 逻辑值,是否将矩阵数值绘制在图中,默认为False """# 开始计时print("\n开始匹配......")start = time.time()# 获取列表l1, l2 = str_to_list(seq1, seq2)# 初始化打分矩阵score_matrix = ini_matrix(l1, l2, gap)# 为矩阵赋分res_matrix= score(score_matrix, l1, l2, match, mismatch, gap)# 回溯path = trace_back(res_matrix, l1, l2, match, mismatch, gap)# 停止计时stop = time.time()print("匹配结束!\n********************")# 打印最终结果print("打印比对结果......")for i in path:i = i[::-1]seq_match1, seq_match2, connect = trace_to_base(i, l1, l2)print(seq_match1)print(connect)print(seq_match2)print("--------------------")print("匹配结果共计%s条, 匹配分数为: %s"%(len(path), np.max(res_matrix)))print("\n计算耗时为: ", stop - start, "秒")if plot:print("\n可视化得分矩阵......")matrix_visual(res_matrix, l1, l2, path, plot_val)print("\n结束!!!")# %% 应用
seq1 = "AACGTACTCAAGTCT"
seq2 = "TCGTACTCTAACGAT"
sw(seq1, seq2, match=9, mismatch=-3, gap=-2, plot=True, plot_val=True)

【序列匹配算法】Smith Waterman算法 v2.0相关推荐

语音识别：时间序列的Smith–Waterman对齐算法
一.基本概念 Smith-Waterman 算法执行局部序列比对:也就是说,用于确定两串核酸序列或蛋白质序列之间的相似区域. Smith-Waterman 算法不是查看整个序列,而是比较所有可能长度的 ...
matlab fft函数说明_【V2.0更新】基于FFT算法的MTALAB傅里叶级数3D可视化
最近这份代码受到很多朋友关注,在此一并感谢!由于之前写1.0版本代码时的时间有限,当时的知识储备也不甚完善,所以只是做出了基本功能. 近期对这份代码进行了更新,优化了代码结构,并加入了FFT等算法让其 ...
【数据结构与算法】二叉堆V2.0的Java实现
更新说明我们在此前已经编写过简单版的二叉大根堆V1.0,这次,换成二叉小根堆,命名为二叉堆V2.0. 大家也知道,堆是完全二叉树,存储方式借助动态数组实现顺序存储,依赖于父子结点之间的index关系 ...
TensorFlow v2.0实现Word2Vec算法
使用TensorFlow v2.0实现Word2Vec算法计算单词的向量表示,这个例子是使用一小部分维基百科文章来训练的. 更多信息请查看论文: Mikolov, Tomas et al. " ...
为什么神经元有数千个突触，一个新皮质中的序列记忆理论(HTM算法基础)
为什么神经元有数千个突触,一个新皮质中的序列记忆理论(HTM算法基础) Jeff Hawkins* and Subutai Ahmad Numenta, Inc., Redwood City, CA, ...
钢七连C2游戏编程具体知识点清单V2.0
以下46项知识的讲解,不再上传.一个学生多会了46项知识.技术.案例项目,是否提高了职业竞争力.青出于蓝,学生成长为第二个老师. 游戏编程的具体知识点清单V2.0 1 腾讯天美工作室游戏资源游戏图书 ...
MMDetection V2.0发布！速度精度全面提升，现有检测框架最优
本文授权转自知乎作者陈恺,https://zhuanlan.zhihu.com/p/145084667.未经作者许可,不得二次转载. MMDetection V1.0 版本发布以来,我们收到了很多用户 ...
MMDetection V2.0：更快更强的通用目标检测平台
MMDetection V1.0 版本发布以来,我们收到了很多用户的反馈,其中有不少有价值的建议,同时也有很多开发者贡献代码,和我们一起不断完善.经过 2 个月的酝酿,再经过 3 个月的开发和打磨(也 ...
AI大事件 | OpenAI员工离职创立机器人新公司，spaCy v2.0.0发布
呜啦啦啦啦啦大家好呀,又到了本周的AI大事件时间了.过去的一周中AI圈都发生了什么?大佬们互撕了哪些问题?研究者们发布了哪些值得一读的论文?又有哪些开源的代码和数据库可以使用了?文摘菌带你盘点过去一周 ...

【序列匹配算法】Smith Waterman算法 v2.0

【序列匹配算法】Smith Waterman算法 v2.0相关推荐

最新文章

热门文章