CN103699817B - 双向多步De Bruijn图的自环双向边识别与去除方法 - Google Patents

Abstract

本发明公开一种双向多步De Bruijn图的自环双向边识别与去除方法，包括步骤，S1、读取测序数据源文件，并构造双向多步De Bruijn图；S2、设定所述双向多步De Bruijn图中的每个顶点u的数据结构，对所述双向多步De Bruijn图自环双向边的识别；S3、对所述双向多步De Bruijn图自环双向边的去除。本发明技术效果包括（1）面向多步双向De Bruijn图；（2）根据自环双向边的目标节点以及自环双向边的目标节点方向来识别区分自环双向边；（3）根据自环双向边的目标节点方向，采取不同的自环双向边的去除方法；（4）可以有效的发现错误的双向边，而且避免对正确的双向边的破坏，从而提高contigs的长度，更重要的是能够同步提高contig的质量。

Description

双向多步De Bruijn图的自环双向边识别与去除方法

技术领域

本发明涉及基因测序领域，具体涉及双向多步De Bruijn图的自环双向边识别与去除方法。

背景技术

基因序列分析以算法与数学模型为核心，包括：基因数据的存储与获取、序列比对、测序与拼接、基因预测、生物进化与系统发育分析、蛋白质结构预测、RNA结构预测、分子设计与药物设计、代谢网络分析、基因芯片、DNA计算等。生物技术和计算机信息处理技术的紧密结合，加快了处理生物信息数据的速度，使得在尽短的时间内对生物学做出准确的诠释，加快生物信息学的发展。

基因序列分析是对海量基因序列数据进行分析，从而提取和挖据新的生物信息知识。涉及到计算机技术中的机器学习、模式识别、书籍分析与挖掘、组合数学、随机模型、字符串、图形算法、分布式计算、高性能计算、并行计算等。

基因是人类最基本的遗传密码，代表着每个人的生命信息。基因序列上存在着遗传位点的细微差异，这些遗传密码的多态性与人类的健康、致病机理、医学治疗有着相当密切的关系。

自1977年Sanger测序技术问世以来，经过三十多年的发展，DNA测序技术发展突飞猛进，以高通量、短序列为特点的第二代测序技术逐渐占领市场，以单分子测序为特点的第三代测序技术也逐渐出现，分别在测序特点上占有不同的优势。传统的基因测序方法的数据提取和分析软件经过近10年来的研究与开发，目前已经较为完善。但是，测序技术的发展，带来了测序数据的变化，使得当前存在的数据处理软件不能满足当前生物医学研究的需求。

新一代高通量测序方法在技术的应用，可以在短时间内完成整个基因组数据的测定。高通量测序方法的日新月异也同时对获取的基因数据的分析处理方法提出了挑战。目前，迫切需要开发能满足高通量测序技术的海量数据处理的生物信息学平台。面对个人基因组计划及未来的个性化医疗前景，高效低成本的测序技术成为必然的趋势。同时，简化高效的一站式完备的生物信息学数据分析平台等完备的测序解决方案，也是极为重要不可或缺的发展方向。

然而新一代的高通量测序方法虽然测序通量高，但是却会引入测序误差，同时测序样本本身由于基因突变，测序不均匀而导致有SNP的出现，而上述测序误差，测序不均将会在基因组组装时构造的多步双向De Bruijn图中引入一些错误的双向边，其中有一部分是自环双向边。而上述错误的自环双向边在De Bruijn图中，能够阻碍图的收缩，contig无法扩展，最终使得contig的长度和质量都很低。

新一代的高通量测序方法产生的短基因片段的组装导致大量的测序错误，加大了组装算法的计算量。大量的测序错误，使得组装错误率增加，严重影响了组装结果。

目前组装算法策略分为两类，一是基于Overlap-Layout-Consensus(OLC)的算法，另一个是基于DeBruijn图的算法。其中，基于OLC组装算法开发的软件，如SSAKE、VCAKE、SHARCGS等，在基因长序列组装中更占有优势，但并不完全适用于新一代的短序列组装。与OLC组装算法不同，DeBruijn算法不再以read为单位组织数据，而是以k-mers为单位进行数据组装，其优点主要有以下几个方面：首先，以k-mers为单位进行序列组装，不影响节点的质量，减少了冗余数据量；其次，在图中重复区域只出现一次，便于识别，可以避免错误的组装，减小出错率；最后，采取将有重叠区域映射到同一条弧上的策略，从而简化了搜索路径。目前，很多短序列组装算法都使用这种框架，如Velvet、IDBA、SOAPdenovo，ABySS。

Velvet有效的利用De Bruijn图，实现了高效的短序列组装。Velvet以k-mer为基本单位构建De Bruijn图，利用图的结构，结合相应的序列特征，简化图的构造，最终找到一条最优路径完成组装过程。Velvet把焦点集中在错误的数据产生的三种结构上，即tip、bubble、以及erroneous connection。依照长度原则和少数性原则，将长度小于2k的均去除；利用Tour Bus算法中的深度优先搜索策略合并bubble，最后利用覆盖度阈值法去除了erroneous connection。该方法也充分利用了paired-end双端信息，进一步解决repeat问题，优化了组装效果。Velvet充分利用图的结构性质，简化了数据冗余，速度较之前的算法有了很大的改进。虽然它没有在预处理阶段对序列进行纠错，但是其对错误的预防机制，很大程度上的弥补了这方面的缺陷。这使得它更好的应用在大型基因组序列的组装中。

IDBA也是基于De Bruijn图，实现简便且高效的短序列组装。IDBA以k-mer为基本单位，采用一个变化的k值域(Kmin-Kmax)，代替使用固定的k值来得到k-mers的长度。由于基因组装以k-mers为单位，通常会形成很多个重叠单元，使得组装面临着错误位置组装、顶点缺失和覆盖度低的问题。正确的选择k值的大小成为组装的一个关键因素。一些错误的reads的产生，也导致产生了大量的branching。K值越小，branching问题越严重，k值越大，则出现的reapt区域则变少，直接影响了组装的质量。IDBA采用不固定的k值进行组装，很好的解决branching问题，提高了组装的质量。另外IDBA通过删除低覆盖率的错误k-mers而使得IDBA的内存使用率明显降低，同时也提升了IDBA的处理速度。

SOAPdenovo能够高效高质量的完成数以亿计的reads的组装。SOAPdenovo继承了OLC算法和De Bruijn图算法的优点，使得其组装质量大为提高。SOAP通过预置k-mer阈值的方法，采取过滤、纠错的方式减少了错误序列的产生。同时，借鉴了Velvet软件的方法成功处理了bubble，使得其平均覆盖度增加。另外，SOAPdenovo利用了双端信息进行进行重叠区域匹配，并合并read生成contig片段，生成基于contig的图结构，从而，SOAPdenovo大大简化了contig图的复杂性。

ABySS引进并行计算的思想，搭建一个linux集群，在集群上建立一个分布式的DeBruijn图结构，将数据分布式存储于每个节点上。其采用MPI通信机制完成节点之间的相互通信。从构建图、纠错处理到后面的定点融合，最后完成整个基因组序列的再现，其在运行时间和内存消耗方面占有很大的优势，并且其错误率极低，在性能方面特别是cluster中单机内存使用上均有很大的提升，正在得到越来越广泛的应用。

然而，上述策略纠错去错的策略都是基于丰度均一的前提假设，但实际测序中部分顶点由于测序丰度比较低而被误认为是错误顶点，而处于重复序列中的SNP或者突变却可能由于去丰度过高，被认为是正确的顶点。另外上述策略都是针对De Bruijn图的顶点进行分类和去错的，针对双向多步De Bruijn图中的自环双向边的识别和去除方法却均未有涉及。

同时，自环双向边的形成原因比较多，其中包括基因组本身的SNP或者突变，同时也可能是由于测序错误，所以并不是所有的自环双向边都有方法去除。

发明内容

本发明目的在于解决现有技术存在的问题，提供一种双向多步De Bruijn图的自环双向边识别与去除方法。

本发明的技术方案包括一种双向多步De Bruijn图的自环双向边识别与去除方法，包括步骤，

S1、读取测序数据源文件，并构造双向多步De Bruijn图；

S2、设定所述双向多步De Bruijn图中的每个顶点u的数据结构，对所述双向多步De Bruijn图自环双向边的识别；

S3、对所述双向多步De Bruijn图自环双向边的去除。

优选地，所述De Bruijn图构造步骤为，

S11、读取一个序列s；

S12、将序列s用滑动窗口切割为多个片段t，选取一片段t其标志数为cur、并标记其前、后的片段的标志数分别为pre、lat；

S13、若t的编码小于其互补片段编码，则交换pre，lat的值；

S14、在cur的正向位置映射表的相应bit位置1来表示指向pre的边；

S15、在cur的反向位置映射表的相应bit位置1来表示指向lat的边；

S16、重复步骤S12-S15，处理序列s的其他片段t，直至完成序列s的全部片段t，执行步骤S17；

S17、读取一个新的序列s，重复步骤S12-S16；直至处理完所有的序列，执行步骤S18；

S18、完成双向多步de Bruijn图的构造。

优选地，所述自环双向边的识别包括，

S21、遍历所述双向多步De Bruijn图中的每个顶点u；

S22、统计所述顶点u中指向所述顶点u的边；

S23、判断所述顶点u是否存在指向所述顶点u的边，是则执行步骤S21；否则判断所述顶点u是否只有一条边或者两条边指向所述顶点u，是则执行步骤S24；否则执行步骤S21；

S24、判断两条自环边指向所述顶点u的方向是否不同，是则，所述自环边无法处理，执行步骤S21；否则执行步骤S25；

S25、如果所述顶点u的出度减去所述顶点u的自环边的个数deg等于1或者2，且两个非自环边的下标一个小于4、另一个大于3,则标记所述顶点u为一可扩展的自环顶点，执行步骤S21。

优选地，所述自环双向边的去除包括，

S31、遍历整个双向多步De Bruijn图，访问所述De Bruijn图中的每个顶点u；

S32、如果所述顶点u被标记为可扩展的自环顶点；则执行步骤S33，否则执行步骤S31；

S33、选择顶点u的一自环边和一非自环边，同时所述自环边指向目的节点的方向和所述非自环边出发的方向相同；

S34、将所述顶点u的所述自环边和非自环边合并，同时将所述顶点u的非自环边删除。

本发明的有益效果包括：

（1）面向多步双向De Bruijn图；

（2）根据自环双向边的目标节点以及自环双向边的目标节点方向来识别区分自环双向边；

（3）根据自环双向边的目标节点方向，采取不同的自环双向边的去除方法；

（4）可以有效的发现错误的双向边，而且避免对正确的双向边的破坏，从而提高contigs的长度，更重要的是能够同步提高contig的质量。

附图说明

图1为本发明一实施例的自环双向边的识别流程图。

图2为本发明一实施例的自环双向边的去除流程图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步详细说明。

本发明提供一种双向多步De Bruijn图的自环双向边识别与去除方法，包括步骤，

S1、读取测序数据源文件，并构造双向多步De Bruijn图；

S2、设定所述双向多步De Bruijn图中的每个顶点u的数据结构，对所述双向多步De Bruijn图自环双向边的识别；

如果所述顶点u的出度减去所述顶点u的自环边的个数deg等于1或者2，且两个非自环边的下标一个小于4、另一个大于3,则标记所述顶点u为一可扩展的自环顶点。

S3、对所述双向多步De Bruijn图自环双向边的去除；

选择顶点u的一自环边和一非自环边，同时所述自环边指向目的节点的方向和所述非自环边出发的方向相同；将所述顶点u的所述自环边和非自环边合并，同时将所述顶点u的非自环边删除。

本发明实施例，（1）面向多步双向De Bruijn图；（2）根据自环双向边的目标节点以及自环双向边的目标节点方向来识别区分自环双向边；（3）根据自环双向边的目标节点方向，采取不同的自环双向边的去除方法；（4）可以有效的发现错误的双向边，而且避免对正确的双向边的破坏，从而提高contigs的长度，更重要的是能够同步提高contig的质量。

优选地，所述De Bruijn图构造步骤为，

S11、读取一个序列s；

S12、将序列s用滑动窗口切割为多个片段t，选取一片段t其标志数为cur、并标记其前、后的片段的标志数分别为pre、lat；

S13、若t的编码小于其互补片段编码，则交换pre，lat的值；

S14、在cur的正向位置映射表的相应bit位置1来表示指向pre的边；

S15、在cur的反向位置映射表的相应bit位置1来表示指向lat的边；

S16、重复步骤S12-S15，处理序列s的其他片段t，直至完成序列s的全部片段t，执行步骤S17；

S17、读取一个新的序列s，重复步骤S12-S16；直至处理完所有的序列，执行步骤S18；

S18、完成双向多步de Bruijn图的构造。

如图1所示，统计整个双向多步De Bruijn图的所有双向边的加权平均覆盖度W。所述自环双向边的识别包括，

S21、遍历所述双向多步De Bruijn图中的每个顶点u；

S22、统计所述顶点u中指向所述顶点u的边；

S23、判断所述顶点u是否存在指向所述顶点u的边，是则执行步骤S21；否则判断所述顶点u是否只有一条边或者两条边指向所述顶点u，是则执行步骤S24；否则执行步骤S21；

S24、判断两条自环边指向所述顶点u的方向是否不同，是则，所述自环边无法处理，执行步骤S21；否则执行步骤S25；

S25、如果所述顶点u的出度减去所述顶点u的自环边的个数deg等于1或者2，且两个非自环边的下标一个小于4、另一个大于3,则标记所述顶点u为一可扩展的自环顶点，执行步骤S21。

如图1所示，所述自环双向边的去除包括，

S31、遍历整个双向多步De Bruijn图，访问所述De Bruijn图中的每个顶点u；

S32、如果所述顶点u被标记为可扩展的自环顶点；则执行步骤S33，否则执行步骤S31；

S33、选择顶点u的一自环边和一非自环边，同时所述自环边指向目的节点的方向和所述非自环边出发的方向相同；

S34、将所述顶点u的所述自环边和非自环边合并，同时将所述顶点u的非自环边删除。

上述去处方法，使De Bruijn图继续收缩，contigs继续扩展；同时不会引入错误，造成contig质量的下降，准确度降低。

以上所述本发明的具体实施方式，并不构成对本发明保护范围的限定。任何根据本发明的技术构思所作出的各种其他相应的改变与变形，均应包含在本发明权利要求的保护范围内。

Claims (3)

1.一种双向多步De Bruijn图的自环双向边识别与去除方法，其特征在于，包括步骤，

S1、读取测序数据源文件，并构造双向多步De Bruijn图；

S2、设定所述双向多步De Bruijn图中的每个顶点u的数据结构，对所述双向多步DeBruijn图自环双向边的识别；

S3、对所述双向多步De Bruijn图自环双向边的去除；

所述自环双向边的识别包括，

S21、遍历所述双向多步De Bruijn图中的每个顶点u；

S22、统计所述顶点u中指向所述顶点u的边；

S23、判断所述顶点u是否存在指向所述顶点u的边，是则执行步骤S21；否则判断所述顶点u是否只有一条边或者两条边指向所述顶点u，是则执行步骤S24；否则执行步骤S21；

S24、判断两条自环边指向所述顶点u的方向是否不同，是则，所述自环边无法处理，执行步骤S21；否则执行步骤S25；

S25、如果所述顶点u的出度减去所述顶点u的自环边的个数deg等于1或者2，且两个非自环边的下标一个小于4、另一个大于3,则标记所述顶点u为一可扩展的自环顶点，执行步骤S21。

2.如权利要求1所述的自环双向边识别与去除方法，其特征在于，所述De Bruijn图构造步骤为，

S11、读取一个序列s；

S12、将序列s用滑动窗口切割为多个片段t，选取一片段t其标志数为cur、并标记其前、后的片段的标志数分别为pre、lat；

S13、若t的编码小于其互补片段编码，则交换pre，lat的值；

S14、在cur的正向位置映射表的相应bit位置1来表示指向pre的边；

S15、在cur的反向位置映射表的相应bit位置1来表示指向lat的边；

S16、重复步骤S12-S15，处理序列s的其他片段t，直至完成序列s的全部片段t，执行步骤S17；

S17、读取一个新的序列s，重复步骤S12-S16；直至处理完所有的序列，执行步骤S18；

S18、完成双向多步de Bruijn图的构造。

3.如权利要求1所述的自环双向边识别与去除方法，其特征在于，所述自环双向边的去除包括，

S31、遍历整个双向多步De Bruijn图，访问所述De Bruijn图中的每个顶点u；

S32、如果所述顶点u被标记为可扩展的自环顶点；则执行步骤S33，否则执行步骤S31；

S33、选择顶点u的一自环边和一非自环边，同时所述自环边指向目的节点的方向和所述非自环边出发的方向相同；

S34、将所述顶点u的所述自环边和非自环边合并，同时将所述顶点u的非自环边删除。