标签：
这是一个山科大的同学给我的一个问题，向我询问一下思路，对于数学建模，我没太多的了解，所以只能用计算机程序的方法来解答。
这是具体的问题：
这个问题来自 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
2015年“深圳杯”数学建模夏令营-B题：DNA序列的k-mer index 问题
标签：
原文地址：http://blog.csdn.net/mayuko2012/article/details/45667215
