序列匹配（五）重复匹配问题的动态规划算法
前言： 蛋白质序列中常有重复的功能域（domain）或模体（motif）拷贝，由此衍生出一个抽象的序列多重匹配的问题，即如何从一个序列中找出另一个序列的某部分（如功能域或模体）的多个无交叠（non-overlapping)拷贝。本文给出了该问题的示例、关键计算公式以及C语言实现代码。
问题及算法描述 更具体地描述上面的问题：有序列x和y，其中y是包含结构域的序列，x是要从中找到多重匹配的序列。将x分割成一段一段的不交叠的子序列，这些子序列要么不参与和y的联配，要么与y的某一段子序列联配且联配的分值不低于一个阈值T。如果将x的某一子序列的联配分值减去T作为其“标准联配分值”，那么最终目标是找到这些参与联配的x的子序列的“标准联配分值”之和的最大值。
在《Biological sequence analysis》这本书中列举了一个蛋白质的例子：
引自《生物序列分析》
上图中显示，在最优联配（即“标准联配分值”之和最大的联配）中，x有两个子序列参与了联配，“标准联配分值”分别是1和8。
那么上图中的“标准联配分值”是如何计算得到的呢？这依然是利用动态规划算法，在《Biological sequence analysis》书中给出了关键的计算公式：
引自《生物序列分析》
其中，F(i, j)假设x(i)参与联配，且对应的联配结束于x(i)与y(j)时的“标准联配分值”之和的最大值；而F(i, 0)指的是子序列x(1,2,…,i)与y（联配可结束于y的任意位置）的“标准联配分值”之和的最大值，假设x(i)不参与联配。
一个困惑 上面计算公式的C代码实现见下文，简称其为alnRepeat，以便和本文另一段代码区分。其中一个示例如下：
没有问题，但是另一个示例的结果让我困惑：
最优分值应该是6啊。回过头再去看上面的计算公式，如果我对书中相应章节的理解无误的话，它似乎没有考虑到x两个子序列都参与联配且这两个子序列紧挨在一起的情况，最简单的就是上面AA的这个例子。理论上，最优联配中，两个连续的A应该都参与了联配，且属于两个不同的“匹配段”。
算法的补充 由此，我重新思考分值的计算公式。F(i, 0)没有问题，而F(i, j)的计算可以分为下面三种情况（注意，F(i, j)的前提条件假设了x(i)“参与”了联配）：
- x(i-1)没有参与联配；
- x(i-1)参与了联配，且与x(i)属于同一个“匹配段”；
- x(i-1)参与了联配，且与x(i)属于不同的“匹配段”。
考虑了上面三种情况，F(i, j)的计算公式变成：
对上述新公式编写代码的过程中，发现输出结果有很多重复，再次检视公式，发现F(i, 1)需要单独处理，否则F(i, 1)计算公式中的F(i, 0)一项会引发歧义。于是，最终的计算公式变成：
其C语言实现代码简称alnRepeat3。 运行alnRepeat以及alnRepeat3比较二者的不同：
alnRepeat3的结果仍有重复，说明代码还要优化。
小结 本文介绍了生物序列重复匹配的问题以及相应的动态规划算法，在代码实现过程中，发现了疑似错误的示例（原计算公式似乎没有考虑到两个“匹配段”紧挨在一起的情况）并补充了计算公式。
对笔者来说，最大的收获在于学习这个新算法的过程中认真地进行过思考，但由于自身能力以及时间精力所限，对这个问题的理解以及代码实现还有很多不足，真切期望能有热心同道能够给出意见和建议！ 最后，谢谢大家看此长文！
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);
}本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2019-04-30，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读