序列比对（18）重复匹配问题的补充说明
- 2019 年 10 月 7 日
- 筆記
前文介绍了重复匹配问题的动态规划算法，但是遗留了重复结果输出的问题。本文对该问题进行了补充说明。
前文《序列匹配(五)——重复匹配问题的动态规划算法》介绍了重复匹配问题的动态规划算法。
但是这个公式在回溯时会出现重复结果输出的问题，比如：
校正公式和代码
这样的公式目前还没有出现重复结果输出的问题：
相应的代码放在了文末。
对比对总长度的估计
实现代码
#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 != ' ') { 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=%dn", 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("nn"); 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=%dn"); 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("nn"); 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: %fn", 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); }
