基于 de Bruijn 图算法概述

1、基于 de Bruijn 图的算法概述de Bruijn 图简介传统的 Sanger 测序的 reads 较长（1000bp），数据量较少，精度较高，所有的组装算法都利用 reads 之间的重叠，通过公共路径的方法解决拼接问题。而新一代测序产生的数据 read 更短、覆盖度更高、序列精度较低，为此这种read 为中心的方法面临海量计算的困境，似乎不可能找到恰当的启发式方法来处理大量的重叠。de Bruijn图框架为处理高覆盖、短序列提供了很好思路，该框架借鉴了 Pevzner和 Waterman 等人针对传统的长 reads 提出的欧拉遍历方法37,38，并在此基础上针对新一代测序数据的特点进

2、行了改进要想以较低的成本快速得到某个新物种的 DNA 分子碱基序列，就要依靠新一代的测序技术和从头测序拼接组装算法。目前新一代测序数据用于从头测序的短序列拼接组装算法普遍采用 de Bruijn图数据结构。在 de Bruijn图上，每一个 k-mer 都构成图的节点，如果两个 k-mer 在某一 read 中相邻，那么这两个节点之间就有一条边。reads 集合中的每个 read 都对它所含的节点和边加权，这样 reads 集合产生一个节点和边都具有权值的 de Bruijn图。在存储每一个 k-mer 时，往往要建一个无冲突的哈希表，以加快查找速度。而建立哈希表可能会消耗更多的内存。但是，由

3、于每个 k-mer 在哈希表中只存储一次，不管该k-mer 在 read 中出现了多少次，所以实际消耗的内存小于存储所有 read 所需要的空间。另外，基因组中的重复片段会在 de Bruijn图中产生环路。环路将在遍历 deBruijn 图时产生障碍。目前的研究主要面临两个问题，一个是基因组中存在大量重复片段，一个是测序错误。这两个问题相互影响，使问题变的更加复杂。本文通过仔细分析这两个问题，来改进以前基于 de Bruijn 图的算法，提出一种新的 deBruijn 图，并且引入了决策表的概念，通过决策表里的信息来选取后继 k-mer，并在适当的时候更新决策表。1 基因组中存在大量重复片段

4、重复片段问题可用如下方法解决：通过比对，可先将重复片段隔离开来，较高的覆盖度有利于重复片段的隔离，但是，较多的测序错误将不利于该过程的进行。因为错误的存在，严格的比对将导致一些重复片段未被发现，而非严格的比对会把一些不是重复片段的区域隔离开来，这不是本文所希望的。如果重复片段比 read 长，可利用 pared end read 来解决；如果重复片段比 read 短，那么该 read又被称为 spanner，一个 spanner 就是一个重复片段两端再加几个碱基组成。利用spanner 解决重复片段问题需要如下两个信息：一是重复片段两端配对的 read，这两个 read 必须不相同；二是重复片

5、段中的一个配对 read，只要知道一个即可，另一个配对 read 可以不在重复片段中2 测序过程中可能出现错误现在主要有两种纠错方法，一种基于多重比对，通过将多个 read 放在一起比对来发现错误，如图1-2 所示。通过图中 4 条 read 比对，可发现 read 3 中的一个碱基错误(read 3 的第 5 个碱基)，该方法在 overlap 过程中比较常用，而在 de Bruijn图中，所使用的纠错方法是：若当前 k-mer 在一条 read 中连续未出现恰好 k 次，可以认为该 read 中存在一个碱基错误。2 基于 de Bruijn 图算法的一般步骤1) 确定 k 值，建立 de

6、Bruijn 图。这时需要扫描所有 read 数据，将每一个长为 L 的 read 拆分成 L-k+1 个 kmer，并用所有 read 的所有 k-mer 来累加，建立节点和边都加权的 de Bruijn图；2) 化简 de Bruijn 图，连续线性延伸节点合并为单一节点，产生一些碱基序列更长的节点；3) 错误校正，删去由于测序错误产生的尖端和泡状结构；4) 通过 read 的配对末端 (pair-end)、环化配对(mate-pair)信息伸展或者删去一些环；5) 依据环上节点和边的权值(覆盖深度信息)进一步伸展或者删去一些环；6) 遍历 de Bruijn 图产生 contig。实际上

7、，de Bruijn图是一种特殊的加权图，不仅图的结点上有权值，而且图的边上也有权值。化简 de Bruijn图是非常关键的一个步骤，通过对 de Bruijn图化简，可降低算法的时间复杂性以及空间复杂性，同时可以保证错误校正顺进行拼接总体思路假设所有满足上述条件（1）的 read 都已经存到了 read 库中，下面就用这些read 来构建 contig。给定 k 值后，长度为 k 的一个 DNA 片段称为一个 k-mer。一般地，k 要小于每条 read 的长度 L，故每条 read 中含有 k-mer 数量为 L-k+1。一个k-mer 的第一个碱基在一个 read 中出现的位置记为 po

8、s，pos 的值从 1 开始，最大为 L-k+1，如图 2-2 所示。选定一个初始 k-mer 后，通过对该 kmer 不断扩展，来得到一条 contig。一个k-mer 上有 k 个碱基，而碱基共有四种，故扩展下一个碱基有四种选择，这样就会形成一个四叉树，如图 2-3 所示。显然，这个四叉树的深度是无限的，任何一个子树的深度也是无限的。算法中要设定终止条件，不能让它无限地扩展下去。实际上，该树中任何一条有限的路径都可能成为一条 contig，每条 contig 都可以使一些 read 成为它的子串，所以可以用 read 库中的 read 来评价 contig的好坏。虽然可以用无信息搜索的方法

9、来拼接 contig，但现在的问题是，有些 contig 长达几万 bp，这样算法要搜索上万层，搜索空间过大，以至于不能在有效的时间内完成。故本文采用启发式搜索，来减少时间开销。在一条 contig开始拼接前，需要根据一定策略，选定树中一个初始 k-mer，接下来就可以在以该 k-mer 为根结点的子树上开始搜索。搜索时采用贪心策略，每一步选择在当前看来最优的后继 k-mer，直到满足事先设定的终止条件，结束一条contig 的拼接，接着开始下一条 contig 的拼接。直到没有合适的初始 k-mer 可供选择，整个拼接过程结束。由于选定初始 k-mer 后，可以向该 k- mer 的两端分别

10、扩展，故初始 k-mer 选取的好坏对拼接结果影响不大。故该问题的关键是选取后继 k-mer。后继 k-mer 如果选择的好，contig 会拼接得较长，会有较多的 read 成功参与拼接；后继 k-mer 如果选择的差，contig会拼接得较短，会有较少的 read 成功参与拼接。基于de Bruijn图的序列拼接技术分析Idba拼接技术Velvet、Soapdenovo等在处理无错误序列、高覆盖度的序列拼接问题时,能够表现很好的性能,但是,由于这些技术在拼接过程中以k-mer为基本单位,就不可避免的会产生很多重叠单元,这使得拼接面临着错误位置拼接、顶点缺失和覆盖度低等问题。由于一些错误re

11、ad的存在,产生了大量的分支区域。k-mer长度越小,分支问题越严重,k-mer长度越大,出现的重叠区域变得越少,这直接影响了拼接的质量。正确的选择k-mer的大小成为影响拼接质量的一个关键因素Idba的主要特点是:按弃了使用固定的k-mer长度构建de Bruijn图的方法采用一个变化的k-mer长度完成read切割过程。首先,它设置k-mer长度的变化区域,其中A为当前k-mer的长度,采用重复迭代算法来计算每一个k-mer的长度;在构建图之前,该技术对低覆盖度的k-mer进行了预处理,去除覆盖度低的k-mer顶点,从而简化了 de Bmijn图的结构,使得其在内存消耗上明显降低,提高了算

12、法的拼接效率。设置k-mer阈值的方法,成功解决了 de Bmijn图中路径多分支问题,提高了DNA序列拼接质量。以往的拼接技术都是基于单线程进行的序列拼接,而Idba采用了多线程技术来进行序列拼接,提高了序列拼接的效率。通过对真实数据进行测试,结果表明,Idba技术能够得到更长的ccmtig长度。对同一组基因组数据进行拼接质量测试时,Velvet得到N50大小为19284,而Idba得到的N50大小为63218,其长度将是Velvet的近6倍;在内存消耗上,Velvet消耗内存为1641M,Idba仅仅消耗内存360M。可见,Idba对de Bmijn图的构建优化效果显著45。, 序列拼接的

13、过程1)假设read集合为>S=将其中每一个read划分为若干个连续碱基组成的k-mer集合尺=為,.人,每一个k-mer为图中的一个顶点。其中,k-rtier的截取规则为:选取一个长度为的read,先以read的左端为起始位置,截取长度为&的碱基,再将起始位置向右移动一位,截取第二个长度为的碱基,直至到达read的右端,这些k-mer组成了 de Bruijn图中的顶点。为了防止DNA序列中互补双链中截取的k-mer相同,在此,取<4:值为奇数;2)对于k-mer集合=,若存在两个k-mer,其中、的后A:-l个碱基与的前A:-l个碱基相同,那么,k、之间必然存在一条由:

14、1指向的一条边。根据k-mer之间的重叠信息,每一条read表示图中的一条路径,由此构建一个以k-mer为基本单位的有向路径;3)根据上述得到的有向路径集合,以及k-mer之间的重叠信息,寻找一条经过所有边一次且仅一次的路径,即将拼接问题转化为图论中寻找欧拉路径求解。根据DNA序列中的pair-end信息,对欧拉路径进行解親,最终找到一条近似的连续的DNA序列46。整个de Bruijn图的构建过程如图2-1描述:更新决策表策略决策表中存了 read 的如下信息：readID，方向 orientation，刚刚进入决策表时当前 k-mer 出现在该 read 中的位置 first_pos，当前

15、 k-mer 出现在该 read 中的位置cur_pos（若当前 k-mer 未出现在该 read 中，则 cur_pos 为 0），最后出现在该 read中的 k-mer 的出现位置 lastappearpos，拼接过程中所用到的 k-mer 在该 read 中的出现次数 k-merappeartimes，未出现次数 k-merunappeartimes，未出现连续块数k-merunapperblocks，及该 read 的状态 status，删除标记 delsign，锁定标记 locked。其中 read 的状态有如下 4 种：拼接状态，终断状态，成功状态，失败状态。以上信息在拼接过程中都

16、会用到，通过更新上述信息，可以知道一个 read 是否可参与k-mer 评分；如果是，该 read 所占权值是多少。contig 的构建过程contig 的构建过程主要有如下几步：步骤一，选取初始 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 上

17、会有初始 k-mer 的 k 个碱基，如图 2-4 所示。这样，就会有至少min_read_num 条 read 参与拼接，这些 read 会影响到初始 k-mer 的扩展过程称图 2-4 为 read 拼接过程图。以后每当有 read 开始参与拼接时，就要将这些read 加入到该图中，每当有 read 结束拼接时，就要将这些 read 从该图中删掉。此时，初始 k-mer 为当前 k-mer，初始 k-mer 出现在某条 read 中的 pos 为该 read 的当前 pos，记为 cur_pos，现在所有 read 的当前 pos 为 1。至此，步骤一结束，接下来进行步骤二。步骤二，选取后

18、继 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 有四个，如图 2-5 所示。接上面的例子，选定的后继 k-mer 是 ACCG，此时，contig上多了一个碱基，

19、当前 k-mer变为 ACCG，read1 至 read4 的当前 pos 变为 2，read5 处于拼接中断状态。由于当前 k-mer 可能出现在其它 read 上，且出现的 pos 为 1，所以要让这些 read 开始参与拼接，如图 2-6 中的 read6，read7，read8。接下来，需要反复重复步骤二，直到找不到合适的后继 k-mer 为止。虽然每次必定能产生四个候选的后继 k-mer，但当出现如下情形时，没有合适的后继 k-mer可供选择：1) 对于任意一个候选后继 k-mer，它不出现在任何当前参与拼接的 read 的合适位置，即不管我们选择哪个后继 k-mer 都将导致 re

20、ad 拼接过程图中所有 read 处于拼接中断状态；2) 对于任意一个候选后继 k-mer，对于任意一条 read 库中的 read，该 k-mer 不出现在该 read 上 pos 为 1 的位置，即不管我们选择哪个后继 k-mer，都不会使新的一批 read 开始参与拼接。当上述情形同时出现时，我们将找不到合适的后继 k-mer，此时一条 contig拼接结束，如果该 contig 足够长（长度不小于 100bp），我们会保存该 contig，否则我们将丢弃该 contig。如果我们丢弃了该 contig，我们要把那些处于成功拼接状态的 read 从新加入到 read 库中，让它们继续参与

21、后面的拼接。如果没有找到合适的后继 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 选取失败，整个拼

22、接过程结束。de Bruijn 图的建立（a）read 的选取。数据文件中不仅保存了 read 的碱基，还保存了 read 的质量值。该质量值能反映 read 碱基的正确率。质量值越大，read 碱基的正确率就越高。故我们选择那些质量较高的 read。显然我们要事先指定一个质量阈值。文件中有的 read 包含未知碱基，这是由于测序过程中有些碱基没有被测出导致的。未知碱基用 N 表示，故我们要过滤掉那些含 N 的 read。另外，有些 read 含碱基 A 较多，甚至全 A。根据测序仪性质，这些 read 的错误率较高，故我们过滤掉了那些 A 的含量超过 90%的 read。最终我们用所有符合上

23、述条件的 read 构建 contig。（b）预读数据文件，初始化 de bruijn图。依次读取每一个 read，判断该 read是否满足参与拼接的条件。若不满足，跳过该 read；若满足，生成该 read 上所有k-mer（共 25 个 k-mer），统计这些 k-mer 的出现次数，填写 k-mer 结构中的 num字段，该过程如图 3-6 所示。统计 k-mer 的出现次数时，必须要找到该 k-mer 所对应的 k-mer 结构。这就需要计算该 k-mer 的哈希值，根据哈希表的首地址找到 k-mer结构的地址，然后将其中的 num 字段加 1，该过程如图 3-7 所示。（c）遍历 d

24、e Bruijn图，根据第（b）步中统计的 k-mer 个数，申请 readID_pos数组所需要的内存空间。实际上，整个系统中所有大型数组都是动态开辟的，因为事先我们不知道数组的大小，而且我们也不想申请一块足够大的内存来使用。另外，函数的栈空间非常有限，大型数组是绝对不能放到栈上的，为避免过多地使用全局变量，和静态变量，我们选择了动态申请内存的方法。再读数据文件，read；若满足，生成该 read 上所有 k-mer，填写每个 k- mer 的 readID_pos 数组，此时获得 k-mer 结构地址的方法与步骤（b）中的一样，如图 3-7 所示。然后填写readID_pos 数组中的第

25、cur 行，填好后更新 cur 的值，即把 cur 的值加 1，该过程如图 3-8 所示。当某个 readID_pos 数组填满后，其 cur 的值应该和 num的值一样，表示当前 readID_pos 数组中所有元素都是有效的。有关 de Bruijn 图的基本操作首先介绍k-mer碱基与k-mer哈希值的相互转换。将k-mer碱基转化成整数时，直接将碱基替换成 2 位二进制整数。例如，k=5，要转换的 k-mer 碱基为 ACTGT，则转换之后结果是 0001111011。显然，这种转换规则是可逆的，我们只要自左至右扫描哈希值的每一位，就可得到该 k-mer 的碱基序列接下来介绍如何查找某

26、 k-mer 出现在某 read 的哪些位置。给定 k-mer 的碱基序列以及一个 readID，可在 de Bruijn图中找到该 k-mer 出现在该 read 的哪些位置。若再给定一个 pos，可以判断该 k-mer 是否出现在该 read 的 pos 置。步骤如下：首先获得该 k-mer 所对应的 k-mer 结构地址（这一步可参考图 3-8），进而获得readID_pos 数组的首地址和数组大小。然后进行二分查找，在 readID_pos 数组中查找给定的 readID。若查找失败，则该 k-mer 不出现在该 read 中。若查找成功，由于给定的 readID 可能在数组中出现多次

27、，故还需向上搜索几步，找到 readID_pos数组中第一个出现 readID 的位置，再从该位置向下遍历，直到遇上的 readID 与给定的 readID 不同。此时我们已经获得了该 k-mer 出现在该 read 的所有位置，同时也可判断该 k-mer 是否出现在该 read 的 pos 位置最后介绍 read 的删除、恢复操作。给定一条 read 的碱基序列以及 readID，要在 de Bruijn图中删除该 read，既要找到所有该 readID 出现的地方，并把那一行的删除标记置 1。看上去好像要遍历整个 de Bruijn图，其实可以在 2O k时间内完成，步骤如下：依次生成该

28、read 的每一个 k-mer，找出该 k-mer 出现在该 read 中的所有位置，修改其删除标记，同时修改该 k-mer 的 cur 字段。后继 k-mer 选取过程后继 k-mer 的选取是一个非常关键的步骤，如果后继 k-mer 选取的合理，就能拼接成质量较高的 contig，其长度相对较长，和参考基因组匹配的较好，并且存在大量 read 拼接成功。反之，如果后继 k-mer 选的不合理，会导致 contig 很快拼接结束，其长度非常短（小于 100bp），并且导致大量 read 拼接失败。由于基因组中存在大量重复片段，故较短的 contig 和参考基因组会匹配的较好，但由于大量 read拼接失败，所以该 contig的质量并不高。理想情况下，在任意选定初始 k-mer 后，在向左右扩展时，会逐渐拼接成整个基因组。但由于测试数据存在错误，并且基因组中存在大量重复片段，导致可能选取了错误的后继 k-mer，进而导致拼接出的contig 长度较短。由此可见，后继 k-mer 的选取是非常关键的。在这一阶段并不需要特殊的数据结构，只需要 4 个存储候选后继 k-mer 碱基的字符数组，以及一个候选后继 k-mer 得分数组。我们用长度为 2 的指针数组存储后继 k-mer。数组的 0 号位置指向正向 k-mer。1 号位置指向反向 k-mer