序列比对(七)——序列比对之线性空间算法

原创: hxj7

一般而言，运用动态规划算法进行序列比对对内存空间的要求是 O(mn) 阶的，本文介绍了一种线性空间要求的序列比对方法。

前文如《序列比对（一）全局比对Needleman-Wunsch算法》所介绍的运用动态规划算法进行序列比对时，对内存空间的要求是 O(mn) 阶的。如果只是要求最大得分而不要求最大得分所对应的序列（也就是不进行回溯）的话，其实只需要保留上一行的得分就够了，这一点可以从计算得分矩阵的迭代公式中看出来。

图片引自https://www.jianshu.com/p/2b99d0d224a2

但是如果要求回溯呢，是否有一种线性空间算法来进行序列比对呢？前人已经给出了多种算法。其中一种在《生物序列分析》一书中给出，描述如下：
图片内容引自《生物序列分析》

如图中所说，关键点就是找到v值，然后通过不断的分划，最终得到全部的比对序列。本文给出了这种算法的一种代码实现。
代码的关键在于终止条件的设置以及必要时巧妙地颠倒行列。与 O(mn) 阶的算法相比，这种算法只能得到其中一种最佳比对方式，而无法得到所有的可能。

代码运行的效果：

具体的代码如下：
（由于代码中运用了“引用（具体地就是指代码中的 int &n 这一用法）”这一方法，所以是以cpp文件编译的）

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MAXSEQ 1000
#define GAP_CHAR '-'
#define max2(a, b) ((a) > (b) ? (a) : (b))
#define exchange(a, b, c) {\
    *(c) = *(a);  \
    *(a) = *(b);  \
    *(b) = *(c);  \
}

// 对空位的罚分是线性的
struct Unit {
    int v;         // 中间折半点(u, v)中的v值
    float M;      // 得分矩阵第(i, j)这个单元的分值，即序列s(1,...,i)与序列r(1,...,j)比对的最高得分
};
typedef struct Unit *pUnit;

void strUpper(char *s);
float getFScore(char a, char b);
float divideAlign(pUnit** a, int ti, int tj, int bi, int bj, char* saln, char* raln, int &n);
void align(void);

char* sraw;
char* rraw;
float gap = -2.5;     // 对空位的罚分
float tm[3];

int main() {
    if ((sraw = (char*) malloc(sizeof(char) * MAXSEQ)) == NULL || \
        (rraw = (char*) malloc(sizeof(char) * MAXSEQ)) == NULL) {
        fputs("Error: out of space!\n", stderr);
        exit(1);
    }
    printf("The 1st seq: ");
    scanf("%s", sraw);
    printf("The 2nd seq: ");
    scanf("%s", rraw);
    align();
    free(sraw);
    free(rraw);
    return 0;
}

void strUpper(char *s) {
    while (*s != '\0') {
        if (*s >= 'a' && *s <= 'z') {
            *s -= 32;
        }
        s++;
    }
}

// 替换矩阵：match分值为5，mismatch分值为-4
// 数组下标是两个字符的ascii码减去65之后的和
float FMatrix[] = {
    5, 0, -4, 0, 5, 0, -4, 0, -4, 0,
    0, 0, 5, 0, 0, 0, 0, 0, 0, -4,
    0, -4, 0, 0, 0, -4, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 5
};

float getFScore(char a, char b) {
    return FMatrix[a + b - 'A' - 'A'];
}

float divideAlign(pUnit** a, int ti, int tj, int bi, int bj, char* saln, char* raln, int &n) {
    // ti, tj, bi, bj 分别是左上角和右下角单元的横纵坐标；
    // 为了统一，在填充序列时不算左上角
    int i, j, k, u, v;
    float maxScore;
    int flag;  // 是否颠倒了行列
    int itmp;
    char* ptmp;
    if (bi - ti == 1 && bj - tj == 1) {
        saln[n] = sraw[bi - 1];
        raln[n] = rraw[bj - 1];
        n++;
        return getFScore(sraw[bi - 1], rraw[bj - 1]);
    }
    if (ti == bi) {
        for (j = tj + 1; j <= bj; j++, n++) {
            saln[n] = GAP_CHAR;
            raln[n] = rraw[j - 1];
        }
        return (bj - tj) * gap;
    }
    if (tj == bj) {
        for (i = ti + 1; i <= bi; i++, n++) {
            saln[n] = sraw[i - 1];
            raln[n] = GAP_CHAR;
        }
        return (bi - ti) * gap;        
    }
    // 为什么有时要颠倒行列？
    // 因为在bi - ti = 1这种情况下有可能会出现u = ti, v = tj从而导致无限循环
    // 如果颠倒行列，可以将分治进行下去，因为此时不可能出现u = ti, v = tj的情况。
    flag = 0;
    if (bi - ti <= 1) { // 注意此时 bj - tj > 1
        exchange(&ti, &tj, &itmp);
        exchange(&bi, &bj, &itmp);
        exchange(&sraw, &rraw, &ptmp);
        flag = 1;        
    }
    u = (ti + bi) / 2;   // 中间行/列
    for (j = tj; j <= bj; j++)
        a[0][j]->M = (j - tj) * gap;
    k = 1;   // 当前行
    for (i = ti + 1; i <= bi; i++, k = 1 - k) {
        // j = tj
        a[k][tj]->M = (i - ti) * gap;
        a[k][tj]->v = tj;
        // 动态规划算法计算每个单元的分值
        for (j = tj + 1; j <= bj; j++) {
            tm[0] = a[k][j - 1]->M + gap;
            tm[1] = a[1 - k][j - 1]->M + getFScore(sraw[i - 1], rraw[j - 1]);
            tm[2] = a[1 - k][j]->M + gap;
            a[k][j]->M = max2(max2(tm[0], tm[1]), tm[2]);
            if (i == u)
                a[k][j]->v = j;
            else if (i > u) {  // 注意这里没有保存所有可能的v，只挑选了其中一种可能的v
                if (tm[0] == a[k][j]->M) a[k][j]->v = a[k][j - 1]->v;
                else if (tm[1] == a[k][j]->M) a[k][j]->v = a[1 - k][j - 1]->v;
                else a[k][j]->v = a[1 - k][j]->v;
            }
        }
    }
    v = a[1 - k][bj]->v;
    maxScore = a[1 - k][bj]->M;
    if (flag) {  // 如果颠倒过一次行列，那么还要颠倒回来
        exchange(&ti, &tj, &itmp);
        exchange(&bi, &bj, &itmp);
        exchange(&u, &v, &itmp);
        exchange(&sraw, &rraw, &ptmp);
    }
    divideAlign(a, ti, tj, u, v, saln, raln, n);
    divideAlign(a, u, v, bi, bj, saln, raln, n);
    return maxScore;
}

void align(void) {
    int i, j, l;
    int m = strlen(sraw);
    int n = strlen(rraw);
    int r = max2(m, n);
    pUnit** aUnit;
    char* salign;
    char* ralign;
    float maxScore;
    // 初始化
    if ((aUnit = (pUnit **) malloc(sizeof(pUnit*) * 2)) == NULL) {
        fputs("Error: Out of space!\n", stderr);
        exit(1);        
    }
    for (i = 0; i <= 1; i++) {
        if ((aUnit[i] = (pUnit *) malloc(sizeof(pUnit) * (r + 1))) == NULL) {
            fputs("Error: Out of space!\n", stderr);
            exit(1);        
        }
        for (j = 0; j <= r; j++) {
            if ((aUnit[i][j] = (pUnit) malloc(sizeof(struct Unit))) == NULL) {
                fputs("Error: Out of space!\n", stderr);
                exit(1);        
            }
        }
    }
    if ((salign = (char*) malloc(sizeof(char) * (m + n))) == NULL || \
        (ralign = (char*) malloc(sizeof(char) * (m + n))) == NULL) {
        fputs("Error: Out of space!\n", stderr);
        exit(1);
    }
    // 将字符串都变成大写
    strUpper(sraw);
    strUpper(rraw); 
    // 递归比对
    l = 0;
    maxScore = divideAlign(aUnit, 0, 0, m, n, salign, ralign, l);
    // 打印比对结果
    printf("max score: %f\n", maxScore);
    for (i = 0; i < l; i++) {
        printf("%c", salign[i]);
    }
    printf("\n");
    for (i = 0; i < l; i++) {
        printf("%c", ralign[i]);
    }
    printf("\n");
    // 释放内存
    free(salign);
    free(ralign);
    for (i = 0; i <= 1; i++) {
        for (j = 0; j <= n; j++)
            free(aUnit[i][j]);
        free(aUnit[i]);
    }
    free(aUnit);
}