Smith-Waterman 算法（C++实现）

本文实现Smith-Waterman 算法案例，用于局部序列比对。该算法是生物信息学中用于寻找两个 DNA、RNA 或蛋白质序列之间最优局部比对的经典算法，广泛应用于序列相似性分析和功能预测。

---

问题描述

给定两个生物序列 seq1 和 seq2，如何找到它们的最优局部比对，使得比对得分最大化？

---

算法思想

Smith-Waterman 算法的核心思想是动态规划，通过构建一个得分矩阵，逐步计算两个序列的比对得分，并回溯找到最优局部比对路径。与 Needleman-Wunsch 算法不同，Smith-Waterman 算法允许比对从任意位置开始和结束，更适合寻找局部相似性。具体步骤如下：

1. 初始化得分矩阵，其中 dp[i][j] 表示 seq1 的前 i 个字符与 seq2 的前 j 个字符的比对得分。
2. 填充得分矩阵，考虑四种可能的比对操作：
   • 匹配或错配：dp[i-1][j-1] + score(seq1[i], seq2[j])
   • 插入空格：dp[i][j-1] + gap_penalty
   • 删除空格：dp[i-1][j] + gap_penalty
   • 比对从当前位置重新开始：0
3. 回溯得分矩阵，找到最优局部比对路径。

---

C++代码实现

#include <iostream>
#include <vector>
#include <algorithm>

using namespace std;

// 定义得分函数
int match_score(char a, char b) {
    return (a == b) ? 1 : -1; // 匹配得分为 1，错配得分为 -1
}

// Smith-Waterman 算法
pair<int, string> smithWaterman(const string& seq1, const string& seq2, int gap_penalty = -1) {
    int m = seq1.size();
    int n = seq2.size();

    // 初始化得分矩阵
    vector<vector<int>> dp(m + 1, vector<int>(n + 1, 0));
    int max_score = 0; // 记录最大得分
    int max_i = 0, max_j = 0; // 记录最大得分的位置

    // 填充得分矩阵
    for (int i = 1; i <= m; i++) {
        for (int j = 1; j <= n; j++) {
            int match = dp[i - 1][j - 1] + match_score(seq1[i - 1], seq2[j - 1]);
            int insert = dp[i][j - 1] + gap_penalty;
            int del = dp[i - 1][j] + gap_penalty;
            dp[i][j] = max({0, match, insert, del});

            // 更新最大得分及其位置
            if (dp[i][j] > max_score) {
                max_score = dp[i][j];
                max_i = i;
                max_j = j;
            }
        }
    }

    // 回溯找到最优局部比对
    string align1, align2;
    int i = max_i, j = max_j;
    while (i > 0 && j > 0 && dp[i][j] != 0) {
        if (dp[i][j] == dp[i - 1][j - 1] + match_score(seq1[i - 1], seq2[j - 1])) {
            align1 = seq1[i - 1] + align1;
            align2 = seq2[j - 1] + align2;
            i--;
            j--;
        } else if (dp[i][j] == dp[i][j - 1] + gap_penalty) {
            align1 = '-' + align1;
            align2 = seq2[j - 1] + align2;
            j--;
        } else {
            align1 = seq1[i - 1] + align1;
            align2 = '-' + align2;
            i--;
        }
    }

    return {max_score, align1 + "\n" + align2};
}

int main() {
    string seq1 = "GATTACA";
    string seq2 = "GCATGCU";

    auto result = smithWaterman(seq1, seq2);
    cout << "最优局部比对得分: " << result.first << endl;
    cout << "最优局部比对结果: " << endl << result.second << endl;

    return 0;
}

---

关键解析

1. 时间复杂度：O(m * n)，其中 m 和 n 分别是两个序列的长度。
2. 空间复杂度：O(m * n)，用于存储得分矩阵。
3. 适用场景：
   • 局部序列比对。
   • 寻找序列中的功能域或保守区域。

---

输出示例

最优局部比对得分: 2
最优局部比对结果: 
AT
AT

---

总结

Smith-Waterman 算法是生物信息学中用于局部序列比对的经典算法，通过动态规划和回溯找到最优局部比对。