DNA序列的kmerindex问题数模论文.doc

重庆交通大学2015年第八届数学建模竞赛参 赛 论 文论 文 选 题 ：B 题学生姓名 学号 所在学院 联 系 电 话 ： 1 E-mail地 址 ： DNA 序列的k-mer index问题摘要 本小组在查阅了相关文献资料后，基于“数据结构”中的“哈希算法26”、“倒排索引12”法及“BKDRHash算法2”，建立相应的数学模型，给出分析和结果，对DNA 序列的k-mer index 问题给出解决方案。 本模型对不同k值采用不同算法建立索引。当k值较小时，利用基因序列其碱基种类较少（仅A,T,G,C四种）的特点，根据哈希算法进制转换的思想，可将k-mer 看成一个四进制的序列数，将其转化为十进制数作为哈希表的关键字2，并采用倒排索引的方法对哈希表关键字分类整理，建立相应的地址存储单元，实现索引；当k值较大时，考虑到内存溢出6的问题，采用“BKDRHash算法”对k-mer进行十进制转化，并结合“倒排索引2”法建立索引，从而对给定的 k-mer片段进行精确查找，最终输出碱基片段所在位置。 此方案将“哈希(Hash)算法”、“BKDRHash算法”和“倒排索引法”相结合，对哈希算法结构进行优化，提升了运算效率，操作简洁、高效。实现了在基因数据库34中对给定的碱基片段的位置进行查找的目的。关键词：倒排索引，哈希(Hash)算法，BKDRHash算法，碱基序列，基因数据库。一问题重述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，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值范围· 建立索引时间二.符号说明及假设1. 符号说明sumFile: 把记录分成sumFile个文件存放sumK-mer1: K-mer的总个数sumK-mer2: K-mer可能出现的不重复的个数maxY: 每个文件的最大行数maxX： 平均每行字符数fileName:文件名fileNumber;文件编号,范围是0 - sumFile-1hashKey1:将K_mer转化后的哈希关键字，是一个十进制整数position:记录K-mer所在的DNA序列及在每个序列中的位置HZJS:文件里平均每行的字符数&&：表示连接，如A && B即表示AB2.假设在题目给出的数据中，利用BKDRHash算法转化后不会出现相同的哈希值三问题分析 题目所给文件中有 100 万行的碱基序列，其中每行序列的长度为 100，给定一个固定的 k 值，则每行序列有（ 100-k+1）个k-mer，K-mer总数有1000000*(100-k+1)个。数据量级达到百万至千万，十分庞大，故考虑采用建立哈希表的方法实现索引。 碱基序列由A、T、C、G四种碱基无序组合，当给定K值后，理论有 种k-mer。比较1000000*（100-k+1）（实际值）和（理论值）的大小：当K小于等于13时， 1000000*（100-k+1）>，即K_mer值会出现重复，K越接近1，重复的越多；当K大于13时， 理论上不会出现重复值。四模型建立与算法分析3.1模型建立3.1.1 当时，即，则实际上出现的k-mer的个数多于理论上可能出现的k-mer的个数，即K-mer有重复。此时k较小，以“哈希算法”和“倒排索引”相结合的方法建立索引。 因为k-mer由四种碱基构成，根据进制转换思想1，将K-mer化为四进制再换算为十进制作为哈希表的关键字。现令碱基A->0，T->1，G->2，C->3，从而这四个碱基可以看成一个四进制的序列数，分别将 A,T,G,C赋值 0，1，2，3，可以根据四进制对十进制的转化方法（四进制化为十进制的方法为：假设取定一个序列 TCAGC，则以四进制13023表示，化为十进制：1*44+3*43+0*42+2*41+3*40=468）可以得到用十进制数表示的碱基序列（即468可以表示序列 TCAGC，为一个哈希关键字）。采用倒排索引的方法对碱基序列进行分类整理： 由于碱基序列号的最大值为1000000*(101-K)，将整数以字符形式保存，最大占用7个字符，每行的K-mer序列最大值为101-K,最大占用3个字符，即每个记录最大占位11个字符，假定一个字符占用一个字节，则索引记录共计大小为1000000*（101-K）*11*1个字节，即11*（101-K）M的字节大小。考虑索引记录保存在一个文件中，会导致文件太大，不利于索引操作，故将记录保存到1000个文件中。 将哈希关键字除1000取余，余数有0,1,2.999,共1000种，将余数作为文件编号fileNumber,文件名用fileNumber &&.tet来表示；商作为记录在每个文件里面的唯一编号即行号Y，同时采用哈希表记录属性的位置，因此在保存记录时将不记录属性值，属性记录采用“碱基序列号 && 每行的K-mer序列”的格式记录K_mer的位置；将记录保存在文件中。 当K_mer出现重复时，按其出现的先后顺序依次放在同一行。 另外，当k<5时，<1000，此时将记录保存到个文件中，较节约内存，其它过程同上。 当K值比较大时，如K取100，根据此种方法算出的哈希值，计算机将无法表示，内存溢出,这种情况参考下述方法。3.1.2 当K大于13时，此时1000000*(101-k)<4k，则实际上出现的K-mer个数少于理论上可能出现的k-mer个数，理论上k-mer序列不会出现重复，但考虑随机事件及内存溢出的问题，采用BKDRHash算法，计算k-mer的哈希值，作为哈希关键字。用“碱基序列号 && 每行的K-mer序列 && k-mer”的格式记录位置信息，其它过程同上。3.1.3检索 输入K_mer按前面的方法转化为不同的十进制整数，按照上面的方法即可确定对应的存储位置，到相应文件直接读出那一行数据即可；但当k的取值大于13，需对取出的数据进行筛选，选出正确的k-mer位置。3.1.4流程图5开始输入k值，建立索引输入k-mer,检索是否存在输出“无该k-mer值”输出位置结束 NY3.1.5建立索引（流程图）将数据写入文件并清除内存文件读取完毕，将最后一份数据写入文件索引创建完成开始输入k判断是否达到内存上限按行读取文件，将每行的k-mer取出后存入内存初始化文件，内存上限YN3.2模型模拟上图表示建立索引过程，程序显示正在读取数据文件的函数上图过程表示建立索引时文件的大小上图表示检索结果3.3分析3.3.1复杂度分析16索引的时间复杂为：O(n)索引时， 需要将待查序列的十进制数值与所有的 k-mer 的十进制数进行比较， 所以时间复杂度为： O(n)3.3.2查询效率及支持的最大K值支持100位的K-mer；该模型在建立索引上方法不同， K大于13时，需要筛选重复值，效率比K<13低3.3.3索引查询速度 K<13时，1秒内可以出结果；K>13时，时间稍长，2到3秒内可以出结果3.3.4索引内存使用2索引时只需要一个两个 int 型变量，内存消耗几乎可以不计3.3.5 G内存下，所能支持的k值范围1-1003.3.6建立索引时间因为要将文件的所有数据读一次，同时将K-mer值进行转化，在达到内存上限后还需将数据写入文件，故数据量越大时间消耗越多五模型评价模型优点1.将碱基序列转化为数值，节省空间；2.采用哈希算法倒排索引相结合的方法，查询速度远快于其他的算法且简单直接；3.索引速度与几乎与K值无关，速度快。模型缺点1. 当K值大于13，K取13、14等较小值时，出现重复的K_mer记录相对较多，既影响了建立索引的时间，又影响查找时间；2.K值较小时，平均到每行的字符数相对较多，则将缓存写入文件的次数增多，增加了时间消耗3. 此方法需在硬盘上创建文件，但在某些情况下用户并没有创建文件的权限，这将导致算法不能工作，且在硬盘上读写文件远没有在内存上操作的速度快，故建立索引和查找时间都受到影响。六参考文献1 吴孟达，成礼智，数据建模的理论与实践，北京出版社，1999 年8 月。2 白其峥，数据建模案例分析，高等教育出版社，2000年。3 张鑫鑫，生物序列数据 K-mer 频次统计与可视化研究， 第一期： 2014 年。4 陈波，何继凌，生物序列数据 K-mer 频次统计问题的算法， 第四期 ：2014 年。5 周建丽，大学计算机基础，中国铁道出版社，2012年8月6 严蔚敏，吴伟民，数据结构M，清华大学出版社，1998年。七 附录使用软件：sublime3.0源代码：<?php /建立索引环境为:ubuntu 1504 php 5.5.12-2,检索环境为lamp 在谷歌浏览器中显示,在值大于13时，建立索引刚开始很慢，后面速度会加快header("Content-Type:text/html;charset=utf8");ini_set("memory_limit","4096M");/试验设置为1G报错set_time_limit(0);$DNAObj = new DNA(8);/echo $DNAObj->strToHash('ATGC');/$DNAObj->run();/create index建立索引$DNAObj->indexK_low('AAGGCAAA');/index class DNA private $K;private $sumFile;private $arrfileNumber;private $HZJS;/平均每行的最大字符数private $startRAM;private $nowRAM;private $romLimit;function _destruct()function DNA($K)$this->startRAM = memory_get_usage();/获取初始内存占用情况/echo $this->startRAM.'<br>'$this->K = $K;switch ($this->K %5) case 0:$this->sumFile = pow(4,$this->K);break;default:$this->sumFile = 1000;break;$this->romLimit = 1024*1;/测试时设置为$this->HZJS= 100000;/function run()for ($i=0; $i < $this->sumFile; $i+) $this->arrfileNumber$i = 10;echo "开始"$this->initFile();if ($this->K<14 & $this->K >0) $this->K_low();elseif($this->K<101)$this->K_high();elseecho "输入出错"function K_low()/值取指较小时建立索引$num = 1;$isWrite = 1;$resFile = fopen("/var/www/html/01.fa", "r");for ($FN=0; $FN < 2; $FN+) if ($resFile) $hang = 1; while (!feof($resFile) if ($num = 1) $num = 2; $DNAstr = fgets($resFile, 150); else $DNAstr = fgets($resFile, 150); for($i=1;$i<=(101-$this->K);$i+) $k_mer = substr($DNAstr, $i,$this->K); $hashKey1 = $this->strToHash($k_mer); $fileNumber = $hashKey1%$this->sumFile; $Y = floor($hashKey1/$this->sumFile); if ($Y >$this->arrfileNumber$fileNumber) $this->fileAdd($fileNumber,$Y-$this->arrfileNumber$fileNumber); $this->arrfileNumber$fileNumber = $Y; $position = ' '.$hang.','.$i; if (isset($arrJL$fileNumber$Y) $arrJL$fileNumber$Y = $arrJL$fileNumber$Y.$position; else $arrJL$fileNumber$Y = $position; $this->nowRAM = memory_get_usage();/获取当前内存占用情况 $ram = ($this->nowRAM - $this->startRAM )/10224/2;/ if($ram > $this->romLimit)/试验时设为1 for ($i=0; $i < $this->sumFile; $i+) $fileName = '/var/www/html/DNAfile/'.$i.'.txt'$NewFile = fopen($fileName,'r');$j = 0;while (!feof($NewFile)if (isset($arrJL$i$j) $arrLS$j = fgets($NewFile, $this->HZJS).$arrJL$i$j;unset($arrJL$i$j);else$arrLS$j = fgets($NewFile, $this->HZJS);$j+;fclose($NewFile);file_put_contents($fileName,"");$NewFile = fopen($fileName,'w');foreach ($arrLS as $key => $value) fputs($NewFile,$value);unset($value);fclose($NewFile); $hang+; $num = 1; echo "2<br>" fclose($resFile); $resFile = fopen("/var/www/html/02.fa", "r"); fclose($resFile); echo "建立索引完成<br>"function K_high()/值取指较大时建立索引$num = 1;$resFile = fopen("/var/www/html/01.fa", "r");for ($FN=0; $FN < 2; $FN+) if ($resFile) $hang = 1; while (!feof($resFile) if ($num = 1) $num = 2; $DNAstr = fgets($resFile, 150); else $DNAstr = fgets($resFile, 150); for($i=1;$i<=(101-$this->K);$i+) $k_mer = substr($DNAstr, $i,$this->K); $hashKey1 = $this->MYHash($k_mer); $fileNumber = $hashKey1%$this->sumFile; $Y = floor($hashKey1/$this->sumFile); if ($Y >$this->arrfileNumber$fileNumber) $this->fileAdd($fileNumber,$Y-$this->arrfileNumber$fileNumber); $this->arrfileNumber$fileNumber = $Y; $position = ' '.$hang.','.$i; if (isset($arrJL$fileNumber$Y) $arrJL$fileNumber$Y = $arrJL$fileNumber$Y.$position.$k_mer; else $arrJL$fileNumber$Y = $position; $this->nowRAM = memory_get_usage();/获取当前内存占用情况 $ram = ($this->nowRAM - $this->startRAM )/10224/2;/ if($ram > $this->romLimit) for ($i=0; $i < $this->sumFile; $i+) $fileName = '/var/www/html/DNAfile/'.$i.'.txt'$NewFile = fopen($fileName,'r');$j = 0;while (!feof($NewFile)$tmpstr2 = fgets($NewFile, $this->HZJS);if (isset($arrJL$i$j) $arrLS$j = $tmpstr2.$arrJL$i$j;unset($arrJL$i$j);else$arrLS$j = $tmpstr2;$j+;fclose($NewFile);$NewFile = fopen($fileName,'w');for ($j=0; $j < $this->sumFile; $j+) if (isset($arrLS$j) fputs($NewFile,$arrLS$j);unset($arrLS$j);elsefclose($NewFile); $hang+; $num = 1; fclose($resFile); $resFile = fopen("/var/www/html/02.fa", "r"); fclose($resFile); echo "建立索引完成<br>"function strToHash($K_mer)$hash = 0;$len = strlen($K_mer);for ($i=0; $i < $len; $i+) $ch = substr($K_mer,$i,1);switch ($ch) case 'A':$hash = $hash+0*pow(4, $len-1-$i);break;case 'T':$hash = $hash+1*pow(4, $len-1-$i);break;case 'G':$hash = $hash+2*pow(4, $len-1-$i);break;case 'C':$hash = $hash+3*pow(4, $len-1-$i);break;default:break;return $hash;function MYHash($str)/BKDRHash哈希算法$seed = 131; / 31 131 1313 13131 131313 etc. $hash = 0; $cnt = strlen($str); for($i = 0; $i < $cnt; $i+) $hash = (floatval($hash * $seed) & 0x7FFFFFFF) + ord($str$i) & 0x7FFFFFFF; return floor($hash & 0x7FFFFFFF)%100000019);/将原算法改进以下，降低值的大小function initFile()for ($i=0; $i < $this->sumFile; $i+) $fileName = '/var/www/html/DNAfile/'.$i.'.txt'$NewFile = fopen($fileName,'w');/初始化文件每行先预留行for ($j=0; $j < 10-1; $j+) /fputs($NewFile,$j);fputs($NewFile,"n");fclose($NewFile);function fileAdd($fileNumber,$add)/在运行中增加行数$fileName = '/var/www/html/DNAfile/'.$fileNumber.'.txt'$NewFile = fopen($fileName,'a');for ($j=0; $j < $add; $j+) /fputs($NewFile,$j);fputs($NewFile,"n");fclose($NewFile);function indexK_low($str)$hashKey1 = $this->strToHash($str);$fileNumber = $hashKey1%$this->sumFile;$Y = floor($hashKey1/$this->sumFile);$fileName = '/var/www/html/DNAfile/'.$fileNumber.'.txt'$NewFile = fopen($fileName,'r');$j = 0;echo $str;echo "<br>"while (!feof($NewFile)$position = fgets($NewFile, $this->HZJS);if ($j=$Y) echo $position.''$j+;fclose($NewFile);function indexK_high($str)$hashKey1 = $this->MYHash($str);$fileNumber = $hashKey1%$this->sumFile;$Y = floor($hashKey1/$this->sumFile);$fileName = '/var/www/html/DNAfile/'.$fileNumber.'.txt'echo $str;echo "<br>"$NewFile = fopen($fileName,'r');for ($j=0; $j <= $Y; $j+) $position = fgets($NewFile, $this->HZJS);fclose($NewFile);$len = strlen($str);preg_match_all("/s+/s",$position,$arr);/ 这里除了匹配 空格，还匹配中文全角的空格 s后面直接加上就是了$arr = $arr0;if ( count($arr) = 1) echo $arr0;exit();foreach ($arr as $key => $value) if (strlen($value) > 20) /判断为含有-mer信息$str2 = substr($value, strlen($value) - $len, $len);if(strcmp($str, $str2) = 0)echo $str2;echo " <br>"elseecho $arr0;?>