2015年“深圳杯”数学建模夏令营-B题：DNA序列的k-mer index 问题

这是一个山科大的同学给我的一个问题，向我询问一下思路，对于数学建模，我没太多的了解，所以只能用计算机程序的方法来解答。

这是具体的问题：

这个问题来自 DNA序列的k-mer index问题。

给定一个DNA序列，这个系列只含有4个字母ATCG，如 S =“CTGTACTGTAT”。给定一个整数值k，从S的第一个位置开始，取一连续k个字母的短串，称之为k-mer（如k= 5，则此短串为CTGTA）， 然后从S的第二个位置， 取另一k-mer（如k= 5，则此短串为TGTAC），这样直至S的末端，就得一个集合，包含全部k-mer 。 如对序列S来说，所有5-mer为

｛CTGTA，TGTAC，GTACT，TACTG，ACTGT，CTGTA，TGTAT｝

通常这些k-mer需一种数据索引方法，可被后面的操作快速访问。例如，对5-mer来说，当查询CTGTA，通过这种数据索引方法，可返回其在DNA序列S中的位置为｛1，6｝。

问题

现在以文件形式给定 100万个 DNA序列，序列编号为1-1000000，每个基因序列长度为100 。

（1）要求对给定k， 给出并实现一种数据索引方法，可返回任意一个k-mer所在的DNA序列编号和相应序列中出现的位置。每次建立索引，只需支持一个k值即可，不需要支持全部k值。

（2）要求索引一旦建立，查询速度尽量快，所用内存尽量小。

（3）给出建立索引所用的计算复杂度，和空间复杂度分析。

（4）给出使用索引查询的计算复杂度，和空间复杂度分析。

（5）假设内存限制为8G，分析所设计索引方法所能支持的最大k值和相应数据查询效率。

（6）按重要性由高到低排列，将依据以下几点，来评价索引方法性能 

· 索引查询速度

· 索引内存使用

· 8G内存下，所能支持的k值范围

· 建立索引时间

 

要想做出来一个索引，就必须知道子序列的位置，要求对子序列编号。

其实还有两个70多M的文档，以.fa结尾，

经过百度查了一下，fa的文件是一个数据文件。全名是：FASTA Formatted Sequence File。

它属于生物信息学中数据文件格式，是一种基于文本用于表示核苷酸序列或氨基酸序列的格式。在这种格式中碱基对或氨基酸用单个字母来编码，且允许在序列前添加序列名及注释。

试着用记事本打开了一下，结果直接卡住了。

那么，这些应该是需要读取的数据了。我们利用二进制方式进行读取。

70多M，1000000个100位的碱基序列，要求在8G内存下建立索引。

既然是建立一个索引，那么对于给定的一个子序列，用字符串的方式来查找，真的很费时间，所以，我建立了一种规定，即：

将A替换成1；

将T替换成2；

将C替换成3；

将G替换成4。

我们可以这样定义代码：

//这是目前只是数字值转化的代码：
#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
int getnum(int,int);
char s[100];
int a[100];//保存数值的数组
int numall=0;
int main()
{
    int len,k;//len:求出有多少个序列 k:k个字符串
    int i,j;
    gets(s);
    cin>>k;//输入整数值k
    len=strlen(s);//获取序列的长度
    for (i=0,j=0; i<(len-k+1); i++,j++)
    {
        a[i]=getnum(i,k);
        cout<<a[i]<<endl;
    }
}
int getnum(int n,int k)//获取单个序列的数字值，并将结果保存在a数组中
{
    int numb,l,m=0;
    l=n;
    numall=0;
    for (; n<k+l; n++,m++)
    {
        numb=1;
        for (int i=0; i<k-m-1; i++)
            numb=numb*10;
        if (s[n]=='A')
            numb*=1;
        else if (s[n]=='T')
            numb*=2;
        else if (s[n]=='C')
            numb*=3;
        else if (s[n]=='G')
            numb*=4;
        numall+=numb;
    }
    return numall;
}

我们将例题中的碱基序列输入进去：

这样，进行数字进行查找，在代码上也会比较方便一些。

另一个是索引的建立，既然查找，我们首先想到的是从头到位一个一个的进行对比，符合了，找到了。

但是，在内存限制下，1000000个序列对于计算机也有些困难。

按照我们的思想，首先进行排序。

冒泡？选择？木桶排序？那个比较好？其实我也在忧郁，其实，最先想到的还是冒泡排序，所以，我们可以选择冒泡排序。

要记住，题目要求还有获得具体的位置，所以我们应该用一个另外一个数组和原数组做对应，也就是进行编号。

排序完成，编号数组不能变，原数组的规则是从小到大。

给定一个子序列，我们利用二分法获得具体位置，再找到编号，从而获得它在原数组的位置。

最后根据规则反译序列。

所有的流程，大概是这样子：

纯属个人见解，如有错误之处还望改正。

@ Mayuko

转载于:https://www.cnblogs.com/mayuko/p/4567500.html