双向多步De Bruijn图的重复双向边识别与去除方法
【专利摘要】本发明公开一种双向多步De?Bruijn图的重复双向边识别与去除方法，包括步骤：S1.读取测序数据源文件，并构造双向多步De?Bruijn图；S2.判断出度为3的顶点u是否指向的同一个顶点v的同一方向或者同一端，指向相同的一端为耦合重复双向边，指向不同的一端为错误重复双向边；S3.对被判定含有耦合重复双向边顶点u和被判定含有错误重复双向边顶点u进行判断删除。本发明根据重复双向边形成，将部分重复双向边分类，并采取正确的措施清除De?Bruijn图中存在的重复双向边；使contigs可以继续延长，De?Bruijn图继续收缩；有效的解耦并去除重复双向边，而且避免对正确的双向边的破坏，从而可以一定程度上提高contigs的长度，提高contig的质量。
【专利说明】 双向多步De Bruijn图的重复双向边识别与去除方法
【技术领域】
[0001]本发明涉及基因测序领域,具体涉及双向多步De Bruijn图的重复双向边识别与去除方法。
【背景技术】
[0002]基因序列分析以算法与数学模型为核心，包括:基因数据的存储与获取、序列比对、测序与拼接、基因预测、生物进化与系统发育分析、蛋白质结构预测、RNA结构预测、分子设计与药物设计、代谢网络分析、基因芯片、DNA计算等。生物技术和计算机信息处理技术的紧密结合，加快了处理生物信息数据的速度，使得在尽短的时间内对生物学做出准确的诠释，加快生物信息学的发展。
[0003]基因序列分析是对海量基因序列数据进行分析，从而提取和挖据新的生物信息知识。涉及到计算机技术中的机器学习、模式识别、书籍分析与挖掘、组合数学、随机模型、字符串、图形算法、分布式计算、高性能计算、并行计算等。
[0004]基因是人类最基本的遗传密码，代表着每个人的生命信息。基因序列上存在着遗传位点的细微差异，这些遗传密码的多态性与人类的健康、致病机理、医学治疗有着相当密切的关系。
[0005]自1977年Sanger测序技术问世以来,经过三十多年的发展,DNA测序技术发展突飞猛进，以高通量、短序列为特点的第二代测序技术逐渐占领市场，以单分子测序为特点的第三代测序技术也逐渐出现，分别在测序特点上占有不同的优势。传统的基因测序方法的数据提取和分析软件经过近10年来的研究与开发，目前已经较为完善。但是，测序技术的发展，带来了测序数据的变化，使得当前存在的数据处理软件不能满足当前生物医学研究的需求。
[0006]新一代高通量测序方法在技术的应用，可以在短时间内完成整个基因组数据的测定。高通量测序方法的日新月异也同时对获取的基因数据的分析处理方法提出了挑战。目前，迫切需要开发能满足高通量测序技术的海量数据处理的生物信息学平台。面对个人基因组计划及未来的个性化医疗前景，高效低成本的测序技术成为必然的趋势。同时，简化高效的一站式完备的生物信息学数据分析平台等完备的测序解决方案，也是极为重要不可或缺的发展方向。
[0007]然而新一代的高通量测序方法虽然测序通量高，但是却会引入测序误差，同时测序样本本身由于基因突变，SNP，测序不均会在基因组组装时构造的双向多步De Bruijn图中引入重复双向边。而这些重复双向边在对整个De Bruijn图中将直接影响图的收缩过程的效果，并阻碍contig的延长。
[0008]新一代的高通量测序方法产生的短基因片段的组装导致大量的测序错误，加大了组装算法的计算量。大量的测序错误，使得组装错误率增加，严重影响了组装结果。
[0009]目前组装算法策略分为两类，一是基于Overlap-Layout-Consensus (OLC)的算法，另一个是基于DeBruijn图的算法。其中，基于OLC组装算法开发的软件，如SSAKE、VCAKE、SHARCGS等，在基因长序列组装中更占有优势，但并不完全适用于新一代的短序列组装。与OLC组装算法不同,De Brui jn算法不再以read为单位组织数据,而是以k-mers为单位进行数据组装，其优点主要有以下几个方面:首先，以k-mers为单位进行序列组装，不影响节点的质量，减少了冗余数据量；其次，在图中重复区域只出现一次，便于识别，可以避免错误的组装，减小出错率；最后，采取将有重叠区域映射到同一条弧上的策略，从而简化了搜索路径。目前，很多短序列组装算法都使用这种框架，如Velvet、IDBA、SOAPdenovo,ABySS0
[0010]Velvet有效的利用De Brui jn图，实现了高效的短序列组装。Velvet以kner为基本单位构建De Bruijn图，利用图的结构，结合相应的序列特征，简化图的构造，最终找到一条最优路径完成组装过程。Velvet把焦点集中在错误的数据产生的三种结构上，即tip、bubble、以及erroneous connection。依照长度原则和少数性原则,将长度小于2k的均去除；利用Tour Bus算法中的深度优先搜索策略合并bubble，最后利用覆盖度阈值法去除了erroneous connection。该方法也充分利用了 paired-end双端信息,进一步解决repeat问题，优化了组装效果。Velvet充分利用图的结构性质，简化了数据冗余，速度较之前的算法有了很大的改进。虽然它没有在预处理阶段对序列进行纠错，但是其对错误的预防机制，很大程度上的弥补了这方面的缺陷。这使得它更好的应用在大型基因组序列的组装中。
[0011]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的处理速度。
[0012]SOAPdenovo能够高效高质量的完成数以亿计的reads的组装。SOAPdenovo继承了OLC算法和De Bruijn图算法的优点，使得其组装质量大为提高。SOAP通过预置k-mer阈值的方法，采取过滤、纠错的方式减少了错误序列的产生。同时，借鉴了 Velvet软件的方法成功处理了 bubble，使得其平均覆盖度增加。另外，SOAPdenovo利用了双端信息进行进行重叠区域匹配，并合并read生成contig片段,生成基于contig的图结构,从而，SOAPdenovo大大简化了 contig图的复杂性。
[0013]ABySS引进并行计算的思想,搭建一个Iinux集群,在集群上建立一个分布式的DeBruijn图结构，将数据分布式存储于每个节点上。其采用MPI通信机制完成节点之间的相互通信。从构建图、纠错处理到后面的定点融合，最后完成整个基因组序列的再现，其在运行时间和内存消耗方面占有很大的优势，并且其错误率极低，在性能方面特别是cluster中单机内存使用上均有很大的提升，正在得到越来越广泛的应用。
[0014]然而以上的纠错去错的策略都是基于丰度均一的前提假设，但实际测序中部分顶点由于测序丰度比较低而被误认为是错误顶点，而那些处于重复序列中的SNP或者突变却可能由于去丰度过高，而被认为是正确的顶点。另外，上述策略都是针对De Bruijn图的顶点进行分类和去错的，而针对双向多不De Bruijn图，特别是基于针对双向多步De Bruijn 图的重复双向边识别以及去除的工作却均未有涉及。
【发明内容】
[0015]本发明目的在于解决现有技术存在的问题，提供一种双向多步De Bruijn图的重复双向边识别与去除方法
[0016]本发明的技术方案包括一种双向多步De Bruijn图的重复双向边识别与去除方法，包括步骤，
[0017]S1、读取测序数据源文件，并构造双向多步De Bruijn图；
[0018]S2、判断出度为3的顶点u是否指向的同一个顶点V的同一方向或者同一端，指向相同的一端为耦合重复双向边，指向不同的一端为错误重复双向边；
[0019]S3、对被判定含有耦合重复双向边顶点u和被判定含有错误重复双向边顶点u进行判断删除。
[0020]优选地,所述De Bruijn图构造步骤为，
[0021]S11、读取一个序列s;
[0022]S12、将序列s用滑动窗口切割为多个片段t，选取一片段t其标志数为cur、并标记其前、后的片段的标志数分别为pre、Iat ；
[0023]S13、若t的编码小于其互补片段编码，则交换pre，Iat的值；
[0024]S14、在cur的正向位置映射表的相应bit位置I来表示指向pre的边；
[0025]S15、在cur的反向位置映射表的相应bit位置I来表示指向Iat的边；
[0026]S16、重复步骤S12-S15，处理序列s的其他片段t，直至完成序列s的全部片段t，执行步骤S17 ;
[0027]S17、读取一个新的序列S，重复步骤S12-S16 ;直至处理完所有的序列，执行步骤S18 ；
[0028]S18、完成双向多步de Bruijn图的构造。
[0029]优选地，所述重复双向边识别包括，
[0030]S21、遍历所述双向多步De Bruijn图中的每个顶点u ；
[0031]S22、如果顶点u中有两条边指向同一个目的顶点V，且所述两条边的下标设为i和j，指向顶点V的方向为Direci和Direcj ；
[0032]S23、统计顶点u的边的出度，如果边的出度是3，则执行步骤S24，否则执行步骤S21 ；
[0033]S24、如果方向Direci不等于方向Direcj,则该边为f禹合重复双向边；
[0034]S25、如果方向Direci等于方向Direc j,则该边为错误重复双向边。
[0035]优选地,所述重复双向边去除包括，
[0036]S31、遍历所述双向多步De Bruijn图中的每个顶点u ；
[0037]S32、如果顶点u中有两条边指向同一个目的顶点V，且所述两条边的下标设为i和j，指向V的方向为Direci和Direcj,另外一条指向其他节点w的边设为k ；
[0038]S33、统计顶点u的边的出度，如果边的出度是3，则执行步骤S34，否则执行步骤S31 ；
[0039]S34、双向边U.arcs [i]和u.arcs [j]的对称边在顶点V的下标记为retlndexi,retlndexj，在 kmer 上的方向为 retDireci，retDirecj ；
[0040]S35、双向边u.arcs [k]的对称边在顶点w的下标记为Retlndex，方向记为RetDirec ；
[0041]S36、如果Direci不等于Direcj，则该边为賴合重复双向边；判断删除；
[0042]S37、如果Direci等于Direcj则该边为错误重复双向边；判断删除。
[0043]优选地，所述步骤S36包括，
[0044]S361、如果I^etIndexa,则执行步骤S362，否则进入S363 ；
[0045]S362、w.arcs [retlndex] +=u.arcs[j]+v.arcs [retlndexi]；
[0046]u.arcs [i]+=v.arcs [retlndexj]+u.arcs [k]；
[0047]删除双向边u.arcs [i]以及双向边的对称双向边v.arcs [retlndexi]:
[0048]u.arcs [i].clear () ；u.multiplicity [i] =0 ；
[0049]v.arcs [retlndexi].clear () ；v.multiplicity [i] =0 ；
[0050]S363、w.arcs [retlndex] =w.arcs [Retlndex] +u.arcs[i]+v.arcs [retlndexj]；
[0051]u.arcs [i] =u.arcs[j]+v.arcs [retlndexi] +u.arcs [k]；
[0052]删除双向边u.arcs [j]以及双向边的对称双向边v.arcs [retlndexj]:
[0053]u.arcs [ j].clear () ；u.multiplicity [j] =0 ；
[0054]v.3JTcs [retlndexjj.cIqsJl 0 ! v.multiplicity [j]—0。
[0055]优选地，所述步骤S37包括，
[0056]S371、如果i〈4，j〈4并且k>3，则执行步骤S372，否则执行步骤S373 ；
[0057]S372、删除双向边u.arcs[j]以及双向边的对称双向边ν.arcs [retlndexj]:
[0058]u.arcs [ j].clear () ；u.multiplicity [j] =0 ；
[0059]v.arcs[retlndexj].clear () ； v.multiplicity[j]=0 ；
[0060]S373、删除双向边u.arcs [i]以及双向边的对称双向边v.arcs [retlndexi]:
[0061]u.arcs [i].clear () ；u.multiplicity [i] =0 ；
[0062]v.arcs [retlndexi].clear () ； v.multiplicity [i] =0。
[0063]本发明的有益效果包括:根据重复双向边形成，将部分重复双向边分类，并采取正确的措施清除De Brui jn图中存在的重复双向边；使contigs可以继续延长,De Bruijn图继续收缩；有效的解耦并去除重复双向边，而且避免对正确的双向边的破坏，从而可以一定程度上提高contigs的长度,提高contig的质量。
【专利附图】
【附图说明】
[0064]图1为本发明一实施例的判定重复双向边的流程图。
[0065]图2为本发明一实施例的错误双向边的识别与去除流程图。
【具体实施方式】
[0066]下面结合附图和具体实施例对本发明作进一步详细说明。
[0067]本发明实施例识别重复双向边，并且将重复双向边去除或者解耦。根据其造成的原因的不同，重复双向边的解决方法也不同。如果重复双向边本身是正确的，只是由于恰巧有一条双向边将两个相邻kmer的链接在一起，这种情况下，只能解耦并合并该双向边到可延长的双向边上；而对于是由测序错误造成的两个kmer同时被两条边链接，其中一条来源于测序错误，那么根据其覆盖度大小而选择正确的一条边，同时去除错误的一条边。以上两种情况都尽力辅助原始的De Bruijn图继续收缩,边继续延长,从而得到更长的contig,并通过尽量去除错误，并不引入新的错误来保障其准确性。
[0068]实施例提供一种双向多步De Bruijn图的重复双向边识别与去除方法，包括步骤，
[0069]S1、读取测序数据源文件，并构造双向多步De Bruijn图；
[0070]S2、判断出度为3的顶点u是否指向的同一个顶点V的同一方向或者同一端，指向相同的一端为耦合重复双向边，指向不同的一端为错误重复双向边；
[0071]S3、对被判定含有耦合重复双向边顶点u和被判定含有错误重复双向边顶点u进行判断删除。
[0072]对一被判定含有I禹合重复双向边顶点u,若指向其他顶点w的边的下标retlndex小于 4,合并 3 条边 w.arcs [retlndex]，u.arcs [ j]和 v.arcs [retlndexi]，否则合并边 w.arcs [retlndex], u.arcs [i]和 v.arcs [retlndexj];同时如果 retlndex〈4 则删除u.arcs [ j], v.arcs [retlndexi];否者删除 u.arcs [i], v.arcs [retlndexj]。
[0073]对一被判定含有错误重复双向边顶点u，假设指向其他顶点其指向同一个顶点的两条边所在下标为i，j ;以及指向其他顶点的边的下标为k，若i〈4，j〈4并且k>3，则删除双向边U.arcs[j]以及双向边的对称双向边V.arcs [retlndexj],否则删除双向边u.arcs[i]以及双向边的对称双向边V.arcs [retlndexi]。
[0074]本发明实施例根据重复双向边形成，将部分重复双向边分类，并采取正确的措施清除De Bruijn图中存在的重复双向边；使contigs可以继续延长，De Bruijn图继续收缩；有效的解耦并去除重复双向边，而且避免对正确的双向边的破坏，从而可以一定程度上提高contigs的长度,提高contig的质量。
[0075]优选地，De Bruijn图构造步骤为，
[0076]S11、读取一个序列s;
[0077]S12、将序列s用滑动窗口切割为多个片段t，选取一片段t其标志数为cur、并标记其前、后的片段的标志数分别为pre、Iat ；
[0078]S13、若t的编码小于其互补片段编码，则交换pre，Iat的值；
[0079]S14、在cur的正向位置映射表的相应bit位置I来表示指向pre的边；
[0080]S15、在cur的反向位置映射表的相应bit位置I来表示指向Iat的边；
[0081]S16、重复步骤S12-S15，处理序列s的其他片段t，直至完成序列s的全部片段t，执行步骤S17 ;
[0082]S17、读取一个新的序列S，重复步骤S12-S16 ;直至处理完所有的序列，执行步骤S18 ；
[0083]S18、完成双向多步de Bruijn图的构造。
[0084]如图1所示，重复双向边识别包括，
[0085]S21、遍历双向多步De Bruijn图中的每个顶点u ；
[0086]S22、如果顶点u中有两条边指向同一个目的顶点V，且两条边的下标设为i和j，指向顶点V的方向为Direci和Direcj ；
[0087]S23、统计顶点u的边的出度(个数)，如果边的出度是3，则执行步骤S24，否则执行步骤S21 ；
[0088]S24、如果方向Direci不等于方向Direcj,则该边为f禹合重复双向边；
[0089]S25、如果方向Direci等于方向Direc j，则该边为错误重复双向边。
[0090]如图2所示，重复双向边去除包括，
[0091]S31、遍历双向多步De Bruijn图中的每个顶点u ；
[0092]S32、如果顶点u中有两条边指向同一个目的顶点V，且两条边的下标设为i和j，指向V的方向为Direci和Direcj,另外一条指向其他节点w的边设为k ；
[0093]S33、统计顶点u的边的出度(个数)，如果边的出度是3，则执行步骤S34，否则执行步骤S31 ；
[0094]S34、双向边U.arcs [i]和u.arcs [j]的对称边在顶点V的下标记为retlndexi,retlndexj,在 kmer 上的方向为 retDireci, retDirecj ；
[0095]S35、双向边u.arcs [k]的对称边在顶点w的下标记为Retlndex,方向记为RetDirec ；
[0096]S36、如果Direci不等于Direcj,则该边为I禹合重复双向边；判断删除；[0097]S37、如果Direci等于Direcj则该边为错误重复双向边；判断删除。
[0098]其中，步骤S36包括，
[0099]S361、如果retlndex〈4，则执行步骤S362，否则进入S363 ；
[0100]S362、W.arcs[retlndex]+=U.arcs[j]+V.arcs[retlndexi]；
[0101]u.arcs[i]+=v.arcs[retlndexj]+u.arcs[k]；
[0102]删除双向边u.arcs [i]以及双向边的对称双向边v.arcs [retlndexi]:
[0103]u.arcs [i].clear () ；u.multiplicity [i] =0 ；
[0104]v.arcs [retlndexi].clear () ； v.multiplicity [i] =0 ；
[0105]S363、w.arcs[retlndex]=w.arcs[Retlndex]+u.arcs[i]+v.arcs[retlndexj]；
[0106]u.arcs [i] =u.arcs [j]+v.arcs [retlndexi]+u.arcs [k]；
[0107]删除双向边u.arcs [j]以及双向边的对称双向边v.arcs [retlndexj]:
[0108]u.arcs [ j].clear () ；u.multiplicity [j] =0 ；
[0109]v.3Jl cs [retlndex jj.c I qsJl 0 ! v.mu ltiplicity [j]—0。
[0110]其中，步骤S37包括，
[0111]S371、如果i〈4，j〈4并且k>3，则执行步骤S372，否则执行步骤S373 ；
[0112]S372、删除双向边u.arcs[j]以及双向边的对称双向边v.arcs [retlndexj]:
[0113]u.arcs [j].clear O ；u.multiplicity [j] =0 ；
[0114]v.arcs [retlndexj].clear () ； v.multiplicity [ j] =0 ；
[0115]S373、删除双向边u.arcs[i]以及双向边的对称双向边v.arcs [retlndexi]:
[0116]u.arcs [i].clear () ；u.multiplicity [i] =0 ；
[0117]v.arcs [retlndexi].clear () ； v.multiplicity [i] =0。
[0118]以上所述本发明的【具体实施方式】，并不构成对本发明保护范围的限定。任何根据本发明的技术构思所作出的各种其他相应的改变与变形，均应包含在本发明权利要求的保护范围内。
【权利要求】
1.一种双向多步De Bruijn图的重复双向边识别与去除方法,其特征在于,包括步骤， 51、读取测序数据源文件，并构造双向多步DeBruijn图； 52、判断出度为3的顶点u是否指向的同一个顶点V的同一方向或者同一端，指向相同的一端为耦合重复双向边，指向不同的一端为错误重复双向边； 53、对被判定含有耦合重复双向边顶点u和被判定含有错误重复双向边顶点u进行判断删除。
2.如权利要求1所述的重复双向边识别与去除方法，其特征在于，所述DeBruijn图构造步骤为， SI 1、读取一个序列S ； 512、将序列s用滑动窗口切割为多个片段t，选取一片段t其标志数为cur、并标记其前、后的片段的标志数分别为pre、Iat ； 513、若t的编码小于其互补片段编码，则交换pre，Iat的值； 514、在cur的正向位置映射表的相应bit位置I来表示指向pre的边； 515、在cur的反向位置映射表的相应bit位置I来表示指向Iat的边； 516、重复步骤S12-S15，处理序列s的其他片段t，直至完成序列s的全部片段t，执行步骤S17 ； 517、读取一个新的序列S，重复步骤S12-S16;直至处理完所有的序列，执行步骤S18 ； 518、完成双向多步deBruijn图的构造。
3.如权利要求1所述的重复双向边识别与去除方法，其特征在于，所述重复双向边识别包括， 521、遍历所述双向多步DeBruijn图中的每个顶点u ； 522、如果顶点u中有两条边指向同一个目的顶点V，且所述两条边的下标设为i和j，指向顶点V的方向为Direci和Direcj ；523、统计顶点u的边的出度，如果边的出度是3，则执行步骤S24，否则执行步骤S21； 524、如果方向Direci不等于方向Direcj,则该边为f禹合重复双向边； 525、如果方向Direci等于方向Direcj,则该边为错误重复双向边。
4.如权利要求3所述的重复双向边识别与去除方法，其特征在于，所述重复双向边去除包括， 531、遍历所述双向多步DeBruijn图中的每个顶点u ； 532、如果顶点u中有两条边指向同一个目的顶点V，且所述两条边的下标设为i和j，指向V的方向为Direci和Direcj,另外一条指向其他节点w的边设为k ；533、统计顶点u的边的出度，如果边的出度是3，则执行步骤S34，否则执行步骤S31； 534、双向边u.arcs [i]和u.arcs [j]的对称边在顶点v的下标记为retlndexi,retlndexj,在 kmer 上的方向为 retDireci, retDirecj ；535、双向边u.arcs [k]的对称边在顶点w的下标记为RetIndex,方向记为RetDirec ； 536、如果Direci不等于Direcj,则该边为稱合重复双向边；判断删除； 537、如果Direci等于Direcj则该边为错误重复双向边；判断删除。
5.如权利要求4所述的重复双向边识别与去除方法，其特征在于，所述步骤S36包括， S361、如果retlndex〈4，则执行步骤S362，否则进入S363 ；.5362、w.arcs [retlndex] +=u.arcs [ j] +v.arcs [retlndexi]；u.arcs[i]+=v.arcs[retlndexj]+u.arcs[k]；删除双向边u.arcs [i]以及双向边的对称双向边v.arcs [retlndexi]:u.arcs [i].clear () ；u.multiplicity [i] =O ；v.arcs [retlndexi].clear () ； v.multiplicity [i] =0 ；.5363、w.arcs[retlndex]=w.arcs[Retlndex]+u.arcs[i]+v.arcs[retlndexj]；u.arcs [i] =u.arcs [j]+v.arcs [retlndexi]+u.arcs [k]；删除双向边u.arcs [j]以及双向边的对称双向边v.arcs [retlndexj]:u.arcs [ j].clear () ；u.multiplicity [j] =0 ； v.arcs [retlndexj].clear () ； v.multiplicity [j] =0。
6.如权利要求4所述的重复双向边识别与去除方法，其特征在于，所述步骤S37包括，.5371、如果i〈4，j〈4并且k>3，则执行步骤S372，否则执行步骤S373；.5372、删除双向边u.arcs[j]以及双向边的对称双向边v.arcs [retlndexj]:u.arcs [ j].clear () ；u.multiplicity [j] =0 ；v.arcs [retlndexj].clear () ； v.multiplicity [ j] =0 ；.5373、删除双向边u.arcs[i]以及双向边的对称双向边v.arcs [retlndexi]:u.arcs [i].clear () ；u.multiplicity [i] =0 ；v.arcs [retlndexi].clear () ； v.multiplicity [i] =0。
【文档编号】G06F19/10GK103699813SQ201310670440
【公开日】2014年4月2日 申请日期:2013年12月10日 优先权日:2013年12月10日
【发明者】孟金涛, 张慧琳, 彭丰斌, 魏彦杰, 冯圣中 申请人:深圳先进技术研究院