基于 de Bruijn 图的算法概述

《基于 de Bruijn 图的算法概述》由会员分享，可在线阅读，更多相关《基于 de Bruijn 图的算法概述（27页珍藏版）》请在装配图网上搜索。

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、之间必然存在一条由:1指向的一条边。根据k-

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

15、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 上会有初始 k-mer 的

17、 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。至此，步骤一结束，接下来进行步骤二。步骤二，选取后继 k-mer。接着要选

18、取当前 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上多了一个碱基，当前 k-mer变为 A

19、CCG，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 都将导致 read 拼接过程图中所有

20、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 选取失败，整个拼接过程结束。de Bru

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

23、 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）遍历 de Bruijn图，根据

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

25、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 的碱基序列接下来介绍如何查找某 k-mer 出现在某

26、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时间内完成，步骤如下：依次生成该 read 的每一个 k-

28、mer，找出该 k-mer 出现在该 read 中的所有位置，修改其删除标记，同时修改该 k-mer 的 cur 字段。后继 k-mer 选取过程后继 k-mer 的选取是一个非常关键的步骤，如果后继 k-mer 选取的合理，就能拼接成质量较高的 contig，其长度相对较长，和参考基因组匹配的较好，并且存在大量 read 拼接成功。反之，如果后继 k-mer 选的不合理，会导致 contig 很快拼接结束，其长度非常短（小于 100bp），并且导致大量 read 拼接失败。由于基因组中存在大量重复片段，故较短的 contig 和参考基因组会匹配的较好，但由于大量 read拼接失败，所以该 c

29、ontig的质量并不高。理想情况下，在任意选定初始 k-mer 后，在向左右扩展时，会逐渐拼接成整个基因组。但由于测试数据存在错误，并且基因组中存在大量重复片段，导致可能选取了错误的后继 k-mer，进而导致拼接出的contig 长度较短。由此可见，后继 k-mer 的选取是非常关键的。在这一阶段并不需要特殊的数据结构，只需要 4 个存储候选后继 k-mer 碱基的字符数组，以及一个候选后继 k-mer 得分数组。我们用长度为 2 的指针数组存储后继 k-mer。数组的 0 号位置指向正向 k-mer。1 号位置指向反向 k-mer。真正的 k-mer都存在 de bruijn图中。存储结构如

30、图 3-15 所示。由于反向 k-mer 的存在仅仅是为了避免生成反向 read 的碱基序列，故我们并不给反向 k-mer 一个得分，但评分过程中还需要反向 k-mer 的参与，因为评分时要遍历决策表，而决策表中有反向的 read，正如前面所述，反向 read 要想发挥作用，必须依靠反向 k-mer。后继 k-mer 要使决策表中较多的 read 继续参与拼接，较少的 read 拼接中断。这样大量的 read 就有成功拼接的希望，而那少量拼接失败的 read 可能存在数据错误。我们应该使正确率高的 read 拼接成功，使正确率较低的 read 拼接失败。这样有利于产生高质量的 contig。其

31、次，如果有些 read 马上拼接成功了，我们应该让这些 read 尽快成功拼接，即我们要是后继 k-mer 尽量出现在这些 read 上。实际上这些 read 就是被锁定的 read，所以有时并不是所有 read 都参与评分。使用锁定 read 的好处是可以保证 contig上每一个碱基都被一个成功拼接的 read 所覆盖。后继 k-mer 的选取过程示意图如图 3-16 所示。首先取当前正向 k-mer 的后 k-1个碱基，分别接上 ACGT，产生四个候选的后继 k-mer，如图 3-17 所示。然后对每一个候选 k-mer，要在 de bruijn图中找到该 k-mer 所对应的 k-me

32、r 结构，这一步详细过程如图 3-7 所示。这时我们也得到了该 k-mer 所对应的 readID_pos 数组。若de bruijn 图中该 k-mer 为空，则不存在相应的 readID_pos 数组，令该 k-mer 得分为0，继续为下一条候选 k-mer 评分。取得得分最高的k-mer作为后继k-mer，若所有k-mer得分为0，则置后继k-mer为空。若有多个 k-mer 得分相同且最高，说明在当前来看，最优路径有多条，但最终结果可能差距较大。我们无法预知哪条路径更优，在此，我们采用的策略是立即停止该条contig的拼接，直接置后继k-mer为空。这样做的好处是可以提高contig的准确度，因为如果选择了错误的路径将使 contig质量变得很差。: