序列比对（24）最长公共子序列

最新推荐文章于 2021-12-26 13:45:36 发布

生信了（公众号同名）

最新推荐文章于 2021-12-26 13:45:36 发布

阅读量434

点赞数

分类专栏： # 序列算法文章标签：序列比对公共子序列

本文链接：https://blog.csdn.net/biocity/article/details/100736953

版权

序列算法专栏收录该内容

32 篇文章 16 订阅

订阅专栏

本文介绍如何求解两个字符串的最长公共子序列。

最长公共子序列问题

前文《序列比对（23）最长公共子字符串》介绍了如何求解两个字符串的最长公共子字符串，本文将介绍如何求解两个字符串的最长公共子序列。二者听起来很像，所以我们首先得说明一下子字符串和子序列的区别。

在本文语境中，子字符串是由原字符串中的连续字符组成；而子序列是从原字符串中选择字符，只需要保持其次序不变（可以说，子字符串都是子序列，但子序列不一定是子字符串）。比如：对于 $\bm{v}=AGCT$ 这个字符串， $A G$ 既是其子字符串也是子序列； $A C$ 是子序列而非子字符串； $C A$ 既非子字符串也非子序列。

给定两个字符串 $\bm{v}$ 和 $\bm{w}$ ，长度分别为 $m$ 和 $n$ ，如何找到这两个字符串的最长公共子序列呢？举例来说：如果 $\bm{v}=AGCT$ ，而 $\bm{w}=TAC$ ，那么二者的最长公共子序列就是 $A C$ 。

与最长公共子字符串问题类似，最长公共子序列问题也是一种序列比对问题，可以用动态规划解决，只是在迭代时允许插入和缺失，而不允许错配而已。如果是匹配，得分为1，否则得分为0。其迭代公式如下：

$\begin{aligned} & \text{$F(i,j)$ is the maximum score of alignments between $x_{1 \ldots i}$ and $y_{1 \ldots j}$.} \\ & F(i, 0) = 0 \quad \text{for $i = 0 \ldots m$.} \\ & F(0, j) = 0 \quad \text{for $j = 1 \ldots n$.} \\ & F(i, j) = \max \begin{cases} F(i - 1, j), \\ F(i, j - 1), \\ F(i - 1, j - 1) + 1 \quad \text{if $x_i = y_j$.} \end{cases} \end{aligned}$

回溯的时候从得分矩阵的右下角开始，一直到值为0的单元。回溯过程中只记录 $x_i = y_j$ 的单元对应的字符。

效果如下：
在这里插入图片描述
（公众号：生信了）

动态规划求解最长公共子序列的代码

具体代码如下：

#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);
}