序列比对（20）基序发现问题的算法及实现代码

最新推荐文章于 2021-05-23 10:36:31 发布

生信了（公众号同名）

最新推荐文章于 2021-05-23 10:36:31 发布

阅读量1k

点赞数 1

分类专栏： # 序列算法文章标签：多序列比对基序

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

版权

序列算法专栏收录该内容

32 篇文章 16 订阅

订阅专栏

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

基序发现问题的简单算法及伪代码

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

简单回顾一下，基序发现问题其实就是要找到使得共有序列得分最大的一组起始位点。
$\underset{\boldsymbol{\vec{s}}}{\mathrm{argmax}} \, Score(\boldsymbol{\vec{s}},DNA)$

由于要遍历所有可能的起始位点，所以一种自然的想法是使用递归。但是为了配合后续的分支定界法，我们采用了树结构，并且进行DFS（深度优先搜索）。既然采用树结构，最简单的算法如下（伪代码）：
在这里插入图片描述

图1到图2引自《生物信息学算法导论》

从上面的代码可以看出，算法构建了一个深度为 $t$ 的搜索树，除了叶顶点外，每个顶点有 $n - l + 1$ 个子顶点。当遍历到叶顶点时，计算共有序列得分，如果该得分比目前的最大得分还高，那么就将该得分记为新的最大得分，并且将此时的起始位点向量记录下来。

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

上述简单算法的效率是非常低的，具体实验效果参见后文。效率低是由于树的顶点数量过大，遍历每一个顶点所需的总时间过长。鉴于此，可以采用一种分支定界策略有选择地跳过一些分支，即以某个顶点为根的分支树中的所有顶点都被跳过，从而可能大大提升效率。

我们可以将分支定界策略这样应用到基序发现问题上：我们发现，一个深度为 $i$ 的顶点，假设其“共有序列得分”记为 $Score(\boldsymbol{\vec{s}}, i, DNA)$ ，那么以该顶点为根的分支树中的所有叶顶点，其共有序列得分的上界是 $Score(\boldsymbol{\vec{s}}, i, DNA) + (t - i) \times l$ 。如果这个上界分值比目前的最大得分还小，那么遍历这个分支树就没有意义了。因此，我们可以跳过这个分支树。这样做与上述简单算法比较的话，多出来的开销就是“至多”要为从根到此深度为 $i$ 的顶点的路径上的共 $i$ 个顶点计算“共有序列得分”，但是一旦跳过分支树，那么至少不用计算 $n - l + 1)^{t-i}$ 个叶顶点的“共有序列得分”。从实际效果来看，往往还是可以节省很多时间的。

用分支定界策略改进上述简单算法的伪代码如下：
在这里插入图片描述

图3到图4引自《生物信息学算法导论》

实验效果

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

CGGGGCTATGCAACTGGGTCGTCACATTCCCCTTTCGATA
TTTGAGGGTGCCCAATAAATGCAACTCCAAAGCGGACAAA
GGATGCAACTGATGCCGTTTGACGACCTAAATCAACGGCC
AAGGATGCAACTCCAGGAGCGCCTTTGCTGGTTCTACCTG
AATTTTCTAAAAAGATTATAATGTCGGTCCATGCAACTTC
CTGCTGTACAACTGAGATCATGCTGCATGCAACTTTCAAC
TACATGATCTTTTGATGCAACTTGGATGAGGGAATGATGC

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

CGGGGCTATcCAgCTGGGTCGTCACATTCCCCTTTCGATA
TTTGAGGGTGCCCAATAAggGCAACTCCAAAGCGGACAAA
GGATGgAtCTGATGCCGTTTGACGACCTAAATCAACGGCC
AAGGAaGCAACcCCAGGAGCGCCTTTGCTGGTTCTACCTG
AATTTTCTAAAAAGATTATAATGTCGGTCCtTGgAACTTC
CTGCTGTACAACTGAGATCATGCTGCATGCcAtTTTCAAC
TACATGATCTTTTGATGgcACTTGGATGAGGGAATGATGC

可惜的是，简单算法的效率过低，当序列增至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(Seq* mulSeq, const int t, const int l, int* a, int d);     /* 得到下一个顶点 */
int maxArray(int* a, const int n);
int getScore(Seq* mulSeq, int* a, int* tmp, const int t, const int l);        /* 计算Profile的得分 */
int indexOfMax(int* a, const int n);
char* getConsensus(Seq* mulSeq, int* a, int* tmp, const int t, const int l);        /* 获取共有序列 */
void findMotif(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);
    findMotif(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(Seq* mulSeq, const int t, const int l, int* a, int d) {
/* 得到下一个顶点 */
    int i;
    if (d < t - 1) {
        a[d + 1] = 0;
        return d + 1;
    }
    for (i = t - 1; i >= 0; i--) {
        if (a[i] < mulSeq[i]->l - l) {
            a[i]++;
            return i;
        }
    }
    return -1;
}

int maxArray(int* a, const int n) {
    int i, max;
    for (i = 1, max = a[0]; i < n; i++)
        if (a[i] > max)
            max = a[i];
    return max;
}

int getScore(Seq* mulSeq, int* a, int* tmp, const int t, const int l) {
/* 计算Profile的得分 */
    int i, j, k;
    int score;
    // 计算每一个位置的得分（即Profile的每一列的得分）
    for (i = 0, score = 0; i < l; i++) {
        for (j = 0; j < nalpha; j++)
            tmp[j] = 0;
        for (j = 0; j < t; j++)
            tmp[mulSeq[j]->a[a[j] + i]]++;
        score += maxArray(tmp, nalpha);
    }
    return score;
}

int indexOfMax(int* a, const int n) {
    int i, m, max;
    for (i = 1, m = 0, max = a[0]; i < n; i++)
        if (a[i] > max) {
            max = a[i];
            m = i;
        }
    return m;
}

char* getConsensus(Seq* mulSeq, int* a, int* tmp, const int t, const int l) {
/* 获取共有序列 */
    int i, j, k;
    char* cs;    // consensus.
    if ((cs = (char*) malloc(sizeof(char) * (l + 1))) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);
    }
    // 计算每一个位置的共有字符
    for (i = 0; i < l; i++) {
        for (j = 0; j < nalpha; j++)
            tmp[j] = 0;
        for (j = 0; j < t; j++)
            tmp[mulSeq[j]->a[a[j] + i]]++;
        cs[i] = sortedAlpha[indexOfMax(tmp, nalpha)];
    }
    cs[l] = '\0';
    return cs;
}

void findMotif(Seq* mulSeq, const int t, const int l) {
/* 寻找基序 */
    int* a;     // 起始位点向量
    int* ma;    // 最大分值对应的起始位点向量
    int i, d, j;
    int score, maxScore;
    int* scoreTmp;
    char* cs;    // consensus.
    if ((a = (int*) malloc(sizeof(int) * t)) == NULL || \
        (ma = (int*) malloc(sizeof(int) * t)) == NULL || \
        (scoreTmp = (int*) malloc(sizeof(int) * nalpha)) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);
    }
    // 遍历所有可能的起始位点，寻找最大分值
    for (i = 0; i < t; i++)
        ma[i] = a[i] = 0;
    maxScore = 0;
    d = 0;
    while (d >= 0) {
        if (d < t - 1)
            d = nextVertex(mulSeq, t, l, a, d);
        else {
            score = getScore(mulSeq, a, scoreTmp, t, l);
            if (score > maxScore) {  // 注意最大值对应的起始位点向量可能有多种，这里只保留了其中一种。其实还可以设置当maxScore = l * t时终止搜索。
                maxScore = score;
                for (i = 0; i < t; i++)
                    ma[i] = a[i];
            }
            d = nextVertex(mulSeq, t, l, a, d);
        }
    }
    // 获取共有序列
    cs = getConsensus(mulSeq, ma, scoreTmp, t, l);
    // 输出结果
    printf("Max score: %d\n", maxScore);
    printf("Starting point vector: ");
    for (i = 0; i < t; i++)
        printf("%d ", ma[i] + 1);
    printf("\n");
    printf("Target seq:\n");
    for (i = 0; i < t; i++) {
        printf("\t");
        for (j = ma[i]; j < ma[i] + l; j++)
            printf("%c", mulSeq[i]->s[j]);
        printf("\n");
    }
    printf("Consensus seq: %s\n", cs);
    // 释放内存
    free(a);
    free(ma);
    free(scoreTmp);
    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(Seq* mulSeq, const int t, const int l, int* a, int d);     /* 得到下一个顶点 */
int byPass(Seq* mulSeq, const int t, const int l, int* a, int d);        /* 跳过分支，得到下一个顶点 */
int maxArray(int* a, const int n);
int getScore(Seq* mulSeq, int* a, int* tmp, const int d, const int l);        /* 计算Profile的得分 */
int indexOfMax(int* a, const int n);
char* getConsensus(Seq* mulSeq, int* a, int* tmp, const int t, const int l);        /* 获取共有序列 */
void findMotif(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);
    findMotif(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(Seq* mulSeq, const int t, const int l, int* a, int d) {
/* 得到下一个顶点 */
    int i;
    if (d < t - 1) {
        a[d + 1] = 0;
        return d + 1;
    }
    for (i = t - 1; i >= 0; i--) {
        if (a[i] < mulSeq[i]->l - l) {
            a[i]++;
            return i;
        }
    }
    return -1;
}

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

int maxArray(int* a, const int n) {
    int i, max;
    for (i = 1, max = a[0]; i < n; i++)
        if (a[i] > max)
            max = a[i];
    return max;
}

int getScore(Seq* mulSeq, int* a, int* tmp, const int d, const int l) {
/* 计算Profile的得分 */
    int i, j, k;
    int score;
    // 计算每一个位置的得分（即Profile的每一列的得分）
    for (i = 0, score = 0; i < l; i++) {
        for (j = 0; j < nalpha; j++)
            tmp[j] = 0;
        for (j = 0; j <= d; j++)
            tmp[mulSeq[j]->a[a[j] + i]]++;
        score += maxArray(tmp, nalpha);
    }
    return score;
}

int indexOfMax(int* a, const int n) {
    int i, m, max;
    for (i = 1, m = 0, max = a[0]; i < n; i++)
        if (a[i] > max) {
            max = a[i];
            m = i;
        }
    return m;
}

char* getConsensus(Seq* mulSeq, int* a, int* tmp, const int t, const int l) {
/* 获取共有序列 */
    int i, j, k;
    char* cs;    // consensus.
    if ((cs = (char*) malloc(sizeof(char) * (l + 1))) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);
    }
    // 计算每一个位置的共有字符
    for (i = 0; i < l; i++) {
        for (j = 0; j < nalpha; j++)
            tmp[j] = 0;
        for (j = 0; j < t; j++)
            tmp[mulSeq[j]->a[a[j] + i]]++;
        cs[i] = sortedAlpha[indexOfMax(tmp, nalpha)];
    }
    cs[l] = '\0';
    return cs;
}

void findMotif(Seq* mulSeq, const int t, const int l) {
/* 寻找基序 */
    int* a;     // 起始位点向量
    int* ma;    // 最大分值对应的起始位点向量
    int i, d, j;
    int score, maxScore;
    int* scoreTmp;
    char* cs;    // consensus.
    if ((a = (int*) malloc(sizeof(int) * t)) == NULL || \
        (ma = (int*) malloc(sizeof(int) * t)) == NULL || \
        (scoreTmp = (int*) malloc(sizeof(int) * nalpha)) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);
    }
    // 遍历所有可能的起始位点，寻找最大分值
    for (i = 0; i < t; i++)
        ma[i] = a[i] = 0;
    maxScore = 0;
    d = 0;
    while (d >= 0) {
        score = getScore(mulSeq, a, scoreTmp, d, l);
        if (d < t - 1) {
            if (score + (t - 1 - d) * l < maxScore)     // 分支定界法
                d = byPass(mulSeq, t, l, a, d);
            else
                d = nextVertex(mulSeq, t, l, a, d);            
        } else {
            if (score > maxScore) {  // 注意最大值对应的起始位点向量可能有多种，这里只保留了其中一种。其实还可以设置当maxScore = l * t时终止搜索。
                maxScore = score;
                for (i = 0; i < t; i++)
                    ma[i] = a[i];
            }
            d = nextVertex(mulSeq, t, l, a, d);
        }
    }
    // 获取共有序列
    cs = getConsensus(mulSeq, ma, scoreTmp, t, l);
    // 输出结果
    printf("Max score: %d\n", maxScore);
    printf("Starting point vector: ");
    for (i = 0; i < t; i++)
        printf("%d ", ma[i] + 1);
    printf("\n");
    printf("Target seq:\n");
    for (i = 0; i < t; i++) {
        printf("\t");
        for (j = ma[i]; j < ma[i] + l; j++)
            printf("%c", mulSeq[i]->s[j]);
        printf("\n");
    }
    printf("Consensus seq: %s\n", cs);
    // 释放内存
    free(a);
    free(ma);
    free(scoreTmp);
    free(cs);
}