CN103699819B - 基于多步双向De Bruijn图的变长kmer查询的顶点扩展方法 - Google Patents

Abstract

本发明涉及基因测序技术领域，提供了一种基于多步双向De Bruijn图的变长kmer查询的顶点扩展方法，包括：步骤A：读取测序数据源文件，构造多步双向De Bruijn图；步骤B：在所述多步双向De Bruijn图中对分叉顶点的变长kmer进行构造和统计；步骤C：在所述多步双向De Bruijn图中基于变长kmer查询的顶点扩展。本发明只选取一些分叉的顶点构建非常少的一些变长kmer，然后对这些分叉顶点进行定向解耦，无需对每种kmer长度都去构建一个De Bruijn图，可以方便快速地解决所有长度小于序列长度的repeat，最大化contig的长度和质量。

Description

基于多步双向De Bruijn图的变长kmer查询的顶点扩展方法

【技术领域】

本发明涉及基因测序技术领域，特别是涉及一种基于多步双向De Bruijn图的变长kmer查询的顶点扩展方法。

【背景技术】

基因序列分析以算法与数学模型为核心，研究内容涉及多个方面，主要包括：基因数据的存储与获取、序列比对、测序与拼接、基因预测、生物进化与系统发育分析、蛋白质结构预测、RNA结构预测、分子设计与药物设计、代谢网络分析、基因芯片、DNA计算等等。现在生物技术和计算机信息处理技术的紧密结合，加快了处理生物信息数据的速度，使得在尽量短的时间内对生物学意义做出尽量准确的诠释，加快了生物信息学的发展。目前，生物信息处理成为当前信息技术领域面临的巨大挑战之一。

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

基因是人类最基本的遗传密码，代表着每个人的生命信息。基因序列上存在着遗传位点的细微差异，这些遗传密码的多态性与人类的健康、致病机理、医学治疗有着相当密切的关系。其中，DNA测序是研究全基因组序列需要完成的基本内容之一。

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

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

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

综上所述，新一代的高通量测序方法产生的短基因片段的组装导致了大量的测序错误，这大大加大了组装算法的计算量。大量的测序错误，使得组装错误率增加，严重影响了组装结果。能否有效地解决这个问题，成为评价一个组装算法优劣的关键。

目前组装算法的策略主要分为两类，一个是基于Overlap-Layout-Consensus(OLC)的算法，另外一个是基于De Bruijn图的算法。其中基于OLC组装算法开发的软件，如SSAKE、VCAKE、SHARCGS等，在基因长序列组装中更占有优势，但并不完全适用于新一代的短序列组装。与OLC组装算法不同，De Bruijn算法不再以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值越大，则出现的repeat区域则变少，这直接影响了组装的质量。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集群，在集群上建立了一个分布式的De Bruijn图结构，将数据分布式存储于每个节点上。其采用MPI通信机制完成节点之间的相互通信。从构建图、纠错处理到后面的定点融合，最后完成整个基因组序列的再现，其在运行时间和内存消耗方面占有很大的优势，并且其错误率极低，在性能方面特别是cluster中单机内存使用上均有很大的提升，正在得到越来越广泛的应用。

现有的主要序列组装软件，例如SOAPdenovo,Velvet,ABySS,Ray等，都是基于给定长度的kmer进行De Bruijn构建，然后进行收缩。其优化的办法也只有去选择最佳的一个kmer长度。这种基于固定长度kmer的组装策略对于所有长度大约kmer长度的重复序列无法解耦。虽然IDBA能够对多种kmer长度进行迭代收缩De Bruijn图，但是它需要对每种kmer长度去将所有的序列进行分解存储和计算，该策略将会耗费巨大的内存和计算时间。

鉴于此，克服该现有技术所存在的缺陷是本技术领域亟待解决的问题。

【发明内容】

本发明要解决的技术问题是提供一种基于多步双向De Bruijn图的变长kmer查询的顶点扩展方法。

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

一种基于多步双向De Bruijn图的变长kmer查询的顶点扩展方法，包括：

步骤A：读取测序数据源文件，构造多步双向De Bruijn图；

步骤B：在所述多步双向De Bruijn图中对分叉顶点的变长kmer进行构造和统计；

步骤C：在所述多步双向De Bruijn图中基于变长kmer查询的顶点扩展。

进一步地，所述步骤B中，对所述多步双向De Bruijn图中的顶点所有可能的分叉合并路径上的k+2长的变长kmer构造权重表，选择权重最高的一组分叉路径组合进行分叉路径上的双向边合并。

进一步地，所述步骤C中，对于一个给定的分叉顶点u，在查询顶点u所有的k+2长的变长kmer权重值之后，选择权重最高的一组分叉路径组合进行分叉路径上的双向边合并，同时删除合并前的被选择的分叉双向边。

进一步地，所述权重为变长kmer出现次数或变长kmer模糊匹配加权次数。

进一步地，所述步骤B进一步包括：

步骤B1：遍历所述多步双向De Bruijn图中的每个顶点u；

步骤B2：统计顶点u中正向边的个数p和反向边的个数q；

步骤B3：若p+q大于等于3且p和q均至少为1，则执行步骤B4，否则返回执行步骤B1；

步骤B4：计算出顶点u的q个反向边的对偶双向边，并将对偶双向边的倒数第k+1个字符取出存到入边字符数组m；

步骤B5：将顶点u的p个正向边的第一个字符存到出边字符数组n；

步骤B6：将（m，顶点u的正向字符串，n）所有的组合构成的k+2长的kmer记录为变长kmer数组。

进一步地，所述步骤C进一步包括：

步骤C1：打开测序序列文件，逐个读取每条序列；

步骤C2：将所述变长kmer数组逐个匹配读入的序列，并对每个变长kmer计数；

步骤C3：遍历所述多步双向De Bruijn图中的每个顶点u；

步骤C4：统计顶点u中正向边的个数p，反向边的个数q；

步骤C5：若p+q大于等于3且p和q均至少为1，则执行步骤C6，否则返回执行步骤C3；

步骤C6：计算出顶点u的q个反向边的对偶双向边，并将对偶双向边的倒数第k+1个字符取出存到入边字符数组m；

步骤C7：将顶点u的p个正向边的第一个字符存到出边字符数组n；

步骤C8：查询由（m，顶点u的正向字符串，n）的所有组合构成的k+2长的kmer的出现次数，选择出现次数最大的一组正向边和反向边进行合并扩展。

进一步地，所述步骤C2中将所述变长kmer数组逐个匹配读入的序列包括对反向序列的匹配。

进一步地，所述步骤A进一步包括：

压缩存储步骤，具体为

A11、读取一个序列s；

A12、将序列s用滑动窗口切割为多个片段t；

A13、对每个片段t，使用核酸编码表进行编码，并表示为一个64位的整数a；

A14、将片段t进行反转，使用对称互补表将反转的片段互补处理，得到互补片段v，并再次使用步骤A13中的核酸编码表将互补片段进行编码，并表示为一个64位的整数b；

A15、取整数a和整数b的最大数，作为片段t和互补片段v的k分子的标志数；

A16、重复步骤A11-A15，直至所有序列完成；

和De Bruijn图构造步骤，具体为

A21、读取一个序列s；

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

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

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

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

A26、重复步骤A22-A25，处理序列s的其他片段t，直至完成序列s的全部片段t，执行步骤A27；

A27、读取一个新的序列s，重复步骤A22-A26；直至处理完所有的序列，执行步骤A28；

A28、完成双向多步De Bruijn图的构造。

进一步地，所述步骤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中的核酸编码表将字符串v进行编码，并表示为一个64位的整数b。

进一步地，所述步骤A22中，若片段t没有之前或之后的片段，则对pre或者lat值赋为空或NULL；

步骤A24中正向位置映射表为{A:0，C:1，G:2，T:3}，位置查询字符为pre的最后一个字符；

步骤A25中反向位置映射表为{A:4，C:5，G:6，T:7}，位置查询字符为lat的第一个字符的互补字符。

与现有技术相比，本发明的有益效果在于：

（1）针对De Bruijn图上存在的分叉顶点构造长度为k+2的变长kmer组合，并在输入序列中统计其出现次数，然后根据其出现次数选择顶点上具有最大权重的双向边合并；而IDBA方法是通过迭代所有的kmer长度，需要将每个kmer长度的所有可能的kmer都构造出来后，再收缩De Bruijn图，其方法将导致更大的内存消耗和计算时间消耗；

（2）选择最优的分叉边的合并，将顶点上的分叉边的错误合并的可能降到最低；

（3）可以显著提高contig的长度，也能够将contig的质量损失降到最小；相比于其他已有方法的提高contig长度就得牺牲contig质量，本发明有了一定程度上的控制和改善。

【附图说明】

图1是本发明实施例基于多步双向De Bruijn图的变长kmer查询的顶点扩展方法流程图；

图2是图1步骤A中的压缩存储步骤流程图；

图3是图1步骤A中De Bruijn图构造步骤流程图；

图4是步骤B在多步双向De Bruijn图中对分叉顶点的变长kmer进行构造和统计的流程图；

图5是步骤C在多步双向De Bruijn图中基于变长kmer查询的顶点扩展的流程图。

【具体实施方式】

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

此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明的目的在于设计一种基于变长kmer查询的分叉顶点扩展方法，它将使De Bruijn图继续收缩，contigs继续扩展，同时不会引入错误，造成contig质量的下降，准确度降低。

本发明提供了一种基于多步双向De Bruijn图的变长kmer查询的顶点扩展方法，如图1所示，该方法包括：

步骤A：读取测序数据源文件，构造多步双向De Bruijn图；

步骤B：在多步双向De Bruijn图中对分叉顶点的变长kmer进行构造和统计；

步骤C：在多步双向De Bruijn图中基于变长kmer查询的顶点扩展。

其中，步骤A可通过如下方式具体实现：

压缩存储步骤，所需原始数据包括第一代，第二代和新一代的测序仪器产生出来的FASTA格式文件，将FASTA文件中的序列逐个切割成k分子并且用二进制编码进行压缩存储为一个64位的长整型k分子的标志数。

如图2所示，具体为

A11、读取一个序列s；其中，序列s取自FASTA格式文件；

A12、将序列s用滑动窗口切割为多个片段t；

A13、对每个片段t，使用核酸编码表进行编码，并表示为一个64位的整数a；

A14、将片段t进行反转，使用对称互补表将反转的片段互补处理，得到互补片段，并再次使用步骤A13中的核酸编码表将互补片段进行编码，并表示为一个64位的整数b；

A15、取整数a和整数b的最大数，作为片段t和互补片段v的k分子的标志数；

A16、重复步骤A11-A15，直至所有序列完成。

通过上述步骤将两个传统的De Brujin图中的kmer，转化为一个64位的k分子的标志数来存储。该步骤可以将其他软件例如velvet、IDBA、SOAPdenovo里的两个压缩kmer存储为一个压缩k分子的标志数，并且在得到k分子的标志数后也可以反过来求出该k分子的长度为k的片段t和它的互补片段v。

步骤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中的核酸编码表将字符串v进行编码，并表示为一个64位的整数b。

和De Bruijn图构造步骤，1、使用上述压缩存储步骤中计算k分子的标志数，2、将每个片段以及和它前后相邻的片段的扩展字符作为该k分子和其前后相邻的片段的对应的k分子的边并初始化k分子数据结构的边；3、将初始化后的k分子数据结构以k分子的标志数为关键值存入hash_map。

如图3所示，具体为

A21、读取一个序列s；

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

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

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

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

A26、重复步骤A22-A25，处理序列s的其他片段t，直至完成序列s的全部片段t，执行步骤S27；

A27、读取一个新的序列s，重复步骤A22-A26；直至处理完所有的序列，执行步骤A28;

A28、完成双向多步De Bruijn图的构造。

步骤A22中，若片段t没有之前或之后的片段，则对pre或者lat值赋为空或NULL；步骤A24中正向位置映射表为{A:0，C:1，G:2，T:3}，位置查询字符为pre的最后一个字符；步骤A25中反向位置映射表为{A:4，C:5，G:6，T:7}，位置查询字符为lat的第一个字符的互补字符。

步骤B中，对多步双向De Bruijn图中的顶点所有可能的分叉合并路径上的k+2长的变长kmer构造权重表，选择权重最高的一组分叉路径组合进行分叉路径上的双向边合并。权重为变长kmer出现次数或变长kmer模糊匹配加权次数。

本方法中设定多步双向De Bruijn图中的每个顶点的数据结构为：

如图4所示，步骤B进一步包括：

步骤B1：遍历多步双向De Bruijn图中的每个顶点u；

步骤B2：统计顶点u中正向边的个数p和反向边的个数q；

步骤B3：若p+q大于等于3且p和q均至少为1，则执行步骤B4，否则返回执行步骤B1；

步骤B4：计算出顶点u的q个反向边的对偶双向边，并将对偶双向边的倒数第k+1个字符取出存到入边字符数组m；

步骤B5：将顶点u的p个正向边的第一个字符存到出边字符数组n；

步骤B6：将（m，顶点u的正向字符串，n）所有的组合构成的k+2长的kmer记录为变长kmer数组。

步骤C中，对于一个给定的分叉顶点u，在查询顶点u所有的k+2长的变长kmer权重值之后，选择权重最高的一组分叉路径组合进行分叉路径上的双向边合并，同时删除合并前的被选择的分叉双向边。权重为变长kmer出现次数或变长kmer模糊匹配加权次数。

如图5所示，步骤C进一步包括：

步骤C1：打开测序序列文件，逐个读取每条序列；

步骤C2：将变长kmer数组逐个匹配读入的序列，并对每个变长kmer计数；其中将变长kmer数组逐个匹配读入的序列包括对反向序列的匹配；

步骤C3：遍历多步双向De Bruijn图中的每个顶点u；

步骤C4：统计顶点u中正向边的个数p，反向边的个数q；

步骤C5：若p+q大于等于3且p和q均至少为1，则执行步骤C6，否则返回执行步骤C3；

步骤C6：计算出顶点u的q个反向边的对偶双向边，并将对偶双向边的倒数第k+1个字符取出存到入边字符数组m；

步骤C7：将顶点u的p个正向边的第一个字符存到出边字符数组n；

步骤C8：查询由（m，顶点u的正向字符串，n）的所有组合构成的k+2长的kmer的出现次数，选择出现次数最大的一组正向边和反向边进行合并扩展。

本发明只选取一些分叉的顶点构建非常少的一些变长kmer，然后对这些分叉顶点进行定向解耦，无需对每种kmer长度都去构建一个De Bruijn图，可以方便快速地解决所有长度小于序列长度的repeat，最大化contig的长度和质量。

与现有技术相比，本发明的有益效果在于：

（1）针对De Bruijn图上存在的分叉顶点构造长度为k+2的变长kmer组合，并在输入序列中统计其出现次数，然后根据其出现次数选择顶点上具有最大权重的双向边合并；而IDBA方法是通过迭代所有的kmer长度，需要将每个kmer长度的所有可能的kmer都构造出来后，再收缩De Bruijn图，其方法将导致更大的内存消耗和计算时间消耗；

（2）选择最优的分叉边的合并，将顶点上的分叉边的错误合并的可能降到最低；

（3）可以显著提高contig的长度，也能够将contig的质量损失降到最小；相比于其他已有方法的提高contig长度就得牺牲contig质量，本发明有了一定程度上的控制和改善。

本领域普通技术人员可以理解实施例的各种方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器（ROM，Read Only Memory）、随机存取存储器（RAM，Random Access Memory）、磁盘或光盘等。

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

Claims (8)

1.一种基于多步双向De Bruijn图的变长kmer查询的顶点扩展方法，其特征在于，包括：

步骤A：读取测序数据源文件，构造多步双向De Bruijn图；

步骤B：在所述多步双向De Bruijn图中对分叉顶点的变长kmer进行构造和统计；

步骤C：在所述多步双向De Bruijn图中基于变长kmer查询的顶点扩展；所述步骤B中，对所述多步双向De Bruijn图中的顶点所有可能的分叉合并路径上的k+2长的变长kmer构造权重表，选择权重最高的一组分叉路径组合进行分叉路径上的双向边合并；所述步骤C中，对于一个给定的分叉顶点u，在查询顶点u所有的k+2长的变长kmer权重值之后，选择权重最高的一组分叉路径组合进行分叉路径上的双向边合并，同时删除合并前的被选择的分叉双向边。

2.如权利要求1所述的方法，其特征在于，所述权重为变长kmer出现次数或变长kmer模糊匹配加权次数。

3.如权利要求1所述的方法，其特征在于，所述步骤B进一步包括：

步骤B1：遍历所述多步双向De Bruijn图中的每个顶点u；

步骤B2：统计顶点u中正向边的个数p和反向边的个数q；

步骤B3：若p+q大于等于3且p和q均至少为1，则执行步骤B4，否则返回执行步骤B1；

步骤B4：计算出顶点u的q个反向边的对偶双向边，并将对偶双向边的倒数第k+1个字符取出存到入边字符数组m；

步骤B5：将顶点u的p个正向边的第一个字符存到出边字符数组n；

步骤B6：将(m，顶点u的正向字符串，n)所有的组合构成的k+2长的kmer记录为变长kmer数组。

4.如权利要求1所述的方法，其特征在于，所述步骤C进一步包括：

步骤C1：打开测序序列文件，逐个读取每条序列；

步骤C2：将所述变长kmer数组逐个匹配读入的序列，并对每个变长kmer计数；

步骤C3：遍历所述多步双向De Bruijn图中的每个顶点u；

步骤C4：统计顶点u中正向边的个数p，反向边的个数q；

步骤C5：若p+q大于等于3且p和q均至少为1，则执行步骤C6，否则返回执行步骤C3；

步骤C6：计算出顶点u的q个反向边的对偶双向边，并将对偶双向边的倒数第k+1个字符取出存到入边字符数组m；

步骤C7：将顶点u的p个正向边的第一个字符存到出边字符数组n；

步骤C8：查询由(m，顶点u的正向字符串，n)的所有组合构成的k+2长的kmer的出现次数，选择出现次数最大的一组正向边和反向边进行合并扩展。

5.如权利要求4所述的方法，其特征在于，所述步骤C2中将所述变长kmer数组逐个匹配读入的序列包括对反向序列的匹配。

6.如权利要求1所述的方法，其特征在于，所述步骤A进一步包括：

压缩存储步骤，具体为

A11、读取一个序列s；

A12、将序列s用滑动窗口切割为多个片段t；

A13、对每个片段t，使用核酸编码表进行编码，并表示为一个64位的整数a；

A14、将片段t进行反转，使用对称互补表将反转的片段互补处理，得到互补片段v，并再次使用步骤A13中的核酸编码表将互补片段进行编码，并表示为一个64位的整数b；

A15、取整数a和整数b的最大数，作为片段t和互补片段v的k分子的标志数；

A16、重复步骤A11-A15，直至所有序列完成；

和De Bruijn图构造步骤，具体为

A21、读取一个序列s；

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

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

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

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

A26、重复步骤A22-A25，处理序列s的其他片段t，直至完成序列s的全部片段t，执行步骤A27；

A27、读取一个新的序列s，重复步骤A22-A26；直至处理完所有的序列，执行步骤A28；

A28、完成双向多步De Bruijn图的构造。

7.如权利要求6所述的方法，其特征在于，所述步骤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中的核酸编码表将字符串v进行编码，并表示为一个64位的整数b。

8.如权利要求6所述的方法，其特征在于，所述步骤A22中，若片段t没有之前或之后的片段，则对pre或者lat值赋为空或NULL；

步骤A24中正向位置映射表为{A:0，C:1，G:2，T:3}，位置查询字符为pre的最后一个字符；

步骤A25中反向位置映射表为{A:4，C:5，G:6，T:7}，位置查询字符为lat的第一个字符的互补字符。