序列比对(六)——交叉匹配问题

原创: hxj7

之前几篇文章介绍了全局匹配以及局部匹配，本文介绍交叉匹配问题并给出代码。

交叉匹配

所谓交叉匹配（overlap alignment 或者叫 glocal alignment），就是两条序列中至少有一条的头部序列要参加比对并且至少有一条的尾部序列要参加比对。
一般而言，就是下面两种情形：
一种是两条序列有重叠的部分，但互不包含。比如x序列的头部与y序列的尾部匹配。

第二种是一条序列包含另一条序列，比如x序列包含y序列。

与全局联配、局部联配的区别

那么如何得到这种匹配的最高匹配得分呢？首先得弄清楚此类匹配与全局联配、局部联配的区别。
不同于全局匹配，交叉匹配中两端的序列可以不参与联配（或者说不罚分）。
不同于局部匹配，交叉匹配中某一条序列的头部必须参与联配且某一条序列的尾部必须参加联配。

交叉匹配的算法

假设x序列和y序列的长度分别是m和n，根据上面的比较可以得到解决交叉匹配问题的关键步骤（依然是利用得分矩阵）：

1. 设置F(0, 0) = 0。F(i, 0) = 0; i = 1, 2, …, m。F(0, j) = 0; j = 1, 2, …, n。
2. 得分矩阵最后一行的最大值以及最后一列的最大值中的较大者为Fmax。
3. 从Fmax所在单元开始回溯，到达i = 0的行（代表x序列的头部）或者j = 0的列（代表y序列的头部）为止。
   得分矩阵的计算同全局联配：
   
   图片引自《生物序列分析》
   具体的代码参见文末的 alnOverlap.c 代码段。

限定x头部与y尾部的交叉匹配

如果在给出交叉匹配问题的同时，限定x头部以及y尾部的序列必须参与联配，又该如何计算最高得分呢？
注意此时F(i,0)=i*gap。也就是说，x序列的头部联配空位要罚分了。
此时比对的最高得分就是最后一列的最大值。至于回溯其实也很简单，这个时候只需要从最后一列的最大值所在单元开始回溯，到i = 0的行结束。
或者，根据以下计算公式得到最高得分：

图片改编自《生物序列分析》

其中X(i)表示序列中x(i)不参与联配的情况下联配的最高得分，X(0) = 0。T是得分阈值。比对的最高得分就是X(m+1)。回溯方式是从X(m)（如果其值等于X(m+1)）或者F(m,n)（如果其值等于X(m+1)）开始，到X(0)或者F(0,k)结束。
这种限定条件下的交叉联配代码参见文末的 alnOverlaoX.c 代码段。

代码

首先是 alnOverlap.c 代码段

    #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  是否往左回溯一格
    // 00000010  是否往左上回溯一格
    // 00000100  是否往上回溯一格
    struct Unit {
   
        int W;     // F(i, j)的回溯指标，不同的bit代表不同的回溯方式
        float M;      // 得分矩阵第(i, j)这个单元的分值，即序列s(1,...,i)与序列r(1,...,j)比对的最高得分
    };
    typedef struct Unit *pUnit;
    
    void strUpper(char *s);
    float maxOfRow(pUnit* a, int n);
    float maxOfCol(pUnit** a, int col, 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);
    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 maxOfRow(pUnit* a, int n) {
   
        int i;
        float max = a[0]->M;
        for (i = 1; i < n; i++)
            if (a[i]->M > max)
                max = a[i]->M;
        return max;
    }
    
    float maxOfCol(pUnit** a, int col, int n) {
   
        int i;
        float max = a[0][col]->M;
        for (i = 1; i < n; i++)
            if (a[i][col]->M > max)
                max = a[i][col]->M;
        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 k;
        pUnit p = a[i][j];
        if (! i) {
     // 到达s序列的头部
            for (k = j; k >= 1; k--, n++) {
    // r序列头部未匹配的部分
                saln[n] = UNMATCH_CHAR;
                raln[n] = r[k - 1];
            }
            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;
        } else if (! j) {
     // 到达r序列的头部
            for (k = i; k >= 1; k--, n++) {
    // s序列头部未匹配的部分
                raln[n] = UNMATCH_CHAR