[发明专利]双向多步deBruijn图的压缩存储和构造方法

说明书

【技术领域】

本发明涉及基因测序领域，尤其涉及一种双向多步deBruijn图的压缩存储和构造方法。

【背景技术】

新一代基因测序技术所产生的序列片段具有序列短、高覆盖率、额外的双端信息等特点，使得原有传统的序列拼接技术无法使用，加大了序列拼接技术的复杂度和难点。面对海量的数据，高效的序列拼接技术成为处理测序数据的关键。

序列拼接技术问题主要转化为图论的知识进行求解。当前，序列拼接技术问题解决途径主要有两条。一种途径是传统的先重叠后扩展方法，即OLC（Overlap-Layout-Consensus）方法。该方法主要是将序列拼接技术转化为哈密尔顿路径进行求解。OLC的拼接过程主要分为三个步骤：1）Overlap：对获得的所有read顶点进行排序，通过序列比对算法，寻找DNA序列片段间的重叠信息；2）Layout：根据read之间的重叠信息，排列所有DNA序列片段，形成新的链接体结构，该步骤主要是为了寻找经过每个顶点一次且仅一次的路径，即寻找哈密尔顿路径；3）Consensus：根据新的链接体结构中原始质量数据，在链接体中寻找质量最重的路径，从而组合成最终的DNA序列。

基于OLC的序列拼接在技术上容易实现，但是该方法仅局限于比对大于某个阈值的read之间的信息，忽略了多个read之间的相互信息，从而使得该方法在处理重复区域问题上受到了很大的限制；另外，该方法进行序列比对时，无论使用贪心算法或者BWA算法，都需要耗费大量的内存，同时在任意两个read序列之间比对，使得该序列拼接的算法复杂度为O(n²)，其中n为序列片段的个数。因此，从时间和空间上考虑，该算法均难以用于拼接长达百万碱基的基因组样本或者由百万条以上的DNA序列片段组成的测序数据。基于OLC的拼接主要适用于基于Sanger测序原理的拼接。测序样本也主要是基因组较小的生物,例如病毒、真菌。目前，PHRAP、TIGR、CAP3、CELERA、ARACHNE、PHUSION、SSAKE、VCAKE、SHARCGS等序列拼接技术均是在此基础上开发的。

DNA序列拼接的另外一种途径，是基于de Bruijn图的序列拼接方法。该方法主要是将DNA序列拼接问题转化为欧拉路径进行求解；通过构造并简化deBruijn图来实现整个拼接过程。基于de Bruijn图的拼接方法理论上需要找到一条欧拉路径来重构源基因组参考序列，但实际上只能通过将构造好的de Bruijn图中的无分支路径进行收缩最终合并为完整的contig进行输出。目前，Velvet、Soapdenovo、Idba、Abyss、Pasha等序列拼接技术均是在此基础上开发的。

与OLC相比，基于de Bruijn图的序列拼接技术有很多的优点。例如，在过滤错误信息、重复区域发现、解耦以及利用双端信息对contig进行扩展等方面。然而，由于大基因组构造的de Bruijn图异常庞大，以至于单个计算机的内存无法释放。例如，对人类基因组的数据进行拼接，构建的de Bruijn图中有30G的顶点，而存储这样的图通常需要消耗大约500G到1T的内存。同时，由于在图的数据结构中顶点关联的随机性，不能使用数据预取技术，导致最终对图的简化也将耗费几周的时间。在拼接大基因组和宏基因组测序数据时，上述拼接技术的性能依然无法满足需求。

【发明内容】

本发明要解决的技术问题在于提高现有技术序列拼接运行速度、降低单机内存消耗。

为此，本发明提供一种双向多步deBruijn图的压缩存储和构造方法，包括

压缩存储步骤，具体为

S11、读取一个序列s；

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

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

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

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

S16、重复步骤S11-S15，直至所有序列完成；

和deBruijn图构造步骤，具体为

S21、读取一个序列s；

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

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