基于 de Bruijn 图的算法概述.docx
《基于 de Bruijn 图的算法概述.docx》由会员分享，可在线阅读，更多相关《基于 de Bruijn 图的算法概述.docx（21页珍藏版）》请在冰点文库上搜索。
基于deBruijn图的算法概述
基于deBruijn图的算法概述
deBruijn图简介
传统的Sanger测序的reads较长（1000bp），数据量较少，精度较高，所有的
组装算法都利用reads之间的重叠，通过公共路径的方法解决拼接问题。
而新一代
测序产生的数据read更短、覆盖度更高、序列精度较低，为此这种―read为中心‖
的方法面临海量计算的困境，似乎不可能找到恰当的启发式方法来处理大量的重
叠。
deBruijn图框架为处理高覆盖、短序列提供了很好思路，该框架借鉴了Pevzner
和Waterman等人针对传统的长reads提出的欧拉遍历方法
[37,38]
，并在此基础上针
对新一代测序数据的特点进行了改进要想以较低的成本快速得到某个新物种的DNA分子碱基序列，就要依靠新一
代的测序技术和从头测序拼接组装算法。
目前新一
代测序数据用于从头测序的短序列拼接组装算法普遍采用deBruijn图数据结构。
在deBruijn图上，每一个k-mer都构成图的节点，如果两个k-mer在某一read中
相邻，那么这两个节点之间就有一条边。
reads集合中的每个read都对它所含的节点和边加权，这样reads集合产生一个节点和边都具有权值的deBruijn图。
在存储
每一个k-mer时，往往要建一个无冲突的哈希表，以加快查找速度。
而建立哈希表
可能会消耗更多的内存。
但是，由于每个k-mer在哈希表中只存储一次，不管该
k-mer在read中出现了多少次，所以实际消耗的内存小于存储所有read所需要的
空间。
另外，基因组中的重复片段会在deBruijn图中产生环路。
环路将在遍历de
Bruijn图时产生障碍。
目前的研究主要面临两个问题，一个是基因组中存在大量重
复片段，一个是测序错误。
这两个问题相互影响，使问题变的更加复杂。
本文通
过仔细分析这两个问题，来改进以前基于deBruijn图的算法，提出一种新的de
Bruijn图，并且引入了决策表的概念，通过决策表里的信息来选取后继k-mer，并在适当的时候更新决策表。
1基因组中存在大量重复片段
重复片段问题可用如下方法解决：
通过比对，可先将重复片段隔离开来，较
高的覆盖度有利于重复片段的隔离，但是，较多的测序错误将不利于该过程的进
行。
因为错误的存在，严格的比对将导致一些重复片段未被发现，而非严格的比
对会把一些不是重复片段的区域隔离开来，这不是本文所希望的。
如果重复片段
比read长，可利用paredendread来解决；如果重复片段比read短，那么该read
又被称为spanner，一个spanner就是一个重复片段两端再加几个碱基组成。
利用
spanner解决重复片段问题需要如下两个信息：
一是重复片段两端配对的read，这
两个read必须不相同；二是重复片段中的一个配对read，只要知道一个即可，另
一个配对read可以不在重复片段中
2测序过程中可能出现错误
现在主要有两种
纠错方法，一种基于多重比对，通过将多个read放在一起比对来发现错误，如图
1-2所示。
通过图中4条read比对，可发现read3中的一个碱基错误（read3的第5个碱
基），该方法在overlap过程中比较常用，而在deBruijn图中，所使用的纠错方法
是：
若当前k-mer在一条read中连续未出现恰好k次，可以认为该read中存在一
个碱基错误。
2基于deBruijn图算法的一般步骤
1）确定k值，建立deBruijn图。
这时需要扫描所有read数据，将每一个长
为L的read拆分成L-k+1个kmer，并用所有read的所有k-mer来累加，建立节点
和边都加权的deBruijn图；
2）化简deBruijn图，连续线性延伸节点合并为单一节点，产生一些碱基序列
更长的节点；
3）错误校正，删去由于测序错误产生的尖端和泡状结构；
4）通过read的配对末端（pair-end）、环化配对（mate-pair）信息伸展或者删去一
些环；
5）依据环上节点和边的权值（覆盖深度信息）进一步伸展或者删去一些环；
6）遍历deBruijn图产生contig。
实际上，deBruijn图是一种特殊的加权
图，不仅图的结点上有权值，而且图的边上也有权值。
化简deBruijn图是非常关
键的一个步骤，通过对deBruijn图化简，可降低算法的时间复杂性以及空间复杂
性，同时可以保证错误校正顺进行
拼接总体思路
假设所有满足上述条件
（1）的read都已经存到了read库中，下面就用这些
read来构建contig。
给定k值后，长度为k的一个DNA片段称为一个k-mer。
一
般地，k要小于每条read的长度L，故每条read中含有k-mer数量为L-k+1。
一个
k-mer的第一个碱基在一个read中出现的位置记为pos，pos的值从1开始，最大
为L-k+1，如图2-2所示。
选定一个初始k-mer后，通过对该kmer不断扩展，来得到一条contig。
一个
k-mer上有k个碱基，而碱基共有四种，故扩展下一个碱基有四种选择，这样就会形成一个四叉树，如图2-3所示。
显然，这个四叉树的深度是无限的，任何一个子树的深度也是无限的。
算法
中要设定终止条件，不能让它无限地扩展下去。
实际上，该树中任何一条有限的
路径都可能成为一条contig，每条contig都可以使一些read成为它的子串，所以
可以用read库中的read来评价contig的好坏。
虽然可以用无信息搜索的方法来拼
接contig，但现在的问题是，有些contig长达几万bp，这样算法要搜索上万层，
搜索空间过大，以至于不能在有效的时间内完成。
故本文采用启发式搜索，来减
少时间开销。
在一条contig开始拼接前，需要根据一定策略，选定树中一个初始k-mer，接
下来就可以在以该k-mer为根结点的子树上开始搜索。
搜索时采用贪心策略，每一
步选择在当前看来最优的后继k-mer，直到满足事先设定的终止条件，结束一条
contig的拼接，接着开始下一条contig的拼接。
直到没有合适的初始k-mer可供选
择，整个拼接过程结束。
由于选定初始k-mer后，可以向该k-mer的两端分别扩展，故初始k-mer选取
的好坏对拼接结果影响不大。
故该问题的关键是选取后继k-mer。
后继k-mer如果
选择的好，contig会拼接得较长，会有较多的read成功参与拼接；后继k-mer如果
选择的差，contig会拼接得较短，会有较少的read成功参与拼接。
基于deBruijn图的序列拼接技术分析
Idba拼接技术
Velvet、Soapdenovo等在处理无错误序列、高覆盖度的序列拼接问题时,能
够表现很好的性能,但是,由于这些技术在拼接过程中以k-mer为基本单位,就
不可避免的会产生很多重叠单元,这使得拼接面临着错误位置拼接、顶点缺失和
覆盖度低等问题。
由于一些错误read的存在,产生了大量的分支区域。
k-mer长
度越小,分支问题越严重,k-mer长度越大,出现的重叠区域变得越少,这直接
影响了拼接的质量。
正确的选择k-mer的大小成为影响拼接质量的一个关键因
素
Idba的主要特点是:
按弃了使用固定的k-mer长度构建deBruijn图的方法采用一个变化的k-mer长度完成read切割过程。
首先,它设置k-mer长度的变化
区域,其中A为当前k-mer的长度,采用重复迭代算法来计算每一个k-mer的长度;
在构建图之前,该技术对低覆盖度的k-mer进行了预处理,去除覆盖度低的k-mer
顶点,从而简化了deBmijn图的结构,使得其在内存消耗上明显降低,提高了算
法的拼接效率。
设置k-mer阈值的方法,成功解决了deBmijn图中路径多分支问题,提高了
DNA序列拼接质量。
以往的拼接技术都是基于单线程进行的序列拼接,而Idba采
用了多线程技术来进行序列拼接,提高了序列拼接的效率。
通过对真实数据进行测试,结果表明,Idba技术能够得到更长的ccmtig长
度。
对同一组基因组数据进行拼接质量测试时,Velvet得到N50大小为19284,而
Idba得到的N50大小为63218,其长度将是Velvet的近6倍;在内存消耗上,
Velvet消耗内存为1641M,Idba仅仅消耗内存360M。
可见,Idba对deBmijn图的
构建优化效果显著[45]。
序列拼接的过程
1）假设read集合为>S=
将其中每一个read划分为若干个连续
碱基组成的k-mer集合尺={^]為,...人},每一个k-mer为图中的一个顶点。
其中,
k-rtier的截取规则为:
选取一个长度为的read,先以read的左端为起始位置,截
取长度为&的碱基,再将起始位置向右移动一位,截取第二个长度为的碱基,直
至到达read的右端,这些k-mer组成了deBruijn图中的顶点。
为了防止DNA序
列中互补双链中截取的k-mer相同,在此,取<4:
值为奇数;
2）对于k-mer集合=
若存在两个k-mer,其中、的后A:
-l个
碱基与的前A:
-l个碱基相同,那么,k、、、之间必然存在一条由^:
1指向的一
条边。
根据k-mer之间的重叠信息,每一条read表示图中的一条路径,由此构建
一个以k-mer为基本单位的有向路径;
3）根据上述得到的有向路径集合,以及k-mer之间的重叠信息,寻找一条经
过所有边一次且仅一次的路径,即将拼接问题转化为图论中寻找欧拉路径求解。
根据DNA序列中的pair-end信息,对欧拉路径进行解親,最终找到一条近似的连
续的DNA序列[46]。
整个deBruijn图的构建过程如图2-1描述:
更新决策表策略
决策表中存了read的如下信息：
readID，方向orientation，刚刚进入决策表时
当前k-mer出现在该read中的位置first_pos，当前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种：
拼接状态，终断状态，成功状态，失败状态。
以上信息在拼接过程中都会用到，通过更新上述信息，可以知道一个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的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。
接着要选取当前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
变为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都将导致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选取失败，整个拼接过程结束。
deBruijn图的建立
（a）read的选取。
数据文件中不仅保存了read的碱基，还保存了read的质量
值。
该质量值能反映read碱基的正确率。
质量值越大，read碱基的正确率就越高。
故我们选择那些质量较高的read。
显然我们要事先指定一个质量阈值。
文件中有
的read包含未知碱基，这是由于测序过程中有些碱基没有被测出导致的。
未知碱
基用N表示，故我们要过滤掉那些含N的read。
另外，有些read含碱基A较多，
甚至全A。
根据测序仪性质，这些read的错误率较高，故我们过滤掉了那些A的
含量超过90%的read。
最终我们用所有符合上述条件的read构建contig。
（b）预读数据文件，初始化debruijn图。
依次读取每一个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）遍历deBruijn图，根据第（b）步中统计的k-mer个数，申请readID_pos
数组所需要的内存空间。
实际上，整个系统中所有大型数组都是动态开辟的，因
为事先我们不知道数组的大小，而且我们也不想申请一块足够大的内存来使用。
另外，函数的栈空间非常有限，大型数组是绝对不能放到栈上的，为避免过多地
使用全局变量，和静态变量，我们选择了动态申请内存的方法。
再读数据文件，read；若满足，生成该read上所有k-mer，填写每个k-mer的readID_pos数组，此
时获得k-mer结构地址的方法与步骤（b）中的一样，如图3-7所示。
然后填写
readID_pos数组中的第cur行，填好后更新cur的值，即把cur的值加1，该过程
如图3-8所示。
当某个readID_pos数组填满后，其cur的值应该和num的值一样，
表示当前readID_pos数组中所有元素都是有效的。
有关deBruijn图的基本操作
首先介绍k-mer碱基与k-mer哈希值的相互转换。
将k-mer碱基转化成整数时，
直接将碱基替换成2位二进制整数。
例如，k=5，要转换的k-mer碱基为ACTGT，
则转换之后结果是0001111011。
显然，这种转换规则是可逆的，我们只要自左至
右扫描哈希值的每一位，就可得到该k-mer的碱基序列
接下来介绍如何查找某k-mer出现在某read的哪些位置。
给定k-mer的碱基
序列以及一个readID，可在deBruijn图中找到该k-mer出现在该read的哪些位置。
若再给定一个pos，可以判断该k-mer是否出现在该read的pos置。
步骤如下：
首
先获得该k-mer所对应的k-mer结构地址（这一步可参考图3-8），进而获得
readID_pos数组的首地址和数组大小。
然后进行二分查找，在readID_pos数组中
查找给定的readID。
若查找失败，则该k-mer不出现在该read中。
若查找成功，
由于给定的readID可能在数组中出现多次，故还需向上搜索几步，找到readID_pos
数组中第一个出现readID的位置，再从该位置向下遍历，直到遇上的readID与给
定的readID不同。
此时我们已经获得了该k-mer出现在该read的所有位置，同时
也可判断该k-mer是否出现在该read的pos位置
最后介绍read的删除、恢复操作。
给定一条read的碱基序列以及readID，要
在deBruijn图中删除该read，既要找到所有该readID出现的地方，并把那一行的
删除标记置1。
看上去好像要遍历整个deBruijn图，其实可以在{}
2
Ok时间内完成，
步骤如下：
依次生成该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。
真正的k-mer
都存在debruijn图中。
存储结构如图3-15所示。
由于反向k-mer的存在仅仅是为了避免生成反向read的碱基序列，故我们并
不给反向k-mer一个得分，但评分过程中还需要反向k-mer的参与，因为评分时要
遍历决策表，而决策表中有反向的read，正如前面所述，反向read要想发挥作用，
必须依靠反向k-mer。
后继k-mer要使决策表中较多
的read继续参与拼接，较少的read拼接中断。
这样大量的read就有成功拼接的希
望，而那少量拼接失败的read可能存在数据错误。
我们应该使正确率高的read拼
接成功，使正确率较低的read拼接失败。
这样有利于产生高质量的contig。
其次，
如果有些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，要在debruijn图中找到该k-mer所对应的k-mer结构，这一步
详细过程如图3-7所示。
这时我们也得到了该k-mer所对应的readID_pos数组。
若
debruijn图中该k-mer为空，则不存在相应的readID_pos数组，令该k-mer得分为
0，继续为下一条候选k-mer评分。
升级会员 