DNA比对

题目描述

脱氧核糖核酸即常说的 DNA，是一类带有遗传信息的生物大分子。它由 4 种主要的脱氧核苷酸( dAMP、dGMP、dCMT 和 dTMP )通过磷酸二酯键连接而成。这 4 种核苷酸可以分别记为：A、G、C、T。

DNA 携带的遗传信息可以用形如：AGGTCGACTCCA······ 的串来表示。DNA 在转录复制的过程中可能会发生随机的偏差，这才最终造就了生物的多样性。

为了简化问题，我们假设，DNA 在复制的时候可能出现的偏差是（理论上，对每个碱基被复制时，都可能出现偏差）：

1. 漏掉某个脱氧核苷酸。例如把 AGGT 复制成为：AGT

2. 错码，例如把 AGGT 复制成了：AGCT

3. 重码，例如把 AGGT 复制成了：AAGGT

如果某 DNA 串 aa，最少要经过 nn 次出错，才能变为 DNA 串 bb，则称这两个 DNA 串的距离为 nn。

例如：AGGTCATATTCC 与 CGGTCATATTC 的距离为 2。

你的任务是：编写程序，找到两个 DNA 串的距离。

输入描述

用户先输入整数 n (n<100)n (n<100)，表示接下来有 2n2n 行数据。

接下来输入的 2n2n 行每 22 行表示一组要比对的 DNA。（每行数据长度< 104104）

输出描述

程序则输出 nn 行，表示这 nn 组 DNA 的距离。

输入输出样例

示例

输入

3
AGCTAAGGCCTT
AGCTAAGGCCT
AGCTAAGGCCTT
AGGCTAAGGCCTT
AGCTAAGGCCTT
AGCTTAAGGCTT

输出

1
1
2

运行限制

• 最大运行时间：1s
• 最大运行内存: 256M

总通过次数: 360  |  总提交次数: 438  |  通过率: 82.2%

难度: 困难   标签: 2012, 国赛, 动态规划

算法思路：动态规划求解编辑距离

编辑距离（Levenshtein Distance）是衡量两个DNA序列差异的核心指标，通过最小编辑操作次数（插入、删除、替换）将一个序列转换为另一个序列。动态规划是解决该问题的经典方法，其核心思想是通过二维状态表记录子问题最优解。

状态转移示意图

    A   G   C   T   A
  ┌───┬───┬───┬───┬───┐
  │ 0 │ 1 │ 2 │ 3 │ 4 │ ← 初始化：空串→"AGCTA"需4次插入
A ├───┼───┼───┼───┼───┤
  │ 1 │ 0 │ 1 │ 2 │ 3 │ ← "A"→"AG": min(1+1,1+1,0+1)=1
G ├───┼───┼───┼───┼───┤
  │ 2 │ 1 │ 0 │ 1 │ 2 │ ← "AG"→"AGC": min(1+1,0+1,0+0)=0
C ├───┼───┼───┼───┼───┤
  │ 3 │ 2 │ 1 │ 0 │ 1 │ ← 状态转移方程见下文
T └───┴───┴───┴───┴───┘

动态规划四步法

1. ​​状态定义​​
   dp[i][j]：将DNA串s1[0..i-1]转换为s2[0..j-1]的最小操作次数

2. ​​初始化边界​​

   • dp[i][0] = i：删除s1的i个字符
   • dp[0][j] = j：插入s2的j个字符

3. 

4. ​​结果获取​​
   dp[m][n]即为最终编辑距离（m、n为两序列长度）

---

C++代码实现（带滚动数组优化）

#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;

int minDistance(const string& s1, const string& s2) {
    int m = s1.length(), n = s2.length();
    vector<int> prev(n + 1), curr(n + 1);
    
    // 初始化第一行（空串→s2）
    for (int j = 0; j <= n; ++j) 
        prev[j] = j;
    
    for (int i = 1; i <= m; ++i) {
        curr[0] = i;  // 初始化第一列（s1→空串）
        
        for (int j = 1; j <= n; ++j) {
            if (s1[i - 1] == s2[j - 1]) {
                curr[j] = prev[j - 1];  // 字符匹配，无需操作
            } else {
                int del = prev[j] + 1;    // 删除s1[i-1]
                int ins = curr[j - 1] + 1; // 在s1插入s2[j-1]
                int rep = prev[j - 1] + 1; // 替换s1[i-1]为s2[j-1]
                curr[j] = min({del, ins, rep});
            }
        }
        swap(prev, curr);  // 滚动数组：仅保留两行
    }
    return prev[n];
}

int main() {
    int n;
    cin >> n;
    string s1, s2;
    
    for (int i = 0; i < n; ++i) {
        cin >> s1 >> s2;
        cout << minDistance(s1, s2) << endl;
    }
    return 0;
}

代码解析

1. ​​滚动数组优化​​

   • 仅用两个数组prev和curr存储状态
   • 空间复杂度从O(mn)优化到O(n)
   • swap(prev, curr)实现状态行轮换

2. ​​核心操作逻辑​​

   if (s1[i-1] == s2[j-1]) 
       curr[j] = prev[j-1];  // 字符相同，直接继承
   else 
       curr[j] = min(prev[j] + 1,    // 删除s1字符
                     curr[j-1] + 1,  // 插入s2字符
                     prev[j-1] + 1); // 替换字符

3. ​​边界处理​​

   • curr[0] = i：每行首列初始化为删除操作次数
   • prev[j] = j：首行初始化为插入操作次数

---

实例验证

测试用例	计算过程	输出	操作说明
s1="AGCTAAGGCCTT"<br>s2="AGCTAAGGCCT"	dp[12][11] = 1	1	删除s1末尾'T'
s1="AGCTAAGGCCTT"<br>s2="AGGCTAAGGCCTT"	dp[12][13] = 1	1	在s1第4位插入'G'
s1="AGCTAAGGCCTT"<br>s2="AGCTTAAGGCTT"	dp[12][12] = 2	2	替换第5位'A'→'T' + 插入'A'

​​操作分解（第三组）​​：

原始： A G C T A A G G C C T T
步骤1：A G C T T A A G G C C T T  // 替换第5位A→T
步骤2：A G C T T A A G G C T T    // 删除第9位C
结果： A G C T T A A G G C T T

---

注意事项

1. ​​边界极端情况​​

   • 空串处理：minDistance("", "ACG") = 3
   • 相同串：minDistance("ATCG", "ATCG") = 0

2. ​​大字符串优化​​

   • 若两串长度差k较大时，可用O(k)空间的改进算法

3. ​​内存限制​​

   • 滚动数组将内存降至256MB内（10^4长度仅需80KB）

4. ​​字符合法性​​

   • 需确保输入仅含A/G/C/T（题目已保证）

---

测试点设计

​​测试类型​​	输入样例	预期输出	验证点
完全相同序列	s1="ATCG" s2="ATCG"	0	基础功能
空串操作	s1="" s2="AGCT"	4	边界处理
单字符差异	s1="A" s2="T"	1	替换操作
长短串混合	s1="AGG" s2="AG"	1	删除操作
复杂混合操作	s1="ATGC" s2="TACG"	3	综合操作路径
最大长度压力测试	s1=s2=10^4个'A'	0	时间/内存边界

---

优化建议

1. ​​空间极致优化​​

   // 仅用单数组（需保存左上角值）
   for (int i=1; i<=m; ++i) {
       int prev_diag = curr[0];  // 保存dp[i-1][j-1]
       curr[0] = i;
       for (int j=1; j<=n; ++j) {
           int temp = curr[j];
           if (s1[i-1] == s2[j-1]) 
               curr[j] = prev_diag;
           else 
               curr[j] = min({curr[j]+1, curr[j-1]+1, prev_diag+1});
           prev_diag = temp;
       }
   }

2. ​​并行计算优化​​

   • 对每组DNA比对启动独立线程（需<thread>库）

   vector<thread> threads;
   vector<int> results(n);
   for (int i=0; i<n; ++i) {
       threads.emplace_back([&,i]{
           results[i] = minDistance(s1_arr[i], s2_arr[i]);
       });
   }
   for (auto& t : threads) t.join();

3. ​​分支预测优化​​

   #define likely(x) __builtin_expect(!!(x), 1)  // GCC专用
   if (likely(s1[i-1] == s2[j-1])) 
       curr[j] = prev[j-1];