原创：hxj7
本文介绍如何求解两个字符串的最长公共子序列。
最长公共子序列问题
前文《序列比对（23）最长公共子字符串》介绍了如何求解两个字符串的最长公共子字符串，本文将介绍如何求解两个字符串的最长公共子序列。二者听起来很像，所以我们首先得说明一下子字符串和子序列的区别。
在本文语境中，子字符串是由原字符串中的连续字符组成；而子序列是从原字符串中选择字符，只需要保持其次序不变（可以说，子字符串都是子序列，但子序列不一定是子字符串）。比如：对于这个字符串，
既是其子字符串也是子序列；
是子序列而非子字符串；
既非子字符串也非子序列。
给定两个字符串和
，长度分别为
和
，如何找到这两个字符串的最长公共子序列呢？举例来说：如果
，而
，那么二者的最长公共子序列就是
。
与最长公共子字符串问题类似，最长公共子序列问题也是一种序列比对问题，可以用动态规划解决，只是在迭代时允许插入和缺失，而不允许错配而已。如果是匹配，得分为1，否则得分为0。其迭代公式如下：
回溯的时候从得分矩阵的右下角开始，一直到值为0的单元。回溯过程中只记录的单元对应的字符。
效果如下：
image
动态规划求解最长公共子序列的代码
具体代码如下：
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MAXSEQ 1000
#define max(a, b) ((a) > (b) ? (a) : (b))
#define min(a, b) ((a) < (b) ? (a) : (b))
struct Unit {
int W1; // 是否往上回溯一格
int W2; // 是否往左上回溯一格
int W3; // 是否往左回溯一格
int M; // 得分矩阵第(i, j)这个单元的分值，即序列s(1,...,i)与序列r(1,...,j)比对的最高得分
};
typedef struct Unit *pUnit;
void strUpper(char *s);
void printAlign(pUnit** a, const int i, const int j, char* s, char* r, int* is, int* ir, 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 != '\0') {
if (*s >= 'a' && *s <= 'z') {
*s -= 32;
}
s++;
}
}
void printAlign(pUnit** a, const int i, const int j, char* s, char* r, int* is, int* ir, int n) {
int k;
pUnit p = a[i][j];
if (p->M == 0) { // 到值为0的矩阵单元就结束
printf("index of common seq in seq1: ");
for (k = n - 1; k >= 0; k--)
printf("%d ", is[k]);
printf("\n");
printf("index of common seq in seq2: ");
for (k = n - 1; k >= 0; k--)
printf("%d ", ir[k]);
printf("\n");
for (k = n - 1; k >= 0; k--)
printf("%c", s[is[k] - 1]);
printf("\n");
for (k = n - 1; k >= 0; k--)
printf("%c", r[ir[k] - 1]);
printf("\n\n");
return;
}
if (p->W1) // 向上回溯一格
printAlign(a, i - 1, j, s, r, is, ir, n);
if (p->W2) { // 向左上回溯一格
is[n] = i;
ir[n] = j;
printAlign(a, i - 1, j - 1, s, r, is, ir, n + 1);
}
if (p->W3) // 向左回溯一格
printAlign(a, i, j - 1, s, r, is, ir, n);
}
void align(char *s, char *r) {
int i, j;
int m = strlen(s);
int n = strlen(r);
int ml = min(m, n);
int m1, m2, m3, maxm;
pUnit **aUnit;
int* sidx;
int* ridx;
// 初始化
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]->W1 = 0;
aUnit[i][j]->W2 = 0;
aUnit[i][j]->W3 = 0;
}
}
for (i = 0; i <= m; i++)
aUnit[i][0]->M = 0;
for (j = 1; j <= n; j++)
aUnit[0][j]->M = 0;
// 将字符串都变成大写
strUpper(s);
strUpper(r);
// 动态规划算法计算得分矩阵每个单元的分值
for (i = 1; i <= m; i++) {
for (j = 1; j <= n; j++) {
m1 = aUnit[i - 1][j]->M;
m2 = s[i - 1] == r[j - 1] ? aUnit[i - 1][j - 1]->M + 1 : -1;
m3 = aUnit[i][j - 1]->M;
maxm = max(max(m1, m2), m3);
aUnit[i][j]->M = maxm;
if (m1 == maxm) aUnit[i][j]->W1 = 1;
if (m2 == maxm) aUnit[i][j]->W2 = 1;
if (m3 == maxm) aUnit[i][j]->W3 = 1;
}
}
/*
// 打印得分矩阵
for (i = 0; i <= m; i++) {
for (j = 0; j <= n; j++)
printf("%d ", aUnit[i][j]->M);
printf("\n");
}
*/
printf("max score: %d\n", aUnit[m][n]->M);
// 打印最优比对结果，如果有多个，全部打印
// 递归法
if (aUnit[m][n]->M == 0) {
fputs("No common seq found.\n", stdout);
} else {
if ((sidx = (int*) malloc(sizeof(int) * ml)) == NULL || \
(ridx = (int*) malloc(sizeof(int) * ml)) == NULL) {
fputs("Error: Out of space!\n", stderr);
exit(1);
}
printAlign(aUnit, m, n, s, r, sidx, ridx, 0);
// 释放内存
free(sidx);
free(ridx);
}
for (i = 0; i <= m; i++) {
for (j = 0; j <= n; j++)
free(aUnit[i][j]);
free(aUnit[i]);
}
free(aUnit);
}
（公众号：生信了）
网友评论