Как правильно осуществить выравнивание последовательности

Question

Я пытаюсь создать программу выравнивания последовательностей и пытаюсь найти какой-то код, который я нашел в Интернете, но я пытаюсь выяснить, почему выходные данные моей диаграммы иногда корректны, особенноближе к концу, как я всегда, кажется, получаю правильную оценку выравнивания, но моя диаграмма всегда отклоняется на несколько пунктов здесь и там, как показано выше. (Требуемый вывод справа, неправильный вывод слева) Я подозреваю, что ошибки находятся где-то в матрице подсчета замещения. Я понимаю, что некоторые функции не нужны, но не решаюсь заменить их из-за боязни испортить то, что у меня уже есть.

#include <iostream>
#include <fstream>
#include <string>
#include <cmath>
#include <algorithm>
#include <time.h>

using namespace std;

// Affine gap function
int gap_affinity(int gap) {
    int gap_aff = gap;

    return gap_aff;
}

// Get maximal score and trace it
int  max_score(int up, int diag, int left, char* ptr) {
    int  max = 0;

    if (diag >= up && diag >= left) {
        max = diag;
        *ptr = '\\';
    }
    else if (up > left) {
        max = up;
        *ptr = '|';
    }
    else {
        max = left;
        *ptr = '-';
    }

    return  max;
}

// Initialise scoring matrix with first row and column
void  init(int** M, char** M_tb, int A_n, int B_n, int gap) {
    M[0][0] = 0;
    M_tb[0][0] = 'n';
    int gapvert = gap;
    int gaphor = gap;

    int i = 0, j = 0;

    for (j = 1; j <= A_n; j++) {
        M[0][j] = -(gaphor);
        M_tb[0][j] = '-';
        gaphor += 2;
    }
    for (i = 1; i <= B_n; i++) {
        M[i][0] = -(gapvert);
        M_tb[i][0] = '|';
        gapvert += 2;
    }
}

// Needleman and Wunsch algorithm
int alignment(int** M, char** M_tb, string A, string B, string& A_al, string& B_al, int A_n, int B_n, int match, int misMatch, int gap) {
    int x = 0, y = 0;
    int scU, scD, scL; // scores for respectively cell above, diagonal and left
    char ptr, nuc;
    int i = 0, j = 0;

    // create substitution scoring matrix
    const int  s[3][3] = { { match, misMatch, misMatch },
                            { misMatch, match, misMatch },
                            { misMatch, misMatch, match } };

    for (i = 1; i <= B_n; i++) {
        for (j = 1; j <= A_n; j++) {
            nuc = A[j - 1];

            switch (nuc)
            {
            case 'C':  x = 0;  break;
            case 'T':  x = 1;  break;
            case 'A':  x = 2;  break;
            case 'G':  x = 3;
            }

            nuc = B[i - 1];

            switch (nuc)
            {
            case 'C':  y = 0;  break;
            case 'T':  y = 1;  break;
            case 'A':  y = 2;  break;
            case 'G':  y = 3;
            }

            scU = M[i - 1][j] - gap_affinity(gap); // get score if trace would go up
            scD = M[i - 1][j - 1] + s[x][y]; // get score if trace would go diagonal
            scL = M[i][j - 1] - gap_affinity(gap); // get score if trace would go left

            M[i][j] = max_score(scU, scD, scL, &ptr); // get max score for current optimal global alignment

            M_tb[i][j] = ptr;
        }
    }
    i--; j--;

    while (i > 0 || j > 0) {
        switch (M_tb[i][j])
        {
        case '|':      A_al += '-';
            B_al += B[i - 1];
            i--;
            break;

        case '\\':      A_al += A[j - 1];
            B_al += B[i - 1];
            i--;  j--;
            break;

        case '-':      A_al += A[j - 1];
            B_al += '-';
            j--;
        }
    }

    reverse(A_al.begin(), A_al.end());
    reverse(B_al.begin(), B_al.end());

    return  0;
}

// Print the scoring matrix
void  print_mtx(int** M, string A, string B, int A_n, int B_n) {
    cout << "        ";
    for (int j = 0; j < A_n; j++) {
        cout << A[j] << "   ";
    }
    cout << "\n  ";

    for (int i = 0; i <= B_n; i++) {
        if (i > 0) {
            cout << B[i - 1] << " ";
        }
        for (int j = 0; j <= A_n; j++) {
            cout.width(3);
            cout << M[i][j] << " ";
        }
        cout << endl;
    }
    cout << endl;
}

// Print the traceback matrix
void  print_tb(char** M_tb, string A, string B, int A_n, int B_n) {
    cout << "        ";
    for (int j = 0; j < A_n; j++) {
        cout << A[j] << "   ";
    }
    cout << "\n  ";

    for (int i = 0; i <= B_n; i++) {
        if (i > 0) {
            cout << B[i - 1] << " ";
        }
        for (int j = 0; j <= A_n; j++) {
            cout.width(3);
            cout << M_tb[i][j] << " ";
        }
        cout << endl;
    }
    cout << endl;
}

// Initiate matrices, align and export
int** NW(string A, string B, string& A_al, string& B_al, int A_n, int B_n, int match, int misMatch, int gap) {
    // Create alignment matrix
    int** M = new int* [B_n + 1];
    for (int i = 0; i <= B_n; i++) {
        M[i] = new int[A_n + 1];
    }

    // Create traceback matrix
    char** M_tb = new char* [B_n + 1];
    for (int i = 0; i <= B_n; i++) {
        M_tb[i] = new char[A_n + 1];
    }

    clock_t t; // for timing execution
    t = clock(); // get time of start

    // Initialize traceback and F matrix (fill in first row and column)
    init(M, M_tb, A_n, B_n, gap);


    // Create alignment
    alignment(M, M_tb, A, B, A_al, B_al, A_n, B_n, match, misMatch, gap);

    int score = M[B_n][A_n]; // get alignment score

    if (true) {
        print_mtx(M, A, B, A_n, B_n);
        print_tb(M_tb, A, B, A_n, B_n);
    }

    if (true) {
        cout << endl << "Alignments: " << endl;
        int start = 0; // start of new line for printing alignments
        int cntr = 0; // iterator for printing alignments
        int Al_n = A_al.length(); // length of alignment
        do {
            cout << start + 1 << " A: ";
            for (cntr = start; cntr < start + 150; cntr++) {
                if (cntr < Al_n) {
                    cout << A_al[cntr];
                }
                else {
                    break;
                }
            }
            cout << " " << endl << start + 1 << " B: ";
            for (cntr = start; cntr < start + 150; cntr++) {
                if (cntr < Al_n) {
                    cout << B_al[cntr];
                }
                else {
                    break;
                }
            }
            cout << " " << endl << endl;
            start += 150;
        } while (start <= Al_n);
    }

    // Show score and runtime
    cout << "Alignment score: " << score << endl;



    // Free memory
    for (int i = 0; i <= B_n; i++)  delete M[i];
    delete[] M;
    for (int i = 0; i <= B_n; i++)  delete M_tb[i];
    delete[] M_tb;

    return 0;
}

int main(int argc, char** argv) {
    string A, B; // sequence to be aligned A and B
    string A_al, B_al = ""; // aligned sequence A and B
    int match = 5; // match
    int misMatch = -3; // mismatch
    int gap = 2; // initial gap penalty. Gap penalty is lower than mismatch: two sequences from same species assumed.
    string line; // reading in data

    // User interface

    cout << "==============================" << endl;
    cout << "Enter the first string ";
    getline(cin, line);
    A = line;
    cout << "Enter the second string ";
    getline(cin, line);
    B = line;

    cout << "Enter score for matching. ";
    cin >> match;
    cout << "Enter score for mismatch. ";
    cin >> misMatch;
    cout << "Enter gap pentalty. ";
    cin >> gap;

    // Read in the alignments
    // Run the alignment script
    int A_n = A.length();
    int B_n = B.length();

    NW(A, B, A_al, B_al, A_n, B_n, match, misMatch, gap);

    // Output parameters
    cout << endl << "Used parameters:" << endl;
    cout << "Match = " << match << endl;
    cout << "Mismatch = " << misMatch << endl;
    cout << "Gap = -" << gap << endl;

    return 0;
}