序列比对（十八）——重复匹配问题的补充说明

2019-11-15 本文已影响0人生信了

原创：hxj7

前文介绍了重复匹配问题的动态规划算法，但是遗留了重复结果输出的问题。本文对该问题进行了补充说明。

前文《序列匹配(五)——重复匹配问题的动态规划算法》介绍了重复匹配问题的动态规划算法。重复匹配问题就是从序列 $x$ 中找到序列 $y$ 的多个拷贝，并且使这多个拷贝的“标准比对分值”之和最大。《生物序列分析》一书中给出的迭代公式是：
$\begin{aligned} & F(i,0) = \max \begin{cases} F(i-1,0), \\ F(i-1,j)-T, \ \ \ \ \ j=1,2,...,n \end{cases} \\ & F(i,j) = \max \begin{cases} F(i,0) \\ F(i-1,j-1) + s(i,j) \\ F(i-1,j) + d \\ F(i, j-1) + d \end{cases} \end{aligned}$

经过笔者分析，在前文中给出的迭代公式是：
$\begin{aligned} & F(i,0) = \max \begin{cases} F(i-1,0), \\ F(i-1,j)-T, \ \ \ \ \ j=1,2,...,n \end{cases} \\ & F(i,j) = \max \begin{cases} F(i,0) \\ F(i-1,j-1) + s(i,j) \\ F(i-1,j) + d \\ F(i, j-1) + d \\ F(i,0) + s(i,j) \end{cases}\ \ \ , \ \ \ \text{when $j > 1$} \\ & F(i,j) = \max \begin{cases} F(i,0) \\ F(i-1,j) + d \\ F(i,0) + s(i,j) \end{cases}\ \ \ , \ \ \ \text{when $j = 1$} \end{aligned}$

但是这个公式在回溯时会出现重复结果输出的问题，比如：

经过笔者近期回顾，发现重复结果输出的可能途径有二（当然可能还有其他途径）：
一是：
$\begin{aligned} & F(i,0) \rightarrow F(i-1,j)-T \rightarrow F(i-1,0) \\ & F(i,0) \ \ \ \ \ \ \ \ \ \ \ \ \ \longrightarrow \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \, F(i-1,0) \end{aligned}$

这两种回溯步骤结果一致，但意义互相矛盾。二者只能取其一。

二是：
$F(i,j) \rightarrow F(i-1,j-1)+s(i,j) \rightarrow F(i-1,0)$

第一个箭头表示的回溯步骤要求 $x_{i-1}$ 参与联配，而第二个箭头表示的回溯步骤要求 $x_{i-1}$ 不参与联配。二者互相矛盾，所以该回溯不成立。

校正公式和代码

笔者认为，上述重复结果输出的问题主要可能是因为 $F(i,j)$ 的定义不清晰，导致 $F(i,j)$ 既能表示 $x_i$ 参与联配，也能表示 $x_i$ 不参与联配。
为此，笔者给出新的更清晰的符号定义以及公式。 $F(i)$ 表示以 $x_i$ 结尾且 $x_i$ 不参与联配的最高得分。 $A(i,j)$ 表示以 $x_i$ 和 $y_j$ 结尾且 $x_i$ 参与联配的最高得分。
那么迭代公式是：
$\begin{aligned} & F(0) = 0, \\ & F(i) = \max \begin{cases} F(i-1), \\ A(i-1,j)-T, \ \ \ \ \ j=1,2,...,n \end{cases} \\ & A(0,0) = d, \\ & A(0,j) = d, \ \ \ \ \ j =1,2,...,n \\ & A(i,0) = F(i) + d, \\ & A(i,j) = \max \begin{cases} F(i) + s(i,j), \\ A(i-1,j-1) + s(i,j), \\ A(i-1,j) + d, \\ A(i,j-1) + d, \end{cases} \end{aligned}$

需要说明的是：

上述公式中令 $A(0,0)=d$ 是为了让 $A(0,0) \neq F(0)$ ，以免产生重复结果输出。
$A(i,0)$ 不光意味着 $x_i$ 和一个空位联配，且代表了一个“旧联配段”的结束和一个“新联配段”的开始，只不过这个新联配是以空位开始的。
$A(i,j)$ 表示 $x_i$ 参与联配。其迭代公式中第一项表示 $x_i$ 开始了一个新的局部联配；而后三项表示 $x_i$ 加入（或者说延续）前一个局部联配。
最终比对的最高得分就是 $F(m+1)$ 。

这样的公式目前还没有出现重复结果输出的问题：

相应的代码放在了文末。

对比对总长度的估计

在写代码的过程中，遇到一个小问题，就是要对最优比对结果的总长度（以下简称比对总长度）进行估计。即 $x$ 序列长度为 $m$ ， $y$ 序列长度为 $n$ 。那么比对总长度 $l$ 是在什么范围内呢？这里说的总长度包括了 $x$ 序列中未参与联配的碱基。
由于 $x$ 序列每个碱基在迭代的过程中都要被遍历到，所以 $l$ 的下界是 $m$ （比如考虑 $x$ 序列是 $GACGACT$ ，而 $y$ 序列是 $GACT$ ，那么 $l=m=7$ ）：
$l \ge m$

假设比对时 $match$ 的得分是 $M$ ，而 $gap$ 的得分是 $G$ 。那么，总长度的一个宽松上界是：
$l \le m + \left[ m \left| \frac{M}{G} \right| \right]$

其中， $\left| a \right|$ 表示 $a$ 的绝对值， $\left[ a \right]$ 表示不大于 $a$ 的最大整数。上述不等式证明如下：

假设将长度为 $m$ 的 $x$ 序列分割成一系列子序列 $S_i$ ，这些子序列构成一个完整分割集 $S$ ，即
$\bigcup_{i} S_i = S \ \ \ \text{and} \ \ \ \bigcap_{i} S_i = \emptyset$
如果将 $S_i$ 的长度记为 $m_i$ ，则有：
$0 < m_i \le m \ \ \ \text{and} \ \ \ \sum_i m_i = m$
并且， $S_i$ 要么是“非联配段”，要么是“联配段”。所谓“非联配段”，即该段序列中的所有碱基都不参与最终的联配；所谓“联配段”，即该段序列中的所有碱基都参与最终的联配。
接着，将每个子序列 $S_i$ 与长度为 $n$ 的序列 $y$ 进行局部比对，比对的长度记为 $L(S_i,y)$ 。
如果 $S_i$ 是“联配段”，由于是局部比对，假设比对时 $match$ 的得分是 $M$ ，而 $gap$ 的得分是 $G$ （这里假设都是简单情形，即所有符号的 $match$ 得分都是一样的，所有符号的 $gap$ 得分也都是一样的，且空位罚分是线性的），
则有：
$L(S_i,y) \le m_i + \left[ \left| \frac{m_{i}M}{G} \right| \right]$
其中， $\left| a \right|$ 表示 $a$ 的绝对值， $\left[ a \right]$ 表示不大于 $a$ 的最大整数。上述不等式是一个宽松上界，证明也很简单，其关键就是局部比对的总分值不能小于0，以至于会限制空位的数量不得高于 $\left[ \left| \frac{m_{i}M}{G} \right| \right]$ 。
如果 $S_i$ 是“非联配段”，那么：
$\begin{aligned} L(S_i,y) & = m_i \\ & \le m_i + \left[ \left| \frac{m_{i}M}{G} \right| \right] \end{aligned}$
所以，不管 $S_i$ 是“非联配段”，还是“联配段”，都有：
$L(S_i,y) \le m_i + \left[ \left| \frac{m_{i}M}{G} \right| \right]$
那么，这些子序列的比对结果拼接起来的长度为：
$\begin{aligned} \sum_i L(S_i, y) & \le \sum_i \left( m_i + \left[ \left| \frac{m_{i}M}{G} \right| \right] \right) \\ & \le \sum_i m_i + \left[ \sum_i \left| \frac{m_{i}M}{G} \right| \right] \\ & = m + \left[ \sum_i m_{i} \left| \frac{M}{G} \right| \right] \\ & = m + \left[ m \left| \frac{M}{G} \right| \right] \end{aligned}$
我们知道，重复匹配问题其实就是将 $x$ 序列分割成一段段互不交叠的子序列，然后每一段子序列都和 $y$ 序列进行局部比对。
假设重复匹配问题最终的比对结果是将 $x$ 序列按上面所述完整分割集分割成一系列子序列 $\hat{S}_i$ ，那么重复匹配问题的比对总长度当然也服从上述不等式，即
$\begin{aligned} l & = \sum_i L(\hat{S}_i, y) \\ & \le m + \left[ m \left| \frac{M}{G} \right| \right] \end{aligned}$

总的来说，重复匹配问题的比对总长度（在上述简单情形下）在以下范围内：
$m \le l \le m + \left[ m \left| \frac{M}{G} \right| \right]$

实现代码

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define MAXSEQ 1000
#define GAP_CHAR '-'
#define UNALIGN_CHAR '.'
#define max2(a, b) ((a) > (b) ? (a) : (b))

// 对空位的罚分是线性的
struct FUnit {
    int W0;   // X{i-1}不参与联配
    int* Wj;      // 跳转到A(i - 1, j)
    int nj;       // Wj数组的大小
    float M;      // F(i,0)的值
};
typedef struct FUnit* pFUnit;

// 00000001  F(i, 0) + s(i, j)
// 00000010  是否往左回溯一格
// 00000100  是否往左上回溯一格
// 00001000  是否往上回溯一格
struct AUnit {
    int Wi;      // 不同的bit代表不同的回溯方式
    float M;
};
typedef struct AUnit* pAUnit;

float gap = -2.5;     // 对空位的罚分
float match = 5;
float unmatch = -4;

void strUpper(char *s);
float maxArray(float *a, int n);
float getFScore(char a, char b);
void printFAlign(pFUnit* f, pAUnit** a, const int i, char* s, char* r, char* saln, char* raln, int n, int flag);
void printAAlign(pFUnit* f, pAUnit** 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 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 printFAlign(pFUnit* f, pAUnit** a, const int i, char* s, char* r, char* saln, char* raln, int n, int flag) {
    // flag: 是否是从F(i+1,0)跳转过来，而不是从F(i, 0) + s(i, j)跳转过来
    int k;
    pFUnit p = f[i];
    //printf("i=%d, n=%d, flag=%d\n", i, n, flag);
    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 (flag) {
        saln[n] = s[i - 1];
        raln[n] = UNALIGN_CHAR;
    }
    if (p->W0) 
        printFAlign(f, a, i - 1, s, r, saln, raln, n + flag, 1);
    for (k = 0; k < p->nj; k++)
        if (p->Wj[k])
            printAAlign(f, a, i - 1, k + 1, s, r, saln, raln, n + flag);
}

void printAAlign(pFUnit* f, pAUnit** a, const int i, const int j, char* s, char* r, char* saln, char* raln, int n) {
    int k;
    pAUnit p = a[i][j];
    //printf("i=%d, j=%d, n=%d\n");
    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) + d
        saln[n] = s[i - 1];
        raln[n] = GAP_CHAR;
        printFAlign(f, a, i, s, r, saln, raln, n + 1, 0);
    } else {
        if (p->Wi & 1) {    // F(i, 0) + s(i, j)
            saln[n] = s[i - 1];
            raln[n] = r[j - 1];
            printFAlign(f, a, i, s, r, saln, raln, n + 1, 0);
        }
        if (p->Wi & 2) {    // 向上回溯一格
            saln[n] = s[i - 1];
            raln[n] = GAP_CHAR;
            printAAlign(f, a, i - 1, j, s, r, saln, raln, n + 1);
        }
        if (p->Wi & 4) {    // 向左上回溯一格
            saln[n] = s[i - 1];
            raln[n] = r[j - 1];
            printAAlign(f, a, i - 1, j - 1, s, r, saln, raln, n + 1);
        }
        if (p->Wi & 8) {    // 向左回溯一格
            saln[n] = GAP_CHAR;
            raln[n] = r[j - 1];
            printAAlign(f, 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 tm[4];
    float em;   // F(m + 1, 0)
    pFUnit* fUnit;
    pAUnit** aUnit;
    char* salign;
    char* ralign;
    int alnSize = m + floor(m * abs(match / gap));
    // 初始化
    if ((fUnit = (pFUnit*) malloc(sizeof(pFUnit) * (m + 1))) == NULL || \
        (aUnit = (pAUnit**) malloc(sizeof(pAUnit*) * (m + 1))) == NULL) {
        fputs("Error: Out of space!\n", stderr);
        exit(1);
    }
    for (i = 0; i <= m; i++) {
        if ((fUnit[i] = (pFUnit) malloc(sizeof(struct FUnit))) == NULL) {
            fputs("Error: Out of space!\n", stderr);
            exit(1);            
        }
        fUnit[i]->W0 = 0;
        fUnit[i]->nj = n;
        // 创建F(i, 0)的跳转数组
        if ((fUnit[i]->Wj = (int*) malloc(sizeof(int) * fUnit[i]->nj)) == NULL) {
            fputs("Error: Out of space!\n", stderr);
            exit(1);     
        }
        for (k = 0; k < fUnit[i]->nj; k++) {
            fUnit[i]->Wj[k] = 0;
        }
    }
    for (i = 0; i <= m; i++) {
        if ((aUnit[i] = (pAUnit *) malloc(sizeof(pAUnit) * (n + 1))) == NULL) {
            fputs("Error: Out of space!\n", stderr);
            exit(1);     
        }
        for (j = 0; j <= n; j++) {
            if ((aUnit[i][j] = (pAUnit) malloc(sizeof(struct AUnit))) == NULL) {
                fputs("Error: Out of space!\n", stderr);
                exit(1);     
            }
            aUnit[i][j]->Wi = 0;
        }
    }
    fUnit[0]->M = 0;
    aUnit[0][0]->M = gap;      // 注意这里设置A(0,0) != 0 是很有必要的，否则A(0,0)=F(0,0)会导致重复结果输出
    for (j = 1; j <= n; j++)
        aUnit[0][j]->M = gap;
    // 将字符串都变成大写
    strUpper(s);
    strUpper(r);
    // 动态规划算法计算得分矩阵每个单元的分值
    for (i = 1; i <= m; i++) {
        // 计算F(i, 0) i >= 1
        fUnit[i]->M = fUnit[i - 1]->M;
        for (j = 1; j <= n; j++)
            if (fUnit[i]->M < aUnit[i - 1][j]->M - t)
                fUnit[i]->M = aUnit[i - 1][j]->M - t;
        if (fUnit[i]->M == fUnit[i - 1]->M)
            fUnit[i]->W0 = 1;
        for (j = 1; j <= n; j++)
            if (fUnit[i]->M == aUnit[i - 1][j]->M - t)
                fUnit[i]->Wj[j - 1] = 1;
        // 计算A(i, 0) i >= 1
        aUnit[i][0]->M = fUnit[i]->M + gap;
        aUnit[i][0]->Wi |= 1;
        // 计算A(i, j) i >= 1 and j >= 1
        for (j = 1; j <= n; j++) {
            tm[0] = fUnit[i]->M + getFScore(s[i - 1], r[j - 1]);
            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;
            aUnit[i][j]->M = maxArray(tm, 4);
            for (k = 0; k < 4; k++)
                if (tm[k] == aUnit[i][j]->M)
                    aUnit[i][j]->Wi |= 1 << k;
        }
    }
    // 计算F(m + 1, 0)
    em = fUnit[m]->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++)
        printf("%f ", fUnit[i]->M);
    printf("\n");
    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) * alnSize)) == NULL) {
            fputs("Error: Out of space!\n", stderr);
            exit(1);
        }
        if ((ralign = (char*) malloc(sizeof(char) * alnSize)) == NULL) {
            fputs("Error: Out of space!\n", stderr);
            exit(1);
        }
        if (em == fUnit[m]->M)
            printFAlign(fUnit, aUnit, m, s, r, salign, ralign, 0, 1);
        for (j = 1; j <= n; j++)
            if (em == aUnit[m][j]->M - t)
                printAAlign(fUnit, 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]);
        free(aUnit[i]);
    }
    free(aUnit);
    for (i = 0; i <= m; i++) {
        free(fUnit[i]->Wj);
        free(fUnit[i]);
    }
    free(fUnit);
}

（公众号：生信了）

序列比对（十八）——重复匹配问题的补充说明

校正公式和代码

对比对总长度的估计

实现代码

猜你喜欢

热点阅读