BWA-mem Smith-Waterman 算法

news2024/12/21 19:49:50

Smith-Waterman 算法简介

Smith-Waterman 算法 是一种用于局部序列比对的动态规划算法，它可以在两个序列中找到最优的局部比对，允许处理序列中的插入、缺失（indels） 和错配。该算法最早由 Temple F. Smith 和 Michael S. Waterman 于 1981 年提出，因此得名 Smith-Waterman。

核心思想

Smith-Waterman 算法通过构建一个得分矩阵，从两个序列中找到得分最高的局部比对，即找到序列中最相似的子序列，哪怕这些序列之间存在错配或 indels（插入和缺失）。

应用场景

该算法通常用于：

局部比对：寻找两个序列中最相似的片段。
比对短序列和长参考序列：如在基因组比对中，BWA-MEM 就使用 Smith-Waterman 算法扩展比对。
序列相似性分析：用来识别基因或蛋白质序列中的相似区域。

与全局比对（如 Needleman-Wunsch 算法）的区别

局部比对：Smith-Waterman 专注于找到两个序列中最相似的一部分，而不是从头到尾的全局比对。这特别有用，当你关心的是局部的相似性，而不是整个序列的相似性。
允许 indels 和错配：在比对过程中，Smith-Waterman 允许插入、缺失和错配，但会根据这些事件给出相应的惩罚分数。

算法步骤

初始化得分矩阵：构建一个二维得分矩阵，序列 A 作为行，序列 B 作为列。矩阵的每个元素表示当前子序列的比对得分。
填充得分矩阵：根据相似性（匹配得分）和惩罚（错配或 indels）来填充得分矩阵。计算方式如下：
- 匹配：如果两个字符匹配，增加匹配得分。
- 错配：如果两个字符不匹配，扣除错配罚分。
- 插入或缺失：如果在一个序列中存在插入或缺失，扣除 indel 罚分。
寻找最大得分：找到得分矩阵中的最大值，这表示局部比对的最优解。
回溯路径：从最大值开始，沿着得分矩阵回溯，找到局部比对的路径（即序列片段）。

示例：

假设有两个序列：

序列 A：GACTTAC
序列 B：CGTGAATTCAT

构建一个得分矩阵，按照上述步骤进行比对。Smith-Waterman 算法将会找到最相似的局部子序列，并返回局部比对结果。

优点：

局部比对：能够准确找到最相似的片段，即使序列中有较长的不匹配部分。
灵活处理 indels：适合处理带有插入和缺失的序列。

缺点：

时间复杂度较高：Smith-Waterman 的时间复杂度为 O(m*n)，其中 m 和 n 分别是两个序列的长度，因此对长序列的比对效率较低。

总结

Smith-Waterman 算法是一种强大的局部序列比对工具，它通过动态规划方法找到最优的局部比对区域，允许插入、缺失和错配。这使得它在基因比对、蛋白质序列分析等领域非常有用。

让我们通过这个矩阵来演示 Smith-Waterman 算法的整个比对过程。我们会逐步填写这个矩阵，并通过这个例子理解算法的流程。

输入序列
序列 A：GACTTAC
序列 B：CGTGAATTCAT
我们将通过填充比对网格，找到它们之间的最优局部比对。

步骤 1：初始化矩阵
首先，我们创建一个大小为 8 x 12 的矩阵（包括初始行和列，行列代表序列的字符）。矩阵的每个格子将存储两个字符之间的得分。

在 Smith-Waterman 算法中，矩阵的第一行和第一列初始化为 0。这是因为比对可以从任何位置开始，不必从头比到尾。
整个 Smith-Waterman 算法的比对过程。接下来，我们将继续基于序列 A (GACTTAC) 和 B (CGTGAATTCAT) 填充矩阵，找出最优的局部比对。

1. 初始化矩阵

首先，初始化矩阵第一行和第一列为 0：

      C   G   T   G   A   A   T   T   C   A   T
    +------------------------------------------
G   0   0   0   0   0   0   0   0   0   0   0
A   0
C   0
T   0
T   0
A   0
C   0

2. 定义打分规则

匹配得分：+1
错配罚分：-1
插入/缺失（indel）罚分：-1

3. 填充矩阵

我们继续按照之前的示例填充剩下的部分。

第 1 行（序列 G 与 C, G, T, G, A, A, T, T, C, A, T）

(G, C)：错配，得分 0。
(G, G)：匹配，得分 1。
(G, T)：错配，得分 0。
(G, G)：匹配，得分 1。
(G, A)：错配，得分 0。
(G, A)：错配，得分 0。
(G, T)：错配，得分 0。
(G, T)：错配，得分 0。
(G, C)：错配，得分 0。
(G, A)：错配，得分 0。
(G, T)：错配，得分 0。

因此，第 1 行填充如下：

      C   G   T   G   A   A   T   T   C   A   T
    +------------------------------------------
G   0   0   1   0   1   0   0   0   0   0   0
A   0
C   0
T   0
T   0
A   0
C   0

第 2 行（序列 A 与 C, G, T, G, A, A, T, T, C, A, T）

(A, C)：错配，得分 0。
(A, G)：错配，得分 0。
(A, T)：错配，得分 0。
(A, G)：错配，得分 0。
(A, A)：匹配，得分 1。
(A, A)：匹配，得分 1。
(A, T)：错配，得分 0。
(A, T)：错配，得分 0。
(A, C)：错配，得分 0。
(A, A)：匹配，得分 1。
(A, T)：错配，得分 0。

因此，第 2 行填充如下：

      C   G   T   G   A   A   T   T   C   A   T
    +------------------------------------------
G   0   0   1   0   1   0   0   0   0   0   0
A   0   0   0   0   0   1   1   0   0   0   1
C   0
T   0
T   0
A   0
C   0

第 3 行（序列 C 与 C, G, T, G, A, A, T, T, C, A, T）

(C, C)：匹配，得分 1。
(C, G)：错配，得分 0。
(C, T)：错配，得分 0。
(C, G)：错配，得分 0。
(C, A)：错配，得分 0。
(C, A)：错配，得分 0。
(C, T)：错配，得分 0。
(C, T)：错配，得分 0。
(C, C)：匹配，得分 1。
(C, A)：错配，得分 0。
(C, T)：错配，得分 0。

因此，第 3 行填充如下：

      C   G   T   G   A   A   T   T   C   A   T
    +------------------------------------------
G   0   0   1   0   1   0   0   0   0   0   0
A   0   0   0   0   0   1   1   0   0   0   1
C   0   1   0   0   0   0   0   0   1   0   0
T   0
T   0
A   0
C   0