蛋白质序列两两比对C++实现回溯环节异常问题求助

问题描述 #

已实现的DNA序列两两比对算法移植到蛋白质序列场景后，回溯（traceback）阶段出现异常。测试用两条相同的蛋白质序列AAGG，生成的得分矩阵如下：

0 -5 -10 -15 -20
-5 4 -1 -6 -11
-10 -1 8 3 -2
-15 -6 3 13 8
-20 -11 -2 8 18

预期回溯过程应输出得分：18、13、8、4、0，但程序仅输出到8即停止，无报错或崩溃现象。

相关C++实现代码：

Prot1_Seq = "AAGG";
Prot2_Seq = "AAGG";

cout << "\nSequence 1: " << Prot1_Seq << endl;
int Prot1_len = Prot1_Seq.length();

cout << "\nSequence 2: " << Prot2_Seq << endl;
int Prot2_len = Prot2_Seq.length();

int matrix[Prot1_len+1][Prot2_len+1] {0};

// applying the gap penalty to the outer border columns
for (int i {1}; i <= Prot1_len; i++)
{
        matrix[i][0] = gap_pen  * i;

}
for (int j {1}; j <= Prot2_len; j++)
{
        matrix[0][j] = gap_pen * j;

}

fstream myfile;
myfile.open("100pam.txt", ios::in);

if (!myfile.is_open())
{
        cerr << "Error can not open file" << endl;
        exit(1);
}


const int len = 20;
int pam_matrix[len][len];

string line;

getline(myfile,line); // get ing rid of the first line of pam matrix file, by using getline to grab it

cout << "\n" << endl;

for (int i = 0; i < len; i++)
{
        char Row_first_letter;
        myfile >> Row_first_letter; // using the >> to grab the first char of each row leaving the pam matrix with just numbers

        for(int j = 0; j < len; j++)
        {
                myfile >> pam_matrix[i][j];
                cout << pam_matrix[i][j] << "\t";
        }
        cout << endl;
}

cout << endl;

int match_score {};
int mismatch {};

char row_labels [20] {'A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G', 'H', 'I', 'L', 'K', 'M', 'F', 'P', 'S', 'T', 'W', 'Y', 'V'};
char col_labels [20] {'A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G', 'H', 'I', 'L', 'K', 'M', 'F', 'P', 'S', 'T', 'W', 'Y', 'V'};

for(int i {0}; i <= Prot1_len; i++)
{
        for(int j {0}; j <= Prot2_len; j++)
        {
                if (i > 0 && j > 0)
                {
                        if (Prot1_Seq[i-1] == Prot2_Seq[j-1])
                        {

                    int row = distance(row_labels, find(row_labels, row_labels + 20, Prot1_Seq[i-1]));
                                        int col = distance(col_labels, find(col_labels, col_labels + 20, Prot2_Seq[j-1]));

                    match_score = pam_matrix[row][col];

                    matrix[i][j] = matrix[i-1][j-1] + match_score;
                        }
                else if (Prot1_Seq[i-1] != Prot2_Seq[j-1])
                {
                    //row_labels[i-1] = Prot1_Seq[i-1];
                    //col_labels[i-1] = Prot2_Seq[i-1];

                    int row = distance(row_labels, find(row_labels, row_labels + 20, Prot1_Seq[i-1]));
                    int col = distance(col_labels, find(col_labels, col_labels + 20, Prot2_Seq[j-1]));

                    mismatch = pam_matrix[row][col];

                    matrix[i][j] = max({matrix[i-1][j-1] + mismatch,
                                                matrix[i-1][j] + gap_pen ,
                                                matrix[i][j-1] + gap_pen});
                }
                }
                cout << matrix[i][j] << "\t";
        }
        cout << endl;
}

int terminal = matrix[Prot1_len+1-1][Prot2_len+1-1]; // rows and colums are length+1 long. and to get the edge value it's col/row length -1.
cout << terminal << endl;

int temp = terminal;
int decrement_X = Prot1_len+1-1;
int decrement_Y = Prot2_len+1-1;

string Prot1_aligned = "";
string Prot2_aligned = "";

int inc_seq1 = Prot1_len-1;
int inc_seq2 = Prot2_len-1;

while (!(decrement_X == 0 || decrement_Y == 0))
{//while
        if (terminal == matrix[decrement_X][decrement_Y-1] + gap_pen) {
                temp = matrix[decrement_X][decrement_Y-1];
                cout << temp << endl;
                terminal = temp;
                decrement_Y--;
                Prot1_aligned += "_";
                Prot2_aligned += Prot2_Seq[inc_seq2];
                inc_seq2--;
        }

        else if (terminal == matrix[decrement_X-1][decrement_Y] + gap_pen) {
                temp = matrix[decrement_X-1][decrement_Y];
            cout << temp << endl;
                terminal = temp;
                decrement_X--;
                Prot1_aligned += Prot1_Seq[inc_seq1];
                Prot2_aligned += "_";
                inc_seq1--;
        }

        else if (terminal == matrix[decrement_X-1][decrement_Y-1] + match_score) {
                temp = matrix[decrement_X-1][decrement_Y-1];
            cout << temp << endl;
                terminal = temp;
                decrement_X--;
                decrement_Y--;
                Prot1_aligned += Prot1_Seq[inc_seq1];
                Prot2_aligned += Prot2_Seq[inc_seq2];
                inc_seq1--;
                inc_seq2--;
        }

        else if (terminal == matrix[decrement_X-1][decrement_Y-1] + mismatch) {
                temp = matrix[decrement_X-1][decrement_Y-1];
            cout << temp << endl;
                terminal = temp;
                decrement_X--;
                decrement_Y--;
                Prot1_aligned += Prot1_Seq[inc_seq1];
                Prot2_aligned += Prot2_Seq[inc_seq2];
                inc_seq1--;
                inc_seq2--;
        }
} //while

reverse(Prot1_aligned.begin(), Prot1_aligned.end());
reverse(Prot2_aligned.begin(), Prot2_aligned.end());

cout << "Performing alignment" << " .................................................. 100%" << endl;
cout << "Alignment complete! Your aligned sequences are:" << endl;

cout << "Sequence 1: " << Prot1_aligned << endl;
cout << "Sequence 2: " << Prot2_aligned << endl;

问题根源 #

1. 回溯阶段依赖失效的变量：match_score和mismatch仅在得分矩阵计算阶段赋值，回溯时使用的是矩阵计算最后一次赋值的旧值，而非当前回溯位置对应的PAM得分。例如，回溯到(2,2)位置（得分8）时，match_score仍为最后一次匹配G-G的得分5，导致8 == 4 + 5判断不成立，四个回溯条件均不满足，循环无法推进。
2. 矩阵计算逻辑不一致：匹配分支直接取对角方向得分加匹配分，未考虑插入/删除的可能，虽不影响当前测试用例，但不符合序列比对算法的规范逻辑。

修复方案 #

1. 回溯阶段实时计算当前位置得分 #

删除对match_score和mismatch的依赖，在回溯循环内部根据当前位置的序列字符，实时从PAM矩阵获取对应得分，统一判断对角方向的来源：

while (!(decrement_X == 0 || decrement_Y == 0))
{
    // 实时获取当前回溯位置对应的PAM矩阵得分
    char curr_prot1 = Prot1_Seq[decrement_X - 1];
    char curr_prot2 = Prot2_Seq[decrement_Y - 1];
    int row_idx = distance(row_labels, find(row_labels, row_labels + 20, curr_prot1));
    int col_idx = distance(col_labels, find(col_labels, col_labels + 20, curr_prot2));
    int current_pam_score = pam_matrix[row_idx][col_idx];

    if (terminal == matrix[decrement_X][decrement_Y-1] + gap_pen) {
        temp = matrix[decrement_X][decrement_Y-1];
        cout << temp << endl;
        terminal = temp;
        decrement_Y--;
        Prot1_aligned += "_";
        Prot2_aligned += Prot2_Seq[inc_seq2];
        inc_seq2--;
    }
    else if (terminal == matrix[decrement_X-1][decrement_Y] + gap_pen) {
        temp = matrix[decrement_X-1][decrement_Y];
        cout << temp << endl;
        terminal = temp;
        decrement_X--;
        Prot1_aligned += Prot1_Seq[inc_seq1];
        Prot2_aligned += "_";
        inc_seq1--;
    }
    else if (terminal == matrix[decrement_X-1][decrement_Y-1] + current_pam_score) {
        temp = matrix[decrement_X-1][decrement_Y-1];
        cout << temp << endl;
        terminal = temp;
        decrement_X--;
        decrement_Y--;
        Prot1_aligned += Prot1_Seq[inc_seq1];
        Prot2_aligned += Prot2_Seq[inc_seq2];
        inc_seq1--;
        inc_seq2--;
    }
}

2. 统一矩阵计算逻辑（可选但规范） #

将匹配分支的逻辑修改为与错配分支一致，取三个方向的最大值，确保算法逻辑统一：

if (Prot1_Seq[i-1] == Prot2_Seq[j-1])
{
    int row = distance(row_labels, find(row_labels, row_labels + 20, Prot1_Seq[i-1]));
    int col = distance(col_labels, find(col_labels, col_labels + 20, Prot2_Seq[j-1]));
    match_score = pam_matrix[row][col];
    // 改为取三个方向的最大值，与错配逻辑保持一致
    matrix[i][j] = max({matrix[i-1][j-1] + match_score,
                        matrix[i-1][j] + gap_pen,
                        matrix[i][j-1] + gap_pen});
}

修复效果 #

修改后，回溯过程将正确输出18、13、8、4、0，并生成正确的比对结果：

Sequence 1: AAGG
Sequence 2: AAGG

内容的提问来源于stack exchange，提问作者TheLordRaj