序列比对（22）中间字符串分支定界方法中更紧的界

前文介绍了中间字符串的算法和代码，但是使用分支定界策略时所使用的界限是很宽松的。本文给出了一个更紧的界限。

对分支定界法的简单回顾

前文《序列比对（21）中间字符串问题的算法及实现代码》介绍了中间字符串的算法和代码，但是使用分支定界策略时所使用的界限是很宽松的。分支定界法的伪代码如下：

引自《生物信息学算法导论》

其中关键的是这几句：
$$\begin{array}{rl}& optimisticDistance\leftarrow TotalDistance(prefix,DNA)\\ & \text{if }optimisticDistance>bestDistance\\ & byPass(...)\end{array}$$

上述给出的界限是$TotalDistance(prefix,DNA)$，这是一个很宽松的界限，为了解释这一点，下面我们将给出一些符号和推导进行详细说明。

对分支定界法界限的详细说明

我们知道，中间字符串问题的分支定界算法的核心就是构建一棵深度为$l$的搜索树，然后在遍历树的非叶顶点的时候，进行预判，如果该顶点在所有可能情况下所能得到的最小总距离仍比目前整体的最小总距离还要大，那么就跳过以该顶点为根的分支树。

对于一个深度为$i$的顶点$D_i$，从根到$D_i$的路径构成一个长度为$i$的$wordu$，以$D_i$为根到任意叶顶点的路径构成一个长度为$l-i$的$wordv$。并且，$u+v$构成一个完整的$l$-元组。那么，如果我们将任意$wordw$的长度记为$len(w)$，则有：
$$\begin{array}{rl}& len(u)=i\\ & len(v)=l-i\end{array}$$

对于一个深度为$i$的特定的顶点$D_i$，当遍历到此顶点时，$wordu$是确定的，就一种可能；而$wordv$则有$4^{l-i}$种可能。如果我们将$TotalDistance$简记为$TD$，那么$u+v$总距离的最小值$optimisticDistance$可以记为：
(1) o p t i m i s t i c D i s t a n c e = min ⁡ { v ∣ l e n ( v ) = l − i }   T D ( u + v , D N A ) optimisticDistance = \underset{\{v|len(v) = l - i\}}{\min} \, TD(u+v, DNA) \tag{1} optimisticDistance={v∣len(v)=l−i}min​TD(u+v,DNA)(1)

我们可以很容易地推断出：
$$TD(u+v,DNA)\ge TD(u,DNA)+TD(v,DNA)\text{\hspace{1em}(2)}$$

如果我们令：
v ∗ = a r g m i n { v ∣ l e n ( v ) = l − i }   T D ( u + v , D N A ) v^* = \underset{\{v|len(v) = l - i\}}{\mathrm{argmin}} \, TD(u + v, DNA) v∗={v∣len(v)=l−i}argmin​TD(u+v,DNA)

那么：
(3) o p t i m i s t i c D i s t a n c e = min ⁡ { v ∣ l e n ( v ) = l − i }   T D ( u + v , D N A ) = T D ( u + v ∗ , D N A ) ≥ T D ( u , D N A ) + T D ( v ∗ , D N A ) ≥ T D ( u , D N A ) + min ⁡ { v ∣ l e n ( v ) = l − i }   T D ( v , D N A ) \begin{aligned} optimisticDistance & = \underset{\{v|len(v) = l - i\}}{\min} \, TD(u + v, DNA) \\ & = TD(u + v^*, DNA) \\ & \ge TD(u, DNA) + TD(v^*, DNA) \\ & \ge TD(u, DNA) + \underset{\{v|len(v) = l - i\}}{\min} \, TD(v, DNA) \end{aligned} \tag{3} optimisticDistance​={v∣len(v)=l−i}min​TD(u+v,DNA)=TD(u+v∗,DNA)≥TD(u,DNA)+TD(v∗,DNA)≥TD(u,DNA)+{v∣len(v)=l−i}min​TD(v,DNA)​(3)

如果要跳过以$D_i$为根的分支树，那么势必$optimisticDistance$的值要比目前的总体最小值要大，即：
$$optimisticDistance>bestDistance\text{\hspace{1em}(4)}$$

假设有：
(5) T D ( u , D N A ) + min ⁡ { v ∣ l e n ( v ) = l − i }   T D ( v , D N A ) > b e s t D i s t a n c e TD(u, DNA) + \underset{\{v|len(v) = l - i\}}{\min} \, TD(v, DNA) > bestDistance \tag{5} TD(u,DNA)+{v∣len(v)=l−i}min​TD(v,DNA)>bestDistance(5)

那么根据式（3），式（4）就一定成立。

也就是说在顶点$D_i$处利用分支定界法进行预判时，$optimisticDistance$的一个比较紧的下界是 T D ( u , D N A ) + min ⁡ { v ∣ l e n ( v ) = l − i }   T D ( v , D N A ) TD(u, DNA) + \underset{\{v|len(v) = l - i\}}{\min} \, TD(v, DNA) TD(u,DNA)+{v∣len(v)=l−i}min​TD(v,DNA)。说到底，所谓$D_i$处分支定界法的下界其实就是通过式（3）得出的。式（3）的不等式越多，说明分支定界法的下界越宽松。

比如，我们可以进一步放宽式（3）中的下界：
(6) o p t i m i s t i c D i s t a n c e ≥ T D ( u , D N A ) + min ⁡ { v ∣ l e n ( v ) = l − i }   T D ( v , D N A ) ≥ T D ( u , D N A ) \begin{aligned} optimisticDistance & \ge TD(u, DNA) + \underset{\{v|len(v) = l - i\}}{\min} \, TD(v, DNA) \\ & \ge TD(u, DNA) \end{aligned} \tag{6} optimisticDistance​≥TD(u,DNA)+{v∣len(v)=l−i}min​TD(v,DNA)≥TD(u,DNA)​(6)

因为这里的$TD(u,DNA)$其实就是伪代码中的$TotalDistance(prefix,DNA)$，所以我们可以看出原先的伪代码中的下界是非常宽松的。

也就是说伪代码中的一段应该改为：
o p t i m i s t i c D i s t a n c e ← T D ( u , D N A ) + min ⁡ { v ∣ l e n ( v ) = l − i }   T D ( v , D N A ) if    o p t i m i s t i c D i s t a n c e > b e s t D i s t a n c e b y P a s s ( . . . ) \begin{aligned} & optimisticDistance \leftarrow TD(u, DNA) + \underset{\{v|len(v) = l - i\}}{\min} \, TD(v, DNA) \\ & \textbf{if \ } optimisticDistance > bestDistance \\ & \qquad byPass(...) \end{aligned} ​optimisticDistance←TD(u,DNA)+{v∣len(v)=l−i}min​TD(v,DNA)if  optimisticDistance>bestDistancebyPass(...)​

进一步讨论

既然上面已经给出一个更紧的下界，关键就是如何计算 min ⁡ { v ∣ l e n ( v ) = l − i }   T D ( v , D N A ) \underset{\{v|len(v) = l - i\}}{\min} \, TD(v, DNA) {v∣len(v)=l−i}min​TD(v,DNA)。如果精确计算的话，可能比较费时。从实现代码的运行效率考虑，所以我们可以试着去稍微放宽其下界。

依照上面的推导，我们可以很容易得出以下结论：
(7) min ⁡ { v ∣ l e n ( v ) = l − i }   T D ( v , D N A ) ≥ min ⁡ { x ∣ l e n ( x ) = j }   T D ( x , D N A ) + min ⁡ { y ∣ l e n ( y ) = l − i − j }   T D ( y , D N A ) , ∀   j ∈ ( 0 , l − i ) \begin{aligned} & \underset{\{v|len(v) = l - i\}}{\min} \, TD(v, DNA) \ge \\ & \underset{\{x|len(x) = j\}}{\min} \, TD(x, DNA) + \underset{\{y|len(y) = l - i - j\}}{\min} \, TD(y, DNA), \quad \forall \, j \in(0, l - i) \end{aligned} \tag{7} ​{v∣len(v)=l−i}min​TD(v,DNA)≥{x∣len(x)=j}min​TD(x,DNA)+{y∣len(y)=l−i−j}min​TD(y,DNA),∀j∈(0,l−i)​(7)

其中：
$$v=x+y$$

所以根据式（7），就有很多种放宽 min ⁡ { v ∣ l e n ( v ) = l − i }   T D ( v , D N A ) \underset{\{v|len(v) = l - i\}}{\min} \, TD(v, DNA) {v∣len(v)=l−i}min​TD(v,DNA)的方法。比如：
$$\begin{array}{rl}\text{let }& f=\left[\frac{l-i}{2}\right]\\ & m=l-i-f\times 2\end{array}$$

其中，$[a]$表示不大于$a$的最大整数。那么：
min ⁡ { v ∣ l e n ( v ) = l − i }   T D ( v , D N A ) ≥ 2 × min ⁡ { x ∣ l e n ( x ) = f }   T D ( x , D N A ) + min ⁡ { y ∣ l e n ( y ) = m }   T D ( y , D N A ) \begin{aligned} & \underset{\{v|len(v) = l - i\}}{\min} \, TD(v, DNA) \ge \\ & \qquad 2 \times \underset{\{x|len(x) = f\}}{\min} \, TD(x, DNA) + \underset{\{y|len(y) = m\}}{\min} \, TD(y, DNA) \end{aligned} ​{v∣len(v)=l−i}min​TD(v,DNA)≥2×{x∣len(x)=f}min​TD(x,DNA)+{y∣len(y)=m}min​TD(y,DNA)​

这样，可以通过减少$word$的长度来减少运算量。当然，这样做的坏处就是式（3）中的界限被放宽了。界限越紧，就可能跳过更早地跳过更大的分支树，从而减少运算量。从这两点来看，最终计算 min ⁡ { v ∣ l e n ( v ) = l − i }   T D ( v , D N A ) \underset{\{v|len(v) = l - i\}}{\min} \, TD(v, DNA) {v∣len(v)=l−i}min​TD(v,DNA)的方案很可能是一个折中的方案。

运行效果

笔者按照上述方案选择了一种更紧的界限及其计算方式，从代码的实际运行效果来看，对效率的提升并不大。这可能和输入数据有关，当然也可能是笔者的实现代码还需优化（比如选择更好的计算界限的方法）。尽管如此，通过本文第一次尝试了为分支定界法估计更紧的界，这也许为以后的学习打下了一个基础。

具体代码

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>
#define MAXSEQ 1000
#define max(a, b) ((a) > (b) ? (a) : (b))

struct Sequence {
    char* s;
    int* a;      // 该向量存储每个字符在字符集中的序号
    int l;
};
typedef struct Sequence* Seq;

const int nalpha = 4;     // 字符种数
const char sortedAlpha[] = {'A', 'C', 'G', 'T'};    // 排过序的字符集

char* sec2time(time_t t);
void strUpper(char *s);
int search(char c);                      /* 在排序后的字符集中寻找字符c，如果找到返回序号；找不到返回-1 */
Seq* readSeq(char* filename, const int t);          /* 从文件中读取多条序列 */
Seq create(char* s, const int n, const int i);  
void destroy(Seq seq);
int nextVertex(int* a, int d, const int l);     /* 得到下一个顶点 */
int byPass(int* a, int d, const int l);           /* 跳过某个分支，得到下一个顶点 */
int getDistance(Seq* mulSeq, const int t, int* a, const int l, const int d);        /* 计算一个l-元组和DNA的最小总距离 */
int* getStartPoints(Seq* mulSeq, const int t, int* a, const int l);    /* 获取起始位点向量 */
int getLowerBound(Seq* mulSeq, const int t, int* a, int* b, const int l);      /* 计算TotalDistance(v, DNA)的下界 */
void findMedianStr(Seq* mulSeq, const int t, const int l);           /* 寻找中间字符串 */

int main(void) {
    time_t start, end;
    char* tstr;
    int i, t, l;
    char filename[30];
    Seq* mulSeq;
    start = time(NULL);  // 开始时间
    printf("t: ");
    scanf(" %d", &t);
    printf("l: ");
    scanf(" %d", &l);
    printf("seq_file path: ");
    scanf(" %s", filename);
    mulSeq = readSeq(filename, t);
    findMedianStr(mulSeq, t, l);
    for (i = 0; i < t; i++)
        destroy(mulSeq[i]);
    free(mulSeq);
    end = time(NULL);   // 结束时间
    tstr = sec2time(end - start);
    printf("time spent: %s\n", tstr);
    free(tstr);
    return 0;
}

char* sec2time(time_t t) {
/* 将秒数转化为时间字符串 */
    int h, m, s;
    char* tstr;
    if ((tstr = malloc(sizeof(char) * 50)) == NULL) {
        fputs("Error: out of space!\n", stderr);
        exit(8);
    }
    h = t / 3600;
    m = (t % 3600) / 60;
    s = t % 60;
    if (sprintf(tstr, "%02d:%02d:%02d", h, m, s) > 0)
        return tstr;
    else {
        sprintf(tstr, "more than %d hours", h);
        return tstr;
    }
}

void strUpper(char *s) {
    while (*s != '\0') {
        if (*s >= 'a' && *s <= 'z')
            *s -= 32;
        s++;
    }
}

int search(char c) {
/* 在排序后的字符集中寻找字符c，如果找到返回序号；找不到返回-1 */
    int i;
    int th = 5;
    int h, l, m;
    if (nalpha < th) {
        for (i = 0; i < nalpha; i++)
            if (sortedAlpha[i] == c)
                return i;
        return -1;
    } else {
        l = 0;
        h = nalpha - 1;
        while (l <= h) {
            m = (l + h) / 2;
            if (sortedAlpha[m] == c)
                return m;
            else if (sortedAlpha[m] < c)
                l = m + 1;
            else 
                h = m - 1;
        }
        return -1;
    }
}

Seq* readSeq(char* filename, const int t) {
/* 从文件中读取多条序列 */
    int i, n;
    char buf[MAXSEQ];
    FILE* fp;
    Seq* mulSeq;
    if ((fp = fopen(filename, "r")) == NULL) {
        fprintf(stderr, "Error: cannot open '%s'\n", filename);
        exit(1);
    }
    if ((mulSeq = (Seq*) malloc(sizeof(Seq) * t)) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);        
    }
    for (i = 0; i < t; i++) {
        if (fgets(buf, MAXSEQ, fp) == NULL) {
            fprintf(stderr, "Error: reading seq from '%s'\n", filename);
            exit(2);
        }
        n = strlen(buf);
        if (buf[n - 1] != '\n' && ! feof(fp)) {
            fprintf(stderr, "Error: line %d is too long.\n", i + 1);
            exit(3);
        }
        if (! feof(fp)) {
            buf[n - 1] = '\0';
            mulSeq[i] = create(buf, n - 1, i);
        } else {
            mulSeq[i] = create(buf, n, i);
        }
    }
    fclose(fp);
    return mulSeq;
}

Seq create(char* s, const int n, const int i) {
/* n: s的长度，不包括末尾的\0 */
    Seq seq;
    int j;
    if ((seq = (Seq) malloc(sizeof(struct Sequence))) == NULL || \
        (seq->s = (char*) malloc(sizeof(char) * (n + 1))) == NULL || \
        (seq->a = (int*) malloc(sizeof(int) * n)) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);
    }
    strcpy(seq->s, s);
    // 将序列中所有字符变为大写
    strUpper(seq->s);
    // 找到每个字符在字符集中的序号。如果是在后期做这个工作，那么就有很多重复，降低效率。
    for (j = 0; j < n; j++)
        if ((seq->a[j] = search(seq->s[j])) < 0) {
            fprintf(stderr, "Error: char %d of seq %d is %c, which is out of range.\n", j + 1, i + 1, seq->s[j]);
            exit(16);
        }
    seq->l = n;
    return seq;
}

void destroy(Seq seq) {
    free(seq->a);
    free(seq->s);
    free(seq);
}

int nextVertex(int* a, int d, const int l) {
/* 得到下一个顶点 */
    int i;
    if (d < l - 1) {
        a[d + 1] = 0;
        return d + 1;
    }
    for (i = l - 1; i >= 0; i--) {
        if (a[i] < nalpha - 1) {
            a[i]++;
            return i;
        }
    }
    return -1;
}

int byPass(int* a, int d, const int l) {
/* 跳过某个分支，得到下一个顶点 */
    int i;
    for (i = d; i >= 0; i--) {
        if (a[i] < nalpha - 1) {
            a[i]++;
            return i;
        }
    }
    return -1;
}        

int getDistance(Seq* mulSeq, const int t, int* a, const int l, const int d) {
/* 计算一个l-元组和DNA的最小总距离 */
    int i, j, k;
    int dist, minDist, totalDist;
    for (i = 0, totalDist = 0; i < t; i++) {
        for (j = 0, minDist = l; j <= mulSeq[i]->l - l; j++) {
            for (k = 0, dist = 0; k <= d; k++)
                dist += a[k] == mulSeq[i]->a[j + k] ? 0 : 1;
            if (dist < minDist)
                minDist = dist;
        }
        totalDist += minDist;
    }
    return totalDist;
}       

int* getStartPoints(Seq* mulSeq, const int t, int* a, const int l) {
/* 获取起始位点向量 */
    int* st;
    int i, j, k;
    int dist, minDist;
    if ((st = (int*) malloc(sizeof(int) * t)) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);
    }
    for (i = 0; i < t; i++) {
        for (j = 0, st[i] = 0, minDist = l; j <= mulSeq[i]->l - l; j++) {
            for (k = 0, dist = 0; k < l; k++)
                dist += a[k] == mulSeq[i]->a[j + k] ? 0 : 1;
            if (dist < minDist) {
                st[i] = j;
                minDist = dist;
            }
        }
    }
    return st;
}

int getLowerBound(Seq* mulSeq, const int t, int* a, int* b, const int l) {
/* 计算TotalDistance(v, DNA)的下界 */
    int i, d;
    int dist, minDist;
    // 遍历所有可能的l-元组，寻找最小距离
    for (i = 0; i < l; i++)
        a[i] = 0;
    minDist = l * t;
    d = 0;
    while (d >= 0) {
        if (d < l - 1) {
            dist = getDistance(mulSeq, t, a, l, d);
            if (dist + b[l - 2 - d] > minDist)
                d = byPass(a, d, l);
            else 
                d = nextVertex(a, d, l);
        } else {
            dist = getDistance(mulSeq, t, a, l, l - 1);
            if (dist < minDist) {  // 注意最小值对应的l-元组可能有多种，这里只保留了其中一种。其实还可以设置如果minDist = 0时终止搜索。
                minDist = dist;
                if (minDist == 0)
                    return 0;
            }
            d = nextVertex(a, d, l);
        }
    }
    return minDist;
}

void findMedianStr(Seq* mulSeq, const int t, const int l) {
/* 寻找中间字符串 */
    int* a;     // l-元组
    int* ma;    // 最小距离对应的l-元组
    int i, d, j;
    int dist, minDist;
    char* cs;    // consensus.
    int* st;     // start point vector
    int* b;     // 各长度的l-元组的最小距离
    int md, lb, ib;
    int ml = 6;     // Minimum length of l-tuple tried.
    if ((a = (int*) malloc(sizeof(int) * l)) == NULL || \
        (ma = (int*) malloc(sizeof(int) * l)) == NULL || \
        (cs = (char*) malloc(sizeof(char) * (l + 1))) == NULL || \
        (b = (int*) malloc(sizeof(int) * (l - 1))) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);
    }
    // 下面一段是为了找出每个长度的L-元组的总距离的下界
    ib = -1;     // 总距离最大的L-元组的长度 - 1。如果小于0，说明所有可能的L-元组其总距离都小于0。
    if (l > 2) {
        md = l - 2 >= ml ? max(ml, (int) sqrt((double) l)) : (int) sqrt((double) l);
        for (i = 0; i < l - 1; i++) { // 获取各个长度l-元组的最小距离
            b[i] = getLowerBound(mulSeq, t, a, b, i + 1);
            if (i >= md && b[i] > 0)
                break;
        }
        if (i < l - 1) 
            ib = i;
    } else if (l == 2) {
        b[0] = getLowerBound(mulSeq, t, a, b, 1);
        if (b[0] > 0)
            ib = 0;
    }
    // 遍历所有可能的l-元组，寻找最小距离
    for (i = 0; i < l; i++)
        ma[i] = a[i] = 0;
    minDist = l * t;
    d = 0;
    while (d >= 0) {
        if (d < l - 1) {   // 如果l == 1的话这一段是不会执行的。
            dist = getDistance(mulSeq, t, a, l, d);
            if (ib >= 0) {     // 如果将word拆分为u+v，此处计算TotalDistance(v, DNA)的下界
                j = (l - 1 - d) % (ib + 1);
                lb = l - 2 - d <= ib ? b[l - 2 - d] : b[ib] * (l - 1 - d) / (ib + 1) + (j == 0 ? 0 : b[j - 1]);
            } else {
                lb = 0;
            }
            if (dist + lb > minDist)
                d = byPass(a, d, l);
            else 
                d = nextVertex(a, d, l);
        } else {
            dist = getDistance(mulSeq, t, a, l, l - 1);
            if (dist < minDist) {  // 注意最小值对应的l-元组可能有多种，这里只保留了其中一种。其实还可以设置如果minDist = 0时终止搜索。
                minDist = dist;
                for (i = 0; i < l; i++)
                    ma[i] = a[i];
            }
            d = nextVertex(a, d, l);
        }
    }
    // 获取共有序列
    for (i = 0; i < l; i++)
        cs[i] = sortedAlpha[ma[i]];
    cs[l] = '\0';
    // 获取最优起始位点向量
    st = getStartPoints(mulSeq, t, ma, l);
    // 输出结果
    printf("Min distance: %d\n", minDist);
    printf("Starting point vector: ");
    for (i = 0; i < t; i++)
        printf("%d ", st[i] + 1);
    printf("\n");
    printf("Target seq:\n");
    for (i = 0; i < t; i++) {
        printf("\t");
        for (j = st[i]; j < st[i] + l; j++)
            printf("%c", mulSeq[i]->s[j]);
        printf("\n");
    }
    printf("Consensus seq: %s\n", cs);
    // 释放内存
    free(a);
    free(ma);
    free(st);
    free(cs);
    free(b);
}