Needleman-Wunsch 算法

2024 年 8 月 29 日 | 阅读 12 分钟

生物信息学序列比对概述

序列比对是生物信息学中的一项基本任务，它涉及比较 DNA、RNA 或蛋白质等生物序列，以识别它们的相似性和差异。此过程对于理解不同物种之间的进化关系、注释基因以及破译遗传物质内的功能元件至关重要。开创序列比对领域的算法之一是 Needleman-Wunsch 算法。

Needleman-Wunsch 算法的起源

Needleman-Wunsch 算法由 Saul B. Needleman 和 Christian D. Wunsch 于 1970 年开发，旨在解决需要一种方法来全面地全局比对两个序列的问题。在此之前，局部比对方法占主导地位，侧重于识别高度相似的子序列。然而，为了更全面地理解序列关系，科学界需要一种能够比对整个序列的工具。

目的和意义

Needleman-Wunsch 算法的主要目的是通过最大化匹配字符的数量并策略性地引入间隙来解释插入或删除，从而确定两个序列的最佳比对。这种全局比对方法确保序列中的每个位置都被考虑在内，从而提供对其相似性和差异的全面分析。该算法的意义在于其在生物信息学中的应用，它是各种分析的基石，包括数据库搜索、系统发育研究和功能注释。

算法基础

Needleman-Wunsch 算法采用动态规划范例，这是一种强大的计算机科学技术，它将复杂的问题分解为更小的重叠子问题。在序列比对的背景下，动态规划使算法能够系统地评估所有可能的比对，并通过逐步过程确定最佳比对。

1. 动态规划矩阵：得分矩阵

Needleman-Wunsch 算法的核心是动态规划矩阵，通常称为“得分矩阵”。该矩阵是一个二维表，其中每个单元格对应于序列中两个位置之间的特定比对状态。矩阵有 n+1 行和 m+1 列，其中 n 是第一个序列的长度，m 是第二个序列的长度。

该矩阵是迭代填充的，每个单元格的值代表特定比对状态的累积得分。值是根据相邻单元格的值计算的，从而全面表示所有可能的比对。

2. 评分系统：匹配、不匹配和间隙罚分

评分系统是 Needleman-Wunsch 算法的关键组成部分。它涉及为匹配、不匹配和间隙罚分分配分数，影响算法在比对过程中的决策。

匹配得分 (S)：当序列中比对位置的字符匹配时分配的正面得分。此得分反映了字符之间的相似性。

不匹配得分 (D)：当序列中比对位置的字符不匹配时分配的罚分。此得分表示不相似性。

间隙罚分 (G)：为在比对中引入间隙而分配的罚分。间隙用于解释序列中的插入或删除。

动态规划矩阵

Needleman-Wunsch 算法的核心是动态规划矩阵，通常称为“得分矩阵”。该矩阵是一个二维表，其中每个单元格对应于序列中两个位置之间的特定比对状态。该矩阵是迭代填充的，每个单元格中的值根据其相邻单元格的值确定。动态规划矩阵作为算法导航比对可能性路径的路线图。

评分系统

Needleman-Wunsch 算法的一个关键方面是评分系统，它涉及为匹配、不匹配和间隙罚分分配分数。这些分数会影响算法在比对过程中的决策。匹配分数对总比对分数产生积极影响，表示字符之间的相似性，而不匹配分数表示不相似性。间隙罚分用于解释在比对中引入间隙，并对它们的出现进行处罚。仔细校准这些分数可确保算法产生具有生物学意义的比对。

动态规划递推关系

动态规划矩阵中分数的计算依赖于递推关系。此关系定义了如何根据相邻单元格的分数来确定单元格的分数。递推关系包含匹配、不匹配和间隙罚分，为评估比对状态提供了一种全面的方法。通过迭代应用递推关系为每个单元格填充分数来构建动态规划矩阵。

算法工作流程

1. 初始化：奠定基础

算法从初始化动态规划矩阵开始，这是一个二维数组，用于评估比对可能性的基础。此阶段涉及设置矩阵的顶部行和最左列的初始分数。

- 顶部行初始化

矩阵的顶部行对应于第一个序列与第二个序列中的间隙的比对。此行中的每个单元格由左侧单元格的分数（表示与第二个序列中的间隙的比对）加上间隙罚分来确定。

Matrix[0][j]=Matrix[0][j-1]+G

此处，

G 代表间隙罚分。此过程继续到顶部行中的每个单元格，为与第二个序列中间隙数量不断增加的比对建立累积分数。

- 最左列初始化

矩阵的最左列对应于第二个序列与第一个序列中的间隙的比对。此列中的每个单元格由上方单元格的分数（表示与第一个序列中的间隙的比对）加上间隙罚分来确定。

Matrix[i][0]=Matrix[i-1][0]+G

与顶部行类似，此过程继续到最左列中的每个单元格，为与第一个序列中间隙数量不断增加的比对建立累积分数。

初始化阶段创建了一个矩阵，其中包含反映比对中引入间隙的累积罚分的得分。它通过为评估比对可能性提供一个起点来为后续步骤奠定基础。

2. 矩阵填充：系统地评估比对状态

矩阵初始化后，算法根据动态规划递推关系继续填充剩余单元格。递推关系涉及在每个单元格处考虑三种潜在的移动：匹配或不匹配（对角线移动）、第一个序列中的间隙（垂直移动）以及第二个序列中的间隙（水平移动）。

- 动态规划递推关系

递推关系的通用形式如下：（第二个序列中的间隙）

Matrix[i][j]=max

Matrix[i-1][j-1]+S （匹配或不匹配），

Matrix[i-1][j]+G （第一个序列中的间隙），

Matrix[i][j-1]+G （第二个序列中的间隙）

此处，

S 代表匹配得分或不匹配罚分，

G 代表间隙罚分。算法选择导致最高得分的移动，反映了当前单元格的最佳比对状态。

- 迭代矩阵填充

算法迭代遍历矩阵的剩余单元格，根据递推关系系统地计算分数。从左上角开始，逐行移动，每个单元格的分数通过考虑三种潜在的移动来确定。该过程一直持续到达到矩阵的右下角。

随着矩阵的填充，每个单元格都包含一个分数，该分数代表与该位置对应的比对状态的累积最优性。迭代矩阵填充确保了对所有可能的比对状态进行全面评估。

3. 回溯：构建最佳比对

动态规划矩阵完成后，算法执行回溯步骤来构建最佳比对。算法从右下角单元格（代表序列的结尾）开始，沿着得分最高的路径在矩阵中回溯。

- 回溯步骤

对角线移动（匹配或不匹配）：如果当前单元格的分数来自对角线移动，则表示匹配或不匹配。序列中的比对字符将添加到比对结果中。

垂直移动（第一个序列中的间隙）：如果分数来自垂直移动，则表示第一个序列中有间隙。第二个序列中的比对字符和间隙符号将添加到比对结果中。

水平移动（第二个序列中的间隙）：如果分数来自水平移动，则表示第二个序列中有间隙。第一个序列中的比对字符和间隙符号将添加到比对结果中。

- 构建比对

随着算法在矩阵中回溯，比对结果逐步构建。过程一直持续到到达矩阵的左上角，代表序列的开头。

最终比对由匹配的字符和引入的间隙组成，反映了两个序列之间的最佳全局比对。

算法工作流的意义

全局优化：Needleman-Wunsch 算法的工作流程通过动态规划矩阵考虑所有可能的比对状态，确保了全局优化。初始化、矩阵填充和回溯步骤共同实现了识别最大化整体得分的最佳比对。

全面比对：算法工作流程确保对所有可能的比对进行彻底评估，使算法能够适应生物序列中遇到的各种字符。通过考虑匹配得分、不匹配罚分和间隙罚分，算法提供了序列关系的整体视图。

适用于各种序列：算法的系统性使其适用于不同长度和组成的序列。它广泛用于生物信息学中比较 DNA、RNA 和蛋白质序列，从而提供对进化关系和功能注释的有价值的见解。

以下是 C++ 中的程序

#include <iostream>
#include <vector>

using namespace std;

const int matchScore = 2;
const int mismatchPenalty = -1;
const int gapPenalty = -1;

void printMatrix(const vector<vector<int>> &matrix) {
    for (const auto &row : matrix) {
        for (int val : row) {
            cout << val << "\t";
        }
        cout << endl;
    }
}

void needlemanWunsch(const string &sequence1, const string &sequence2) {
    int m = sequence1.size();
    int n = sequence2.size();

    // Initialize the dynamic programming matrix
    vector<vector<int>> dpMatrix(m + 1, vector<int>(n + 1, 0));

    // Initialize the first row and column with gap penalties
    for (int i = 1; i <= m; ++i) {
        dpMatrix[i][0] = dpMatrix[i - 1][0] + gapPenalty;
    }

    for (int j = 1; j <= n; ++j) {
        dpMatrix[0][j] = dpMatrix[0][j - 1] + gapPenalty;
    }

    // Fill the dynamic programming matrix
    for (int i = 1; i <= m; ++i) {
        for (int j = 1; j <= n; ++j) {
            int match = dpMatrix[i - 1][j - 1] + (sequence1[i - 1] == sequence2[j - 1] ? matchScore : mismatchPenalty);
            int gap1 = dpMatrix[i - 1][j] + gapPenalty;
            int gap2 = dpMatrix[i][j - 1] + gapPenalty;

            dpMatrix[i][j] = max({match, gap1, gap2});
        }
    }

    // Print the dynamic programming matrix
    cout << "Dynamic Programming Matrix:" << endl;
    printMatrix(dpMatrix);
    cout << endl;

    // Backtrack to find the optimal alignment
    int i = m, j = n;
    string alignedSequence1, alignedSequence2;

    while (i > 0 || j > 0) {
        if (i > 0 && j > 0 && dpMatrix[i][j] == dpMatrix[i - 1][j - 1] + (sequence1[i - 1] == sequence2[j - 1] ? matchScore : mismatchPenalty)) {
            alignedSequence1 = sequence1[i - 1] + alignedSequence1;
            alignedSequence2 = sequence2[j - 1] + alignedSequence2;
            --i;
            --j;
        } else if (i > 0 && dpMatrix[i][j] == dpMatrix[i - 1][j] + gapPenalty) {
            alignedSequence1 = sequence1[i - 1] + alignedSequence1;
            alignedSequence2 = '-' + alignedSequence2;
            --i;
        } else {
            alignedSequence1 = '-' + alignedSequence1;
            alignedSequence2 = sequence2[j - 1] + alignedSequence2;
            --j;
        }
    }

    // Print the optimal alignment
    cout << "Optimal Alignment:" << endl;
    cout << alignedSequence1 << endl;
    cout << alignedSequence2 << endl;
}

int main() {
    string sequence1 = "AGTAC";
    string sequence2 = "GAGC";

    cout << "Sequence 1: " << sequence1 << endl;
    cout << "Sequence 2: " << sequence2 << endl << endl;

    needlemanWunsch(sequence1, sequence2);

    return 0;
}

输出

Sequence 1: AGTAC
Sequence 2: GAGC

Dynamic Programming Matrix:
0       -1      -2      -3      -4      -5
-1      2       1       0       -1     -2
-2      1       4       3       2

说明

1. 包含语句

#include <iostream>：包含用于输入输出操作的输入输出流库。
#include <vector>：包含用于动态数组实现的向量库。

2. 常量声明

const int matchScore = 2;：定义匹配得分常量。
const int mismatchPenalty = -1;：定义不匹配罚分常量。
const int gapPenalty = -1;：定义间隙罚分常量。

3. 矩阵打印函数

void printMatrix(const vector<vector<int>> &matrix)：声明一个函数来打印整数二维矩阵。
for (const auto &row : matrix)：遍历矩阵中的每一行。
for (int val : row)：遍历行中的每个元素。
cout << val << "\t";：打印每个矩阵元素，后跟一个制表符。
cout << endl;：打印完每一行后换行。

4. Needleman-Wunsch 函数初始化

void needlemanWunsch(const string &sequence1, const string &sequence2)：声明 Needleman-Wunsch 函数，该函数接受两个输入序列。
int m = sequence1.size();：获取第一个序列的长度。
int n = sequence2.size();：获取第二个序列的长度。

5. 矩阵初始化

初始化一个尺寸为 (m + 1) x (n + 1) 的二维向量 dpMatrix，并将所有元素初始化为 0。

6. 使用间隙罚分初始化矩阵

使用累积间隙罚分初始化第一行。
使用累积间隙罚分初始化第一列。

7. 矩阵填充循环

嵌套循环遍历动态规划矩阵，根据 Needleman-Wunsch 递推关系计算得分。

8. 打印动态规划矩阵

使用先前定义的 printMatrix 函数打印动态规划矩阵。

9. 回溯循环

回溯遍历动态规划矩阵以查找最佳比对。

10. 打印最佳比对

打印回溯步骤中找到的最佳比对。

11. 主函数

初始化两个示例序列。
打印序列。
使用输入序列调用 needlemanWunsch 函数。
返回 0 表示程序成功执行。

时间复杂度分析

Needleman-Wunsch 算法的时间复杂度主要受其执行的两个主要步骤的影响：矩阵填充和回溯。

矩阵填充

负责矩阵填充的嵌套循环迭代动态规划矩阵中的每个单元格，在每一步考虑三种潜在的移动（匹配/不匹配、第一个序列中的间隙以及第二个序列中的间隙）。循环以 m 行和 n 列运行，其中 m 和 n 是输入序列的长度。

在每次迭代中，都会执行常数时间操作来根据递推关系计算得分，从而导致矩阵填充步骤的总体时间复杂度为 O(m⋅n)。

回溯

回溯步骤涉及遍历动态规划矩阵以追踪最佳比对路径。在最坏的情况下，算法可能需要从右下角回溯到左上角，需要 m+n 步。每次回溯步骤都涉及常数时间操作，为总体时间复杂度贡献 O(m+n)。

因此，时间复杂度中的主导因素是矩阵填充步骤，Needleman-Wunsch 算法的时间复杂度为 O(m⋅n)。

空间复杂度

Needleman-Wunsch 算法的空间复杂度主要由动态规划矩阵的存储需求决定。此外，还为索引和分数等变量使用了恒定的空间。

动态规划矩阵

动态规划矩阵是一个二维数组，其尺寸为 (m+1)×(n+1)，其中 m 和 n 是输入序列的长度。矩阵中的每个单元格都存储一个整数值，表示特定比对状态的累积得分。因此，归因于动态规划矩阵的空间复杂度为 O(m⋅n)。

附加空间

除了矩阵之外，该算法还为变量（如 matchScore、mismatchPenalty、gapPenalty、循环索引和临时变量）使用了恒定的附加空间。这种恒定空间复杂度为 O(1)。

总体空间复杂度

Needleman-Wunsch 算法的总体空间复杂度是动态规划矩阵的空间需求和恒定附加空间的总和。因此，该算法的总体空间复杂度为 O(m⋅n)+O(1)，简化为 O(m⋅n)。总之，Needleman-Wunsch 算法具有二次时间复杂度 (O(m⋅n))，主要由矩阵填充步骤驱动，以及线性空间复杂度 (O(m⋅n))，这归因于动态规划矩阵的存储需求。此分析提供了对算法的可伸缩性和效率的见解，使其适用于生物信息学应用中不同长度序列的比较。

在生物信息学中的应用

Needleman-Wunsch 算法在生物信息学和计算生物学中有着广泛的应用。

数据库搜索：在基因组学中，研究人员使用该算法将查询序列与已知序列数据库进行比较，从而识别同源区域并推断功能相似性。

系统发育分析：该算法通过比对来自不同物种的序列来帮助构建系统发育树。这有助于研究进化关系和分化。

功能注释：理解基因和蛋白质的功能是分子生物学的一个关键方面。Needleman-Wunsch 算法通过揭示序列中的保守区域来帮助注释功能元件。

结论

Needleman-Wunsch 算法是生物信息学领域的开创性成就，为序列比对提供了一种强大而通用的工具。其动态规划方法与精心校准的评分系统相结合，能够识别生物序列之间的最佳全局比对。该算法的影响超出了简单的比对；它构成了各种分析的基础，这些分析有助于我们理解基因密码中错综复杂的相互关系。随着技术的发展和生物数据的量不断增加，Needleman-Wunsch 算法仍然是生物信息学家工具包中的基石，在揭示生命序列中编码的奥秘方面发挥着关键作用。

下一个主题C++ 中的二元 GCD 算法

← 上一个下一个 →

Needleman-Wunsch 算法

生物信息学序列比对概述