2015年深圳杯B题-DNA 序列的k-mer index 问题.doc

封一答卷编号（参赛学校填写）：答卷编号（竞赛组委会填写）：论文题目： b题 组 别： 本科生 参赛队员信息(必填)： 姓 名年级专业学 号联系电话参赛队员114机制二班参赛队员214机制二班参赛队员314机制二班 指导教师：_ 参赛学校： 沈阳农业大学 封二答卷编号（参赛学校填写）：答卷编号（竞赛组委会填写）：评阅情况（学校评阅专家填写）：学校评阅1.学校评阅2.学校评阅3_ _评阅情况（联赛评阅专家填写）：联赛评阅1.联赛评阅2.联赛评阅3.dna 序列的k-mer index 问题摘要： 随着生物技术和计算机技术交叉，合作的范围不断地扩大，深度不断地加深，字符串匹配在生物科学的信息检索领域有着广泛的应用，生物研究对dna序列的碱基信息快速搜索的要求也与日俱增，dna 序列的k-mer index 问题研究日益突出。本文在查阅了相关文献资料后，基于“数据结构”中的“sunday算法7”，我们给出了一种函数覆盖优化sunday算法，对dna 序列的k-mer index 问题中给出固定值k ，进行碱基片段查找。函数覆盖优化sunday算法对所给出的大量数据的碱基序列先运用函数覆盖，同时把dna序列进行分组预处理。当数据链接时用函数覆盖，将整个dna序列数据划分成若干个小的区域，然后再针对若干个小区域进行碱基片段的匹配。在划分小区域时，首先将每个小区域的末端多划分四个属于下个区域的碱基，其次在匹配的过程中需要按给定k值对小区域进行分组，之后再用算法对碱基片段进行匹配，最后再统计所有小区域与所要查找的碱基片段的匹配信息。匹配时经过优化后的算法能够实现同时对划分的所有区域进行查找，直到与所有dna序列都进行了匹配。当匹配成功时，首先会自动将char *格式强制转换成void *格式，然后通过dna碱基片段中第一个字符的指针访问dna序列，最后输出dna碱基片段中第一个字符的（常量）地址。 函数覆盖优化sunday算法后的主要特点：查找速度简洁、高效、执行速度快。函数覆盖优化sunday算法，因为快速跳跃匹配和多区域同步搜索的特征使它可以快速的生成数据和信息。实现了对dna序列数据库的查找，在dna序列中可以得到固定k 长度的碱基片段的位置。关键词：sunday算法，函数覆盖，字符串匹配，碱基序列，数据库。1问题重述给定一个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 值范围，建立索引时间。2符号说明iostream类(文件输入流和输出流类)2a：构造方法： #include创建文件输入流类对象和已存在的文件相关联。文件不存在的话，并创建。如：#include(c:usersdengxindesktopdna)；b：#include 创建文件输出流类对象和已存在的文件相关联，并设置该该流对文件的操作是否为续写。如：#include(c:usersdengxindesktopdna，at+)； c: 表示在 #include 对文件再次写入时，会在该文件的结尾续写，并不会覆盖掉。主要方法：#include写入字符串。当执行完此方法后，字符数据还并没有写入到目的文件中去。此时字符数据会保存在缓冲区中。此时在使用memcpy函数就可以使数据将目标数组地址增加到要追加数据的地址，从而保存到目的文件中去。d:定义函数：void exit(int status)exit()用来正常终结目前进程的执行， 并把参数status 返回给父进程， 而进程所有的缓冲区数据会自动写回并关闭未关闭的文件.。在关闭后，再写入或者刷新的话，会显示eof异常。字符串读写函数a：定义一个头文件，此文件是输入类对象，头文件是关联于源文件。 #include= #includeb：主要方法读字符函数fgets(字符数组名，n，文件指针)； 其中的n是一个正整数。表示从文件中读出的字符串不超过 n-1个字符。在读入的最后一个字符后加上串结束标志0。例如：fgets(str，n，fp)；的意义是从fp所指的文件中读出n-1个字符送入字符数组str中。写字符函数fputc(字符量， 文件指针 )；fputc函数的功能是把一个字符写入指定的文件中。其中，待写入的字符量可以是字符常量或变量。例如：fputc(a，fp)；其意义是把字符a写入fp所指向的文件中。exit()用来正常终结目前进程的执行， 释放与之关联的所有资源。数据块读写函数fread和fwrite；c语言还提供了用于整块数据的读写函数。可用来读写一组数据，如一个数组元素，一个结构变量的值等。 读数据块函数fread(buffer，size，count，fp)；6 写数据块函数fwrite(buffer，size，count，fp)；其中： buffer：是一个指针，在fread函数中，它表示存放输入数据的首地址。在fwrite函数中，它表示存放输出数据的首地址。 size：表示数据块的字节数。 count：表示要读写的数据块块数。 fp：表示文件指针。假设在发生碱基不匹配时sitj，1in，1jm。此时已经匹配的碱基片段为u，并假设字符串u的长度为k，如图1。8明显的，sl+i+1肯定要参加下一轮的匹配，并且tm至少要移动到这个位置(即模式串t至少向右移动一个字符的位置)。每组dna碱基序列中，当第一个k长度中没有子串（碱基片段）出现时，这个时候碱基片段的第一个碱基移动到k长度碱基之后的碱基的位置。图1 sunday算法不匹配的情况图1 分如下两种情况： (1) sl+i+1在模式串t中没有出现。这个时候模式串t0移动到sl+i+1之后的字符的位置，如图2。图2 sunday算法移动的第1种情况图2 (2)sl+i+1在模式串中出现。这里sl+i+1从模式串t的右侧，即按tm-1、tm-2、t0的次序查找。如果发现sl+i+1和t中的某个碱基相同，则记下这个位置，记为k，1km，且tk=sl+i+1。此时，应该把模式串t向右移动m-k个字符的位置，即移动到tk和sl+i+1对齐的位置，如图3。图3 图3 sunday算法移动的第2种情况 依次类推，如果碱基完全匹配了，则碱基片段匹配成功；否则，再进行下一轮的移动，直到主串s的最右端结束。对于碱基片段的匹配问题，该算法执行速度较快。 3模型建立与求解3.1函数覆盖分组 首先根据题目所给的文件中有10 0 万行的碱基序列，1其中每行序列的长度为10 0，我们用二维数组a 1 0 0 0 0 0 0 1 0 0 存储所给的碱基序列， 函数覆盖对dna序列进行分组，用指针数组将整个dna序列分成若干小区域，我们要做的就是取定一个固定的k 值，将每行序列分成10 0 - k + 1 个长度为k 的序列， 查询某序列在碱基序列中的编号和位置。3.2算法匹配 sunday算法是从前往后匹配，4在匹配失败时关注的是碱基片段中参加匹配的最末位碱基的下一位碱基，如果该碱基没有在匹配串中出现则直接跳过，即移动步长= 匹配串长度+ 1，如果dna序列中匹配字符串的右侧一个碱基没在子串中，碱基片段移动步长=整个k长度碱基的距离+1，如果dna序列中匹配范围内的右侧一个字符在子串中，碱基片段移动距离=子串长度-这个字符在子串中的位置。3.3输入k，s将输入的特定的k 个碱基片段在10 0 万行分成不同几个小序组中，5以跳跃匹配的形式进行匹配。3.4位置输出从左向右匹配过程中，当遇到不匹配的时候，看dan序列中匹配范围之外的右侧第一个字符在子串中的最右位置 ，根据事先计算好的移动步长移动dan序列指针，直到匹配，程序会输出所有要查找碱基片段的位置。（如果遇到重复的碱基片段和查找的长度碱基片段相同的情况，我们在文章后面进行了分析和处理。） 用函数覆盖优化sunday算法来解决dna 序列的k-mer index问题,数据已提前加载，划分成若干个小区域，输入固定值k，碱基序s即可匹配碱基片段 输入k， 碱基序s 开始将给定序列所有k-mer 列出再将所有dna序列用函数覆盖分组 用sunday算法对划分的小区域进行同步搜索匹配成功，强制char *转换成void *，输出碱基片段地址（行号，列数） 结束匹配不成功，无此碱基序列图43.5程序流程 dna序列数据.fasunday算法匹配中，字符串是以空终止符（0）结尾的字符数组，k-mer 当通过dna碱基片段中第一个字符的指针访问dna序列，输出的dna碱基片段的值是字符串中第一个字符的（常量）地址图5建立索引（数据流图）：匹配完所有dna序列后，对匹配成功的碱基片段输出相应地址。3.6模型模拟3.61用函数覆盖优化sunday算法对数据进行检验，先在100行序列dna序列中试运行此算法，当k=5，序列s是aacag时。图6检验序列aacag的结果：图7图83.62 用优化后的算法对所有dna进行检索，在100wdna序列划分成的若干个小区域同时搜索，当k=10时，程序运行后的结果和生成的文件。3.63给出固定值，序列后，对所有dna序列数据进行匹配后，运行出的文件。图9由于k-mer【10】长度如果过长，导致字符串没有被分配的空间，会造成内存不够溢出的现象。3.64复杂度分析索引的时间复杂为：o(n*m)索引时，需要将准备搜索的k长度碱基片段与所有的dna序列进行匹配，所以该算法最坏情况下的时间复杂度为o(n*m)。 4模型检验4.1 贪婪算法的概念贪婪算法是一种对某些求最优解问题的更简单、更迅速的设计技术。3用贪婪法设计算法的特点是一步步地进行，常以当前情况为基础根据某个优化测度作最优选择，而不考虑各种可能的整体情况，它省去了为找最优解要穷尽所有可能而必须消耗的大量时间，它采用自顶向下，以迭代的方法做出相继的贪心选择，每做一次贪心选择就将所求问题简化为一个规模更小的子问题，通过每一步贪心选择，可得到问题的一个最优解。4.2 贪婪算法的求解 开始链接数据库，用函数覆盖分组，k-mer序列读取匹配串加载，固定值长度的模式串将模式串与匹配串匹配，如果成功循环的在各组中进行匹配，找到目标串用强制转换输出匹配成功的碱基片段的位置 先把碱基片段（子串）与dna序列（字符串）左边对齐，当在第二个碱基处发现不匹配，要把子串往后移动，第二个碱基是必定要参加下一步的比较的，需要移动子串，使子串中的最右边的这个碱基与字符串对齐。现在k长度的子串不存在，则说明可以直接跳过一大片，从与第二个碱基一样的碱基之后的那个碱基开始作下一步的比较，直到第k个碱基匹配成功。查找出这组中完全匹配成功的碱基片段，如果没有则继续向下一组中循环查找，直到查找到它的位置，再用强制转换输出它的地址。流程图如下：图10贪婪算法流程图5模型评价5.1 模型的优缺点5.1.1 优点 1. 函数覆盖优化sunday算法后的方法通过对数据先进行划分，划分成若干个小区域后再进行同步搜索，查询速度远快于sunday算法并且简洁有效。 2. 函数覆盖使用时，因为dna序列地址在基类中被说明为虚函数，11那么我们所查找的碱基片段也自动成为虚函数，查找时函数覆盖会根据指针所指向的碱基调用对应的虚函数.实现动态联编。 3. 使用函数覆盖时指针的类型是可以强行转换的，便于输出匹配碱基的地址。5.1.2 缺点 虽然对sunday算法用函数覆盖进行了优化，但是当查找模式串越长时，此算法的下一个字符在模式串中出现的概率就越大，无效匹配的次数就越多，导致内存的浪费。 6参考文献1.严蔚敏，吴伟民；数据结构m；清华大学出版社；1998.2.谭浩强；程序设计m；清华大学出版社；1991.3.白其峥；数据建模案例分析；高等教育出版社；2000.4.朱永强，秦志光，江雪；基于sunday算法的改良单模式匹配算法计算机应用，2014， 34(1)；doi：10.11772/j.issn.1001-9081.2014.01.0208.5.巫喜红；分块法的模式匹配算法的研究；重庆邮电大学学报（自然科学版）；2014， 26(4)；doi：10.3979/j.issn.1673-825x.2014.04.022.6.袁晓洲；c语言中函数之间地址传递方式浅析；2014， (20)；journal：computer cd software and applications.7.潘冠桦，张兴忠；sunday算法效率分析；计算机应用；2012， 32(11)；doi：10.3724/sp.j.1087.2012.03082.8.鲁宏伟，魏凯，孔华锋；一种改进的kmp高效模式匹配算法；华中科技大学学报（自然科学版）；2006，34(10)；i10.3321/j.issn:1671-4512.2006.10.013.9.吴孟达，成礼智；数据建模的理论与实践；北京出版社1999 年8 月.10.张鑫鑫，陈波，何继凌，徐云；生物序列数据k-mer频次统计问题的算法；计算机系统应用；2014， (4)；doi：10.3969/j.issn.1003-3254.2014.04.024.11.夏承遗，董玉涛，赵德新，唐树刚，陈国章，张桦；c + +中虚函数的实现机制；天津理工学院学报；2004，20(3)；doi：10.3969/j.issn.1673-095x.2004.03.017.7附录源代码：#include #include #include #include cspreadsheet.h using std:string; #pragma warning(disable:4146) #pragma warning(disable:4786) #import c:usersdengxindesktopdna(eof,adoeof) /插入到数据库 bool insertexcel(cstring str1,cstring str2) try cdatabase m_db; if (!m_db.isopen() m_db.openex(c:usersdengxindesktopdna新建文件夹,0); cstring sql(insert into students(myname,age) values(+ str1+,+str2+); m_db.executesql(sql); if(m_db.isopen() m_db.close(); return true; catch(_com_error e) string errormessage(数据库连接关闭失败:),description,source; description=e.description(); source=e.source(); errormessage+=e.errormessage(); errormessage=errormessage+rn+source+rn+description; :messagebox(null,errormessage.c_str(),错误,mb_ok); return false; /获取路径 cstring getaddr() cstring sfile,spath; /获取主程序所在路径,存在spath中 / 检索已安装的驱动是否有excel. do if (strstr(pszbuf,excel) != 0) /发现 ! sdriver = cstring(pszbuf); break; pszbuf = strchr(pszbuf, 0) + 1; while (pszbuf1 != 0); return sdriver; /读取excel void readfromexcel() try cstring str=getaddr(); if(str.isempty() :messagebox(null,无法获取当前路径,null,mb_ok); else cspreadsheet ss(str,students); cstringarray rows, column; cstring strcontents = ; cstring sitem3=0; for (int i = 2; i ss.gettotalrows()+1; i+) / 读取一行 ss.readrow(rows, i); strcontents.empty(); for (int j = 0; j m_strerror); end_catch; /写excel void writefromexcel(int num,cstring str1,cstring str2,cstring str3) cstring path=getaddr(); if(path.isempty() :messagebox(null,获取路径错误,null,mb_ok); else / 新建excel文件名及路径，testsheet为内部表名 bool selectexcel() _connectionptr m_pconnection; /connection objects pointer _commandptr m_pcommand; /command objects pointer _parameterptr m_pparameter; /parameter objects pointer _recordsetptr m_precordset; hresult hr; try / 创建连接对象 hr=m_pconnection.createinstance(_uuidof(connection); m_precordset.createinstance(_uuidof(recordset); if(!succeeded(hr) return false; / 连接数据库 while (!m_precordset-adoeof) #include #include #include / 预处理子串 void kmp_prepare(char *target, char *pattern, int *overlay) memset(overlay, 0, sizeof(overlay); int i = -1, j = 0, preoverlay; overlay0 = -1; for(i=1; i= 0 & patterni != patternpreoverlay+1) preoverlay = overlaypreoverlay; if(patterni = patternpreoverlay+1) overlayi = preoverlay + 1; else overlayi = -1; for(i=0; istrlen(pattern); i+) printf(overlay%d: %dn, i, overlayi); / 查询子串 public class sunday public sunday() / / todo: 在此处添加构造函数逻辑 / public int qfindchr(string str, string sfind) int str_length = 0; int fin_length = 0; int find_count = 0; int start = 0; int movenum = 0; if (str.length sfind.length) return find_count; str_length = str.length; fin_length = sfind.length; while (start + fin_length = str_length) movenum+; bool isfind = false;/是否在这次移动中找到 string s_temp = str.substring(start, fin_length); if (s_temp = sfind) find_count+; start = start + fin_length; isfind = true; if (isfind = false)/如果没找到计算下次移动位置 int forwardpos = qfindpos(str, sfind, start, fin_length); start = forwardpos; return find_count; /找字符在字符串(不算最后一个字符)的位置(倒数)