malloc 空间 是否 连续_序列比对(七)——序列比对之线性空间算法

原创: 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: %fn", 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);
}

（公众号：生信了）