前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >序列比对(七)序列比对之线性空间算法

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

作者头像
一只羊
发布2019-07-27 18:50:53
1.5K0
发布2019-07-27 18:50:53
举报
文章被收录于专栏:生信了

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

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

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

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

图片内容引自《生物序列分析》

如图中所说,关键点就是找到v值,然后通过不断的分划,最终得到全部的比对序列。本文给出了这种算法的一种代码实现。

代码的关键在于终止条件的设置以及必要时巧妙地颠倒行列。与 O(mn) 阶的算法相比,这种算法只能得到其中一种最佳比对方式,而无法得到所有的可能。

代码运行的效果:

具体的代码如下:

(由于代码中运用了“引用(具体地就是指代码中的 int &n 这一用法)”这一方法,所以是以cpp文件编译的)

代码语言:javascript
复制
#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);
}
本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2019-06-29,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信了 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档