CN103258145B - 一种基于De Bruijn图的并行基因拼接方法 - Google Patents

Abstract

本发明涉及基因测序技术领域，提供了一种基于De Bruijn图的并行基因拼接方法，包括：S1、并行构建分布式De Bruijn图；S2、剔除错误路径；S3、基于深度图遍历方法对De Bruijn图进行化简；S4、合并contig，生成scaffold；S5、输出scaffold。本发明基于集群***，并行构造De Bruijn图，解决了大基因组拼接时由于其数据量太大，传统单机串行的基因拼接算法无法构图和无法进一步处理的问题；同时，在化简过程中，进行基于深度图遍历的并行化简，图化简过程简单，并行度高，拼接速度快。

Description

一种基于De Bruijn图的并行基因拼接方法

技术领域

本发明涉及基因测序技术领域，特别是涉及一种基于DeBruijn图的并行基因拼接方法。

背景技术

基因测序是现代生物信息领域中最重要的问题之一。随着现代生物的发展，基因测序已经越来越广泛地应用于社会的各个领域中，比如基因诊断、基因治疗、药物设计等。基因测序中很重要的一步就是基因拼接。

随着基因测序的广泛应用，一方面，它需要对大量大基因组的生物进行测序。在对大基因组进行测序时，其数据量非常大；另一方面，这也要求基因拼接的算法越来越快。

基因序列拼接算法主要有两类。第一类是基于overlap图的算法。在overlap图中，每一个短序列read（DNA分子被随机打断成很多小片段，每个小片段被称为read）都被当作一个顶点，如果两个read之间存在重叠且重叠的长度超过一定阈值，那么就有一条有向边相连。因此，序列拼接问题转换为在overlap图中寻找一条经过每个顶点的Hamilton路径，这是NP-Hard问题。第二类是基于DeBruijn图的算法。在DeBruijn图中，每一个read被切分成长度为k的小片段，称为k-mer。每一个k-mer为一个顶点。如果存在read，使两个k-mers相邻且重叠k-1个字符，那么它们之间存在一条有向边。这样每个read被映射成图中的一条路径。因此，序列拼接问题变成了在DeBruijn图中寻找一条包含所有read的路径。测序仪在测序时会引入错误，下一代高通量测序仪错误率在1%左右，同时原始序列中存在不同长度的重复片段，这两个问题使得序列拼接问题更加复杂。

早期基于Sanger测序法（自动化桑格测序法）得到的序列片段，长度可以达到1000BP(碱基对)。基于overlap图的拼接算法比较有效，但测序成本比较高。比如通过第一代测序技术完成的人类基因组计划，花费了30亿美元，耗时三年。当第二代测序技术（又称为下一代测序技术）如Solexa、454、SOLID技术出现后，基因测序才开始真正进入大规模应用。第二代测序技术有三个显著的特点，高通量，短序列，高覆盖。高通量，是测序仪一次可以同时测定大量的read序列，极大降低了测序成本。短序列，是序列长度一般在25-500base之间。高覆盖，是因为序列短，为了保证信息的完整性，需要极大的提高DNA的覆盖度（即coverage）。但随着覆盖度的提高，read数量成倍的增加，如果继续采用基于overlap图的算法，图的规模也会成倍的增长。如果采用基于DeBruijn图的算法，图的规模与DNA长度呈线性关系。对同一基因组而言，其规模几乎不变。因此，面对第二代基因测序技术的大规模应用，DeBruijn图的基因拼接算法有很大的优势。其相关的算法有Euler、ALLPATHS、Velvet、IDBA、SOAPdenovo、分布式的ABySS和分布式的YAGA。其中Euler、ALLPATHS、Velvet、IDBA算法是串行算法，适用于小数据集的拼接，SOAPdenovo是基于SMP大型机的多线程拼接算法，可以拼接大型数据集如人类基因组，但最快拼接时间也需要40多个小时。

在基于DeBruijn图的拼接过程中，最为消耗内存的步骤是构建DeBruijn图，尤其是对大基因组生物，其构造的DeBruijn图非常大，现有的单机串行的拼接算法，无法完成构图。同时，对DeBruijn图的化简需要占用最多的处理时间。现有的并行拼接算法，其拼接速度仍然无法达到大规模应用的要求，主要的难点就是化简过程的并行度不高。

综上所述，现有技术中，传统单机串行的基因拼接算法无法对大基因组的海量数据进行拼接，而现行的并行基因拼接算法不能快速对大基因组进行拼接。

发明内容

本发明提供一种基于DeBruijn图的并行基因拼接方法，旨在解决现有技术中传统单机串行的基因拼接算法无法对大基因组的海量数据进行拼接，现行的并行基因拼接算法不能快速对大基因组进行拼接的技术问题。

本发明采用如下技术方案：

一种基于DeBruijn图的并行基因拼接方法，包括：

S1、并行构建分布式DeBruijn图；

S2、剔除错误路径；

S3、基于深度图遍历方法对DeBruijn图进行化简；

S4、合并contig，生成scaffold；

S5、输出scaffold。

优选地，所述步骤S1具体包括：

S11、所有的处理器并行读取原始的短序列文件，每个处理器读取短序列文件的一部分；

S12、每个处理器对读取的短序列进行切分，每个短序列被切分成长度为k的小片段k-mer，获取所述k-mer，把每个k-mer构造成一个节点；

S13、按照哈希规则，对所有的节点进行重分布。

优选地，所述步骤S2具体为：

对步骤S1构建的所述分布式DeBruijn图进行并行去错。

优选地，在所述步骤S2中，对步骤S1构建的所述分布式DeBruijn图中的Tips、Bubble和SpuriousLink三种结构进行剔除。

优选地，所述步骤S3具体包括：

S31、初始化，将DeBruijn图分布存储在每个处理器的locationMap中，同时新建一个subGraphMap；

S32、所有处理器并行从本地的locationMap中的端节点出发遍历单链；

S33、将所有的单链进行化简；

S34、判断是否有未访问的端节点，若有，返回执行步骤S32，否则，执行步骤S35；

S35、图化简完成，生成contig。

优选地，所述步骤S32具体为：

处理器从端节点出发访问单链节点；如果所述链节点是第一次被访问，则记录本次访问的端节点ID号；如果处理器访问了单链的所有链节点直到另一个端节点都没有遇到冲突，则将所述单链分配给所述处理器；

如果处理器在访问第一个链节点时发现该单链已被访问过，则回退，继续访问下一条单链；

如果处理器访问的非第一个单链节点的链节点同时也被另外一个处理器访问，则比较自己端节点ID号与对方端节点ID号的大小，如果自己端节点ID号小，则继续访问，反之，则放弃，继续访问下一条单链。

优选地，所述步骤S33具体为：

将单链节点全部删除，节点信息合并到端节点中，并将端节点***subGraphMap中。

优选地，所述步骤S4具体为：

利用测序获得的pair-end信息，合并步骤S3得到的contig，生成scaffold，所述scaffold为更长的contig。

本发明基于集群***，并行构造DeBruijn图，解决了大基因组拼接时由于其数据量太大，传统单机串行的基因拼接算法无法构图和无法进一步处理的问题；同时，在化简过程中，进行基于深度图遍历的并行化简，图化简过程简单，并行度高，拼接速度快。

附图说明

图1为本发明实施例基于DeBruijn图的并行基因拼接方法流程图；

图2为本发明实施例中并行构建分布式DeBruijn图的方法流程图；

图3为本发明实施例中基于深度图遍历方法对DeBruijn图进行化简的方法流程图；

图4为本发明实施例中深度图遍历寻找单链过程的示意图；

图5为本发明实施例中Yeast数据集的时间开销图；

图6为本发明实施例中C.elegans数据集的时间开销图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

如图1所示，本发明实施例提供了一种基于DeBruijn图的并行基因拼接方法，该方法包括下述步骤：

步骤S1：并行构建分布式DeBruijn图。

本实施例中，通过如下步骤进行分布式DeBruijn图的并行构建，如图2所示：

步骤S11、初始化，所有的处理器并行读取原始的read短序列文件，每个处理器读取read短序列文件的一部分。根据处理器的编号，根据read文件的大小，每个处理器读取文件的不同块。假设文件的长度为L，有n个处理器，那么第i个处理器读取从((i-1)*L)/n到(i*L)/n。

步骤S12、每个处理器对读取的read短序列进行切分，每个read短序列被切分成长度为k的小片段k-mer，获取k-mer，把每个k-mer构造成一个节点。

步骤S13、按照哈希规则（处理器号=hash(ID)），对所有的节点进行重分布。该步骤完成后，任何处理器都能根据相同规则直接定位任何一个节点，这为后续处理中单链遍历和图的化简提供了便利。

步骤S2：剔除错误路径。

由于测序时引入了错误，DeBruijn图构建完成后，图中有大量错误路径。错误有三类：Tips、Bubble和SpuriousLink。因此，本步骤对步骤S1构造的分布式DeBruijn图进行并行去错，主要是对构造的分布式DeBruijn图中的Tips、Bubble和SpuriousLink三种结构进行剔除。去错后，图可看作由error-free数据构成的图，这时的图非常稀疏，存在大量的单链。下一步需进行单链化简。

步骤S3：基于深度图遍历方法对DeBruijn图进行化简。

化简的目的就是找出所有的单链。本实施例中，通过如下步骤对DeBruijn图进行化简，如图3所示：

步骤S31、初始化，将DeBruijn图分布存储在每个处理器的locationMap中，同时新建一个subGraphMap。

每个处理器上有两个Map。一个是locationMap，存储图化简前的节点。另一个是subGraphMap，存储图化简后的节点。初始化时，subGraphMap为空，在化简过程中，会不断地***单链合并后剩下的节点。

步骤S32、所有处理器并行从本地的locationMap中的端节点出发遍历单链。

本步骤中，如果不加限制，那么遍历单链时存在如下三种情况：同一条链，同一个处理器先后访问；同一条链，两个处理器先后访问；同一条链，两个处理器同时访问。为了减少节点的重复访问，节点中增加一个endNodeID变量，存放访问过它的端节点的ID号。根据这个端节点的ID号来判断当前遍历是否继续或者退出。

下面先对几个概念进行定义：1度节点是指度为1的节点。链节点是指度为2且能通过的节点。分叉节点是指度大于2的节点或度为2但不能通过的节点。单链是指由链节点构成的没有分叉的链。端节点是指单链两端的节点。端节点是分叉节点或1度节点，虽然不属于单链，但是可以通过其遍历得到整个单链。由链节点构成的单链是化简的对象。

（1）、处理器从端节点出发访问单链节点，如果链节点是第一次被访问，则记录本次访问的端节点ID号；如果处理器访问了单链的所有链节点直到另一个端节点都没有遇到冲突，则将单链分配给该处理器（参见图4中的1、2）；

（2）、如果处理器在访问第一个链节点时发现该单链已被访问过（查看endNodeID值），则回退，继续访问下一条单链（参见图4中的3）；

（3）、如果处理器访问的链节点（非第一个单链节点）同时也被另外一个处理器访问，则比较自己端节点ID号与对方端节点ID号的大小，如果自己端节点ID号小，则继续访问，反之，则放弃，继续访问下一条单链（参见图4中的4、5）。

步骤S33、将所有的单链进行化简。

本步骤中，将单链节点全部删除，节点信息合并到端节点中，并将端节点***subGraphMap中。

步骤S34、判断是否有未访问的端节点，若有，返回执行步骤S32，否则，执行步骤S35。

步骤S35、图化简完成，生成contig。

当没有未访问的端节点时，化简完成。

步骤S4：合并contig（重叠群），生成scaffold（基因组支架）。

利用测序获得的pair-end信息，合并步骤S3得到的contig，生成更长的contig(即scaffold)。

步骤S5：输出scaffold。

本步骤输出scaffold到指定文件。

本发明实施例基于集群***，并行构造DeBruijn图，解决了大基因组拼接时由于其数据量太大，传统单机串行的基因拼接算法无法构图和无法进一步处理的问题；同时，在化简过程中，进行基于深度图遍历的并行化简，图化简过程简单，并行度高，拼接速度快。

实验：

本发明已经经过实验验证。

用C++、MPI实现了该基因拼接算法。实验平台包括10台曙光5000的高性能计算机，用InfiniBand网络互联。每台机器的配置为16核、32G共享内存。实验选取了Yeast，C.elegans两个基因组，使用Perl脚本自动生成了50倍Yeast测试数据，该数据包括17,007,362read；50倍C.elegans测试数据，该数据包括140,396,108read。read的长度从36bp到50bp不等，错误率为0。实验的目标是验证DeBruijn图的构建和化简算法的扩展性。

首先，测试了Yeast数据集，运行过程被分为三个阶段：并行文件I/O，构图，图的化简。对三个阶段分别计时。第一阶段，从一个分布式文件***读入文件。第二阶段，构建一个DeBruijn图。第三阶段，对DeBruijn图进行化简。第三阶段的时间开销有很大一部分是网络传输的时间开销。该算法中大部分节点只需移动一次，这对第三阶段的优化效果明显。如图5所示，其实验结果表明该算法有比较好的扩展性。当处理器的数目从8增加到128时，总时间从最开始的490s减少为63s，平均加速比为8倍。

然后，测试了C.elegans数据集。C.elegans数据集是Yeast数据集的10倍，运行时间见图6。当处理器的数目从8增加到128时，总时间从最开始的3844s减少为375s，平均加速比为10倍。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于DeBruijn图的并行基因拼接方法，其特征在于，包括：

S1、并行构建分布式DeBruijn图；

S2、剔除错误路径；

S3、基于深度图遍历方法对DeBruijn图进行化简；所述步骤S3具体包括：S31、初始化，将DeBruijn图分布存储在每个处理器的locationMap中，同时新建一个subGraphMap；S32、所有处理器并行从本地的locationMap中的端节点出发遍历单链；S33、将所有的单链进行化简；S34、判断是否有未访问的端节点，若有，返回执行步骤S32，否则，执行步骤S35；S35、图化简完成，生成contig；

S4、合并contig，生成scaffold；

S5、输出scaffold。

2.根据权利要求1所述的方法，其特征在于，所述步骤S1具体包括：

S11、所有的处理器并行读取原始的短序列文件，每个处理器读取短序列文件的一部分；

S12、每个处理器对读取的短序列进行切分，每个短序列被切分成长度为k的小片段k-mer，获取所述k-mer，把每个k-mer构造成一个节点；

S13、按照哈希规则，对所有的节点进行重分布。

3.根据权利要求1所述的方法，其特征在于，所述步骤S2具体为：

对步骤S1构建的所述分布式DeBruijn图进行并行去错。

4.根据权利要求3所述的方法，其特征在于，在所述步骤S2中，对步骤S1构建的所述分布式DeBruijn图中的Tips、Bubble和SpuriousLink三种结构进行剔除。

5.根据权利要求1所述的方法，其特征在于，所述步骤S32具体为：

处理器从端节点出发访问单链节点；如果所述链节点是第一次被访问，则记录本次访问的端节点ID号；如果处理器访问了单链的所有链节点直到另一个端节点都没有遇到冲突，则将所述单链分配给所述处理器；

如果处理器在访问第一个链节点时发现该单链已被访问过，则回退，继续访问下一条单链；

如果处理器访问的非第一个单链节点的链节点同时也被另外一个处理器访问，则比较自己端节点ID号与对方端节点ID号的大小，如果自己端节点ID号小，则继续访问，反之，则放弃，继续访问下一条单链。

6.根据权利要求1所述的方法，其特征在于，所述步骤S33具体为：

将单链节点全部删除，节点信息合并到端节点中，并将端节点***subGraphMap中。

7.根据权利要求1所述的方法，其特征在于，所述步骤S4具体为：

利用测序获得的pair-end信息，合并步骤S3得到的contig，生成scaffold，所述scaffold为更长的contig。