专栏首页生信了序列比对(24)最长公共子序列

序列比对(24)最长公共子序列

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

最长公共子序列问题

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

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

效果如下:

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

具体代码如下:

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

本文分享自微信公众号 - 生信了(gh_ed36a29a9a9d),作者:hxj7

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-09-11

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 序列比对(九)从掷骰子说起HMM

    前面八篇文章介绍了动态规划在序列比对中的基础应用。从本文开始,开始介绍HMM(隐马尔可夫模型)。为什么要介绍它呢?因为该模型在发现Motif、预测CpG岛等多方...

    一只羊
  • 序列比对(14)viterbi算法和后验解码的比较

    前文《序列比对(十)viterbi算法求解最可能路径》介绍了用viterbi算法求解最可能路径:在符号序列已知而状态序列未知时,最可能路径是:

    一只羊
  • 序列比对(22)中间字符串分支定界方法中更紧的界

    前文《序列比对(21)中间字符串问题的算法及实现代码》介绍了中间字符串的算法和代码,但是使用分支定界策略时所使用的界限是很宽松的。分支定界法的伪代码如下:

    一只羊
  • Spring 复盘 | AOP

    Spring AOP 基础 Java 动态代理实现,阅读文章之前,你最好有以下基础:

    一个优秀的废人
  • 【Basic algorithm学习笔记】 排序

    原理: 其实原理很简单 大白话来讲 就是 先定义很多 容器 即这里形象的表达 为桶,之后每次输入的数字则为选择哪个桶,将对应的 桶中扔一个标志物表示有这个 桶...

    Python进击者
  • Fence Repair (PKU 3253)

    思路:用朴素的方法实现的话,时间复杂度为O(n^n)。因为只需要从板的集合中取出最短的两块,并且把长度为两块长度之和的板加入集合中即可,所有使用优先队列就可以高...

    用户7727433
  • Google APAC 2015 University Graduates Test Round D

    Google APAC 2015 University Graduates Test Round D

    宅男潇涧
  • 聊聊artemis的NetworkHealthCheck

    activemq-artemis-2.11.0/artemis-commons/src/main/java/org/apache/activemq/artemi...

    codecraft
  • LeetCode 1390. 四因数

    来源:力扣(LeetCode) 链接:https://leetcode-cn.com/problems/four-divisors 著作权归领扣网络所有。商...

    Michael阿明
  • 聊聊artemis的NetworkHealthCheck

    activemq-artemis-2.11.0/artemis-commons/src/main/java/org/apache/activemq/artemi...

    codecraft

扫码关注云+社区

领取腾讯云代金券