基于be_Bruij图的基因组装算法 (1) (1)

摘要

快速和准确地或其生物提的遗传信息对生命科学研究具有重要的意义。测序技术从第一代到现在普遍应用的第二代以及正在兴起的第三代，能直接读取的碱基对序列长度远小于基因组长度。所以测序之前DNA分子要经过复制若干份、随机打断成短片段。要获得整个DNA片段，需要把这些片段利用重合部分信息组装连接。如何在保证组装序列的连续性、完整性和准确性的同时设计耗时短、内存小的组装算法是本题的关键。

新型测序技术使以往的基于重叠图的拼接算法不能胜任，本文中，提出了一种新的重叠群生成算法。该算法基于de bruijn图，将从多头测序转化成在de bruijn图的欧拉路径问题，并采用启发式搜索，能够快速地处理海量测序数据，而且能得到质量较高的重叠群。

本文详细叙述了算法的逻辑原理以及实现过程。确定k-mer长度后，将这些k-mer存入de bruijn图中。de bruijn图用哈希表储存，发现重叠关系式并不需要所有read之间进行两两比对，只要寻找de bruijn图或子图中的一条欧拉路径就可以找到contig。以初始k-mer为节点，采用贪婪策略获得质量较高的后继k-mer，保证了contig的高质量拼接，从而还原基因组。

本算法较为成功的弥补了新一代测序方法带来的一些弊端，在有限时间内对大数据的处理存在较大优势。但由于一些客观原因，对一些测序误差没有做到有效控制。最终在第二问的实践中也获得了质量较高的contig序列。

关键词：de Bruijn图 贪婪图方法 启发式搜索

一、问题的重述与问题的分析

快速和准确的获取生物体的遗传信息对生命科学研究具有重要意义。随着测序技术的不断发展，新一代测序技术产生的在高通量、低成本的同时也带来了错误率略有则加、读长较短等缺点。本题要求利用数学模型，设计算法解决如下几个问题：

（1）测序过程中可能出现的个别碱基对识别错误； （2）基因组中存在重复片段； （3）快速的处理海量的序列比对。

本题是基于新一代测序技术的基因组装算法问题，要求设计算法针对性的解决新一代测序技术带来的一些弊端。

1.1 read长度较短，数量较多——de bruijn图

新一代测序技术所得的read长度较短，数量较多，不易发现read之间的重叠关系。可以将read转化成定长的k-mer，然后寻找k-mer之间的重叠关系。然后建立de bruijn图，把短序列拼接问题转化为de bruijn图中的欧拉路径问题。

1.2 个别碱基对识别错误——多重对比纠错

通过将多个read放在一起比对来发现错误，如图2.1-1所示 。 图中通过途中4条read比对，可发现read3中的一个碱基错误（read3的第五个碱基）

AACA TGCA TGCT TGAC …… read1 TGCA TGCT TGAC ACAG …… reda2 TGCT CGAC ACAG CGTT …… read3 TGAC ACAG CGTT …… read4 图2.1-1

1.3基因组中存在大量重复片段

重复片段可能导致拼接错误，或者导致不连续的较短contig出现。重

叠片段类型主要有以下几种，如图1.3-1所示

重复片段问题可以用如下问题解决：通过对比，可先将重复片段隔离开来，较高的覆盖度有利于重复片段的隔离，但是，较多的测序错误将不利于该过程的进行。如果重复片段比read 长，可利用pared end read 来解决；如果重复片段比read 短，那么该read又被称为 spanner，一个spanner 就是一个重复片段两端再加几个碱基组成。利用spanner 解决重复片段问题需要如下两个信息：一是重复片段两端配对的read ，这两个read 必须不相同；二是重复片段中的一个配对read ，只要知道一个即可，另一个配对read 可以不在重复片段中。通过分析已知的基因组，可获得有关重复片段的更多信息，如，重复片段的长度，重复片段的模式等。

图1.3-1

二、问题假设

(1)假设测序过程中没有其他因素的干扰；

(2)假设题目所给定的序列相对位置的碱基全部遵循GU-AC法则； (3)假设题目中所有的序列都是正常可判别的序列，没有出现序列的基因突变等情况；

(4)假设一个完整基因组，打断成500bp的片段是随机的； (5)假设基因组每个位置被测到的几率是等可能的； (6)所有片段上的碱基都已经被识别出来，不存在未知碱基。

三、符号假设

符号 read

定义

利用现有的测序技术，可按一定的测序策略获得长度约为50–100个碱基对的序列，称为读长；

由read经过一定算法拼接产生3kb～10Mb以内

的一些基因组片段;

使用contig作为参考序列延伸，并进行合并得

到更长的contig,即supercontig; 根据本题数据，每一个read都含有一个质量

contig（C）：

superconti（S）：

quality（Q）: 值，该值能反映该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存入be bruijn图中，以备之后查找使用。此时要设定的一个重要参数是k-mer的长度。选定k值之后，要将长度L的read拆成L-K+1个k-mer。

根据一定策略，选定一个初始k-mer，接下来就可以在该k-mer为结点开始搜索后继的k-mer。搜索时采用贪婪图策略，每一步选择在当时看来最优的后继k-mer，直到满足事先设定的终止条件，结束一条contig的拼接，接着开始下一条contig的拼接。直到没有合适的初始k-mer可供选择，整个拼接过程结束。

4.2.3简述de Bruijn图的建立过程

基于de Bruijn图数据结构的read之间对比拼接算法可概括一下几个步骤。

（1）把筛选过得序列集合S作为参与比对的read库；

（2）确定k值，建立de bruijn图。这是需要扫描所有read数据，将每一个长为L的read拆分成L-k+1个k-mer,并用所有read的所有k-mer来累加，建立节点和边都加权的de bruijn图；

（3）化简de bruijn图，连续线性延伸节点合并成为单一节点，产生一些碱基序列更长的节点；

（4）遍历de bruijn图产生contig。

以上是对基于de bruijn图的算法做了简单介绍。是加上，de bruijn图是一种特殊的加权图，不仅图的及节点上有权值，而且边上也有权值。使用de bruijn图只能将较短的read拼接成较长的conig，要得到全基因组，还需要contig的组装过程。

五、建模与求解

5.1问题建模

根据新一代测序过程，本文建立了如下模型：设G是一个字符串，称之为基因组。取G的50-100个副本，然后将每个副本在随机的位置断开，形成一些小的字符串，每个片段长度小于100个字符。接着出去一些稍长的和烧断的字符串片段，将剩余字符串片段切成定长的字符串保存在多重

集U中，称U中的字符创为read。字符串上的每一个字符都来自字符集{A,C,G,T},每个字符就是一个碱基。为了方便运算，字符集与二进制字符集一一映射，这些字符串的长度都是L，考虑到新一代测序技术L的范围是10＜L＜100。

5.2模型求解

5.2.1基于k-mer对比法的conyig拼接过程

用满足S集合的所有read来构建contig。给定k值后，长度为k 的一个碱基片段k-mer。一般地，k要小于每条read的长度L，故每条read中含有k-mer数量L-k+1。一个k-mer的第一个碱基在一个read中出现的位置记为pos，pos的值从1开始，最大为L-k+1。

步骤一，选取初始k-mer。

拼接contig时，时候先要选定一个初始k-mer。初始k-mer要满足以条件：

（1）给定阈值min_read_num，要是该k-mer至少出现在min_read_num条read上；

（2）该k-mer出现在每条read上的pos为1。

只要有k-mer出现在某条read上的pos为1，该read就可以开始参与拼接。这时，contig上会有初始k-mer的k个碱基，如图所示。这样，就会有至少min_read_num条read参与拼接，这些read会影响到初始k_mer的扩展过程。

称图为read的拼接过程图。以后每当有read开始参与拼接时，就要将这些read加入到该图中，每当有read结束拼接时，就要将这些read从该图中删掉。此时，初始k-mer为当前k-mer，初始k-mer出现在某条read中的pos为该当前pos，记为cur_pos,现在所有read的当前pos为1。至此第一步骤结束，进入第二步骤。

步骤二，选取后继k-mer。

接着要选取当前k-mer的后继k-mer。后继k-mer至少要满足如下条件：

（1）后继k-mer的前k-1碱基与当前k-mer的后k-1个碱基相同； （2）后继k-mer要尽可能多的出现在正在参与拼接的read中，且出现的位置为read的当前pos+1，这时称该k-mer出现在该read的合适位

置上；

（3）后继k-mer要是尽可能多的read开始参与拼接。

显然，给定当前k-mer后，候选的后继k-mer有四个，如图6.2.1-1所示。接上面的例子，选定的后继k-mer是ACCG，此时，contig上多了一个碱基，当前k-mer变成ACCG，read1至read4当前的pos变为2，read5处于拼接中断状态。由于当前k-mer可能出现在其它read上，且出现的pos为1，所以要让这些read开始参与拼接，如图6.2.1-2所示的read6，read7，read8。

初始k-mer A A C C

候选后继k-mer1 A C C A 候选后继k-mer2 A C C G 候选后继K-mer3 A C C G 候选后继k-mer4 A C C T

图2.初始k-mer的后继k-mer

read1 A A C C G C T A C T A T read2 A A C C G C A A C T A T read3 A A C C G C T C C T A T read4 A A C C G G T A C T A T read5 A A C C T C T A G T A T read6 A C C G C A A C T A T C read7 A C C G T T A C T A T C read8 A C C G C T G C T A T C 初始kmer A A C C 后继k-mer A C C G contig A A C C G 图3. 新的read 开始参与拼接（k=4 ，L=12 ，min_read_num=5 ）

接下来，需要反复重复步骤二，知道找不到合适的后继k-mer位置。虽然每次必能产生四个候选k-mer，当出现如下情形时，没有合适的后继k-mer可供选择：

（1）对于任意一个候选后继k-mer，它不出现在任何当前参与拼接的read的合适位置，即不管我们选择哪个后继k-mer都将导致read拼接过程图中所有read处于拼接中断状态；

（2）对于任意一个候选后继k-mer，对于任意一条read库中的read，

该k-mer不出现在该read上pos为1的位置，即不管我们选择哪个后继k-mer，都不会使新的一批read开始参与拼接。

当上述情形同时出现时，我们将找不到合适的后继k-mer，此时一条contig拼接结束，如果该contig足够长（长度不小于100bp），我们会保存该contig，否则我们将丢弃该contig。如果我们丢弃了该contig，我们要把那些处于成功拼接状态的read从新加入到read库中，让它们继续参与后面的拼接。

如果没有找到合适的后继k-mer，则转到步骤一，开始下一条contig的拼接过程。直到某时刻初始k-mer选取失败时，整个拼接过程结束。当出现如下情形时，初始k-mer选取失败：

（1）read库中所剩read已经不多了，任意一个k-mer都出现在至多min_read_num-1条read上；

（2）有些k-mer虽然出现在至少min_read_num条read上，但这些k-mer中的任意一个都至多出现在min_read_num-1条read上pos为1的位置。

如果有上述情形之一出现时，将导致初始k-mer选取失败，整个拼接过程结束。拼接过程中，有如下情形需要特殊处理：

a）若某条read的当前pos达到了L-k+1，那么这条read上的所有碱基都出现在contig上了，即这是一条成功拼接read，上述模型中的read集合S就是由这些read组成的。我们要记录该read并将其在read库中删除掉，以保证在以后的拼接过程中不在使用该read。该过程如图6.2.1-3所示。

b） 若当前k-mer没有出现在某条正在参与拼接read的合适位置上，则该read处于拼接中断状态（如图2.6中的read5）。设x是一个k-mer，y是x的后继k-mer，x出现在该read的合适位置，而y没有出现在该read的合适位置，我们记录x出现在该read上的pos，称为last_appear_pos，并将该read的当前pos置0，这就是中断read的处理过程。

5.3各模块算法实现 5.3.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。

集合S S 集合

更新

readread库库

更新

read状态拼接图 read_i当前pos为L-k+1(最大pos) 集合S read库

read状态拼接图

图4

图5

readID pos del_flag 图6. readID_pos 数组

我们将所有的k-mer存在一个哈希表中，哈希表实际上就是一个超大数组，以每个k-mer的哈希值作为下表就可找到该k-mer的所有信息。数组中每个元素都是只想如下所示的数据结构的指针，该数据结构就是k-mer结构，如图所示

Kmer_seq addr

num cur 图7.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 btuijn图建好后如图所示

6.3.2 de bruijn图的建立过程

下面将详细叙述de bruijn图的建立过程。该过程主要分两步。第一步先预读一遍数据文件，统计每个readID_pos数组的大小，记录在k-mer结构的num字段中。然后根据统计结构开辟相应的内存空间；第二步再读

一遍数据文件，填写每个k-mer的readID_pos数组，该过程如图8所示。

图8. 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字段，该过程如图10所示。统计k-mer的出现次数时，必须要找到该k-mer所对应的k-mer结构。这就需要计算该k-mer的哈希值，根据哈希表的首地址找到k-mer结构的地址，然后将其中的num字段加1，该过程如图9所示。

read碱基 k-mer碱基 k-mer哈希值 哈希表首地址

相加 k-mer结构地址 图9

开始 读下一条read 数据已读完？该read符合条件？ 计算每个Kmer结构的地址 把Kmer结构中的num字段加1 开始

图10

（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，该过程如图11所示。当某个readID_pos数组填满后，其cur的值应该和num的值一样，表示当前readID_pos数组中所有元素都是有效的。

开始 读下一条read 数据已读该read符合条生成25个kmer 计算每个kmer结构的地址 填写该kmer的readID-pos数组中第cur行 将cur值加结束 图11

下面我们来分析一下建立de Bruijn图的时间复杂度和空间复杂度。设数据文件里有N条read，每条read长度为L，k-mer长度为k。首先，每条read都需要通过过滤器，过滤时要遍历read上每一个碱基，统计质量值。故过滤的时间与read的长度L成正比，为O（L）。设有n条read通过了过滤器。对这些read，每一条都要生成L –k+1个k-mer，生成每个k-mer时间与k成正比，为O（k）。这时我们仅仅得到了k-mer的碱基，我们还要计算k-mer的哈希值，该计算时间依然与k成正比，为O（k）。由于哈希表是无冲突的，故接下来可在常数时间内找到该k-mer所对应的k-mer结构以及readID_pos数组，并更新相应字段，故该时间为O（1）。从而，建立de Bruijn图的总时间为：

2NOLnLK1OKO1 (6.2.1-a)

一般地，L与k成正比，即L=C*k，C是大于1的常数。N与n近似相等，在我们的实验中，L=3*k，n=0.82*N。代入公式3.1，经过化简可得建立de Bruijn图的时间复杂度为O(N*L^2)。

接下来分析空间复杂度。上述的k-mer结构大小是16字节，其中kmer_seq是无符号整型，占4字节，addr是指针，在位操作系统上占8字节，num是无符号短整型，占2字节，cur也是无符号短整型，占2字节。readID_pos数组中每一行是4字节，其中readID占25位，pos占5位，del_flag占1位，保留位占1位。故建立de Bruijn图所需的空间总共为：

LK1n44k816 (6.2.1-b)

考虑到N与n近似相等，并将代入L=C*k建立de Bruijn图的空间复杂度为O(L*N+4^k).

接下来介绍如何查找某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图，其实可以在O{K^2}时间内完成，步骤如下：依次生成该read的每一个k-mer，找出该k-mer出现在该read中的所有位置，修改其删除标记，同时修改该k-mer的cur字段。

至此，有关de Bruijn图的内容介绍完毕。我们先给出了有关de Bruijn图的数据结构，然后介绍了这些结构初始化及更新的方法，最后给出了有关de Bruijn图的基本操作。同时，我们还分析了建立de Bruijn图的时间复杂度和空间复杂度，以及一些基本操作的时间复杂度。由于de Bruijn图建立完成后，那些基本操作消耗极少内存，故我们并不关心它们的空间复杂度。

3-2，化简可得

5.3.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。如果de Bruijn图中所剩的read较多，后继k-mer是不会为空的，只有当大量的read被删除，contig拼接的较长时，后继k-mer才有可能空。一旦选定了非空的后继k-mer，就用后继k-mer替换掉当前k-mer。

6.3.4contig 构建过程的其它模块

拼接contig的其他模块主要有如下几个：（1）初始k-mer的选取与contig的初始化。

这两个模块已经在第2章中介绍过了，参见图6.3.4-1。可在常数时间，常数空间内完成初始k-mer的选取，在O(k)时间，O(k)空间内完成contig的初始化。毕竟它们不是主要模块，故在此不在详细描述。

（1）contig的更新与保存。

若后继k-mer非空，则要更新contig，即往contig上加一个碱基。contig本身是一个带头尾指针的单链表，该操作可在常数时间内完成，所需空间为O(max)，其中max是contig的最大长度。若后继k-mer为空，则一条contig拼接结束，此时要将其保存在文件中，并释放所占空间。此操作可在O(max)时间，O(max)空间内完成。

（2）de bruijn图中成功read的删除。

read1 A A C C G C T A C T A T read2 A A C C G C A A C T A T read3 A A C C G C T C C T A T read4 A A C C G G T A C T A T read5 A A C C T C T A G T A T 初始kmer A A C

contig A A C C

C

图12. 初始k-mer （k=4 ，L=12 ，min_read_num=5 ）

这是de bruijn图的基本操作可在O(k*k)时间，O(k)空间内完成。 这一章介绍了整个系统流程。一开始给出了系统流程图，接下来详细叙述了主要模块。其中，该系统有两大主要模块，分别是de bruijn图的建立模块和contig构建模块。在contig构建模块中又有两个主要过程，分别是决策表更新过程与后继k-mer选取过程。本章详细分析了这些模块，并给出了时间复杂度与空间复杂度。

六、模型二

6.1. 问题分析

本题题目提供全长约为120,000个碱基对的细菌人工染色体，采用新一代的Hiseq2000测序仪进行测序。附件提供了筛选好的定长read数据文件。使用第一题设计的基于de bruijn图的组装算法对read数据进行组装，并对结果进行误差分析。 6.2.数据分析

由测序策略可知，read1和read2为相互补的两条单链，故选用read1的数据带入算法进行组装，read2可以作为校准链备用。分析read1文件可知，本题数据已满足如下条件：

（1）序列片段被切成固定长度L=88bp； （2）经过复制，原基因至少有6个副本；

（3）所有片段上的碱基都已经被识别出来，不存在未知碱基； （4）由于技术，本文不对质量数进行讨论，假设read1中的所有

片段满足正确率要求。 6.3.带入模型求解 6.3.1建立de bruijn图

（1）将k值定为4。把上述read1文件中的序列存入库中，开始建立read条目的数据结构和k-mer条目的数据结构。预读数据，逐条读取read数据，每条readID进行升序保存；生成该read上所有k-mer（共85个k-mer），统计这些k-mer出现的次数，填写k-mer结构中的num字段。如图所示，为相关代码片段。相关数据录入程序源代码见附录。

图13

（2）遍历de bruijn图，根据上一步统计的k-mer数量，申请readID_pos数组所需要的内存空间。依次读取每一个read，填写readID_pos数组中的第cur行，填好之后把cur值加1。

（3）将碱基替换成2位二进制数。{A=00，C=01，G=10，T=11}。 6.3.2模型求解

由于数据非常庞大，演算拼接过程不能完整的展示，接下来将列举一段算法拼接的过程：

（1）初始k-mer定为 AAAC即 00000001，该k-mer出现在4条read

上，且k-mer出现在每条read上的pos为1.这四条read开始参与拼接。如图为k-mer比对拼接相关代码

图14

（2）此时num 为4，addr为1，cur为-1；

图15

（3）当前候选后继k-mer情况如下图：

图16

（4）选定后继k-mer为 AACT,即00000111；进行下一段拼接，此时num为5，前驱结点cur 为1，addr为2；此时contig增加了一个碱基；

read1 read2 read3 read4 read5 初始kmer 后

继kmer contig

0 0

00 0 0 0

0 00 00 00 0 00 0 00

00 00 00 00 00 00 00 00

01 01 01 01 01 01 01 01

01 01

101 00 01 01 00 0

01 00 10 11 01 1

11 01 01 10 11

11 10 00 10 11

01 10 00 01 10

00 00 10 10 01

00 00 10 10 11

00 10 11 10 11

11 10 10 10 10

10 11 11 00 00

11 00 11 00 00

11 00 11 01 0

00 0

1101

1

（5）重复（4）步骤，直到无法找到符合条件的后继k-mer则该条contig

拼接结束；若该条contig的长度大于100bp则保存下来，否则删除；程序继续运行以上案列，得到一个长度为360bp的contig如下：

{ 110111110010010100001000100011000111111000111010010110010011100011011010111011000000010000111000010111110100001011110100010000111011000011001101111000000111111111010000001010010100100010010011011001011001111011001110100010010110001001110010010110010100111110110101111101010001010001011100001110100000110111101011011011011011011011010101111000101111010111110001001000111101101001101000110010000000110010000101001010001011100010011000111011101100010000010011100110101111100110101110000010100100000101111011001100110110001010101010110010001101001111010111001011111011001101001111111111101111110100111011110111111111001100001111001110111111001101001001001101110010011100101111010011101101110011110111001111001101001101111011 }

（6）到此一条contig拼接结束，开始重复（1）步骤；最终得到若干条长度大于100bp 的contig。

{ 101100010011100100101100101001111101101011111010100010100010111000011101000001101111010110110110110110110110101011110001011110101111100010010001111011010011010001100100000001100100001010010100010111000100110001110111011000100000100111001100 } 6.3.3 模型二总结

通过基于de bruijn图的基因组装算法，成功拼接2条长度在100bp以上congtig；中途算法自动放弃长度不够的contig 3条。

contig拼接成功率为0.54945 %。理想情况下，S中的read 也能完 整覆盖原基因组，即使是这样，也不能保证 g与原基因组完全一样，主要原因是基因组中存在大量重复片段。另外，测序过程中的错误也不容忽视。设ra，rb是S中的两条完全一样的read ，则有如下几种可能：

（1）rb 是ra 的复本（或ra 是rb 的复本）；

（2 ）rb 不是ra 的复本，他们是来自重复片段的两条read ； （3 ）rb 中包含测序错误，ra 中无错误，导致rb 碱基与ra 相同； （4 ）ra 中包含测序错误，rb 中无错误，导致rb 碱基与ra 相同； （5 ）ra和rb中都有错误，它们的碱基相同纯属巧合。所以，现在能做的是使碱基序列 g 尽可能与原始基因组一致，存在误差在所难免，本文要使该误差尽可能的小。

七、模型评估

7.1模型评价

（1）本模型针对新一代测序技术出现的问题逐一进行了解决。新一代读长较短，不易发现read之间的重叠关系，本文放弃了传统的重叠图算法。把基因重组问题成功的转化为de bruijn图中的欧拉路径的问题，配合二分法、启发式搜索法实现了在较短的时间内对海量read数据进行比对处理。

（2）由于比赛时间，本模型没有对组装中可能出现的误差进行较好的规避。如上述的基因组中出现重复片段干绕问题；对于海量数据比对

过程中应有的内存优化问题。

参考文献

【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日）。 【8】《基因组组装新技术研讨》，http://www.novogene.com,访问日期（2014年5月11日）。

附录：

源代码文件

（1）DataInput.java package 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{

File file=new File(\"E:\\\\基因数据\ /*File file2=new File(\"D:\\\\基因数组\ String []a = new String[60000]; // TODO Auto-generated method stub

String []b = new String[60000]; String []c = new String[60000]; DataInput myData = new DataInput(); try{

FileReader 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){ a[j] = s; j++;

}else if(i%4 == 2){

b[k] = myData.StringToByte(s); k++;

}else if(i%4 == 0){ c[p] = 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;i<100;i++)

System.out.println(\"a[\"+i+\"]\"+a[i]); for(int i = 0;i<100;i++){ System.out.println(b[i]);

System.out.println(\"B length:\"+b.length );

}

System.out.println(\"a[40000]\"+a[40000]); try {

Class.forName(\"com.mysql.jdbc.Driver\"); String

url

=

\"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;i=

\"insert

into

datasave(dataId,data,quality)

values('\"+a[i]+\"','\"+b[i]+\"','\"+c[i]+\"')\";

pstmt.setString(1,a[i]); pstmt.setString(2,b[i]); pstmt.setString(3,c[i]);

pstmt.executeUpdate();

//stmt.executeUpdate(sql); }

//stmt.close(); pstmt.close(); conn.close();

} catch (ClassNotFoundException e1) { // TODO Auto-generated catch block e1.printStackTrace(); }

System.out.println(\"a1:\"+a[1]); }

public String StringToByte(String str){ String reString=\"\";

char c[] = str.toCharArray(); for(int i=0;i}else if(c[i]=='G'){ reString +=\"10\"; }else if(c[i]=='T'){ reString +=\"11\"; } }

return reString; }

}

(2) DataToKmer.java

package 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 contig = \"\";

readBank []readBank =new readBank [46847]; Read read[] = new Read[46847]; HashMap hsToRead = new HashMap();

List bankList = new ArrayList(); String getRead = \"\"; //初始Kmer try {

Class.forName(\"com.mysql.jdbc.Driver\"); String

url

=

\"jdbc:mysql://localhost/genmdata?user=cyh&password=1011325\";

Connection conn =DriverManager.getConnection(url); ResultSet rs1;

Statement stmt = conn.createStatement(); PreparedStatement pstmt;

String sql = \"select count(*) from datasave\"; int datalength = 0;

rs1 = stmt.executeQuery(sql); if(rs1.next())

datalength

Integer.parseInt(rs1.getString(\"count(*)\")); System.out.println(\"datalength=\"+datalength); sql = \"select * from datasave\"; rs1 = stmt.executeQuery(sql); //该数组用来判定是否为该序列 /*String[] readId ;*/ String readById = \"\";

//该数组用来存取第一个data的值（data） String[] readFirst = new String [386]; int i = 0; //readBank的录入

//第一个序列的data数据已经拿到 while(rs1.next()){

readBank readbank = new readBank();

readbank.setDataIndex(i);

=

readbank.setReadData(rs1.getString(\"data\")); bankList.add(readbank); String[]

readId

=

rs1.getString(\"dataId\").split(\":\");

/*System.out.println(\"readID:\"+readId[2]);*/ readById = readId[2]; if(readById.equals(\"1101\")){

readFirst[i] =rs1.getString(\"data\"); /*System.out.println(\"第\"+i+\"个data为

\"+readFirst[i]);*/ /*

} i++;

/*System.out.println(rs1.getString(\"data\"));*/ }

System.out.println(\"list46847=\"+bankList.get(46846).re

adData);*/

//开始选取初始Kmer(firstKmer) 先坐1101序列 /*1.先取出第一个data找与它前22个相同的read,用一

个dataReConunt记录，若是后面的比当前的大就覆盖

*2.并且返回这些相同read的数组序号

*/

String dataCut1 = \"\"; String dataCut2 = \"\"; ArrayList

finalList=

new

ArrayList();

ArrayList

middleList1=

new

ArrayList();

ArrayList

middleList2=

new

ArrayList();

int finalCount = 0; int middleCount1 = 0; int middleCount = 0; int index[] = new int[6];

Reading firstReading = new Reading(); Reading middleReading = new Reading(); for(int j = 0;jfor(int k = 0 ;k/*

if(dataCut2.equals(dataCut1)){ if(middleList1.size()==0){ firstReading.setData(dataCut1);

middleList1.add(middleCount,firstReading); index[middleCount] = j; middleCount++; }*/

middleReading.setData(readFirst[k]); middleReading.setPos(\"1\"); middleReading.setDel_flag(\"0\");

middleList1.add(middleCount,middleReading); index[middleCount] = k; middleCount++; } }

middleCount1 = middleCount; middleCount = 0; middleList2

=

(ArrayList)

middleList1.clone();

middleList1.clear(); if(finalCountfinalList

(ArrayList)middleList2.clone(); } }

System.out.println(\"

初

Kmer:\"+finalList.get(0).getData()); //初始Kmer设定,并且定义Kmer库

ArrayList

kmerBanking

=

ArrayList(); Kmer firstKmer = new Kmer();

firstKmer.setKmer_seq(finalList.get(0).getData()); firstKmer.setNum(String.valueOf(finalCount)); firstKmer.setCur(\"-1\"); firstKmer.setAddr(\"0\"); kmerBanking.add(firstKmer);

=

始

new

ArrayList readBanking = new

ArrayList();

for(int l = 0;lread1.setPos(String.valueOf(index[l])); //0表示正在参与拼接1表示使用完毕已删除 System.out.println(index[l]); read1.setDel_flag(\"0\");

read1.setReadData(readFirst[l]); readBanking.add(read1); } /*

* 开始选择后继的Kmer * */

//此处为判定read库中是否为空 Kmer kmering = firstKmer; String checkNull = \"\"; String checkCount = \"\"; int KmerCount ; /*;*/

int count[] = new int[4]; int countMax = 0; int myIndex = 0; int Pos = 2;

for(int l = 0;lSystem.out.println(\"contig=\"+contig); for(int m = 0;mcheckNull = \"lose\"; } }

for(int m = 0;mcheckNull = \"lose\"; }

}

KmerCount = readFirst.length-4;

System.out.println(\"Kmer共有:\"+KmerCount); for(int m = 0;msubString

=

kmerBanking.get(Integer.parseInt((kmering.getAddr()))).getKmer_seq().substring(2,8);

if(readFirst[m].substring(0,6).equals(subString)){

if(readFirst[m].substring(7,8).equals(\"00\")){ count[0]++; }else

if(readFirst[m].substring(7,8).equals(\"01\")){

count[1]++; }else

if(readFirst[m].substring(7,8).equals(\"10\")){

count[2]++; }else

if(readFirst[m].substring(7,8).equals(\"11\")){

count[3]++; }

countMax=Math.max(Math.max(count[0],count[1]),

Math.max(count[2],count[3]));

for(int k = 0;kcheckCount = \"00\"; }else if(k==1){ checkCount = \"01\"; }else if(k==2){ checkCount = \"10\"; }else if(k==3){ checkCount = \"11\"; } } }

kmering.setKmer_seq(subString+checkCount);

String addString = subString+checkCount; for(int k=0;kindex1[myIndex] = k; myIndex++; middleCount++; } }

/* for(int l=0;lfor(int p = 0;pmiddleReading.setPos(String.valueOf(Pos)); middleReading.setDel_flag(\"0\"); middleList1.add(middleReading); for(int q=0;qif(88-Integer.parseInt(finalList.get(q).getPos())==1){

finalList.get(q).setDel_flag(\"1\");

} }

finalList.add(middleReading); } } }

for(int m = 0;mcheckNull = \"lose\"; } }

/*System.out.println(finalList.get(3).getData());*/

/*System.out.println(finalList.get(1).getData());*/ //将正在参与拼接的read从read库中清除 for(int l=0;l/* int dataSaveArray [] = new int[100];*/

/*

String dataCut1 = \"\"; String dataCut2 = \"\";

Reading firstReading = new Reading(); List

indexSaveArray

=

new

ArrayList();

for(Integer j = 0;j

indexSaveArray2

=

new

ArrayList();

int dataReCount2 = 0;

dataCut1 = readFirst[j].substring(0,43); firstKmer.setKmer_seq(dataCut1); firstKmer.setAddr(\"0\"); firstReading.setData(dataCut1); for(Integer

k=0;k&&

j!=k;k++ ){

dataCut2 = readFirst[k].substring(0,43); if(dataCut1.equals(dataCut2)

&&

dataReCount2>dataReCount){

Reading reading = new Reading();

reading.setData(dataCut2); if(indexSaveArray2.get(0)==null){

indexSaveArray2.add(index,firstReading); index++; }

dataReCount2++;

indexSaveArray2.add(index, reading); index++;

System.out.println(dataCut2); } }

if(dataReCount2>dataReCount){ dataReCount = dataReCount2; }*/

/* indexSaveArray = indexSaveArray2;*/ /*}*/

/* System.out.println(indexSaveArray.get(0));*/ rs1.close(); stmt.close(); } catch (Exception e) {

// TODO Auto-generated catch block

e.printStackTrace(); } }

}

(3)Kmer.java

package cn.data;

public class Kmer { String Kmer_seq; String addr; //当前正在拼接的read的个数 String num; String cur; public String getKmer_seq() { return Kmer_seq; } public void setKmer_seq(String kmer_seq) { Kmer_seq = kmer_seq; } public String getAddr() { return addr; } public void setAddr(String addr) { this.addr = addr; } public String getNum() { return num; } public void setNum(String num) { this.num = num; } public String getCur() { return cur;

}

public void setCur(String cur) { this.cur = cur; }

}

(4) Read.java

package cn.data;

public class Read { String readData; String pos; String del_flag;

public String getReadData() { return readData; }

public void setReadData(String readData) { this.readData = readData; }

public String getPos() { return pos; }

public void setPos(String pos) { this.pos = pos; }

public String getDel_flag() { return del_flag; }

public void setDel_flag(String del_flag) { this.del_flag = del_flag; }

}

(5)ReadBank.java

package cn.data;

public class readBank { String readData; int dataIndex;

public String getReadData() { return readData; }

public void setReadData(String readData) { this.readData = readData;

}

public int getDataIndex() { return dataIndex; }

public void setDataIndex(int dataIndex) { this.dataIndex = dataIndex;

}

}

(6) Reading.java

package cn.data;

public class Reading { String data; String pos; String dataId; String del_flag;

public String getDel_flag() { return del_flag; }

public void setDel_flag(String del_flag) { this.del_flag = del_flag;

}

public String getData() { return data; }

public void setData(String data) { this.data = data; }

public String getPos() { return pos; }

public void setPos(String pos) { this.pos = pos; }

public String getDataId() { return dataId; }

public void setDataId(String dataId) { this.dataId = dataId;

}

}