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

一只羊

发布于 2019-08-29 16:38:13

8970

发布于 2019-08-29 16:38:13

文章被收录于专栏：生信了生信了

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

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

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

由于要遍历所有可能的起始位点，如前文《序列比对（20）基序发现问题的算法及实现代码》一样，我们采用树结构以及DFS（深度优先搜索）。其简单算法如下（伪代码）：

引自《生物信息学算法导论》

引自《生物信息学算法导论》

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

用分支定界策略改进上述简单算法的伪代码如下：

引自《生物信息学算法导论》

引自《生物信息学算法导论》

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

实验效果

我们用《生物信息学算法导论》书中的两组序列进行实验，分别用简单算法（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);
}

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2019-08-29，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法

存储

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

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

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

编程算法

存储

登录后参与评论

0 条评论

热度