基于多步双向De Bruijn图的变长kmer查询的双向边扩展方法
【专利摘要】本发明涉及基因测序【技术领域】，提供了一种基于多步双向De?Bruijn图的变长kmer查询的顶点扩展方法，包括：步骤A：读取测序数据源文件，构造多步双向De?Bruijn图；步骤B：在所述多步双向De?Bruijn图中对交叉双向边进行统计；步骤C：在所述多步双向De?Bruijn图中基于变长kmer查询的双向边扩展。本发明只针对De?Bruijn图上已经存在的分叉边构造变长的kmer组合，并在输入序列中查询其出现次数，然后根据其次数选择双向边的合并；选择最优的分叉边的合并，将双向边的错误合并的可能降到最低；可以显著提高contigs的长度，也能够将contig的质量损失降到最小。
【专利说明】基于多步双向De Bru i jn图的变长kmer查询的双向边扩展方法
【【技术领域】】
[0001]本发明涉及基因测序【技术领域】，特别是涉及一种基于多步双向De Bruijn图的变长kmer查询的双向边扩展方法。
【【背景技术】】
[0002]基因序列分析以算法与数学模型为核心，研究内容涉及多个方面，主要包括:基因数据的存储与获取、序列比对、测序与拼接、基因预测、生物进化与系统发育分析、蛋白质结构预测、RNA结构预测、分子设计与药物设计、代谢网络分析、基因芯片、DNA计算等等。现在生物技术和计算机信息处理技术的紧密结合，加快了处理生物信息数据的速度，使得在尽量短的时间内对生物学意义做出尽量准确的诠释，加快了生物信息学的发展。目前，生物信息处理成为当前信息【技术领域】面临的巨大挑战之一。
[0003]基因序列分析是对海量基因序列数据进行分析，从而提取和挖掘新的生物信息知识。其中，涉及到计算机技术中的机器学习、模式识别、书籍分析与挖掘、组合数学、随机模型、字符串、图形算法、分布式计算、高性能计算、并行计算等知识。其中，全基因组学的研究是当前生物信息学研究的核心之一。
[0004]基因是人类最基本的遗传密码，代表着每个人的生命信息。基因序列上存在着遗传位点的细微差异，这些遗传密码的多态性与人类的健康、致病机理、医学治疗有着相当密切的关系。其中，DNA测序 是研究全基因组序列需要完成的基本内容之一。
[0005]自1977年Sanger测序技术问世以来,经过三十多年的发展,DNA测序技术发展突飞猛进，以高通量、短序列为特点的第二代测序技术逐渐占领市场，以单分子测序为特点的第三代测序技术也逐渐出现，它们分别在测序特点上占有不同的优势。传统的基因测序方法的数据提取和分析软件经过近10年来的研究与开发，目前已经较为完善。但是，测序技术的发展，带来了测序数据的变化，使得当前存在的数据处理软件不能满足当前生物医学研究的需求。
[0006]新一代高通量测序方法的应用，可以在短时间内完成整个基因组数据的测定。高通量测序方法的日新月异也同时对获取的基因数据的分析处理方法提出了挑战。在这个目前炙手可热的研究领域，迫切需要开发能满足高通量测序技术的海量数据处理的生物信息学平台。面对个人基因组计划及未来的个性化医疗前景，高效低成本的测序技术成为必然的趋势。同时，简化高效的一站式完备的生物信息学数据分析平台等完备的测序解决方案，也是极为重要不可或缺的发展方向。
[0007]然而新一代的高通量测序方法虽然测序通量高，但是却会引入测序误差，同时测序样本本身的测序错误，测序不均匀，SNP,以及基因组本身的重复序列R印eat，而这些测序误差、SNP、重复序列将会在基因组组装时构造的多步双向De Bruijn图中引入一些错误双向边或者顶点，而使得很多双向边无法继续扩展。同时由于固定kmer长度，使得测序序列的有效信息损失，无法解耦所有长度超过kmer长度的重复序列。以上这些情况使得DeBruijn图无法继续收缩，contig无法扩展,最终使得contig的长度和质量都很低。
[0008]新一代的高通量测序方法产生的短基因片段的组装导致了大量的测序错误，这大大加大了组装算法的计算量。大量的测序错误，使得组装错误率增加，严重影响了组装结果。能否有效地解决这个问题，成为评价一个组装算法优劣的关键。
[0009]目前组装算法的策略主要分为两类，一个就是基于Overlap-Layout-Consensus (OLC)的算法，另外一个就是基于De Bruijn图的算法。其中基于OLC组装算法开发的软件，如SSAKE、VCAKE, SHARCGS等，在基因长序列组装中更占有优势，但并不完全适用于新一代的短序列组装。与OLC组装算法不同，De Bruijn算法不再以read为单位组织数据，而是以k-mers为单位进行数据组装，其优点主要有以下几个方面:首先，以k-mers为单位进行序列组装，不影响节点的质量，减少了冗余数据量。其次，在图中重复区域只出现一次，便于识别，可以避免错误的组装，减小出错率。最后，采取将有重叠区域映射到同一条弧上的策略，从而简化了搜索路径。目前，很多短序列组装算法都使用这种框架，如 Velvet、IDBA、SOAPdenovo, ABySS 等。
[0010]其中Velvet有效地利用了 De Brui jn图，实现了高效的短序列组装。Velvet以k-mer为基本单位构建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值越大,贝U出现的repeat区域则变少,这直接影响了组装的质量。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集群，在集群上建立了一个分布式的De Bruijn图结构，将数据分布式存储于每个节点上。其采用MPI通信机制完成节点之间的相互通信。从构建图、纠错处理到后面的定点融合，最后完成整个基因组序列的再现，其在运行时间和内存消耗方面占有很大的优势，并且其错误率极低，在性能方面特别是cluster中单机内存使用上均有很大的提升，正在得到越来越广泛的应用。
[0014]现有的主要序列组装软件，例如SOAPdenovo, Velvet, ABySS, Ray等，都是基于给定长度的kmer进行De Brui jn构建,然后进行收缩。其优化的办法也只有去选择最佳的一个kmer长度。这种基于固定长度kmer的组装策略对于所有长度大约kmer长度的重复序列无法解耦。虽然IDBA能够对多种kmer长度进行迭代收缩De Bruijn图，但是它需要对每种kmer长度去将所有的序列进行分解存储和计算，该策略将会耗费巨大的内存和计算时间。
[0015]鉴于此，克服该现有技术所存在的缺陷是本【技术领域】亟待解决的问题。
【
【发明内容】
】
[0016]本发明要解决的技术问题是提供一种基于多步双向De Bruijn图的变长kmer查询的双向边扩展方法。
[0017]本发明采用如下技术方案:
[0018]一种基于多步双向De Bruijn图的变长kmer查询的双向边扩展方法,包括:
[0019]步骤A:读取测序数据源文件，构造多步双向De Bruijn图；
[0020]步骤B:在所述多步双向De Bruijn图中对交叉双向边进行统计；
[0021]步骤C:在所述多步双向De Bruijn图中基于变长kmer查询的双向边扩展。
[0022]进一步地,所述步骤B中，对所述多步双向De Brui jn图中的所有左右两边顶点分叉的双向边与所述双向边的分叉边合并的所有可能情况，组装为一个更长的kmer，并查询相应地更长的kmer的出现次数，以所述出现次数的大小为权重，进行分叉边走向的选取和合并操作。
[0023]进一步地，所述步骤C中，在对每个分叉的双向边的所有可能的分叉转换为更长的kmer后，依据所述出现次数的高低，将出现次数最高的分叉路径作为最终的分叉合并方案，并将相应的分叉路径上的双向边在合并后一起删除。
[0024]进一步地,所述步骤B进一步包括:
[0025]步骤B1:遍历De Bruijn图中的每个双向边e ；
[0026]步骤B2:统计所有的输入文件的测序序列的最大长度，并赋值给Lmax ；
[0027]步骤B3:若双向边的长度大于Lmax-kmer的长度,则返回执行步骤BI,否则执行步骤B4 ；
[0028]步骤B4:双向边e的起始顶点为u,起始方向为Direc_u,终止顶点为V,终止方向为 Direc_v ；
[0029]步骤B5:所有进入顶点U、终止方向为Direc_u的边取其双向边的最后k+Ι个字符放入字符串数组m,所有离开顶点V、起始方向为Direc_v的双向边取其最后一个字符放入字符串数组η ;
[0030]步骤Β6:将双向边左右分别和m和η数组中的字符串的构造所有组合，构造变长kmer数组,其中变长kmer数组中的每个元素为变长kmer=m+e+n。[0031]进一步地,所述步骤C进一步包括:
[0032]步骤Cl:打开测序序列文件，逐个读取每条序列；
[0033]步骤C2:将变长kmer集合组各匹配读入的序列(包括对其反向序列的匹配)，并对每个变长kmer计数；
[0034]步骤C3:遍历De Bruijn图中的每个双向边e ；
[0035]步骤C4:统计所有的输入文件的测序序列的最大长度，并赋值给Lmax ；
[0036]步骤C5:若双向边的长度大于Lmax-kmer的长度,则返回执行步骤C3,否则执行步骤C6 ；
[0037]步骤C6:双向边e的起始顶点为u,起始方向为Direc_u,终止顶点为V,终止方向为 Direc_v ；
[0038]步骤C7:所有进入顶点U、终止方向为Direc_u的边取其最后双向边e_i的最后k+Ι个字符放入字符串数组m,所有离开顶点V、起始方向为Direc_v的双向边e_j取其最后一个字符放入字符串数组η ;
[0039]步骤C8:将双向边e左右分别和m和η数组中的字符串的构造所有组合，构造变长kmer数组,其中变长kmer数组中的每个元素为变长kmer=m+e+n。在计数的变长kmer中查询这些变长kmer，并取出计数值最大的一个，并将相应的e_i，e，e_j三条边进行合并，然后删除 e_i, e, e_j。
[0040]进一步地，所述步骤C2中将变长kmer集合组各匹配读入的序列包括对反向序列的匹配。
[0041]进一步地,所述步骤A进一步包括:
[0042]压缩存储步骤，具体为
[0043]All、读取一个序列s ;
[0044]A12、将序列s用滑动窗口切割为多个片段t ；
[0045]A13、对每个片段t，使用核酸编码表进行编码，并表示为一个64位的整数a ；
[0046]A14、将片段t进行反转，使用对称互补表将反转的片段互补处理，得到互补片段V,并再次使用步骤A13中的核酸编码表将互补片段进行编码，并表示为一个64位的整数b ；
[0047]A15、取整数a和整数b的最大数，作为片段t和互补片段ν的k分子的标志数；
[0048]A16、重复步骤A11-A15，直至所有序列完成；
[0049]和De Brui jn图构造步骤,具体为
[0050]A21、读取一个序列s;
[0051]A22、将序列s用滑动窗口切割为多个片段t，选取一片段t其标志数为cur、并标记其前、后的片段的标志数分别为pre、Iat ；
[0052]A23、若t的编码小于其互补片段编码，则交换pre，Iat的值；
[0053]A24、在cur的正向位置映射表的相应bit位置I来表示指向pre的边；
[0054]A25、在cur的反向位置映射表的相应bit位置I来表示指向Iat的边；
[0055]A26、重复步骤A22-A25，处理序列s的其他片段t，直至完成序列s的全部片段t，执行步骤A27 ；
[0056]A27、读取一个新的序列S，重复步骤A22-A26 ;直至处理完所有的序列，执行步骤A28 ；
[0057]A28、完成双向多步De Bruijn图的构造。
[0058]进一步地，所述步骤A12、A22中的滑动窗口为长度为k的滑动窗口，其中0〈k〈32且k为奇数；
[0059]所述步骤A13中的核酸编码表为{A:00，C:01，G:10，T:11}；
[0060]所述步骤A14中的对称互补表为{A->T，C->G, G->C, T_>A}；
[0061]所述步骤A14具体为，将片段t的字符串进行反转，使用对称互补表将反转的字符串中每个字符变为其互补字符，得到互补字符的字符串V，并再次使用步骤A13中的核酸编码表将字符串ν进行编码，并表示为一个64位的整数b。
[0062]进一步地，所述步骤A22中，若片段t没有之前或之后的片段，则对pre或者Iat值赋为空或NULL ；
[0063]步骤A24中正向位置映射表为{A:0，C:l，G:2，T:3}，位置查询字符为pre的最后一个字符；
[0064]步骤A25中反向位置映射表为{A:4，C:5, G:6, T:7}，位置查询字符为Iat的第一个字符的互补字符。
[0065]与现有技术相比， 本发明的有益效果在于:
[0066](I)只针对De Bruijn图上已经存在的分叉边构造变长的kmer组合，并在输入序列中查询其出现次数，然后根据其次数选择双向边的合并；而IDBA则是通过迭代所有的kmer长度,需要将每个kmer长度的所有可能的kmer都构造出来后,再收缩De Bruijn图，其方法将导致更大的内存消耗和计算时间消耗；
[0067](2)选择最优的分叉边的合并，将双向边的错误合并的可能降到最低；
[0068](3)可以显著提高contigs的长度,也能够将contig的质量损失降到最小；相比于其他已有方法的提高contig长度就要牺牲contig质量，本发明有了一定程度上的控制和改善。
【【专利附图】
【附图说明】】
[0069]图1是本发明实施例基于多步双向De Bruijn图的变长kmer查询的双向边扩展方法流程图；
[0070]图2是图1步骤A中的压缩存储步骤流程图；
[0071]图3是图1步骤A中De Bruijn图构造步骤流程图；
[0072]图4是步骤B在所述多步双向De Bruijn图中对交叉双向边进行统计的流程图；
[0073]图5是步骤C在所述多步双向De Brui jn图中基于变长kmer查询的双向边扩展的流程图。
【【具体实施方式】】
[0074]为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0075]此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
[0076]本发明的目的在于设计一种基于变长kmer查询的交叉双向边的扩展方法，它将使De Bruijn图继续收缩,contigs继续扩展，同时不会引入错误,造成contig质量的下降，准确度降低。
[0077]本发明提供了一种基于多步双向De Bruijn图的变长kmer查询的双向边扩展方法,如图1所示，该方法包括:
[0078]步骤A:读取测序数据源文件，构造多步双向De Bruijn图；
[0079]步骤B:在多步双向De Bruijn图中对交叉双向边进行统计；
[0080]步骤C:在多步双向De Bruijn图中基于变长kmer查询的双向边扩展。
[0081]其中，步骤A可通过如下方式具体实现:
[0082]压缩存储步骤，所需原始数据包括第一代，第二代和新一代的测序仪器产生出来的FASTA格式文件，将FASTA文件中的序列逐个切割成k分子并且用二进制编码进行压缩存储为一个64位的长整型k分子的标志数。
[0083]如图2所示，具体为
[0084]All、读取一个序列s ;其中，序列s取自FASTA格式文件；
[0085]A12、将序列s用滑动窗口切割为多个片段t ；
[0086]A13、对每个片段t，使用核酸编码表进行编码，并表示为一个64位的整数a ；
[0087]A14、将片段t进行反转，使用对称互补表将反转的片段互补处理，得到互补片段，并再次使用步骤A13中的核酸编码表将互补片段进行编码，并表示为一个64位的整数b ;
[0088]A15、取整数a和整数b的最大数，作为片段t和互补片段ν的k分子的标志数；
[0089]A16、重复步骤A11-A15，直至所有序列完成。
[0090]通过上述步骤将两个传统的De Brujin图中的kmer,转化为一个64位的k分子的标志数来存储。该步骤可以将其他软件例如velvet、IDBA、S0APdenovo里的两个压缩kmer存储为一个压缩k分子的标志数，并且在得到k分子的标志数后也可以反过来求出该k分子的长度为k的片段t和它的互补片段V。
[0091]步骤A12、A22中的滑动窗口为长度为k的滑动窗口，其中0〈k〈32且k为奇数；步骤A13中的核酸编码表为{A:00，C:01，G:10，T:11};步骤A14中的对称互补表为{A->T, C->G, G->C, T->A};步骤A14具体为，将片段t的字符串进行反转，使用对称互补表将反转的字符串中每个字符变为其互补字符，得到互补字符的字符串V，并再次使用步骤A13中的核酸编码表将字符串ν进行编码，并表示为一个64位的整数b。
[0092]和De Bruijn图构造步骤，1、使用上述压缩存储步骤中计算k分子的标志数,2、将每个片段以及和它前后相邻的片段的扩展字符作为该k分子和其前后相邻的片段的对应的k分子的边并初始化k分子数据结构的边；3、将初始化后的k分子数据结构以k分子的标志数为关键值存入hash_map。
[0093]如图3所示，具体为
[0094]A21、读取一个序列s;
[0095]A22、将序列s用滑动窗口切割为多个片段t，选取一片段t其标志数为cur、并标记其前、后的片段的标志数分别为pre、Iat ；
[0096]A23、若t的编码小于其互补片段编码，则交换pre，Iat的值；[0097]A24、在cur的正向位置映射表的相应bit位置I来表示指向pre的边；
[0098]A25、在cur的反向位置映射表的相应bit位置I来表示指向Iat的边；
[0099]A26、重复步骤A22-A25，处理序列s的其他片段t，直至完成序列s的全部片段t，执行步骤S27 ；
[0100]A27、读取一个新的序列S，重复步骤A22-A26 ;直至处理完所有的序列，执行步骤A28;
[0101]A28、完成双向多步De Bruijn图的构造。
[0102]步骤A22中，若片段t没有之前或之后的片段，则对pre或者Iat值赋为空或NULL ；步骤A24中正向位置映射表为{A:0，C:1，G:2，T:3}，位置查询字符为pre的最后一个字符；步骤A25中反向位置映射表为{A:4，C:5, G:6, T:7}，位置查询字符为Iat的第一个字符的
互补字符。
[0103]步骤B中，对所述多步双向De Brui jn图中的所有左右两边顶点分叉的双向边与所述双向边的分叉边合并的所有可能情况，组装为一个更长的kmer，并查询相应地更长的kmer的出现次数，以所述出现次数的大小为权重，进行分叉边走向的选取和合并操作。
[0104]本方法中设定 多步双向De Bruijn图中的每个顶点的数据结构为:
[0105]
【权利要求】
1.一种基于多步双向De Bruijn图的变长kmer查询的双向边扩展方法,其特征在于，包括: 步骤A:读取测序数据源文件，构造多步双向De Bruijn图； 步骤B:在所述多步双向De Bruijn图中对交叉双向边进行统计； 步骤C:在所述多步双向De Bruijn图中基于变长kmer查询的双向边扩展。
2.如权利要求1所述的方法，其特征在于，所述步骤B中，对所述多步双向DeBruijn图中的所有左右两边顶点分叉的双向边与所述双向边的分叉边合并的所有可能情况，组装为一个更长的kmer，并查询相应地更长的kmer的出现次数，以所述出现次数的大小为权重，进行分叉边走向的选取和合并操作。
3.如权利要求2所述的方法，其特征在于，所述步骤C中，在对每个分叉的双向边的所有可能的分叉转换为更长的kmer后，依据所述出现次数的高低，将出现次数最高的分叉路径作为最终的分叉合并方案，并将相应的分叉路径上的双向边在合并后一起删除。
4.如权利要求1所述的方法，其特征在于，所述步骤B进一步包括: 步骤B1:遍历De Brui jn图中的每个双向边e ; 步骤B2:统计所有的输入文件的测序序列的最大长度，并赋值给Lmax ； 步骤B3:若双向边的长度大于Lmax-kmer的长度，则返回执行步骤BI，否则执行步骤B4 ； 步骤B4:双向边e的起始顶点为u,起始方向为Direc_u,终止顶点为V,终止方向为Direc_v ； 步骤B5:所有进入顶点U、终止方向为Direc_u的边取其双向边的最后k+Ι个字符放入字符串数组m,所有离开顶点V、起始方向为Direc_v的双向边取其最后一个字符放入字符串数组η ； 步骤Β6:将双向边左右分别和m和η数组中的字符串的构造所有组合，构造变长kmer数组，其中变长kmer数组中的每个元素为变长kmer=m+e+n。
5.如权利要求1所述的方法，其特征在于，所述步骤C进一步包括: 步骤Cl:打开测序序列文件，逐个读取每条序列； 步骤C2:将变长kmer集合组各匹配读入的序列(包括对其反向序列的匹配)，并对每个变长kmer计数； 步骤C3:遍历De Bruijn图中的每个双向边e ； 步骤C4:统计所有的输入文件的测序序列的最大长度，并赋值给Lmax ； 步骤C5:若双向边的长度大于Lmax-kmer的长度，则返回执行步骤C3，否则执行步骤C6 ； 步骤C6:双向边e的起始顶点为u,起始方向为Direc_u,终止顶点为V,终止方向为Direc_v ； 步骤C7:所有进入顶点U、终止方向为Direc_u的边取其最后双向边e_i的最后k+Ι个字符放入字符串数组m,所有离开顶点V、起始方向为Direc_v的双向边e_j取其最后一个字符放入字符串数组η ; 步骤C8:将双向边左右分别和m和η数组中的字符串的构造所有组合，构造变长kmer数组，其中变长kmer数组中的每个元素为变长kmer=m+e+n,在计数的变长kmer中查询这些变长kmer，并取出计数值最大的一个，并将相应的e_i，e，e_j三条边进行合并，然后删除e_i，e，e_j。
6.如权利要求5所述的方法，其特征在于，所述步骤C2中将变长kmer集合组各匹配读入的序列包括对反向序列的匹配。
7.如权利要求1所述的方法，其特征在于，所述步骤A进一步包括: 压缩存储步骤，具体为 Al 1、读取一个序列s ； A12、将序列s用滑动窗口切割为多个片段t ; A13、对每个片段t，使用核酸编码表进行编码，并表示为一个64位的整数a ； A14、将片段t进行反转，使用对称互补表将反转的片段互补处理，得到互补片段V，并再次使用步骤A13中的核酸编码表将互补片段进行编码，并表示为一个64位的整数b ;A15、取整数a和整数b的最大数，作为片段t和互补片段ν的k分子的标志数； A16、重复步骤A11-A15,直至所有序列完成； 和De Brui jn图构造步骤,具体为 A21、读取一 个序列s ； A22、将序列s用滑动窗口切割为多个片段t，选取一片段t其标志数为cur、并标记其前、后的片段的标志数分别为pre、Iat ； A23、若t的编码小于其互补片段编码，则交换pre，Iat的值； A24、在cur的正向位置映射表的相应bit位置I来表示指向pre的边； A25、在cur的反向位置映射表的相应bit位置I来表示指向Iat的边； A26、重复步骤A22-A25，处理序列s的其他片段t，直至完成序列s的全部片段t，执行步骤A27 ； A27、读取一个新的序列S，重复步骤A22-A26 ;直至处理完所有的序列，执行步骤A28 ； A28、完成双向多步De Bruijn图的构造。
8.如权利要求7所述的方法，其特征在于，所述步骤A12、A22中的滑动窗口为长度为k的滑动窗口，其中0〈k〈32且k为奇数； 所述步骤A13中的核酸编码表为{A:00，C:01，G:10，T:11}； 所述步骤A14中的对称互补表为{A->T，C->G，G->C，T->A}； 所述步骤A14具体为，将片段t的字符串进行反转，使用对称互补表将反转的字符串中每个字符变为其互补字符，得到互补字符的字符串V，并再次使用步骤A13中的核酸编码表将字符串ν进行编码，并表示为一个64位的整数b。
9.如权利要求7所述的方法，其特征在于，所述步骤A22中，若片段t没有之前或之后的片段，则对pre或者Iat值赋为空或NULL ； 步骤A24中正向位置映射表为{A:0，C:l, G:2, Τ: 3}，位置查询字符为pre的最后一个字符； 步骤A25中反向位置映射表为{A:4，C:5, G:6, T:7}，位置查询字符为Iat的第一个字符的互补字符。
【文档编号】G06F19/22GK103699818SQ201310670740
【公开日】2014年4月2日 申请日期:2013年12月10日 优先权日:2013年12月10日
【发明者】孟金涛, 张慧琳, 彭丰斌, 魏彦杰, 冯圣中 申请人:深圳先进技术研究院