背景

这两天给神经科学研究所的学长们写一个算法，就是要在隐球菌的染色体上提取出mRNA和外显子，并用已知的84个ESE去匹配，看各个ESE分别出现多少次，问题的实质其实就是个KMP算法的字符串匹配，希望我的程序对他们的研究有帮助。

步骤

1.爬虫

爬虫的部分不是我写的，是让一个同学用python爬的，我没有接触过多少python，但是听说不是很难，所以最近有时间的话去学一学爬虫，这样以后如果再有类似的工作我自己都可以做了，因为两个人做的话格式要沟通，而且两个人的时间也很难统一，不见得一定会更有效率。

爬取染色体的全部序列

这个过程比较简单，不详细写了，爬到的结果如下：

GAATTCTAAAACAGTTGCATTCTATAATTACAAAATAATTGAAACACTTCCATTCTCTTGACTAATATTATTTAATTAGACACCAACTCGACATTCTGTCTTCGACCTATCTTTCTCCTCTCCCAGTCATCGCCCAGTAGAATTACCAGGCAATGAACCACGGCCTTTCATCCCAACGGCACAGCATATGGGTTCACTCCAACAGTGAACCATTCCAAAAGGCCTTGCCTGCGTGGCCATCTCCTCACAAACCCACCATCCTGCATCATCTCAGGTATCATACCTTCGATTCCTTCATCACCCCAAAGATCTTTCTGTTTGCCTTTGCGATTTGAGTGATACCAAACAAGAGAATACAGAAAGTTAGGGACAAAGGGGGAGGGCTCTCTGAAGCGGCATGCTCCTCTTTCAACATGAAAATGGAAAGCCCCATCATCCATGTGGTATATGGGCAGGCAGCAAACCTTACAAACACCACTTTTCAGATTTTCCAGAGCATCCAACCAGT…

一个文件是一个染色体的全部序列，所以一共会有14个文件，分别命名为string1~string14。

爬取外显子的起始位置

首先在NCBI的网站上（https://www.ncbi.nlm.nih.gov/genome/61?genome_assembly_id=52487）找到隐球菌的染色体，14条。

然后点第一条染色体（NC_026745.1）搜索mRNA。

每个mRNA后面的xxx…xxx代表一个外显子，所以爬虫的结果就是要爬取这些数对，但是过程中有的是单个的一个数字，学长说不用考虑这些单个的数字。

爬虫的结果如下。

1 c
100 1497,1556 1767,1823 2274,2322 2782,2846 3126,3206 3890,3958 5263,5322 5422,5494 5645 -1
2 c
5928 7626,7685 7780,7857 7982 -1
3 c
5928 7626,7685 7785,7857 7982 -1
4 c
8766 8986,9041 9050,9102 9246,9311 9535,9603 -1
5 c
10835 10923,11007 11307,11423 11573,11638 11664 -1

…

格式：同样是每个文件代表一个染色体，第一行第一个数是这个染色体的mRNA序号，后面的字母如果是j就正常处理（稍后说什么是正常处理），如果是c的话就要先把序列反向互补后再正常处理。

接下来以逗号作为分割，比如：100 1497,1556 1767,代表两个字符串，分别是染色体中序号100~1497和1556~1767的两个序列，用逗号分割的话可以判断哪个数字是单个的，不成对，我就不用考虑这个数字了（其实这种数据爬虫时删掉最好，我处理的时候就很方便了。但是爬虫的同学貌似不会…），他这爬虫还有个问题就是-1前面没有逗号，当时他在上课所以我还要特殊处理。。。每个mRNA以-1为终止。这样最后会有14个标记每个染色体的每个mRNA的每个外显子的起始位置的文件，命名为1.txt~14.txt。

2.处理出ESE所在区间

爬虫得到的两个文件作为我第一个程序的输入文件，输出的是ESE所在的区间，这个区间是这么定义的：还是以上面数据100 1497为例，把这个外显子后面3个碱基删掉，取后67个碱基，这就是ESE所在的区间，也就是我上面说的正常处理。当然，如果标记为c的话就要先把100到1497反向互补，然后在这样截取。

值得一提的是，会有外显子长度在70以内，这样的数据我们只删掉后三个碱基，然后输出剩下的碱基。

所有这些碱基是我第一个程序的输出文件（分别命名为output1~output14），也就是说我要在这些字符串中找那84个ESE。

输出格式（flag标记每个mRNA的结束位置，方便后面处理）

1
CTTTCCGGCACCAAGTAACGGACCATTAAGATGACCCGCTACTGACCCTCTCCTCTTTCTATCCAGC
CTAACGACTCCCGTACGACCTAATCACTTCACTGGCTTTATTCCCACGGTCACTGTAACTTAAACCG
GAAAAACAGTCATCTTCTGCAAGGAGCTTTTCTTAAGGTTCGAACCAATCGTCGGTTACGAGTGGAT
TACGATGGTAGTGGACCTTGATGAGTTCTGGTTCGTCACTATCCTAGTCGGTTGCGTCTCCTCGCAG
GATACCCTAAGTGATAGGTTTGAACCGCCGTACGTTTACAACCATATCACGCGGTTATGTCTTACTT
CCAGACCCTACCAGAACACTGCGATGATAGAGTAATGGTGTTTCCGGTTCGGAATTCACCCTTCTCT
CTCGACACCAGTAGTTCTTTTAAGACCAACTAAACCCGTGCTATTCTGTCTGTGGAGTTTGCAATGT
CTTCCTTCTCCACACGACGTGGTCTACGTTCCCTTCTTAAGGTCTCCTCTAAAAATGTAGTACCGTC
GTAGGTAAGTAAAGTAGGTAGGTGGGTAGGTAGGTAGGTAGATAGGTAGGTAGGTAGATGAATAGTA
flag
2
AGCTGCCAGGGGATTCGAGTGAACCACCTAAAGGTGAACGTCAACACGTACTTCTAGTAAAAAAAAA
CTTCCTTCTCCACACGACGTGGTCCACGTTCCCTTCTTACGATCTCCTCTAAAAATGTAGTACCGTC
GGTAGGTAGGTAGGTAGATAGGTAGGTAGGTAGGTAGATAGGTAGGTAGGTAGGTAGATGAATAGTA
flag
3
AGCTGCCAGGGGATTCGAGTGAACCACCTAAAGGTGAACGTCAACACGTACTTCTAGTAAAAAAAAA
CTTCCTTCTCCACACGACGTGGTCCACGTTCCCTTCTTACGATCTCCTCTAAAAATGTAGTACCGTC
GGTAGGTAGGTAGGTAGATAGGTAGGTAGGTAGGTAGATAGGTAGGTAGGTAGGTAGATGAATAGTA
flag

…

这段程序的代码：

#include<cstdio>
#include<string>
#include<cstring>
#include<iostream>
using namespace std;
FILE *fp0, *fp1, *fp3; //, *fp2
char all[9000005];
int a, b; //每对区间端点
char c; //mRNA是否带c
char ch; //string中字符
int n; //输入的mRNA序号
int cnt = 1;  //输出的mRNA序号
int pos;  //string中位置
void withC(){
//cout << “test” << endl; if(b-a+1 >= 70){
for(int i = a+3+66; i >= a+3; i—){
ch = all[i];
if(ch == ‘A’) ch = ‘T’;
else if(ch == ‘T’) ch = ‘A’;
else if(ch == ‘C’) ch = ‘G’;
else if(ch == ‘G’) ch = ‘C’;
        fprintf(fp3, "%c", ch);
    }
}else{
    for(int i = b; i &gt;= a+3; i--){
        ch = all[i];
        if(ch == 'A') ch = 'T';
        else if(ch == 'T') ch = 'A';
        else if(ch == 'C') ch = 'G';
        else if(ch == 'G') ch = 'C';

        fprintf(fp3, "%c", ch);
    }
}
}
void withoutC(){
if(b-a+1 >= 70){
for(int i = b-69; i <= b-3; i++)
fprintf(fp3, “%c”, all[i]);
}else{

    for(int i = a; i &lt;= b-3; i++)
        fprintf(fp3, "%c", all[i]);

}
}
int main(){
fp0 = fopen(“string12.txt”, “r”);
fp1 = fopen(“12.txt”, “r”);
//fp2 = fopen(“output11.txt”, “w”);//原始串
fp3 = fopen(“output12.txt”, “w”);//-3
pos = 1;
    all[0] = ' ';
    fscanf(fp0, "%s", all+1);
    //cout &lt;&lt; all &lt;&lt; endl;



while(~fscanf(fp1, "%d %c", &n, &c)){
    //cout &lt;&lt; n &lt;&lt; endl;
    //cout &lt;&lt; c &lt;&lt; endl;
    //输出mRNA序号

    //fprintf(fp2, "%d\n", cnt);
    fprintf(fp3, "%d\n", n);
    cnt++;

    while(~fscanf(fp1, "%d", &a)){

        if(a == -1) break;
        char flag;
        fscanf(fp1, "%c", &flag);
        if(flag == ',') continue;

        fscanf(fp1, "%d", &b);
        if(b == -1) break;
        fscanf(fp1, "%c", &flag);



        if(c == 'c')
            withC();
        else withoutC();
//cout << a << ” ” <<b <<endl;
//fprintf(fp2, “\n”);
        fprintf(fp3, "\n");

    }
        fprintf(fp3, "flag\n");
}

fclose(fp0);
fclose(fp1);
//fclose(fp2);
fclose(fp3);
return 0;
}

3.KMP算法匹配出每个外显子中的每个ESE的个数

ese.txt是已知的84个ESE碱基，如下：

CCTGGA
CGAGGA
CGAAGA
CTGAAG
CAGAAG
CAAGGA
CAAGAT
CAAGAA
……

然后就是算出output文件中的每个字符串内分别有ese.txt中的每个字符串的次数，还有一些总和的信息，很明显的KMP的匹配问题。KMP算法我之前转载的一个文章内有介绍，这里不细说。

最后输出的是一个可以用excel打开的txt文档（用\t分隔），命名为A1.txt~ A14.txt：

Exon CCTGGA CGAGGA CGAAGA CTGAAG CAGAAG CAAGGA CAAGAT CAAGAA CAAAGA GCAGAA GCAAGA GGAGCA GGAGGA GGAGAT GGAGAA GGAAGA GTGAAG GTTGGA GACCTG GACCAG GACGAA GACTTC GACATC GAGGAG GAGGAT GAGGAA GAGAAG GAGAAA GATGCA GATGGA GATGAA GAAGCA GAAGGA GAAGTT GAAGTA GAAGAC GAAGAG GAAGAT GAAGAA GAAACT GAAAGA TCCTGG TCAAGA TGTGGA TGAGAA TGAAGC TGAAGG TGAAGA TTGGAT TACCTG ACTGGA ACAGAA ACAAGA AGGAAC AGTGAC AGACGA AGAGGA AGAGAA AGAAGC AGAAGT AGAAGA AGAAAC ATCTGC ATTGGA AACCTG AACCAG AACTGG AACTTC AACTAC AACAGA AACAAC AAGGAC AAGGAA AAGACA AAGAGA AAGATC AAGAAC AAGAAG AAGAAT AAGAAA AATGGA AATGAC AAAGGA AAAGAA Total Sequence Sum_mRNA_ESE_count AVG_mRNA_ESE_count Sum_染色体_ESE_count
1-1-1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 CTTTCCGGCACCAAGTAACGGACCATTAAGATGACCCGCTACTGACCCTCTCCTCTTTCTATCCAGC
1-1-2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 CTAACGACTCCCGTACGACCTAATCACTTCACTGGCTTTATTCCCACGGTCACTGTAACTTAAACCG

看起来比较乱，用excel打开它：

成果不错吧！

这段的代码：

#include<cstring>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<vector>
#include<string>
#include<iostream>
#include<algorithm>
using namespace std;
FILE *fp0, *fp1, *fp2;
int n;
int nextt[100][10];
char ese[100][10];
char x[50], y[100];//x:ESE模式串 y:主串（长度67）
void kmp_pre(char x[],int m,int w){
int i, j;
j = nextt[w][0] = -1;
i = 0;
while(i < m){
while(-1!=j && x[i]!=x[j])
j = nextt[w][j];
nextt[w][++i] = ++j;
}
}
void init(){
fp1 = fopen(“ese.txt”, “r”);
int i = 0;
while(~fscanf(fp1, “%s”, x)){
strcpy(ese[i], x);
kmp_pre(x, 6, i);
i++;
}
fprintf(fp2, “Exon\t”);
for(int i = 0; i < 84; i++)
fprintf(fp2, “%s\t”, ese[i]);
fprintf(fp2, "Total\tSequence\tSum_mRNA_ESE_count\tAVG_mRNA_ESE_count\tSum_染色体_ESE_count\n");
}
/*

返回 x 在 y 中出现的次数，可以重叠
*/

int KMPCount(char x[], int m, char y[], int n, int w){//x 是模式串， y 是主串
int i, j;
int ans = 0;
//preKMP(x,m,nextt);
//kmppre(x, m, nextt);
i = j = 0;
while(i<n){ while(-1!=j && y[i]!=x[j]) j = nextt[w][j]; i++; j++; if(j>=m){
int p = i-m+1;
//fprintf(fp2, “%d ”, p); //输出起点位置
ans++;
j = nextt[w][j];
}
}
//fprintf(fp2, “\n”);
return ans;
}
int main(){
fp0 = fopen(“output14.txt”, “r”);
fp2 = fopen(“A14.txt”, “w”);
init();
int sum1 = 0; //染色体
int sum2 = 0; //mRNA
int sum3 = 0; //外显子
while(~fscanf(fp0, "%d", &n)){
    //fprintf(fp2, "mRNA序号：%d\n", n);

    int qq = 1;

    while(1){
        fscanf(fp0, "%s", y);
        if(!strcmp(y, "flag")){
            fprintf(fp2, "%d\n", sum2);
            sum2 = 0;
            break;
        }else{
            if(qq != 1)
                fprintf(fp2, " \n");
        }
        fprintf(fp2, "14-%d-%d\t", n, qq++);
        for(int i = 0; i &lt; 84; i++){
            int ans = KMP_Count(ese[i], 6, y, strlen(y), i);
            fprintf(fp2, "%d\t", ans);
            sum3 += ans;

        //if(ans != 0)
        //    fprintf(fp2, "\n%s %d\n", ese[i], ans);
        }
        sum1 += sum3;
        sum2 += sum3;
        //fprintf(fp2, "该mRNA 外显子%d ESE和为：%d\n", qq++, sum3);
        fprintf(fp2, "%d\t%s\t", sum3, y);
        sum3 = 0;
    }
    //fprintf(fp2, "该mRNA ESE和为：%d\n", sum2);
    //sum2 = 0;

}
for(int i = 0; i &lt; 88; i++)
    fprintf(fp2, " \t");
fprintf(fp2, "%f\t%d\t", 1.0*sum1/n, sum1);

//fprintf(fp2, "该染色体 ESC和为：%d\n", sum1);
//fprintf(fp2, "该染色体 mRNA上ESE的平均数为：%f\n", );
/*
for(int i = 0; i &lt; 84; i++){
    for(int j = 1; j &lt;= 6; j++)
        cout &lt;&lt; nextt[i][j];
    cout &lt;&lt; endl;
}*/
fclose(fp0);
fclose(fp1);
fclose(fp2);
return 0;
}

感想

这是我第一次意识到ACM在实际工程和研究中的重要性，以前只是为了刷题，没有太多的体会。听说现在很多找工作和实习的学长都在各种OJ平台刷算法题，我这种搞ACM的会有些优势吧。