1、重庆交通大学 2015 年第八届数学建模竞赛参 赛 论 文论 文 选 题 ：B 题学生姓名 学号 所在学院联 系 电 话 ： 1 E-mail 地 址 ： 0DNA 序列的 k-mer index 问题摘要本小组在查阅了相关文献资料后，基于“数据结构”中的“哈希算法 26”、“倒排索引 12”法及“BKDRHash 算法 2”，建立相应的数学模型，给出分析和结果，对 DNA 序列的 k-mer index 问题给出解决方案。本模型对不同 k 值采用不同算法建立索引。当 k 值较小时，利用基因序列其碱基种类较少（仅 A,T,G,C 四种）的特点，根据哈希算法进制转换的思想，可将 k-mer 看成
2、一个四进制的序列数，将其转化为十进制数作为哈希表的关键字 2，并采用倒排索引的方法对哈希表关键字分类整理，建立相应的地址存储单元，实现索引；当 k 值较大时，考虑到内存溢出 6的问题，采用“BKDRHash 算法”对 k-mer 进行十进制转化，并结合“倒排索引 2”法建立索引，从而对给定的 k-mer 片段进行精确查找，最终输出碱基片段所在位置。此方案将“哈希(Hash)算法”、“BKDRHash 算法”和“倒排索引法”相结合，对哈希算法结构进行优化，提升了运算效率，操作简洁、高效。实现了在基因数据库 34中对给定的碱基片段的位置进行查找的目的。关键词：倒排索引，哈希(Hash)算法，BKD
3、RHash 算法，碱基序列，基因数据库。1一问题重述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通常这些
4、 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、出使用索引查询的计算复杂度，和空间复杂度分析。（5）假设内存限制为 8G，分析所设计索引方法所能支持的最大 k 值和相应数据查询效率。（6）按重要性由高到低排列，将依据以下几点，来评价索引方法性能 索引查询速度 索引内存使用 8G 内存下，所能支持的 k 值范围 建立索引时间2二.符号说明及假设1. 符号说明sumFile: 把记录分成 sumFile 个文件存放sumK-mer1: K-mer 的总个数sumK-mer2: K-mer 可能出现的不重复的个数maxY: 每个文件的最大行数maxX： 平均每行字符数fileName:文件名fileNumber;文件编号,范围是 0 - sumF
6、ile-1hashKey1:将 K_mer 转化后的哈希关键字，是一个十进制整数position:记录 K-mer 所在的 DNA 序列及在每个序列中的位置HZJS:文件里平均每行的字符数/$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
7、;function _destruct()function DNA($K)$this-startRAM = memory_get_usage();/获取初始内存占用情况11/echo $this-startRAM.;$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 s
8、umFile; $i+) $this-arrfileNumber$i = 10;echo “开始“;$this-initFile();if ($this-KK 0) $this-K_low();elseif($this-KK_high();elseecho “输入出错“;function K_low()/值取指较小时建立索引$num = 1;$isWrite = 1;$resFile = fopen(“/var/www/html/01.fa“, “r“);for ($FN=0; $FN K);$i+)$k_mer = substr($DNAstr, $i,$this-K);$hashKey1
9、= $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$fi
10、leNumber$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)/试验时设为 1for ($i=0; $i sumFile; $i+)$fileName = /var/www/html/DNAfile/.$i.txt;$NewFile = fopen($fileName,
11、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);13$j+;fclose($NewFile);file_put_contents($fileName,“);$NewFile = fopen($fileName,w);foreach ($arrLS as $key = $value) fputs($NewFile,$v
12、alue);unset($value);fclose($NewFile);$hang+;$num = 1;echo “2“;fclose($resFile);$resFile = fopen(“/var/www/html/02.fa“, “r“);fclose($resFile);echo “建立索引完成“;function K_high()/值取指较大时建立索引$num = 1;$resFile = fopen(“/var/www/html/01.fa“, “r“);for ($FN=0; $FN K);$i+)$k_mer = substr($DNAstr, $i,$this-K);14$
13、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) $ar
14、rJL$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 sumFile; $i+)$fileName = /var/www/html/DNAfile/.$i.txt;$NewFile = fopen($file
15、Name,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 sumFile; $j+) if (isset($arrLS$j) fputs($NewFile,$arrLS$j);15uns
16、et($arrLS$j);elsefclose($NewFile);$hang+;$num = 1;fclose($resFile);$resFile = fopen(“/var/www/html/02.fa“, “r“); fclose($resFile);echo “建立索引完成“;function strToHash($K_mer)$hash = 0;$len = strlen($K_mer);for ($i=0; $i sumFile; $i+) $fileName = /var/www/html/DNAfile/.$i.txt;$NewFile = fopen($fileName,w
17、);/初始化文件每行先预留行for ($j=0; $j strToHash($str);$fileNumber = $hashKey1%$this-sumFile;$Y = floor($hashKey1/$this-sumFile);$fileName = /var/www/html/DNAfile/.$fileNumber.txt;17$NewFile = fopen($fileName,r);$j = 0;echo $str;echo “;while (!feof($NewFile)$position = fgets($NewFile, $this-HZJS);if ($j=$Y) ec
18、ho $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 “;$NewFile = fopen($fileName,r);for ($j=0; $j HZJS);fclose($NewFile)
19、;$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 “ “;18elseecho $arr0;?
