美文网首页
序列匹配(五)——重复匹配问题的动态规划算法

序列匹配(五)——重复匹配问题的动态规划算法

作者: 生信了 | 来源:发表于2019-11-09 10:06 被阅读0次

    原创: hxj7

    前言:
    蛋白质序列中常有重复的功能域(domain)或模体(motif)拷贝,由此衍生出一个抽象的序列多重匹配的问题,即如何从一个序列中找出另一个序列的某部分(如功能域或模体)的多个无交叠(non-overlapping)拷贝。本文给出了该问题的示例、关键计算公式以及C语言实现代码。

    问题及算法描述

    更具体地描述上面的问题:有序列x和y,其中y是包含结构域的序列,x是要从中找到多重匹配的序列。将x分割成一段一段的不交叠的子序列,这些子序列要么不参与和y的联配,要么与y的某一段子序列联配且联配的分值不低于一个阈值T。如果将x的某一子序列的联配分值减去T作为其“标准联配分值”,那么最终目标是找到这些参与联配的x的子序列的“标准联配分值”之和的最大值。

    在《Biological sequence analysis》这本书中列举了一个蛋白质的例子:


    image

    引自《生物序列分析》

    上图中显示,在最优联配(即“标准联配分值”之和最大的联配)中,x有两个子序列参与了联配,“标准联配分值”分别是1和8。
    那么上图中的“标准联配分值”是如何计算得到的呢?这依然是利用动态规划算法,在《Biological sequence analysis》书中给出了关键的计算公式:


    image

    引自《生物序列分析》

    其中,F(i, j)假设x(i)参与联配,且对应的联配结束于x(i)与y(j)时的“标准联配分值”之和的最大值;而F(i, 0)指的是子序列x(1,2,…,i)与y(联配可结束于y的任意位置)的“标准联配分值”之和的最大值,假设x(i)不参与联配。

    一个困惑

    上面计算公式的C代码实现见下文,简称其为alnRepeat,以便和本文另一段代码区分。其中一个示例如下:


    image

    没有问题,但是另一个示例的结果让我困惑:


    image

    最优分值应该是6啊。回过头再去看上面的计算公式,如果我对书中相应章节的理解无误的话,它似乎没有考虑到x两个子序列都参与联配且这两个子序列紧挨在一起的情况,最简单的就是上面AA的这个例子。理论上,最优联配中,两个连续的A应该都参与了联配,且属于两个不同的“匹配段”。

    算法的补充

    由此,我重新思考分值的计算公式。F(i, 0)没有问题,而F(i, j)的计算可以分为下面三种情况(注意,F(i, j)的前提条件假设了x(i)“参与”了联配):

    1. x(i-1)没有参与联配;
    2. x(i-1)参与了联配,且与x(i)属于同一个“匹配段”;
    3. x(i-1)参与了联配,且与x(i)属于不同的“匹配段”。

    考虑了上面三种情况,F(i, j)的计算公式变成:


    image

    对上述新公式编写代码的过程中,发现输出结果有很多重复,再次检视公式,发现F(i, 1)需要单独处理,否则F(i, 1)计算公式中的F(i, 0)一项会引发歧义。于是,最终的计算公式变成:


    image

    其C语言实现代码简称alnRepeat3。
    运行alnRepeat以及alnRepeat3比较二者的不同:


    image

    alnRepeat3的结果仍有重复,说明代码还要优化。


    image

    小结

    本文介绍了生物序列重复匹配的问题以及相应的动态规划算法,在代码实现过程中,发现了疑似错误的示例(原计算公式似乎没有考虑到两个“匹配段”紧挨在一起的情况)并补充了计算公式。

    对笔者来说,最大的收获在于学习这个新算法的过程中认真地进行过思考,但由于自身能力以及时间精力所限,对这个问题的理解以及代码实现还有很多不足,真切期望能有热心同道能够给出意见和建议!

    最后,谢谢大家看此长文!

    C代码

    首先是alnRepeat

    #include <stdio.h>
    #include <stdlib.h>
    #include <string.h>
    #define MAXSEQ 1000
    #define GAP_CHAR '-'
    #define UNMATCH_CHAR '.'
    
    // 对空位的罚分是线性的
    struct Unit {
        int W0;   // 不参与联配,为F(i, 0)和F(i, j)共用
        int *Wj;      // 跳转到F(i - 1, j)
        int nj;       // Wj数组的大小
        int W1;   // 是否往左回溯一格
        int W2;   // 是否往左上回溯一格
        int W3;   // 是否往上回溯一格
        float M;      // 得分矩阵第(i, j)这个单元的分值,即序列s(1,...,i)与序列r(1,...,j)比对的最高得分
    };
    typedef struct Unit *pUnit;
    
    void strUpper(char *s);
    float max2(float a, float b);
    float max3(float a, float b, float c);
    float getFScore(char a, char b);
    void printAlign(pUnit** a, const int i, const int j, char* s, char* r, char* saln, char* raln, int n);
    void align(char *s, char *r, float t);
    
    int main() {
        char s[MAXSEQ];
        char r[MAXSEQ];
        float t;
        printf("The 1st seq: ");
        scanf("%s", s);
        printf("The 2nd seq: ");
        scanf("%s", r);
        printf("T (threshold): ");
        scanf("%f", &t);
        align(s, r, t);
        return 0;
    }
    
    void strUpper(char *s) {
        while (*s != '\0') {
            if (*s >= 'a' && *s <= 'z') {
                *s -= 32;
            }
            s++;
        }
    }
    
    float max2(float a, float b) {
        return a > b ? a : b;
    }
    
    float max3(float a, float b, float c) {
        float f = a > b ? a : b;
        return f > c ? f : c;
    }
    
    // 替换矩阵: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'];
    }
    
    void printAlign(pUnit** a, const int i, const int j, char* s, char* r, char* saln, char* raln, int n) {
        int k, l;
        pUnit p = a[i][j];
        if (! i) {  // 保证序列s的每个字符都比对上
            for (k = n - 1; k >= 0; k--)
                printf("%c", saln[k]);
            printf("\n");
            for (k = n - 1; k >= 0; k--)
                printf("%c", raln[k]);
            printf("\n\n");
            return;
        }
        if (! j) {  // F(i, 0)
            saln[n] = s[i - 1];
            raln[n] = UNMATCH_CHAR;
            if (p->W0)
                printAlign(a, i - 1, 0, s, r, saln, raln, n + 1);
            for (k = 0; k < p->nj; k++)
                if (p->Wj[k])
                    printAlign(a, i - 1, k + 1, s, r, saln, raln, n + 1);
        } else {
            if (p->W0) {
                printAlign(a, i, 0, s, r, saln, raln, n);
            }
            if (p->W1) {    // 向上回溯一格
                saln[n] = s[i - 1];
                raln[n] = GAP_CHAR;
                printAlign(a, i - 1, j, s, r, saln, raln, n + 1);
            }
            if (p->W2) {    // 向左上回溯一格
                saln[n] = s[i - 1];
                raln[n] = r[j - 1];
                printAlign(a, i - 1, j - 1, s, r, saln, raln, n + 1);
            }
            if (p->W3) {    // 向左回溯一格
                saln[n] = GAP_CHAR;
                raln[n] = r[j - 1];
                printAlign(a, i, j - 1, s, r, saln, raln, n + 1);
            }
        }
    }
    
    void align(char *s, char *r, float t) {
        int i, j, k;
        int m = strlen(s);
        int n = strlen(r);
        float gap = -2.5;     // 对空位的罚分
        float m1, m2, m3, maxm;
        float em;   // F(m + 1, 0)
        pUnit **aUnit;
        char* salign;
        char* ralign;
        // 初始化
        if ((aUnit = (pUnit **) malloc(sizeof(pUnit*) * (m + 1))) == NULL) {
            fputs("Error: Out of space!\n", stderr);
            exit(1);
        }
        for (i = 0; i <= m; i++) {
            if ((aUnit[i] = (pUnit *) malloc(sizeof(pUnit) * (n + 1))) == NULL) {
                fputs("Error: Out of space!\n", stderr);
                exit(1);     
            }
            for (j = 0; j <= n; j++) {
                if ((aUnit[i][j] = (pUnit) malloc(sizeof(struct Unit))) == NULL) {
                    fputs("Error: Out of space!\n", stderr);
                    exit(1);     
                }
                aUnit[i][j]->W0 = 0;
                aUnit[i][j]->W1 = 0;
                aUnit[i][j]->W2 = 0;
                aUnit[i][j]->W3 = 0;
                // 创建F(i, 0)的跳转数组
                aUnit[i][j]->nj = n;
                if ((aUnit[i][j]->Wj = (int*) malloc(sizeof(int) * aUnit[i][j]->nj)) == NULL) {
                    fputs("Error: Out of space!\n", stderr);
                    exit(1);     
                }
                for (k = 0; k < aUnit[i][j]->nj; k++) {
                    aUnit[i][j]->Wj[k] = 0;
                }
            }
        }
        for (j = 0; j <= n; j++) {
            aUnit[0][j]->M = 0;
        }
        // 将字符串都变成大写
        strUpper(s);
        strUpper(r);
        // 动态规划算法计算得分矩阵每个单元的分值
        for (i = 1; i <= m; i++) {
            // 计算F(i, 0)
            aUnit[i][0]->M = aUnit[i - 1][0]->M;
            for (j = 1; j <= n; j++)
                if (aUnit[i][0]->M < aUnit[i - 1][j]->M - t)
                    aUnit[i][0]->M = aUnit[i - 1][j]->M - t;
            if (aUnit[i][0]->M == aUnit[i - 1][0]->M)
                aUnit[i][0]->W0 = 1;
            for (j = 1; j <= n; j++)
                if (aUnit[i][0]->M == aUnit[i - 1][j]->M - t)
                    aUnit[i][0]->Wj[j - 1] = 1;
            // 计算F(i, j), j>=1
            for (j = 1; j <= n; j++) {
                m1 = aUnit[i - 1][j]->M + gap;
                m2 = aUnit[i - 1][j - 1]->M + getFScore(s[i - 1], r[j - 1]);
                m3 = aUnit[i][j - 1]->M + gap;
                maxm = max2(max3(m1, m2, m3), aUnit[i][0]->M);
                aUnit[i][j]->M = maxm;
                if (aUnit[i][0]->M == maxm) aUnit[i][j]->W0 = 1;
                if (m1 == maxm) aUnit[i][j]->W1 = 1;
                if (m2 == maxm) aUnit[i][j]->W2 = 1;
                if (m3 == maxm) aUnit[i][j]->W3 = 1;
            }
        }
        // 计算F(m + 1, 0)
        em = aUnit[m][0]->M;
        for (j = 1; j <= n; j++)
            if (em < aUnit[m][j]->M - t)
                em = aUnit[m][j]->M - t;
    /*
        // 打印得分矩阵
        for (i = 0; i <= m; i++) {
            for (j = 0; j <= n; j++)
                printf("%f ", aUnit[i][j]->M);
            printf("\n");
        }
    */
        printf("max score: %f\n", em);
        // 打印最优比对结果,如果有多个,全部打印
        // 递归法
        if (em == 0) {
            fputs("No seq aligned.\n", stdout);
        } else {
            if ((salign = (char*) malloc(sizeof(char) * (3 * m))) == NULL) {  // 3m的大小是根据match以及gap的分值来估算的。
                fputs("Error: Out of space!\n", stderr);
                exit(1);
            }
            if ((ralign = (char*) malloc(sizeof(char) * (3 * m))) == NULL) {
                fputs("Error: Out of space!\n", stderr);
                exit(1);
            }
            if (em == aUnit[m][0]->M)
                printAlign(aUnit, m, 0, s, r, salign, ralign, 0);
            for (j = 1; j <= n; j++)
                if (em == aUnit[m][j]->M - t)
                    printAlign(aUnit, m, j, s, r, salign, ralign, 0);
            // 释放内存
            free(salign);
            free(ralign);
        }
        for (i = 0; i <= m; i++) {
            for (j = 0; j <= n; j++) {
                free(aUnit[i][j]->Wj);
                free(aUnit[i][j]);
            }
            free(aUnit[i]);
        }
        free(aUnit);
    }
    

    然后是alnRepeat3

    #include <stdio.h>
    #include <stdlib.h>
    #include <string.h>
    #define MAXSEQ 1000
    #define GAP_CHAR '-'
    #define UNMATCH_CHAR '.'
    #define max2(a, b) ((a) > (b) ? (a) : (b))
    
    // 对空位的罚分是线性的
    // 00000001  F(i, 0)
    // 00000010  是否往左回溯一格
    // 00000100  是否往左上回溯一格
    // 00001000  是否往上回溯一格
    // 00010000  F(i, 0) + s(i, j)
    struct Unit {
        int W0;   // Xi不参与联配
        int *Wj;      // 跳转到F(i - 1, j)
        int nj;       // Wj数组的大小
        int Wi;     // F(i, j)的回溯指标,不同的bit代表不同的回溯方式
        float M;      // 得分矩阵第(i, j)这个单元的分值,即序列s(1,...,i)与序列r(1,...,j)比对的最高得分
    };
    typedef struct Unit *pUnit;
    
    void strUpper(char *s);
    float maxArray(float *a, int n);
    float getFScore(char a, char b);
    void printAlign(pUnit** a, const int i, const int j, char* s, char* r, char* saln, char* raln, int n, int flag);
    void align(char *s, char *r, float t);
    
    int main() {
        char s[MAXSEQ];
        char r[MAXSEQ];
        float t;
        printf("The 1st seq: ");
        scanf("%s", s);
        printf("The 2nd seq: ");
        scanf("%s", r);
        printf("T (threshold): ");
        scanf("%f", &t);
        align(s, r, t);
        return 0;
    }
    
    void strUpper(char *s) {
        while (*s != '\0') {
            if (*s >= 'a' && *s <= 'z') {
                *s -= 32;
            }
            s++;
        }
    }
    
    float maxArray(float *a, int n) {
        float max = a[0];
        int i;
        for (i = 1; i < n; i++) {
            if (a[i] > max)
                max = a[i];
        }
        return max;
    }
    
    // 替换矩阵: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'];
    }
    
    void printAlign(pUnit** a, const int i, const int j, char* s, char* r, char* saln, char* raln, int n, int flag) {
        // flag: 是否打印F(i, 0)所代表的(Xi, '-')
        int k;
        pUnit p = a[i][j];
        if (! i) {  // 保证序列s的每个字符都比对上
            for (k = n - 1; k >= 0; k--)
                printf("%c", saln[k]);
            printf("\n");
            for (k = n - 1; k >= 0; k--)
                printf("%c", raln[k]);
            printf("\n\n");
            return;
        }
        if (! j) {  // F(i, 0)
            if (flag) {
                saln[n] = s[i - 1];
                raln[n] = UNMATCH_CHAR;
            }
            if (p->W0)
                printAlign(a, i - 1, 0, s, r, saln, raln, n + flag, 1);
            for (k = 0; k < p->nj; k++)
                if (p->Wj[k])
                    printAlign(a, i - 1, k + 1, s, r, saln, raln, n + flag, 1);
        } else {
            if (p->Wi & 1) {    // F(i, 0)
                printAlign(a, i, 0, s, r, saln, raln, n, 1);
            }
            if (p->Wi & 2) {    // 向上回溯一格
                saln[n] = s[i - 1];
                raln[n] = GAP_CHAR;
                printAlign(a, i - 1, j, s, r, saln, raln, n + 1, 1);
            }
            if (p->Wi & 4) {    // 向左上回溯一格
                saln[n] = s[i - 1];
                raln[n] = r[j - 1];
                printAlign(a, i - 1, j - 1, s, r, saln, raln, n + 1, 1);
            }
            if (p->Wi & 8) {    // 向左回溯一格
                saln[n] = GAP_CHAR;
                raln[n] = r[j - 1];
                printAlign(a, i, j - 1, s, r, saln, raln, n + 1, 1);
            }
            if (p->Wi & 16) {   // F(i, 0) + s(i, j)
                saln[n] = s[i - 1];
                raln[n] = r[j - 1];
                printAlign(a, i, 0, s, r, saln, raln, n + 1, 0);
            }
        }
    }
    
    void align(char *s, char *r, float t) {
        int i, j, k;
        int m = strlen(s);
        int n = strlen(r);
        float gap = -2.5;     // 对空位的罚分
        float tm[5];
        float em;   // F(m + 1, 0)
        pUnit **aUnit;
        char* salign;
        char* ralign;
        // 初始化
        if ((aUnit = (pUnit **) malloc(sizeof(pUnit*) * (m + 1))) == NULL) {
            fputs("Error: Out of space!\n", stderr);
            exit(1);
        }
        for (i = 0; i <= m; i++) {
            if ((aUnit[i] = (pUnit *) malloc(sizeof(pUnit) * (n + 1))) == NULL) {
                fputs("Error: Out of space!\n", stderr);
                exit(1);     
            }
            for (j = 0; j <= n; j++) {
                if ((aUnit[i][j] = (pUnit) malloc(sizeof(struct Unit))) == NULL) {
                    fputs("Error: Out of space!\n", stderr);
                    exit(1);     
                }
                aUnit[i][j]->Wi = 0;
                aUnit[i][j]->W0 = 0;
                // 创建F(i, 0)的跳转数组
                aUnit[i][j]->nj = n;
                if ((aUnit[i][j]->Wj = (int*) malloc(sizeof(int) * aUnit[i][j]->nj)) == NULL) {
                    fputs("Error: Out of space!\n", stderr);
                    exit(1);     
                }
                for (k = 0; k < aUnit[i][j]->nj; k++) {
                    aUnit[i][j]->Wj[k] = 0;
                }
            }
        }
        for (j = 0; j <= n; j++) {
            aUnit[0][j]->M = 0;
        }
        // 将字符串都变成大写
        strUpper(s);
        strUpper(r);
        // 动态规划算法计算得分矩阵每个单元的分值
        for (i = 1; i <= m; i++) {
            // 计算F(i, 0)
            aUnit[i][0]->M = aUnit[i - 1][0]->M;
            for (j = 1; j <= n; j++)
                if (aUnit[i][0]->M < aUnit[i - 1][j]->M - t)
                    aUnit[i][0]->M = aUnit[i - 1][j]->M - t;
            if (aUnit[i][0]->M == aUnit[i - 1][0]->M)
                aUnit[i][0]->W0 = 1;
            for (j = 1; j <= n; j++)
                if (aUnit[i][0]->M == aUnit[i - 1][j]->M - t)
                    aUnit[i][0]->Wj[j - 1] = 1;
            // 计算F(i, 1)
            tm[0] = aUnit[i][0]->M;
            tm[3] = aUnit[i][0]->M + gap;
            tm[4] = aUnit[i][0]->M + getFScore(s[i - 1], r[0]);
            aUnit[i][1]->M = max2(max2(tm[0], tm[3]), tm[4]);
            if (tm[0] == aUnit[i][1]->M)  aUnit[i][1]->Wi |= 1;
            if (tm[3] == aUnit[i][1]->M)  aUnit[i][1]->Wi |= 8;
            if (tm[4] == aUnit[i][1]->M)  aUnit[i][1]->Wi |= 16;
            // 计算F(i, j), j>=1
            for (j = 2; j <= n; j++) {
                tm[0] = aUnit[i][0]->M;
                tm[1] = aUnit[i - 1][j]->M + gap;
                tm[2] = aUnit[i - 1][j - 1]->M + getFScore(s[i - 1], r[j - 1]);
                tm[3] = aUnit[i][j - 1]->M + gap;
                tm[4] = aUnit[i][0]->M + getFScore(s[i - 1], r[j - 1]);
                aUnit[i][j]->M = maxArray(tm, 5);
                for (k = 0; k < 5; k++)
                    if (tm[k] == aUnit[i][j]->M)
                        aUnit[i][j]->Wi |= 1 << k;
            }
        }
        // 计算F(m + 1, 0)
        em = aUnit[m][0]->M;
        for (j = 1; j <= n; j++)
            if (em < aUnit[m][j]->M - t)
                em = aUnit[m][j]->M - t;
    /*
        // 打印得分矩阵
        for (i = 0; i <= m; i++) {
            for (j = 0; j <= n; j++)
                printf("%f ", aUnit[i][j]->M);
            printf("\n");
        }
    */
        printf("max score: %f\n", em);
        // 打印最优比对结果,如果有多个,全部打印
        // 递归法
        if (em == 0) {
            fputs("No seq aligned.\n", stdout);
        } else {
            if ((salign = (char*) malloc(sizeof(char) * (3 * m))) == NULL) {  // 3m的大小是根据match以及gap的分值来估算的。
                fputs("Error: Out of space!\n", stderr);
                exit(1);
            }
            if ((ralign = (char*) malloc(sizeof(char) * (3 * m))) == NULL) {
                fputs("Error: Out of space!\n", stderr);
                exit(1);
            }
            if (em == aUnit[m][0]->M)
                printAlign(aUnit, m, 0, s, r, salign, ralign, 0, 1);
            for (j = 1; j <= n; j++)
                if (em == aUnit[m][j]->M - t)
                    printAlign(aUnit, m, j, s, r, salign, ralign, 0, 1);
            // 释放内存
            free(salign);
            free(ralign);
        }
        for (i = 0; i <= m; i++) {
            for (j = 0; j <= n; j++) {
                free(aUnit[i][j]->Wj);
                free(aUnit[i][j]);
            }
            free(aUnit[i]);
        }
        free(aUnit);
    }
    

    (公众号:生信了)

    相关文章

      网友评论

          本文标题:序列匹配(五)——重复匹配问题的动态规划算法

          本文链接:https://www.haomeiwen.com/subject/djlwbctx.html