Sequence Alignment

</pre><p class="p1"><pre name="code" class="cpp">#include <iostream>#include <vector>#include <string>#include <algorithm>using namespace std;// Edit penalty:// gap penalty: theta// mismatch penalty: alphaint optimalAlignment(string s1, string s2, int theta, int alpha, vector<string> &align){    int m = s1.length();    int n = s2.length();    vector<vector<int>> M(m + 1, vector<int>(n + 1, 0));    for (int i = 0; i <= m; ++i)        M[i][0] = i * theta;    for (int j = 0; j <= n; ++j)        M[0][j] = j * theta;        for (int i = 1; i <= m; ++i)    {        for (int j = 1; j <= n; ++j)        {            int cost = 0; int kind = 1;            if (s1[i - 1] == s2[j - 1])                cost = M[i - 1][j - 1];            else                cost = M[i - 1][j - 1] + alpha;            cost = min(cost, theta + M[i - 1][j]);            cost = min(cost, theta + M[i][j - 1]);            M[i][j] = cost;                    }    }        int i = m, j = n;    // backtracking to get the solution    while (i >= 1 || j >= 1)    {        if ((i > 0 && M[i][j] == theta + M[i - 1][j]) || j == 0)        {            align[0].push_back(s1[i - 1]);            align[1].push_back('-');            --i;        }        else if ((j > 0 && M[i][j] == theta + M[i][j - 1]) || i == 0)        {            align[0].push_back('-');            align[1].push_back(s2[j - 1]);            --j;        }        else if (M[i][j] == M[i - 1][j - 1] + alpha || M[i][j] == M[i - 1][j - 1] )        {            align[0].push_back(s1[i - 1]);            align[1].push_back(s2[j - 1]);            --i, --j;        }    }    reverse(align[0].begin(), align[0].end());    reverse(align[1].begin(), align[1].end());        return M[m][n];}// x   y  cost// ------------// A   T   1// A   A   0// C   -   2// A   A   0// G   G   0// T   G   1// T   T   0// A   -   2// C   C   0// C   A   1//        ---//   7int main(){    string x = "";    string y = "TAAGGTCA";    string x1 = "CTACCG";    string y1 = "TACATG";    vector<string> ret;    ret.resize(2);    cout << "s1: " << x << endl;    cout << "s2: " << y << endl;    cout << "The minimum cost of alignment: " << optimalAlignment(x, y, 2, 1, ret) << endl;    cout << ret[0] << endl;    cout << ret[1] << endl;    ret.clear();    ret.resize(2);    cout << "s1: " << x1 << endl;    cout << "s2: " << y1 << endl;    cout << "The minimum cost of alignment: " << optimalAlignment(x1, y1, 2, 1, ret) << endl;    cout << ret[0] << endl;    cout << ret[1] << endl;    return 0;}

Reference:

http://www.cs.princeton.edu/courses/archive/fall14/cos126/assignments/sequence.html