Smith Waterman算法是用于两条序列局部匹配的经典算法。与Needleman Wunsch算法不同的是,Smith Waterman算法侧重于寻找局部最优匹配而非全局匹配。
阅读本文前强烈建议先阅读我的上一篇文章,因为重复的地方笔者没有再写一遍了:
那么在本文中,我只介绍与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)