前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >序列比对(23)最长公共子字符串

序列比对(23)最长公共子字符串

作者头像
一只羊
发布2019-09-17 14:07:51
5600
发布2019-09-17 14:07:51
举报
文章被收录于专栏:生信了

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

其实这个问题可以放在序列比对专题的最开始,只是笔者是个新手,所以当初只是照《生物序列分析》教材的进度写的,教材是直接从全局比对开始讲的。Anyway,我们在本文介绍也不迟。

刚开始接受编程训练时,很容易想到利用三层循环求解。在此就不赘述了。

回溯的时候从得分矩阵的最大值所在单元开始,一直到值为0的单元。

效果如下:

当然,笔者还想过如果是用多层循环的话,可以考虑结合KMP算法。当然,这只是一个想法,没有去实现。

点击此处,等你留言

动态规划解法的代码

具体代码如下: (代码是在《序列比对(一)——全局比对Needleman-Wunsch算法》一文代码的基础上修改,没有优化,但足以说明本文问题了。)

代码语言:javascript
复制
#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 != '\0') {
        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 %d\n", i - n + 1, i);
        printf("start and end position of common seq in seq2: %d %d\n", 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("\n\n");
        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: %d\n", 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);
}
本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2019-09-11,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信了 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 动态规划解法的代码
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档