序列比对（21）中间字符串问题的算法及实现代码

最新推荐文章于 2021-05-23 10:36:31 发布
最新推荐文章于 2021-05-23 10:36:31 发布
阅读量488 收藏 3
点赞数
分类专栏： # 序列算法 文章标签： 多序列比对 中间字符串
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/biocity/article/details/100126851
版权

前文介绍了基序发现问题和中间字符串问题。本文给出了中间字符串的算法和实现代码。

中间字符串问题的简单算法及伪代码

前文《序列比对（19）基序发现和中间字符串问题》介绍了基序发现问题和中间字符串问题；《序列比对（20）基序发现问题的算法及实现代码》给出了基序问题的算法和实现代码。本文将介绍中间字符串问题的算法，并给出实现代码。

简单回顾一下，中间字符串问题其实就是要找到使得总距离最小的一个 l l l-元组。即
a r g m i n v   min ⁡ s ⃗   d H ( v , s ⃗ ) \underset{v}{\mathrm{argmin}} \, \underset{\boldsymbol{\vec{s}}}{\min} \, d_H(v,\boldsymbol{\vec{s}}) vargmin​s min​dH​(v,s )

或者：
a r g m i n v   T o t a l D i s t a n c e ( v , D N A ) \underset{v}{\mathrm{argmin}} \, TotalDistance(v,DNA) vargmin​TotalDistance(v,DNA)

由于要遍历所有可能的起始位点，如前文《序列比对（20）基序发现问题的算法及实现代码》一样，我们采用树结构以及DFS（深度优先搜索）。其简单算法如下（伪代码）：
在这里插入图片描述
在这里插入图片描述
图1到图2引自《生物信息学算法导论》

从上面的代码可以看出，算法构建了一个深度为 l l l的搜索树，除了叶顶点外，每个顶点有 4 4 4个子顶点（表示 l l l-元组的每个位置都有可能有4种碱基）。当遍历到叶顶点时，计算总距离 T o t a l D i s t a n c e ( v , D N A ) TotalDistance(v,DNA) TotalDistance(v,DNA)，如果该总距离比目前的最小总距离还小，那么就将该总距离记为新的最小总距离，并且将此时的 l l l-元组记录下来。

分支定界策略改进简单算法

分支定界法已经在前文《序列比对（20）基序发现问题的算法及实现代码》中介绍过。一个深度为 i i i的顶点，如果其距离 T o t a l D i s t a n c e ( v , i , D N A ) TotalDistance(v,i,DNA) TotalDistance(v,i,DNA)比目前的最小总距离还大，那么以这个顶点为根的分支树就没有继续遍历的意义了。

用分支定界策略改进上述简单算法的伪代码如下：
在这里插入图片描述
在这里插入图片描述
图3到图4引自《生物信息学算法导论》

注意，上述顶点的分支界限是很宽松的，后续文章我们会分析更紧的界。

实验效果

我们用《生物信息学算法导论》书中的两组序列进行实验，分别用简单算法（medianStr_BF.exe)和分支定界法（medianStr_BB.exe）去计算中间字符串。
第一组序列（存储在identity.txt文件中，下面示例结果中用到）：

CGGGGCTATGCAACTGGGTCGTCACATTCCCCTTTCGATA
TTTGAGGGTGCCCAATAAATGCAACTCCAAAGCGGACAAA
GGATGCAACTGATGCCGTTTGACGACCTAAATCAACGGCC
AAGGATGCAACTCCAGGAGCGCCTTTGCTGGTTCTACCTG
AATTTTCTAAAAAGATTATAATGTCGGTCCATGCAACTTC
CTGCTGTACAACTGAGATCATGCTGCATGCAACTTTCAAC
TACATGATCTTTTGATGCAACTTGGATGAGGGAATGATGC

第二组序列（存储在mutated.txt文件中，下面示例结果中用到）：

CGGGGCTATcCAgCTGGGTCGTCACATTCCCCTTTCGATA
TTTGAGGGTGCCCAATAAggGCAACTCCAAAGCGGACAAA
GGATGgAtCTGATGCCGTTTGACGACCTAAATCAACGGCC
AAGGAaGCAACcCCAGGAGCGCCTTTGCTGGTTCTACCTG
AATTTTCTAAAAAGATTATAATGTCGGTCCtTGgAACTTC
CTGCTGTACAACTGAGATCATGCTGCATGCcAtTTTCAAC
TACATGATCTTTTGATGgcACTTGGATGAGGGAATGATGC

令人惊喜的是简单算法的效率也非常高：
（只要对基序发现问题和中间字符串问题的简单算法的运行时间做简单分析）
在这里插入图片描述
为identity.txt文件中的7条序列计算中间字符串

在这里插入图片描述
为mutated.txt文件中的7条序列计算中间字符串

分支定界法的结果如下：
在这里插入图片描述
为identity.txt文件中的7条序列计算中间字符串

在这里插入图片描述
为mutated.txt文件中的7条序列计算中间字符串

具体代码

上文及前文都假定多条序列的长度是一样的，但是实际情况并不总是如此。代码实现过程中考虑到这一点，做了改进，使得多条序列长度不一致的情况下也可以用此代码来计算中间字符串。

（公众号：生信了）

首先是简单算法

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#define MAXSEQ 1000

struct Sequence {
    char* s;
    int* a;      // 该向量存储每个字符在字符集中的序号
    int l;
};
typedef struct Sequence* Seq;

const int nalpha = 4;     // 字符种数
const char sortedAlpha[] = {'A', 'C', 'G', 'T'};    // 排过序的字符集

char* sec2time(time_t t);
void strUpper(char *s);
int search(char c);                      /* 在排序后的字符集中寻找字符c，如果找到返回序号；找不到返回-1 */
Seq* readSeq(char* filename, const int t);          /* 从文件中读取多条序列 */
Seq create(char* s, const int n, const int i);  
void destroy(Seq seq);
int nextVertex(int* a, int d, const int l);     /* 得到下一个顶点 */
int getDistance(Seq* mulSeq, const int t, int* a, const int l);        /* 计算一个l-元组和DNA的最小总距离 */
int* getStartPoints(Seq* mulSeq, const int t, int* a, const int l);    /* 获取起始位点向量 */
void findMedianStr(Seq* mulSeq, const int t, const int l);           /* 寻找中间字符串 */

int main(void) {
    time_t start, end;
    char* tstr;
    int i, t, l;
    char filename[30];
    Seq* mulSeq;
    start = time(NULL);  // 开始时间
    printf("t: ");
    scanf(" %d", &t);
    printf("l: ");
    scanf(" %d", &l);
    printf("seq_file path: ");
    scanf(" %s", filename);
    mulSeq = readSeq(filename, t);
    findMedianStr(mulSeq, t, l);
    for (i = 0; i < t; i++)
        destroy(mulSeq[i]);
    free(mulSeq);
    end = time(NULL);   // 结束时间
    tstr = sec2time(end - start);
    printf("time spent: %s\n", tstr);
    free(tstr);
    return 0;
}

char* sec2time(time_t t) {
/* 将秒数转化为时间字符串 */
    int h, m, s;
    char* tstr;
    if ((tstr = malloc(sizeof(char) * 50)) == NULL) {
        fputs("Error: out of space!\n", stderr);
        exit(8);
    }
    h = t / 3600;
    m = (t % 3600) / 60;
    s = t % 60;
    if (sprintf(tstr, "%02d:%02d:%02d", h, m, s) > 0)
        return tstr;
    else {
        sprintf(tstr, "more than %d hours", h);
        return tstr;
    }
}

void strUpper(char *s) {
    while (*s != '\0') {
        if (*s >= 'a' && *s <= 'z')
            *s -= 32;
        s++;
    }
}

int search(char c) {
/* 在排序后的字符集中寻找字符c，如果找到返回序号；找不到返回-1 */
    int i;
    int th = 5;
    int h, l, m;
    if (nalpha < th) {
        for (i = 0; i < nalpha; i++)
            if (sortedAlpha[i] == c)
                return i;
        return -1;
    } else {
        l = 0;
        h = nalpha - 1;
        while (l <= h) {
            m = (l + h) / 2;
            if (sortedAlpha[m] == c)
                return m;
            else if (sortedAlpha[m] < c)
                l = m + 1;
            else 
                h = m - 1;
        }
        return -1;
    }
}

Seq* readSeq(char* filename, const int t) {
/* 从文件中读取多条序列 */
    int i, n;
    char buf[MAXSEQ];
    FILE* fp;
    Seq* mulSeq;
    if ((fp = fopen(filename, "r")) == NULL) {
        fprintf(stderr, "Error: cannot open '%s'\n", filename);
        exit(1);
    }
    if ((mulSeq = (Seq*) malloc(sizeof(Seq) * t)) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);        
    }
    for (i = 0; i < t; i++) {
        if (fgets(buf, MAXSEQ, fp) == NULL) {
            fprintf(stderr, "Error: reading seq from '%s'\n", filename);
            exit(2);
        }
        n = strlen(buf);
        if (buf[n - 1] != '\n' && ! feof(fp)) {
            fprintf(stderr, "Error: line %d is too long.\n", i + 1);
            exit(3);
        }
        if (! feof(fp)) {
            buf[n - 1] = '\0';
            mulSeq[i] = create(buf, n - 1, i);
        } else {
            mulSeq[i] = create(buf, n, i);
        }
    }
    fclose(fp);
    return mulSeq;
}

Seq create(char* s, const int n, const int i) {
/* n: s的长度，不包括末尾的\0 */
    Seq seq;
    int j;
    if ((seq = (Seq) malloc(sizeof(struct Sequence))) == NULL || \
        (seq->s = (char*) malloc(sizeof(char) * (n + 1))) == NULL || \
        (seq->a = (int*) malloc(sizeof(int) * n)) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);
    }
    strcpy(seq->s, s);
    // 将序列中所有字符变为大写
    strUpper(seq->s);
    // 找到每个字符在字符集中的序号。如果是在后期做这个工作，那么就有很多重复，降低效率。
    for (j = 0; j < n; j++)
        if ((seq->a[j] = search(seq->s[j])) < 0) {
            fprintf(stderr, "Error: char %d of seq %d is %c, which is out of range.\n", j + 1, i + 1, seq->s[j]);
            exit(16);
        }
    seq->l = n;
    return seq;
}

void destroy(Seq seq) {
    free(seq->a);
    free(seq->s);
    free(seq);
}

int nextVertex(int* a, int d, const int l) {
/* 得到下一个顶点 */
    int i;
    if (d < l - 1) {
        a[d + 1] = 0;
        return d + 1;
    }
    for (i = l - 1; i >= 0; i--) {
        if (a[i] < nalpha - 1) {
            a[i]++;
            return i;
        }
    }
    return -1;
}

int getDistance(Seq* mulSeq, const int t, int* a, const int l) {
/* 计算一个l-元组和DNA的最小总距离 */
    int i, j, k;
    int dist, minDist, totalDist;
    for (i = 0, totalDist = 0; i < t; i++) {
        for (j = 0, minDist = l; j <= mulSeq[i]->l - l; j++) {
            for (k = 0, dist = 0; k < l; k++)
                dist += a[k] == mulSeq[i]->a[j + k] ? 0 : 1;
            if (dist < minDist)
                minDist = dist;
        }
        totalDist += minDist;
    }
    return totalDist;
}       

int* getStartPoints(Seq* mulSeq, const int t, int* a, const int l) {
/* 获取起始位点向量 */
    int* st;
    int i, j, k;
    int dist, minDist;
    if ((st = (int*) malloc(sizeof(int) * t)) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);
    }
    for (i = 0; i < t; i++) {
        for (j = 0, st[i] = 0, minDist = l; j <= mulSeq[i]->l - l; j++) {
            for (k = 0, dist = 0; k < l; k++)
                dist += a[k] == mulSeq[i]->a[j + k] ? 0 : 1;
            if (dist < minDist) {
                st[i] = j;
                minDist = dist;
            }
        }
    }
    return st;
}

void findMedianStr(Seq* mulSeq, const int t, const int l) {
/* 寻找中间字符串 */
    int* a;     // l-元组
    int* ma;    // 最小距离对应的l-元组
    int i, d, j;
    int dist, minDist;
    char* cs;    // consensus.
    int* st;     // start point vector
    if ((a = (int*) malloc(sizeof(int) * l)) == NULL || \
        (ma = (int*) malloc(sizeof(int) * l)) == NULL || \
        (cs = (char*) malloc(sizeof(char) * (l + 1))) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);
    }
    // 遍历所有可能的l-元组，寻找最小距离
    for (i = 0; i < l; i++)
        ma[i] = a[i] = 0;
    minDist = l * t;
    d = 0;
    while (d >= 0) {
        if (d < l - 1)
            d = nextVertex(a, d, l);
        else {
            dist = getDistance(mulSeq, t, a, l);
            if (dist < minDist) {  // 注意最小值对应的l-元组可能有多种，这里只保留了其中一种。其实还可以设置当minDist = 0时终止搜索。
                minDist = dist;
                for (i = 0; i < l; i++)
                    ma[i] = a[i];
            }
            d = nextVertex(a, d, l);
        }
    }
    // 获取共有序列
    for (i = 0; i < l; i++)
        cs[i] = sortedAlpha[ma[i]];
    cs[l] = '\0';
    // 获取最优起始位点向量
    st = getStartPoints(mulSeq, t, ma, l);
    // 输出结果
    printf("Min distance: %d\n", minDist);
    printf("Starting point vector: ");
    for (i = 0; i < t; i++)
        printf("%d ", st[i] + 1);
    printf("\n");
    printf("Target seq:\n");
    for (i = 0; i < t; i++) {
        printf("\t");
        for (j = st[i]; j < st[i] + l; j++)
            printf("%c", mulSeq[i]->s[j]);
        printf("\n");
    }
    printf("Consensus seq: %s\n", cs);
    // 释放内存
    free(a);
    free(ma);
    free(st);
    free(cs);
}

然后是分支定界法

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#define MAXSEQ 1000

struct Sequence {
    char* s;
    int* a;      // 该向量存储每个字符在字符集中的序号
    int l;
};
typedef struct Sequence* Seq;

const int nalpha = 4;     // 字符种数
const char sortedAlpha[] = {'A', 'C', 'G', 'T'};    // 排过序的字符集

char* sec2time(time_t t);
void strUpper(char *s);
int search(char c);                      /* 在排序后的字符集中寻找字符c，如果找到返回序号；找不到返回-1 */
Seq* readSeq(char* filename, const int t);          /* 从文件中读取多条序列 */
Seq create(char* s, const int n, const int i);  
void destroy(Seq seq);
int nextVertex(int* a, int d, const int l);     /* 得到下一个顶点 */
int byPass(int* a, int d, const int l);           /* 跳过某个分支，得到下一个顶点 */
int getDistance(Seq* mulSeq, const int t, int* a, const int l, const int d);        /* 计算一个l-元组和DNA的最小总距离 */
int* getStartPoints(Seq* mulSeq, const int t, int* a, const int l);    /* 获取起始位点向量 */
void findMedianStr(Seq* mulSeq, const int t, const int l);           /* 寻找中间字符串 */

int main(void) {
    time_t start, end;
    char* tstr;
    int i, t, l;
    char filename[30];
    Seq* mulSeq;
    start = time(NULL);  // 开始时间
    printf("t: ");
    scanf(" %d", &t);
    printf("l: ");
    scanf(" %d", &l);
    printf("seq_file path: ");
    scanf(" %s", filename);
    mulSeq = readSeq(filename, t);
    findMedianStr(mulSeq, t, l);
    for (i = 0; i < t; i++)
        destroy(mulSeq[i]);
    free(mulSeq);
    end = time(NULL);   // 结束时间
    tstr = sec2time(end - start);
    printf("time spent: %s\n", tstr);
    free(tstr);
    return 0;
}

char* sec2time(time_t t) {
/* 将秒数转化为时间字符串 */
    int h, m, s;
    char* tstr;
    if ((tstr = malloc(sizeof(char) * 50)) == NULL) {
        fputs("Error: out of space!\n", stderr);
        exit(8);
    }
    h = t / 3600;
    m = (t % 3600) / 60;
    s = t % 60;
    if (sprintf(tstr, "%02d:%02d:%02d", h, m, s) > 0)
        return tstr;
    else {
        sprintf(tstr, "more than %d hours", h);
        return tstr;
    }
}

void strUpper(char *s) {
    while (*s != '\0') {
        if (*s >= 'a' && *s <= 'z')
            *s -= 32;
        s++;
    }
}

int search(char c) {
/* 在排序后的字符集中寻找字符c，如果找到返回序号；找不到返回-1 */
    int i;
    int th = 5;
    int h, l, m;
    if (nalpha < th) {
        for (i = 0; i < nalpha; i++)
            if (sortedAlpha[i] == c)
                return i;
        return -1;
    } else {
        l = 0;
        h = nalpha - 1;
        while (l <= h) {
            m = (l + h) / 2;
            if (sortedAlpha[m] == c)
                return m;
            else if (sortedAlpha[m] < c)
                l = m + 1;
            else 
                h = m - 1;
        }
        return -1;
    }
}

Seq* readSeq(char* filename, const int t) {
/* 从文件中读取多条序列 */
    int i, n;
    char buf[MAXSEQ];
    FILE* fp;
    Seq* mulSeq;
    if ((fp = fopen(filename, "r")) == NULL) {
        fprintf(stderr, "Error: cannot open '%s'\n", filename);
        exit(1);
    }
    if ((mulSeq = (Seq*) malloc(sizeof(Seq) * t)) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);        
    }
    for (i = 0; i < t; i++) {
        if (fgets(buf, MAXSEQ, fp) == NULL) {
            fprintf(stderr, "Error: reading seq from '%s'\n", filename);
            exit(2);
        }
        n = strlen(buf);
        if (buf[n - 1] != '\n' && ! feof(fp)) {
            fprintf(stderr, "Error: line %d is too long.\n", i + 1);
            exit(3);
        }
        if (! feof(fp)) {
            buf[n - 1] = '\0';
            mulSeq[i] = create(buf, n - 1, i);
        } else {
            mulSeq[i] = create(buf, n, i);
        }
    }
    fclose(fp);
    return mulSeq;
}

Seq create(char* s, const int n, const int i) {
/* n: s的长度，不包括末尾的\0 */
    Seq seq;
    int j;
    if ((seq = (Seq) malloc(sizeof(struct Sequence))) == NULL || \
        (seq->s = (char*) malloc(sizeof(char) * (n + 1))) == NULL || \
        (seq->a = (int*) malloc(sizeof(int) * n)) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);
    }
    strcpy(seq->s, s);
    // 将序列中所有字符变为大写
    strUpper(seq->s);
    // 找到每个字符在字符集中的序号。如果是在后期做这个工作，那么就有很多重复，降低效率。
    for (j = 0; j < n; j++)
        if ((seq->a[j] = search(seq->s[j])) < 0) {
            fprintf(stderr, "Error: char %d of seq %d is %c, which is out of range.\n", j + 1, i + 1, seq->s[j]);
            exit(16);
        }
    seq->l = n;
    return seq;
}

void destroy(Seq seq) {
    free(seq->a);
    free(seq->s);
    free(seq);
}

int nextVertex(int* a, int d, const int l) {
/* 得到下一个顶点 */
    int i;
    if (d < l - 1) {
        a[d + 1] = 0;
        return d + 1;
    }
    for (i = l - 1; i >= 0; i--) {
        if (a[i] < nalpha - 1) {
            a[i]++;
            return i;
        }
    }
    return -1;
}

int byPass(int* a, int d, const int l) {
/* 跳过某个分支，得到下一个顶点 */
    int i;
    for (i = d; i >= 0; i--) {
        if (a[i] < nalpha - 1) {
            a[i]++;
            return i;
        }
    }
    return -1;
}           


int getDistance(Seq* mulSeq, const int t, int* a, const int l, const int d) {
/* 计算一个l-元组和DNA的最小总距离 */
    int i, j, k;
    int dist, minDist, totalDist;
    for (i = 0, totalDist = 0; i < t; i++) {
        for (j = 0, minDist = l; j <= mulSeq[i]->l - l; j++) {
            for (k = 0, dist = 0; k <= d; k++)
                dist += a[k] == mulSeq[i]->a[j + k] ? 0 : 1;
            if (dist < minDist)
                minDist = dist;
        }
        totalDist += minDist;
    }
    return totalDist;
}       

int* getStartPoints(Seq* mulSeq, const int t, int* a, const int l) {
/* 获取起始位点向量 */
    int* st;
    int i, j, k;
    int dist, minDist;
    if ((st = (int*) malloc(sizeof(int) * t)) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);
    }
    for (i = 0; i < t; i++) {
        for (j = 0, st[i] = 0, minDist = l; j <= mulSeq[i]->l - l; j++) {
            for (k = 0, dist = 0; k < l; k++)
                dist += a[k] == mulSeq[i]->a[j + k] ? 0 : 1;
            if (dist < minDist) {
                st[i] = j;
                minDist = dist;
            }
        }
    }
    return st;
}

void findMedianStr(Seq* mulSeq, const int t, const int l) {
/* 寻找中间字符串 */
    int* a;     // l-元组
    int* ma;    // 最小距离对应的l-元组
    int i, d, j;
    int dist, minDist;
    char* cs;    // consensus.
    int* st;     // start point vector
    if ((a = (int*) malloc(sizeof(int) * l)) == NULL || \
        (ma = (int*) malloc(sizeof(int) * l)) == NULL || \
        (cs = (char*) malloc(sizeof(char) * (l + 1))) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);
    }
    // 遍历所有可能的l-元组，寻找最小距离
    for (i = 0; i < l; i++)
        ma[i] = a[i] = 0;
    minDist = l * t;
    d = 0;
    while (d >= 0) {
        if (d < l - 1) {
            dist = getDistance(mulSeq, t, a, l, d);
            if (dist > minDist)
                d = byPass(a, d, l);
            else
                d = nextVertex(a, d, l);
        } else {
            dist = getDistance(mulSeq, t, a, l, l - 1);
            if (dist < minDist) {  // 注意最小值对应的l-元组可能有多种，这里只保留了其中一种。其实还可以设置当minDist = 0时终止搜索。
                minDist = dist;
                for (i = 0; i < l; i++)
                    ma[i] = a[i];
            }
            d = nextVertex(a, d, l);
        }
    }
    // 获取共有序列
    for (i = 0; i < l; i++)
        cs[i] = sortedAlpha[ma[i]];
    cs[l] = '\0';
    // 获取最优起始位点向量
    st = getStartPoints(mulSeq, t, ma, l);
    // 输出结果
    printf("Min distance: %d\n", minDist);
    printf("Starting point vector: ");
    for (i = 0; i < t; i++)
        printf("%d ", st[i] + 1);
    printf("\n");
    printf("Target seq:\n");
    for (i = 0; i < t; i++) {
        printf("\t");
        for (j = st[i]; j < st[i] + l; j++)
            printf("%c", mulSeq[i]->s[j]);
        printf("\n");
    }
    printf("Consensus seq: %s\n", cs);
    // 释放内存
    free(a);
    free(ma);
    free(st);
    free(cs);
}
生信了（公众号同名）
关注
专栏目录
全局序列比对 c语言实现,序列比对(一)——全局比对Needleman-Wunsch算法
weixin_28755447的博客
05-22 1340
原创：hxj7前言序列比对是生信领域的一个古老课题，在这一波NGS的浪潮中重新引起大家的广泛关注。由于生物序列的特殊性，在比对的时候允许插入缺失，所以往往是一种不精确匹配。从本文开始，我们陆续学习前人开发出的流行算法。全局比对算法所谓全局比对算法，就是根据一个打分矩阵(替换矩阵)计算出两个序列比对最高得分的算法。关于它的介绍网上已经非常多了，我们只需看看其中的关键点及实现代码。关键点1. 打...
详解基于python的全局与局部序列比对的实现(DNA)
12-16
程序能实现什么 a.完成gap值的自定义输入以及两条需比对序列的输入 b.完成得分矩阵的计算及输出 c.输出序列比对结果 d.使用matplotlib对得分矩阵路径的绘制 一、实现步骤 1.用户输入步骤 a.输入自定义的gap值 b.输入需要比对的碱基序列1（A,T,C,G）换行表示输入完成 b.输入需要比对的碱基序列2（A,T,C,G）换行表示输入完成 输入（示例）： 2.代码实现步骤 1.获取到用户输入的gap，s以及t 2.调用构建得分矩阵函数，得到得分矩阵以及方向矩阵 3.将得到的得分矩阵及方向矩阵作为参数传到回溯函数中开始回溯得到路径，路径存储使用的是全局变量，存的仍然是方向而
基序发现问题和中间字符串问题
深未来技术
12-31 585
1、如果P(s)表示对应于起始位点s的剖面矩阵，我们用MP(s)(j)来表示P（s）第j列中的最大计数。定义共有序列得分Score(s,DNA)为： l∑ MP(s)(j)，可用来度量起始点为s的序列谱的强度。j=12、给定一组DNA序列，从每条序列中发现一组L元组片段，使得共有序列得分为最大。输入：一个t*n阶的DNA矩阵，以及试图发现的模式的长度L输出：含t个起始位点的数组S=(S
字符串和序列的比较
weixin_43455016的博客
10-24 1548
字符串和序列的比较 字符串是根据字母的字符排列顺序进行比较的 &gt;&gt;&gt; "alpha"&lt;"beta" True 涉及大写字母时，排列顺序可能不同 &gt;&gt;&gt; "a"&lt;"B" False &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
DNA序列动态规划法比对的C代码
12-29
DNA序列动态规划法比对的C代码,生物信息学的内容，动态规划法
序列比对（19）基序发现和中间字符串问题
生信了
08-19 917
本文介绍了基序发现问题和中间字符串问题。 引言：DNA调控元件 我们知道，DNA调控元件往往是一段相似的DNA序列。理想情况下这些序列完全一致，比如下面这样： 图片引自《生物信息学算法导论》 但实际上，这些序列不会完全一样，总会有若干位点发生“变异”，从而不同，比如下面这样： 图片引自《生物信息学算法导论》 如果给定一组DNA序列（暂且假定它们长度相等），那么如何找出这些相似的序列呢？由此...
Java基于动态规划法实现求最长公共子序列及最长公共子字符串示例
08-27
这两个概念在计算机科学中尤其是在字符串处理和序列比对领域非常重要。 动态规划是一种解决优化问题的策略，它通过将大问题分解为相互重叠的子问题来逐步求解。在求解最长公共子序列时，我们不是简单地将两个字符串...
使用Java实现的计算两字符串相似度+最长公共子序列.zip
最新发布
05-10
在计算机科学领域，字符串处理是常见的一类任务，其中计算两个字符串的相似度以及找到它们的最长公共子序列（Longest Common Subsequence, LCS）是一个经典的问题。本篇将深入探讨这个问题，以及如何使用Java来解决...
DELPHI 计算两个字符串相似度 LCS算法（附源代码）
11-30
- LCSUnit.pas: 这是DELPHI的单元文件，包含了LCS算法的实现代码，可能包括字符串处理和动态规划数组的声明及填充逻辑。 - LCSProject.res: 资源文件，可能包含应用程序的图标和其他资源。 通过这些文件，你可以...
生物序列两两比对算法的研究和实现.doc
05-26
**后缀树**是一种用于存储字符串所有后缀的特殊树形数据结构，因其在处理字符串匹配问题上的高效性而在生物信息学中被广泛应用。尤其是**Ukkonen算法**，这是一种经典的线性时间内构造后缀树的方法，它不仅速度快，...
序列比对软件
08-19
生物软件 序列比对软件 源代码没有
基因序列比对
01-19
通过编程实现课程中介绍的全局比对和局部比对的动态规划算法，并对“data.txt”中的2 条序列（每行代表1 条序列）进行测试，打分矩阵使用BLOSUM62 矩阵（BLOSUM62.txt）。
序列比对(一)——全局比对Needleman-Wunsch算法
热门推荐
生信了
04-08 2万+
原创：hxj7 前言 序列比对是生信领域的一个古老课题，在这一波NGS的浪潮中重新引起大家的广泛关注。由于生物序列的特殊性，在比对的时候允许插入缺失，所以往往是一种不精确匹配。从本文开始，我们陆续学习前人开发出的流行算法。 全局比对算法 所谓全局比对算法，就是根据一个打分矩阵（替换矩阵）计算出两个序列比对最高得分的算法。关于它的介绍网上已经非常多了，我们只需看看其中的关键点及实现代码。 关键点 1...
【数据结构与算法】子串与子序列问题
qq_23269761的博客
09-07 966
1.替换最多k字符，找到最长的xxx字符串 问题描述： 有一个01字符串，有之多k次将0替换成1的机会，求最长的全1子串 输入：k=2 s=’1001010101’ 输出：5 （即‘10101’变为‘11111’） 方法1： 本质思想就是：不定长滑动窗口法，但是要保证窗口内0的个数是相同的 def deal(s,k): change = 0 left = 0 ...
算法——字符串、序列问题
上官宏竹的专栏
07-06 754
1、求一个字符串中连续出现次数最多的子串《程序员面试宝典》 P237 字符串中连续出现次数最多的子串，例如字符串“abababc”,最多连续出现的为ab，连续出现三次。 【算法思路】 首先 获得后缀数组 例如：abcbcbcabc的后缀数组为： abcbcbcabc bcbcbcabc cbcbcabc bcbcabc cbcabc
截取特殊字符中间的字符串
夜、醉梦红尘
11-24 1809
方法一： 采用函数strtok的方式截取，示例如下： [html] view plain copy #include stdlib.h>   #include iostream>      using namespace std;      int main(int argc, char* argv[])   {       
序列比对（待续）
Rys_Ben_Blogs
07-11 1129
一、双序列比对 动态规划（Dynamic programming） 动态规划 (Dynamic programming, DP) 是通过把原问题分解为相对简单的子问题的方式求解复杂问题的方法。动态规划常常适用于有重叠子问题和最优子结构性质的问题。 https://zh.wikipedia.org/wiki/动态规划 重叠子问题：自底向上 下面三个函数分别使用了递归、递归+列表、迭代方法求解斐波那契数列。后两种方法仅对重叠子问题求解一次，因而时间复杂度急剧下降。其中迭代呈现的自底向上策略就是动态规划对重
算法设计二（5）—— 字符串比较问题
weixin_43973089的博客
04-28 1226
代码链接：pan.baidu.com/s/1w-VSMWmr9ntEWWdxZoD4Yw 码：jnlh 算法分析与设计第 2 次实验 时间 2020.4.23 实验名称 字符串比较问题 实验目的 通过上机实验，要求掌握动态规划算法的问题描述、算法设计思想、程序设计。......
多序列比对 c语言代码,如何在R中执行基本的多序列比对？
weixin_35687240的博客
05-23 305
(我已经尝试过在biostar上问过这个问题了,但是由于文本挖掘中的某个人认为有更好的解决方案的可能性很小,我也在这里重新发布)我试图实现的任务是对齐几个序列.我没有匹配的基本模式.我所知道的只是“真实”模式的长度应为“30”,并且我在随机点引入了缺失值的序列.这是一个这样的序列的例子,在左边我们看到了缺失值的真实位置,在右边我们看到了我们能够观察到的序列.我的目标是仅使用我在右栏中获得的序列重建...
动态规划解决最长公共子序列问题及代码实现
这段代码提供了实现最长公共子序列算法的一个实例，适合学习和理解动态规划在解决字符串匹配问题中的应用。由于代码已经包含了详细的注释，初学者可以通过阅读和运行代码来深入理解算法的工作原理。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值