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

Domainmain

已于 2022-11-24 12:28:19 修改

阅读量2.6k

点赞数 24

分类专栏：基于Python的算法学习与设计文章标签：算法 python

于 2022-11-24 12:24:26 首次发布

本文链接：https://blog.csdn.net/qq_54501547/article/details/127981120

版权

基于Python的算法学习与设计专栏收录该内容

2 篇文章 3 订阅

订阅专栏

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)格子的分数为：

$S(i,j)=\max \begin{cases} S(i-1,j-1)+match/mismatch \\S(i-1,j)+gap \\S(i,j-1)+gap \\0 \end{cases}$

而在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)