基因组组装 (1.doc

基因组组装摘要快速和准确地获取生物体的遗传信息对生命科学研究具有重要的意义。本文建立了基因组组装的数学模型及算法，解决了基因组的组装问题。如何在保证组装序列的连续性、完整性和准确性的同时，设计出耗时短、内存小的组装算法是本文的关键。针对问题一，本文主要是对重叠群（congtig)拼接问题进行了数学模型的建立，设计出一种新的重叠群（congtig）生成算法。该算法基于deBruijn图，将从多头测序转化为deBruijn图的欧拉路径问题，并采用启发式搜索，能快速的处理大量测序数据，而且能得到质量较高的重叠群，提高了基因组组装的准确性。针对问题二，由于数据量庞大，所以本文先对原始数据进行分析与处理。首先对原始数据中重复片段进行过滤；另外，通过对原始数据中质量序列进行分析和研究，剔除质量偏差较大错误数据，使有效的读长进一步减少。其次，利用问题一的deBruijn图法对读长（read）进行拼接；最后，本文依据组装序列的连续性、完整性和准确性对组装效果进行评价，得到组装序列的匹配度和准确度都较高，达到较为理想的组装效果。关键词：基因组组装 deBruijn图 congtig拼接 贪婪图方法 启发式搜索一 问题重述快速和准确的获取生物体的遗传信息对生命科学研究具有重要意义。随着测序技术的不断发展，新一代测序技术产生的在高通量、低成本的同时也带来了错误率略有则加、读长较短等缺点。本题要求利用数学模型，设计算法解决如下几个问题：（1）测序过程中可能出现的个别碱基对识别错误；（2）基因组中存在重复片段；（3）快速的处理海量的序列比对。二 问题分析问题一是在新一代测序技术的测序策略的基础上，分析了基因组组装所面临的主要挑战，要求设计算法解决新一代测序技术带来的一些弊端。1.新一代测序技术所得的reads长度较短，数量较多，不易发现reads之间的重叠关系。因此我们可以将reads转化成定长的k-mer,然后寻找k-mer之间的重叠关系。最后建立de Bruijin图，把短序列拼接问题转化为de Bruijin图中的欧拉路径问题。2.个别碱基对识别错误，我们可以通过将多个reads放在一起对比来发现错误。3.基因组中存在重复片段，重复片段可能导致拼接错误，或者导致不连续的较短contig出现。重叠片段类型主要有以下几种， 重复片段问题可以用如下问题解决：通过对比，可先将重复片段隔离开来，较高的覆盖度有利于重复片段的隔离，但是，较多的测序错误将不利于该过程的进行。如果重复片段比read 长，可利用pared end read 来解决；如果重复片段比read 短，那么该read又被称为 spanner，一个spanner 就是一个重复片段两端再加几个碱基组成。利用spanner 解决重复片段问题需要如下两个信息：一是重复片段两端配对的read ，这两个read 必须不相同；二是重复片段中的一个配对read ，只要知道一个即可，另一个配对read 可以不在重复片段中。通过分析已知的基因组，可获得有关重复片段的更多信息，如，重复片段的长度，重复片段的模式等。问题二是提供全长约为120,000个碱基对的细菌人工染色体，采用新一代的Hiseq2000测序仪进行测序。附件提供了筛选好的定长read数据文件。使用第一题设计的基于de Bruijn图的组装算法对read数据进行组装，并对结果进行误差分析。三 模型假设与符号说明3.1 模型假设3.2 符号说明符 号定 义read利用现有的测序技术，可按一定的测序策略获得长度约为50100个碱基对的序列，称为读长；contig（C）由read经过一定算法拼接产生3kb10Mb以内的一些基因组片段;supercontig（S）使用contig作为参考序列延伸，并进行合并得到更长的contig,即supercontig;quality（Q）根据本题数据，每一个read都含有一个质量值，该值能反映该read的正确率。质量值越高，read的正确率越高。k-mer长度为k的一段DNA片段四 数据分析与模型原理4.1 数据分析本题中，采用HiSeq2000测序技术产生的数据。HiSeq2000是目前通量最高的测序仪器，但产生的读长较短一般为100bp(本题read长度为88bp)，使拼接问题变得更加复杂。HiSeq2000测序仪测出的数据有如下特征：（1）read的副本较多，约为50-100；（2）基因组中有些位置被较多的read所覆盖，有些位置被较少的read覆盖，这 些位置是随机的，不可预测；（3） 每一个read都含有一个质量值，该值能反映该read的正确率。质量值越高， read的正确率越高；（4）有些read上存在个别碱基对识别错误，无法在比对前甄别出来。4.2模型原理由于该问题比较复杂，直接由read得到整个基因组想到困难，甚至是不可能的，故需要降温题划分为几个子问题，分如下几个阶段解决：（1）由read集合构建contig集合；（2）由contig集合构建supercontig集合；（3）由supercontig集合构建整个基因组。supercontig（更长）contig（100bp）read（50-100bp） 把拼接问题转化成de Bruijn图中的欧拉路径问题，不断迭代得到尽可能长的序列 图1 4.2.1初始序列筛选根据本题提供的数据，本文通过两个步骤的操作筛选序列集合：操作1：设G是一个序列集合，成为基因组。由HiSeq2000测序仪处理形成一些小的片段，每个片段长度在50-100bp。接着去掉一些稍长的和稍短的片段，将剩余的片段切成定长的序列（本题提供定长为88bp）保存在多重集合U中，称U中的序列为read。本题操作1 的具体内容不需要考虑。操作2：找一个U的子集S，满足下列条件：（1）给定阈值quality，S中的read上的碱基质量平均值不小于quality；（2）S中所有read都是成功的；（3）S中的read要尽可能的多。4.2.2 congtig拼接总体思路将read转化成定长的k-mer，并将这些k-mer存入de Bruijn图中，以备之后查找使用。此时要设定的一个重要参数是k-mer的长度。选定值之后，要将长度的read拆成个k-mer。根据一定策略，选定一个初始k-mer，接下来就可以在该k-mer为结点开始搜索后继的k-mer。搜索时采用贪婪图策略，每一步选择在当时看来最优的后继k-mer，直到满足事先设定的终止条件，结束一条contig的拼接，接着开始下一条contig的拼接。直到没有合适的初始k-mer可供选择，整个拼接过程结束。五 模型建立与求解5.1 模型一的建立与求解5.1.1 模型一的建立根据基因的测序过程，建立了如下模型：设 G 是一个字符串，称之为基因 组。取 G 的 10 个复本，然后将每个复本在随机的位置断开，形成一些小的字符串片段，每个片段长度小于 100 个字符。接着扔掉一些稍长的和稍短的字符串片段，将剩余字符串片段切成定长的字符串保存在多重集 U 中，称 U 中的字符串为 read。字符串上的每个字符都来自字符集A，C，G，T，每个字符就是一个碱基。这些字符串长度都是，考虑到新一代测序技术，的范围是 。同时，每个字符串上的每个字符都有一个正确率，并且这些字符串能完整的覆盖原始基因组 G。1. congtig拼接步骤： 将read转化成定长的k-mer，并将这些k-mer存入de Bruijn图中，以备之后查找使用。此时要设定的一个重要参数是k-mer的长度。选定值之后，要将长度的read拆成个k-mer。根据一定策略，选定一个初始k-mer，接下来就可以在该k-mer为结点开始搜索后继的k-mer。搜索时采用贪婪图策略，每一步选择在当时看来最优的后继k-mer，直到满足事先设定的终止条件，结束一条contig的拼接，接着开始下一条contig的组装。直到没有合适的初始k-mer可供选择，整个组装过程结束。2. 基于de Bruijin图算法的一般步骤：1） 从 read 库中筛选 read。由于一些 read 错误率较高，不适合直接拿来建立 de Bruijn 图，故需要现将这些 read 过滤掉。 2） 确定值，建立 de Bruijn 图。这时需要扫描所有 read 数据，将每一个长为 的 read 拆分成个k-mer，并用所有 read 的所有k-mer来累加， 建立节点和边都加权的 de Bruijn 图；3） 化简de Bruijn 图，连续线性延伸节点合并为单一节点，产生一些碱基序列更 长的节点；4）遍历de Bruijn图产生contig。5.1.2 模型一的求解1. de Bruijn的数据结构 我们要将测序仪得到的read合理的保存在内存中。并且给定k-mer后，要快速的找到该k-mer出现在哪些read的哪些位置上。由于read很多，我们不可能将read的碱基保存在内存中，故要给每个read一个编号，称为readID，我们只需保存readID即可。如果确实需要保存碱基时，我们将碱基按照如下规则转化成整数：A-00，C-01，G-10，T-11（其中整数都是二进制的）。我们设计了如下的数据结构，称为readID_pos数组，每个k-mer都关联一个readID_pos数组，如图6.2.1-5所示。其中readID是我们给没条read的一个编号，在de Bruijn图中，一个readID能唯一识别一条read。readID是升序保存的，便于以后的二分查找操作。pos是该k-mer出现在编号为readID的read中的位置。另外，我们要经常删除拼接成功的read,故在结构中保存了删除标记。由于我们可能误删了某个read,所以使用删除标记的好处是便于我们及时恢复误删的read。更新read库更新集合Sread_i当前pos为L-k+1(最大pos)read状态拼接图read状态拼接图图6.2.1-3图6.2.1-4readID posdel_flag图6.2.1-5 readID_pos 数组我们将所有的k-mer存在一个哈希表中，哈希表实际上就是一个超大数组，以每个k-mer的哈希值作为下表就可找到该k-mer的所有信息。数组中每个元素都是只想如下所示的数据结构的指针，该数据结构就是k-mer结构，如图所示Kmer_seq addr num cur图6.2.1-6 k-mer 结构其中，kmer_seq是转化成整数之后的k-mer碱基，实际上该整数就是该k-mer的哈希值。addr是该k-mer所对应的readID_pos数组中元素的个数。addr和num课唯一标识一个数组。cur记录了readID_pos数组中下一个空闲的位置，因为填写readID_pos数组是不是一次性从头到尾填写，在遍历数据文件时。其中的数据会随机的落到某个readID_pos数组中，故每一个数组都要对应一个cur。de Bruijn图建立完成后，cur会有其它用途。由于我们要对read进行删除、恢复操作，故我们需要用cur字段来记录当前某个readID_pos数组中有效元素的个数。 整个 de Bruijn 图建好后如图 3-4 所示。 图 3-4 de Bruijn 图 2. de Bruijn图的建立过程建立de Bruijn图过程如下:（a） read的选取。数据文件中不仅保存了read的碱基，还保存了read的质量值。该质量值能反映read碱基的正确率。质量值越大，read碱基的正确率就越高。故我们选择那些质量较高的read。显然我们要事先指定一个质量阈值。文件中有的read包含未知碱基，这是由于测序过程中有些碱基没有被测出导致的。（b）预读数据文件，初始化de Bruijn图。依次读取每一个read，判断该read是否满足参与拼接的条件。若不满足，跳过该read；若满足，生成该read上所有k-mer（共85个k-mer），统计这些k-mer的出现次数，填写k-mer结构中的num字段，该过程如图6.2.1-9所示。统计k-mer的出现次数时，必须要找到该k-mer所对应的k-mer结构。这就需要计算该k-mer的哈希值，根据哈希表的首地址找到k-mer结构的地址，然后将其中的num字段加1，该过程如图6.2.1-8所示。k-mer碱基read碱基k-mer哈希值哈希表首地址相加k-mer结构地址图6.2.1-8开始读下一条read数据已读完？该read符合条件？计算每个Kmer结构的地址把Kmer结构中的num字段加1开始图6.2.1-9（c）遍历de Bruijn图，根据第（b）步中统计的k-mer个数，申请readID_pos 数组所需要的内存空间。实际上，整个系统中所有大型数组都是动态开辟的，因为事先我们不知道数组的大小，而且我们也不想申请一块足够大的内存来使用。另外，函数的栈空间非常有限，大型数组是绝对不能放到栈上的，为避免过多地使用全局变量，和静态变量，我们选择了动态申请内存的方法。再读数据文件，依次读取每一个read，判断该read是否满足参与拼接的条件。若不满足，跳过该read；若满足，生成该read上所有k-mer，填写每个k-mer的readID_pos数组，此时获得k-mer结构地址的方法与步骤（b）中的一样，如图3-7所示。然后填写readID_pos数组中的第cur行，填好后更新cur的值，即把cur的值加1，该过程如图6.2.1-10所示。当某个readID_pos数组填满后，其cur的值应该和num的值一样，表示当前readID_pos数组中所有元素都是有效的。开始读下一条read数据已读完？该read符合条件？生成25个kmer 填写该kmer的readID-pos数组中第cur行将cur值加1计算每个kmer结构的地址结束图6.2.1-10下面我们来分析一下建立deBruijn图的时间复杂度和空间复杂度。设数据文件里有条read，每条read长度为，k-mer长度为。首先，每条read都需要通过过滤器，过滤时要遍历read上每一个碱基，统计质量值。故过滤的时间与read的长度成正比，为。设有条read通过了过滤器。对这些read，每一条都要生成 个k-mer，生成每个k-mer时间与成正比，为。这时我们仅仅得到了k-mer的碱基，我们还要计算k-mer的哈希值，该计算时间依然与成正比，为。由于哈希表是无冲突的，故接下来可在常数时间内找到该k-mer所对应的k-mer结构以及readID_pos数组，并更新相应字段，故该时间为 。从而，建立de Bruijn图的总时间为： +-+ (3.1)一般地，与成正比，即，是大于1的常数。与近似相等，在我们的实验中，。代入公式3.1，经过化简可得建立deBruijn图的时间复杂度为。接下来分析空间复杂度。上述的k-mer结构大小是16字节，其中kmer_seq是无符号整型，占4字节，addr是指针，在64位操作系统上占8字节，num是无符号短整型，占2字节，cur也是无符号短整型，占2字节。readID_pos数组中每一行是4字节，其中readID占25位，pos占5位，del_flag占1位，保留位占1位。故建立de Bruijn图所需的空间总共为： (6.2.1-b)考虑到N与n近似相等，并将L=C*k代入公式3-2，化简可得建立de Bruijn图的空间复杂度为。接下来介绍如何查找某k-mer出现在某read的哪些位置。给定k-mer的碱基序列以及一个readID，可在de Bruijn图中找到该k-mer出现在该read的哪些位置。若再给定一个pos，可以判断该k-mer是否出现在该read的pos置。步骤如下：首先获得该k-mer所对应的k-mer结构地址（这一步可参考图6.2.1-10），进而获得readID_pos数组的首地址和数组大小。最后介绍read的删除、恢复操作。给定一条read的碱基序列以及readID，要在de Bruijn图中删除该read，既要找到所有该readID出现的地方，并把那一行的删除标记置1。看上去好像要遍历整个de Bruijn图，其实可以在OK2时间内完成，步骤如下：依次生成该read的每一个k-mer，找出该k-mer出现在该read中的所有位置，修改其删除标记，同时修改该k-mer的cur字段。至此，有关de Bruijn图的内容介绍完毕。我们先给出了有关de Bruijn图的数据结构，然后介绍了这些结构初始化及更新的方法，最后给出了有关de Bruijn图的基本操作。同时，我们还分析了建立de Bruijn图的时间复杂度和空间复杂度，以及一些基本操作的时间复杂度。由于de Bruijn图建立完成后，那些基本操作消耗极少内存，故我们并不关心它们的空间复杂度。3. 后继k-mer选取过程后继k-mer的选取是一个非常关键的步骤，如果后继k-mer选取的合理，就能拼接成质量较高的contig，其长度相对较长，和参考基因组匹配的较好，并且存在大量read拼接成功。反之，如果后继k-mer选的不合理，会导致contig很快拼接结束，其长度非常短（小于100bp），并且导致大量read拼接失败。由于基因组中存在大量重复片段，故较短的contig和参考基因组会匹配的较好，但由于大量read拼接失败，所以该contig的质量并不高。理想情况下，在任意选定初始k-mer后，在向左右扩展时，会逐渐拼接成整个基因组。但由于测试数据存在错误，并且基因组中存在大量重复片段，导致可能选取了错误的后继k-mer，进而导致拼接出的contig长度较短。由此可见，后继k-mer的选取是非常关键的。值得说明的是，并不是每次都能选取到非空的后继k-mer。如果deBruijn图中所剩的read较多，后继k-mer是不会为空的，只有当大量的read被删除，contig拼接的较长时，后继k-mer才有可能空。一旦选定了非空的后继k-mer，就用后继k-mer替换掉当前k-mer。4. contig 构建过程 （1）初始k-mer的选取与contig的初始化。 （2）contig的更新与保存。若后继k-mer非空，则要更新contig，即往contig上加一个碱基。contig本身是一个带头尾指针的单链表，该操作可在常数时间内完成，所需空间为，其中max是contig的最大长度。若后继k-mer为空，则一条contig拼接结束，此时要将其保存在文件中，并释放所占空间。此操作可在时间，空间内完成。（3）de Bruijn图中成功read的删除。这是de Bruijn图的基本操作可在时间,空间内完成。5.2 模型二的建立与求解5.2.1 模型二的建立模型二是建立在模型一的基础上，首先建立de Bruijin图，将k值定为4。把上述read1文件中的序列存入库中，开始建立read条目的数据结构和k-mer条目的数据结构。预读数据，逐条读取read数据，每条readID进行升序保存；生成该read上所有k-mer（共85个k-mer），统计这些k-mer出现的次数，填写k-mer结构中的num字段。相关数据录入程序源代码见附录。然后遍历de bruijn图，根据统计的k-mer数量，申请readID_pos数组所需要的内存空间。依次读取每一个read，填写readID_pos数组中的第cur行，填好之后把cur值加1。最后将碱基替换成2位二进制数。A=00，C=01，G=10，T=11。5.2.2 模型二的求解由于数据非常庞大，演算拼接过程不能完整的展示，接下来将列举一段算法拼接的过程： （1）初始k-mer定为 AAAC即 00000001，该k-mer出现在4条read上，且k-mer出现在每条read上的pos为1.这四条read开始参与拼接。如表为k-mer比对拼接相关代码；read100000001011111010100001011101111000100001001110100001000read200000001000011100000101010010000011000111000110011010100read300000001110001001010101110111111110101101101000010010010read400000001111110001110101000000011010001011111111111000010后继kmercontig00000001（2）此时num 为4，addr为1，cur为-1；（3）当前候选后继k-mer情况如下表：初始kmer 00 00 00 016候选后继kmer100 00 01 01候选后继kmer200 00 01 00候选后继kmer300 00 01 11候选后继kmer400 00 01 11 （4）选定后继k-mer为 AACT,即00000111；进行下一段拼接，此时num为5，前驱结点cur 为1，addr为2；此时contig增加了一个碱基；read100000001011111010100001011101111read200000001000011100000101010010000read300000001110001001010101110111111read400000001111110001110101000000011read500000100111111001111110000000010初始kmer00000001后继kmer00000111contig0000000111（5）重复（4）步骤，直到无法找到符合条件的后继k-mer则该条contig拼接结束；若该条contig的长度大于100bp则保存下来，否则删除；程序继续运行以上案列，得到一个长度为360bp的contig如下：110111110010010100001000100011000111111000111010010110010011100011011010111011000000010000111000010111110100001011110100010000111011000011001101111000000111111111010000001010010100100010010011011001011001111011001110100010010110001001110010010110010100111110110101111101010001010001011100001110100000110111101011011011011011011011010101111000101111010111110001001000111101101001101000110010000000110010000101001010001011100010011000111011101100010000010011100110101111100110101110000010100100000101111011001100110110001010101010110010001101001111010111001011111011001101001111111111101111110100111011110111111111001100001111001110111111001101001001001101110010011100101111010011101101110011110111001111001101001101111011 （6）到此一条contig拼接结束，开始重复（1）步骤；最终得到若干条长度大于100bp 的contig。 101100010011100100101100101001111101101011111010100010100010111000011101000001101111010110110110110110110110101011110001011110101111100010010001111011010011010001100100000001100100001010010100010111000100110001110111011000100000100111001100 六 模型评价与改进6.1模型的优点（1）本模型针对新一代测序技术出现的问题逐一进行了解决，把基因重组问题成功的转化为de bruijn图中的欧拉路径的问题，配合二分法、启发式搜索法实现了在较短的时间内对海量read数据进行比对处理。（2）本模型的算法高效，容易推广到实际的基因组组装中，具有一定的实际应用价值。6.2模型的缺点（1）由于比赛时间限制，本模型没有对组装中可能出现的误差进行较好的规避。如上述的基因组中出现重复片段干绕问题；对于海量数据比对过程中也应有的内存优化问题。七 参考文献1 王勇献，王正华 ，生物信息学算法导论，北京：清华大学出版社， （2011-06）.2 胡松年，薛庆中，基因组数据分析手册，浙江：浙江大学出版社.3 逯雯雯，卢志远，王亚旭，面向新一代基因组测序顺序的序列拼接算法， 文章编号1627-5565（2010）-03-248-06.4 牛北方，张西广，刘涛，郎显宇，陆忠华，迟学斌，基于新测序技术的比对 与组长算法.5 张法，刘志勇，乔香珍，刘玮，生物序列拼接算法phrap的并行化研 究.6 曾培龙，基于reads引导的基因组序列拼接，分类号：TP3919，（2012 年12月25日）.7 徐静怡，基因组序列拼接算法及ncRNA新基因的发现，（2004年6月 30日）.八 附录（1）DataInput.javapackage cn.data;import java.io.*;import java.sql.*;public class DataInput /* *A=00 *C=01 *G=10 *T=11 */Public static void main(String args) throws SQLException,ClassNotFoundException/ TODO Auto-generated method stubFile file=new File(E:基因数据,data1.txt);/*File file2=new File(D:基因数组,oo.txt);*/String a = new String60000;String b = new String60000;String c = new String60000;DataInput myData = new DataInput();tryFileReader inOne=new FileReader (file);BufferedReader inTwo=new BufferedReader(inOne);/*FileWriter outOne=new FileWriter(file2);*/*BufferedWriter outTwo=new BufferedWriter(outOne);*/String s=null;int i;int j;int k;int p;for(j=0,k=0,p=0,i=1;(s=inTwo.readLine()!=null;i+)if(i%4 = 1)aj = s;j+;else if(i%4 = 2)bk = myData.StringToByte(s);k+;else if(i%4 = 0)cp = s;p+;/*System.out.println(i+:+s);outTwo.write(i+:+s);/write()函数用法outTwo.newLine();*/inOne.close();inTwo.close();/*outTwo.flush();/outOne.close();outTwo.close();*/catch(Exception e)System.out.println(e);for(int i = 0;i100;i+)System.out.println(a+i+ai);for(int i = 0;i100;i+)System.out.println(bi);System.out.println(B length:+b.length );System.out.println(a40000+a40000);try Class.forName(com.mysql.jdbc.Driver);Stringurl= jdbc:mysql:/localhost/genmdata?user=cyh&password=1011325;Connection conn =DriverManager.getConnection(url);/Statement stmt = conn.createStatement();String sql = ;sql =insert into datasave(dataId,data,quality) values(?,?,?);PreparedStatement pstmt = conn.prepareStatement(sql, Statement.RETURN_GENERATED_KEYS);for(int i = 0;ia.length;i+)/sql = insert into datasave(dataId,data,quality) values(+ai+,+bi+,+ci+);pstmt.setString(1,ai);pstmt.setString(2,bi);pstmt.setString(3,ci);pstmt.executeUpdate();/stmt.executeUpdate(sql);/stmt.close();pstmt.close();conn.close(); catch (ClassNotFoundException e1) / TODO Auto-generated catch blocke1.printStackTrace();System.out.println(a1:+a1);public String StringToByte(String str)String reString=;char c = str.toCharArray();for(int i=0;ic.length;i+)if(ci=A)reString +=00;else if(ci=C)reString +=01;else if(ci=G)reString +=10;else if(ci=T)reString +=11;return reString;(2) DataToKmer.javapackage cn.data;import java.util.*;import java.sql.*;public class DataToKmer /* * param args */public static void main(String args) throws SQLException,ClassNotFoundException/ TODO Auto-generated method stub/*HashMap hs =new HashMap();*/* * 建立两个集合一个村存储read 一个存储Kmer * * */*String sql1 = ;*/String