プログラミングと絵と音楽

コンピューター科学を専攻し、絵と音楽を趣味とするエンジニアのブログです。

動的計画法によるDNA配列の最長共通部分列

動的計画法を学習してから数年、全然使っていないと思って DNA 配列の LCS を求めるプログラムを書きました。
長さ n の文字列 s と長さ m の文字列 t が与えられたとして、部分文字列がより長く一致する空白(ギャップ)の与え方を求めます。
s の i 番目まで s[1..i] と t の j 番目まで t[1..j] の最長の部分文字列の長さ max(LCS) は次のいずれかの最大値になります。

  • s[1..(i-1)] と t[1..j] の max(LCS)
  • s[1..i] と t[1..(j-1)] の max(LCS)
  • s[1..(i-1)] と t[1..(j-1)] の max(LCS) + (if s[i] == t[i] then 1 else 0)

したがってこれを for などのループで繰り返して求め、 dp[n][m] を求めることで結果がわかります。
どれが最大値だったかを記録しておくことで、空白の与え方も求めることが出来ます。

#include <iostream>
#include <string>

struct Cell{
    int score;
    int pre_i;
    int pre_j;

    Cell(int score = 0, int pre_i = 0, int pre_j = 0){
        this->score = score;
        this->pre_i = pre_i;
        this->pre_j = pre_j;
    }
};

int main(){
    std::string s, t;
    Cell **dp;

    std::cin >> s >> t;

    int n = s.size();
    int m = t.size();
    
    dp = new Cell*[n + 1];
    for(int i = 0; i < n + 1; i++){
        dp[i] = new Cell[m + 1];
    }

    for(int i = 0; i < n + 1; i++){
        for(int j = 0; j < m + 1; j++){
            if(i == 0 || j == 0){
                continue;
            }else{
                int s1 = dp[i - 1][j].score;
                int s2 = dp[i][j - 1].score;
                int s3 = (s[i - 1] == t[j - 1] ? 1 : 0) + dp[i - 1][j - 1].score;

                if(s3 >= s1 && s3 >= s2){
                    dp[i][j] = Cell(s3, i - 1, j - 1);
                }else if(s1 >= s2){
                    dp[i][j] = Cell(s1, i - 1, j);
                }else{
                    dp[i][j] = Cell(s2, i, j - 1);
                }
            }
        }
    }

    std::cout << dp[n][m].score << std::endl;

    std::string gs;
    std::string gt;

    while(n * m > 0){
        int i = dp[n][m].pre_i;
        int j = dp[n][m].pre_j;

        if(n != i && m != j){
            gs = s[i] + gs;
            gt = t[j] + gt;
        }else if(n != i){
            gs = s[i] + gs;
            gt = '-' + gt;
        }else if(m != j){
            gs = '-' + gs;
            gt = t[j] + gt;
        }

        n = i;
        m = j;
    }

    if(n == 0){
        gs = std::string(m, '-') + gs;
        gt = t.substr(0, m) + gt;
    }
    if(m == 0){
        gs = s.substr(0, n) + gs;
        gt = std::string(n, '-') + gt;
    }

    std::cout << gs << std::endl;
    std::cout << gt << std::endl;

    for(int i = 0; i < n + 1; i++){
        delete [] dp[i];
    }
    delete [] dp;

    return 0;
}

dp の変数に値を代入するところが動的計画法です。上記の 3 つの最大値を求めて、最大値とそれを与えた i,j の値を代入しています。
その後に文字列を空白付きで復元しています。現在の i,j と前の i,j を比較し、どちらかだけが増えているならばもう片方に空白を、両方増えている場合は両方の文字を足していくことで文字列を復元できます。最後に片方の文字列が残ってしまった場合は、そちら側に残りの文字列を、もう片方にその分の空白を与えることで完全に復元できます。
入力例

AATCGGCCTAGGCAACTCTTTCATCGTTCAGATCAA
AGGGTCGGCCTCGACTCAGCCACCGTGCAA
22
A--ATCGGCCTAGGCAACTCTTTCATCGTTCAGATCAA
AGGGTCGGCCT---CGACTCAGCCACCG-T--G--CAA