基于deBruijn图的算法概述Word格式.docx

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

2、段中的一个配对 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

3、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。实际上

4、，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

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

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

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

8、支区域。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图的结构,使得其在内存消耗上明显降低,提高了算法的拼接效率。设置k-me

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

10、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-mer之间的重叠信息

11、,每一条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 出现在该 read 中的位置c

12、ur_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 种：拼接状态，终断状态，成功状态，失败状态。以上信息在拼接过程中都会用到，通过更新上述信息，可以知道一个 re

13、ad 是否可参与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 的 k 个碱基，如图

14、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。接着要选取当前 k-mer

15、的后继 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变为 ACCG，read1

16、至 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 拼接过程图中所有 read 处于拼接中

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

18、r，则转到步骤一，开始下一条 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 Bruijn 图的建立（a

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

20、）预读数据文件，初始化 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图，根据第（b）步中统计的

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

22、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 出现在某 read 的哪些位置

23、。给定 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 可能在数组中出现多次，故还需向上搜索几步，找到 readID_p

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

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

26、况下，在任意选定初始 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图中。存储结构如图 3-15 所示。由于反向

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

28、接成功了，我们应该让这些 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-mer 结构，这一步详细过程如图 3-7 所示。这时我们也得到了该 k-mer 所对应的 readID_pos 数组。若de bruijn 图中该 k-mer 为空，则不存在相应的 readID_pos 数组，令该 k-mer 得分为0，继续为下一条候选 k-mer 评分。取得得