基于C程序的DNA序列的kmerindex数据查找

1、论文题目：基于c程序的dna序列的k-mer index数据查找姓名 学院 年级 专业 学号 联系电话 数学分析 高等代数 高等数学 线性代数 概率统计 数学实验 数学模型 cet4 cet6 电气工程学院2013 级电气工程及其自动化918996560电气工程学院2013 级电气工程及其自动化918885578540电气工程学院2013 级电气工程及其自动化858989554基于c程序的dna序列的k-mer index数据查找 摘要dna 是生命体的基本遗传物质，其组成和序列变化创造了形形色色的生命世界。快速、准确地获取生物体的遗传信息对于生命科学的研究具有重要意义1。现需要给定一种数据索

2、引方法 ，利用一种查询算法查询百万条序列中是否存在相应的片段，如果存在，则输出相应片段所在的位置。针对问题一，运用karp-rabin算法，在c程序环境下编写字符串匹配算法。具体做法是将碱基序列映射成四进制的数串，对给定的k，构造合适的哈希函数，将四进制数串内每个长度为k的子数串译为唯一的十进制数，按顺序放进索引数组（哈希表）。查找相同的字符串等价于判断相应的hash值是否相同。此法可以大大提高建立索引和查询的时间。针对问题二，对不同k值经过大量多次的尝试，一般来说建立索引约10秒，查询约0.3秒（visual c+ 6.0运行环境）。针对问题三和问题四，本算法使用for循环函数，先计算出建立

3、索引与使用索引的算法中每一个语句的执行次数，然后再相加，最后依据去低阶项，去掉常数项，去掉高阶项的常参的原则得到时间复杂度。然后根据算法临时存储的空间大小来计算空间复杂度，观察有无临时存储大小以及临时存储大小与输入变量的关系。经分析，本算法在建立索引时的时间复杂度为o（n*m），空间复杂度为o（n）；在使用索引时的时间复杂度为o（n*m），空间复杂度为o（1）。针对问题五，首先分析了整形数据最大值的限制，k最大只能取14。为了扩展k的最大值，我们提出将字符串分成几个不相交的完备的字串，其长度不超过14，分别比较每一个字串，然后对结果取交集。并给出了这种方法下所需要的内存与k的关系，得出理论上可

4、以在8g的内存下支持所有k值的查找。关键词：dna 数据索引 查询算法 算法复杂度 karp-rabin 哈希函数c程序 字符串匹配一、 问题重述这个问题来自 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，gt

5、act，tactg，actgt，tgtat通常这些k-mer需一种数据索引方法，可被后面的操作快速访问。例如，对5-mer来说，当查询ctgta，通过这种数据索引方法，可返回其在dna序列s中的位置为1，6。问题现在以文件形式给定 100万个 dna序列，序列编号为1-1000000，每个基因序列长度为100 。（1）要求对给定k， 给出并实现一种数据索引方法，可返回任意一个k-mer所在的dna序列编号和相应序列中出现的位置。每次建立索引，只需支持一个k值即可，不需要支持全部k值。（2）要求索引一旦建立，查询速度尽量快，所用内存尽量小。（3）给出建立索引所用的计算复杂度，和空间复杂度分析。（

6、4）给出使用索引查询的计算复杂度，和空间复杂度分析。（5）假设内存限制为8g，分析所设计索引方法所能支持的最大k值和相应数据查询效率。（6）按重要性由高到低排列，将依据以下几点，来评价索引方法性能 索引查询速度索引内存使用8g内存下，所能支持的k值范围建立索引时间二、问题分析针对问题一，要求对给定k，给出并实现一种数据索引方法，可返回任意一个k-mer所在的dna序列编号和相应序列中出现的位置。以及每次建立索引，只需支持一个k值即可，不需要支持全部k值。首先利用c程序的中的fp文件读取函数将目标文件中的一亿个单字符数据读出，用于建立数据库，由于给定的是100万个dna序列，每个基因序列长度为1

7、00。而c程序中静态数组的范围限制远远小于100万，于是定义全局变量建立动态数组，动态分配储存空间，即建立行为100万，列为100的1000000100的动态数组，也就是索引的建立过程。采用karp-rabin算法2，可令a=0,c=1,g=2,t=3，将碱基序列映射成一个四进制数串。对于给定的k，可以构造这样一个集合，即这个集合囊括了四进制数串中所有长度为k的子数串。构造相应的哈希函数，将每个四进制子数串转换成唯一的十进制数字。待查kmer也做相同的处理，然后遍历整个集合里的十进制数，记录下与待查kmer对应十进制数相同的位置。由于转换的是唯一的十进制数字，所以只要相等，便找到了相应的k-m

8、er所在的dna序列编号和相应序列中出现的位置，然后以txt文件形式输出所有找到的dna序列编号以及相应的位置。针对问题二，要求索引一旦建立，查询速度尽量快，所用内存尽量小。首先要求查询速度非常快，所以对算法的要求就比较高，要求算法的时间复杂度与空间复杂度都应该比较低，也就是在c语言的程序中尽量少一些嵌入式循环，本次的索引算法中最多只有2个for循环嵌入使用，所以能够保证时间复杂度较低，故查询速度较快。本算法采用动态数组，并且最后的索引使用也与k值的大小有关，k值越大，使用的内存越少，所以能够在一定程度上保证所用的内存较小。针对问题三，此处给出建立索引的时间复杂度的分析。一个算法的空间复杂度只

9、考虑在运行过程中为局部变量分配的存储空间的大小，它包括为参数表中形参变量分配的存储空间和为在函数体中定义的局部变量分配的存储空间两个部分。所以本算法的空间复杂度第一要考虑所使用的堆栈空间的大小，它等于一次调用所分配的临时存储空间的大小乘以被调用的次数。本算法使用for循环函数，先计算出建立索引的算法中每一个语句的执行次数，然后再相加，最后依据去低阶项，去掉常数项，去掉高阶项的常参的原则得到时间复杂度。然后根据算法临时存储的空间大小来计算空间复杂度，观察有无临时存储大小以及临时存储大小与输入变量的关系。针对问题四，与建立索引的复杂度的分析相似，先找到 for循环函数，计算出使用索引的算法中每一个

10、语句的执行次数，然后再相加，最后依据去低阶项，去掉常数项，去掉高阶项的常参的原则得到时间复杂度。然后根据算法临时存储的空间大小来计算空间复杂度，观察有无临时存储大小以及临时存储大小与输入变量的关系。针对问题五，假设内存限制为8g，分析所设计索引方法所能支持的最大k值和相应数据查询效率。内存消耗最大的应该是程序当中的两个超大数组，一个是存放文件的数组，另一个是存放索引表（hash表）的数组。通过对两个数组所占内存的分析，我们得出了算法所需要的内存与k的函数关系。显然，k若越大，则对应的十进制数也越大，在整型数据有最大值的条件下，我们通过理论分析给出了k最大能取14。为了扩展k的值，我们考虑将长度

11、超过14的分成多个部分，每个部分限制最长为14。这样，判断两个字符串相等的条件成为要每个部分都相等，对应在算法里就是对每个部分的比较结果取交集。显然，这种方法会导致使用更多的大数组，相应会占更多的内存。三、符号说明符号说明fp文件d:bdata放置dna序列文件的位置solexa_100_170_1.fa前半部分dna序列文件solexa_100_170_2.fa后半部分dna序列文件shuju.txt存放查找出的kmerndna总序列数m每个序列的长度n宏定义1000000total 将4进制转换为10进制后的存储数组p读取文件时的二维存储数组k需查找的kmer序列的长度内存函数扩展内存函数

12、四、模型建立与求解4.1 问题一的模型建立与求解4.1.1文件处理对一百万条dna序列的读取采用c语言的文件操作语句循环读入。在每条序列循环读取的时候，首先判断是否为“eof”，若为空，则表示文件读取完毕，结束，若不为“eof”，则继续读取dna序列，若读取的是“atcg”中的任意一个，则以具体的碱基序列存入二维数组中，若不是“atcg”任意一个，继续从文件中读取下一个字符。算法流程图如下：将该字符放进一个二维字符数组结束是该字符是否为“atcg”中的一个否否从文件中读取一个字符，判断是否为eof开始打开文件是4.1.2 建立数据索引在建立索引时，首先从存放字符的数组中读取一个字符，判断是否为

13、该行（100个）的最后一个字母，如果是，则跳到下一行，如果不是，则分别令可令a=0,c=1,g=2,t=3，即转化为四进制数串。然后将长度为k的字串按递推的顺序放进索引数组。递推的方式为将整个长度为k的字符串舍去第一个字符，再加上下一个字符，成为一个新的长度为k的字符串。这种递推直到这一行的最后一个字符为止。如 s =“ctgtactgtat”。则字串为ctgta，tgtac，gtact，tactg，actgt，tgtat。对每一个字串，可以构造这样的哈希函数：其中sj指的是字串里从右往左数第j个字母对应的数值（对应规则为a=0,c=1,g=2,t=3）。每一个dna链只用计算一次哈希函数，后

14、面的k-mer对应的hash值（十进制数值）可以由以下的递推公式给出：其中sj+k是下一个字符对应的数值，sj是当前第一个字符对应的数值，这两个字符的距离刚好是k。 如此一来，就已经将每个k-mer与一个唯一的十进制数对应起来了。然后将这个十进制数存放到一个二维数组total1000000100里面去。每个k-mer在数组中的位置由如下规则确定：行号是所在dna链的编号，列号是这个k-mer的第一个字符在该列dna链中的位置（从左往右数）。这样就建立好了索引数组，程序中用total表示。算法流程图如下：读取下一个字符 结束是否所有字符转换为了相应的hash值是否是是否是第a个字符 （a=k）跳

15、到下一行从存放字符的数组中读取一字符开始将这个四进制数转化为十进制；存入索引数组是否是这一行最后一个字符是否a=0t=3g=2c=1用递推公式计算下一个hash值否是否是第k个字符 是否4.1.3 使用数据索引计算输入字符串的kmer值，与索引数组的total中的字符串依次比较，如果匹配，则说明找到了相应字符串在100万个序列中的位置，然后将位置输出到txt文件中，利用循环遍历整个索引数组，找出所有匹配的字符串的位置。算法流程图如下：否j=0，i+结束是i是否大于1000000j是否大于（100-）否是kmer值与索引数组totalij是否匹配开始i=0,j=0 计算输入字符串的kmer值是否

16、输出（i,j）到一个txt文件中j+是4.2问题三的模型建立与求解4.2.1建立索引所用的计算复杂度，和空间复杂度分析4.2.1.1建立索引时的时间复杂度分析：建立索引部分程序如下：该部分个语句需要运行的次数如下图所示：其中n1=1000000，m1=100，代表dna的1000000个序列，m代表每个序列的长度，h为输入的kmer值，该部分算法的时间复杂度等于各个语句的执行次数相加，即o（n1+n1+n1*h+n1*h+n1*h+n1*h+n1*（m1-h）+n1*（m1-h）+n1+n1）=o（2n1*m1+2n1*h+2n1*h+2n1）o表示去低阶项，去掉常数项，去掉高阶项的常参得到时

17、间复杂度o,即o（n1*m1）.4.2.1.2建立索引时的空间复杂度分析：临时存储空间大小的部分为：空间复杂度是对一个算法在运行过程中临时占用存储空间大小的量度，建立索引过程中数组hash是作为临时占用储存空间的数组，由于在for循环中有a=1）于是（k=1）（单位：mb）即k每增大1，所占用的内存就会减少3.8147mb。n（k）max=476.8374mb4.3.2 k的最大取值分析限制的因素：限制k的主要原因是整型数据的最大值只有231-1；而在建立哈希表的过程中要用到指数函数pow（4，k-1） ；k的最大值可以由以下的表达式确定：即，可以求得。4.3.3对 k15时索引情况的讨论事实

18、上，对于pow函数（指数函数）来讲，如果计算结果溢出，其输出将自动归为0，在这里，最多能计算pow（4,15） 。因此，计算机在运行时，当建立k15的索引时（哈希表），计算机并不会报错，甚至可以实现查找,不过这种查找可能不会是精确的。这种不精确性体现在下面两方面：1、由于pow函数溢出后输出为0引起的。2、由于在建立哈希表时使用了递推公式，可能会使数据溢出。4.3.4假设内存限制为8g，分析所设计索引方法所能支持的最大k值和相应数据查询效率。4.3.4.1改进后的最大k值分析我们可以扩展k的取值。我们定义从一个字符串中从第i个位置到第j个位置的字串为“字串（i，j）”显然，比较两个字符串是否相

19、等，等价于他们每个字串都分别向相等，特别地，也等价于”字串（1，i）“和”字串（i+1，k）“这两个字串分别相等。 于是，我们采用这样的方法，将字符串分成”字串（1,14）“、”字串（15,28）“.”字串（85，98）“、”字串（99,100）“这样的8个字串，即建立一共8个索引数组，然后根据不同的k值决定要使用几个索引数组。对每个待查的k-mer也采用相同的办法计算hash值，将其与索引数组进行对比。然后对最后的结果取交集，可以正确找出匹配项的位置。显然，对给定的k，一共需要（ ）个索引数组，因此，这样看来，k越大，对内存的需求也越大。最后，所需要的内存与k的关系可以由以下的公式给出：将n

20、（k）代入上式可以得3.8147*（101-k）+95.3674对上式求导数，得:故这是一个单调递增的函数，最大值在k=100处取得，于是 n*（k）max=n*（100）=793.4568mb4.3.4.2效率分析：为了扩展k的最大值，我们无疑是牺牲了效率，因为索引数组被建立了多次，但如果不是因为整形数据溢出，这些操作是完全没有必要的。如若原来时的效率为a，则当时的效率为。五、模型一般性测试1. 查找k=5，k-mer=atcga的dna序列程序运行图在原始fa文件中查找对比来验证正确性，选取其中2个位置查找验证：在文件1（solexa_100_170_1）中，第7条dna序列中read_1

21、70_7_1 random_genome_10000000 1211051 100 gaggactcagcccctcctagcgataggattgtatgtagccagcttggaacagaatatcgaggctaacgcgatgtcccgcactgagcccatggcactgtat在文件2（solexa_100_170_2）中，第499987条dna序列中（所有序列中的第999987条）read_170_499987_2 random_genome_10000000 7315838 100 aatccgattccttatctacgatatgctatcgatacgacaaattgtcaatca

22、aagggcgtggacgacctagtgcgtgaagtgcccatgaacgaaggttact1. 查找k=10，k-mer=atcgtgcata的dna序列程序运行图：在文件1（solexa_100_170_1）中，第12727条dna序列中read_170_12727_1 random_genome_10000000 8805313 100 gaaagggcccatcatcgtgcatatccgcaaaagaacccaacatagtacaaaattcgatgaggatagtccaaatgttatgaaccatactagaaaaagttat在文件2（solexa_100_170_2）中，

23、第486940条dna序列中（所有序列中的第986840条）read_170_486940_2 random_genome_10000000 9507192 100 tgaacgcatctttggaggaaatatgctccttaagagtggccataccatgacgtcccagcgtggaatacgatcgtgcataatgatcgtcctgatgtagtga3.查找k=15，k-mer=的dna序列在文件1（solexa_100_170_1）中，第21823条dna序列中read_170_21823_1 random_genome_10000000 3021063 100 ctgtggg

24、tagacggatacaactacaagatcggaatgcccttcgtgaggcacaggctgtttcagtgtgcggtcaacgattgaggattgccggcttcggaat在文件2（solexa_100_170_2）中，第407297条dna序列中（所有序列中的第907297条）read_170_407297_2 random_genome_10000000 3021120 100 ctggcgtcgcatagttaggctccgcggcaatactccttgcatgagtgtatgaagtacctgtgggtagacggatacaactacaagatcggaatgcccttcg

25、六 模型的优缺点及改进方向6.1模型的优缺点6.1.1模型的优点：优点：1.在程序算法中使用了小的数据类型，并且通过减少运算的强度使得索引建立时间短，检索速度快。2.算法使用索引所占内存空间小，空间复杂度低，因为只使用了o(1)。6.1.2缺点：1.该算法是基于4进制与10进制间的转换，当k14时，数据溢出，会覆盖掉其他内存地址的数据，导致检索缺失。参考文献1 郝甜甜 李强飞 李国治 陈禹翰 邓卫东 畜牧与饲料科学 2014年 第3期 云南农业大学动物科学技术学院 云南昆明 6502012 何建强 对karp-rabin串匹配随机算法的改进 广西科学院学报 第18卷第4期 2002年11月附录

26、实际使用的软件名称：visual c+ 6.0附录dna序列的k-mer index数据查找程序#include #include #include #include #define n 1000000double totaln100;char pn100;int main()printf(准备运行程序.n);int o,n; int i,j; char l100;o=1000100;n=105;file *fp1;file *fp2;if(fp1=fopen(d:bdatasolexa_100_170_1.fa , r)=null)/文件路径一定要正确，请将数据fa文件放置在d盘bdata文

27、件夹中。printf(打开文件失败，请确认文件所在路径正确);getchar();return(1);if(fp2=fopen(d:bdatasolexa_100_170_2.fa , r)=null) /文件路径一定要正确，请将数据fa文件放置在d盘bdata文件夹中。printf(打开文件失败，请确认文件所在路径正确);getchar();return(1);j=0;i=0;int k;char ch1,ch2;printf(- 读取文件中. -n);while( ch1!=eof )ch1=fgetc(fp1);switch(ch1)caset:k=0;break;casea:k=0;b

28、reak;caseg:k=0;break;casec:k=0;break;default:k=1;if(k=1)continue;if(i=100)i=0;j+;if(j=500000)break;pji=ch1 ;i+;fclose(fp1);while( ch2!=eof )ch2=fgetc(fp2);switch(ch2)caset:k=0;break;casea:k=0;break;caseg:k=0;break;casec:k=0;break;default:k=1;if(k=1)continue;if(i=100)i=0;j+;pji=ch2 ;i+;if(j=1000000)b

29、reak;fclose(fp2);int h,m;printf(请输入您想要查找的k值:);scanf(%d,&h);printf(- 正在建立索引. -n);clock_t begin, end;double cost;begin = clock();int f=1000100,g=105; int b,q; int a,r=0;double hash100,s100;for(q=0;q=1000000;q+)totalq0=0;for(b=0;bh;b+)switch(pqb)casea:hashb=0;break;caseg:hashb=2;break;caset:hashb=3;break;casec:hashb=1;break; sb=hashb;totalq0=totalq0+hashb*pow(4,(h-b-1);for(a=0;a100-h;a+,b+)switch(pqb)casea:hasha=0;break;caseg:hasha=2;break;caset:hasha=3;break;casec:hasha=1;break;sa+h=hasha;totalqa+1=hasha+4*(totalqa-sa*pow(4,(h-1);end = clock(