序列比对（23）最长公共子字符串
- 2019 年 10 月 4 日
- 筆記
本文介绍如何求解两个字符串的最长公共子字符串。
其实这个问题可以放在序列比对专题的最开始，只是笔者是个新手，所以当初只是照《生物序列分析》教材的进度写的，教材是直接从全局比对开始讲的。Anyway，我们在本文介绍也不迟。
刚开始接受编程训练时，很容易想到利用三层循环求解。在此就不赘述了。
回溯的时候从得分矩阵的最大值所在单元开始，一直到值为0的单元。
效果如下：
当然，笔者还想过如果是用多层循环的话，可以考虑结合KMP算法。当然，这只是一个想法，没有去实现。
动态规划解法的代码
具体代码如下： （代码是在《序列比对(一)——全局比对Needleman-Wunsch算法》一文代码的基础上修改，没有优化，但足以说明本文问题了。）
#include <stdio.h> #include <stdlib.h> #include <string.h> #define MAXSEQ 1000 void strUpper(char *s); void printAlign(int** a, const int i, const int j, char* s, char* r, int n); void align(char *s, char *r); int main() { char s[MAXSEQ]; char r[MAXSEQ]; printf("The 1st seq: "); scanf("%s", s); printf("The 2nd seq: "); scanf("%s", r); align(s, r); return 0; } void strUpper(char *s) { while (*s != ' ') { if (*s >= 'a' && *s <= 'z') { *s -= 32; } s++; } } void printAlign(int** a, const int i, const int j, char* s, char* r, int n) { int k; int p = a[i - n][j - n]; if (p == 0) { // 到值为0的矩阵单元就结束 printf("start and end position of common seq in seq1: %d %dn", i - n + 1, i); printf("start and end position of common seq in seq2: %d %dn", j - n + 1, j); for (k = 0; k < n; k++) printf("%c", s[i - n + k]); printf("n"); for (k = 0; k < n; k++) printf("%c", r[j - n + k]); printf("nn"); return; } printAlign(a, i, j, s, r, n + 1); } void align(char *s, char *r) { int i, j; int m = strlen(s); int n = strlen(r); int **aUnit; int max; // 初始化 if ((aUnit = (int**) malloc(sizeof(int*) * (m + 1))) == NULL) { fputs("Error: Out of space!n", stderr); exit(1); } for (i = 0; i <= m; i++) { if ((aUnit[i] = (int*) malloc(sizeof(int) * (n + 1))) == NULL) { fputs("Error: Out of space!n", stderr); exit(1); } } for (i = 0; i <= m; i++) aUnit[i][0] = 0; for (j = 1; j <= n; j++) aUnit[0][j] = 0; // 将字符串都变成大写 strUpper(s); strUpper(r); // 动态规划算法计算得分矩阵每个单元的分值 for (i = 1; i <= m; i++) for (j = 1; j <= n; j++) aUnit[i][j] = s[i - 1] == r[j - 1] ? aUnit[i - 1][j - 1] + 1 : 0; /* // 打印得分矩阵 for (i = 0; i <= m; i++) { for (j = 0; j <= n; j++) printf("%d ", aUnit[i][j]); printf("n"); } */ for (i = 1, max = 0; i <= m; i++) for (j = 1; j <= n; j++) if (aUnit[i][j] > max) max = aUnit[i][j]; printf("max score: %dn", max); // 打印最优比对结果，如果有多个，全部打印 // 递归法 if (max == 0) { fputs("No common sub str found.n", stdout); } else { for (i = 1; i <= m; i++) for (j = 1; j <= n; j++) if (aUnit[i][j] == max) printAlign(aUnit, i, j, s, r, 0); } for (i = 0; i <= m; i++) free(aUnit[i]); free(aUnit); }
