本发明涉及测序技术领域，尤其涉及一种二代序列的短序列纠错的方法和装置。
背景技术：
目前，基因组组装项目以全基因组鸟枪法测序(whole-genomeshotgunsequencing,wgs)为主流设计方案，它主要根据基因组的重复序列的具体特点，搭配不同长度的dna插入片段进行双末端测序，在全基因组的平均测序深度足够的情况下可保证单碱基的准确性和基因组的完整性。随着第二代测序技术(next-generationsequencing,ngs)的成熟和普及，测序成本大大降低，基于第二代测序技术的全基因组鸟枪法测序成为各种基因组项目的测序的主流方案。
基因组第二代序列组装主要基于德布鲁恩方法进行组装，组装思路主要是将二代测序序列依次截取出长度为k的短序列k-mer；将k-mer存储到散列表中，形成德布鲁因图的顶点；在测序序列上前后相继的k-mer相连，形成德布鲁因图的边；将所有测序序列都处理完得到整个德布鲁因图；去除德布鲁因图中由测序错误、杂合位点引起的路径；将线性的k-mer路径连接起来形成的重叠群。由于第二代测序技术所得的序列一般存在1％错误的问题，基因组越大，测序的数据量越多，这些错误会大大增加k-mer的种类，进而大大增加将k-mer存储到散列表中的内存，所需要的内存峰值就越大，内存峰值可能会超过几百gb，甚至1t以上，这就对计算机的内存要求非常高。为了降低内存峰值，有必要在基因组组装之前对二代序列进行纠错。
技术实现要素：
本发明提供一种二代序列的短序列纠错的方法和装置，能够大大降低二代序列基于德布鲁因组装所需要的内存峰值，从而能够降低基因组组装的难度。
根据本发明的第一方面，本发明提供一种二代序列的短序列纠错的方法，包括：获取二代序列的k-mer序列；对上述k-mer序列进行分类，得到多个分类文件；对每个上述分类文件，统计上述k-mer序列的频数；筛选并获得低于预定频数的k-mer序列；用上述预定频数的k-mer序列对上述二代序列进行纠错；其中，上述纠错包括：将上述预定频数的k-mer序列比对回上述二代序列，若比对上，则从该k-mer头尾两端对应在上述二代序列上的位置将上述二代序列截断；保留截断后的长度大于预定长度的序列，舍弃截断后的长度小于上述预定长度的序列。
进一步地，上述方法在上述分类和上述统计步骤之间，还包括：压缩上述k-mer序列。
进一步地，上述压缩上述k-mer序列包括：对上述k-mer序列中的每三个碱基用1位ascii字符替代。
进一步地，上述k-mer序列的长度范围是17至75个碱基长度。
进一步地，上述k-mer序列的长度是3的整数倍，优选39个碱基长度。
进一步地，上述对上述k-mer序列进行分类包括：截取上述k-mer序列的连续n位碱基序列，其中含n的k-mer序列去掉，将上述k-mer序列分类到4的n次方份文件中。
进一步地，上述预定频数是4。
进一步地，上述预定长度是80-120个碱基长度，优选100个碱基长度。
根据本发明的第二方面，本发明提供一种二代序列的短序列纠错的装置，包括：获取单元，用于获取二代序列的k-mer序列；分类单元，用于对上述k-mer序列进行分类，得到多个分类文件；统计单元，用于对每个上述分类文件，统计上述k-mer序列的频数；筛选单元，用于筛选并获得低于预定频数的k-mer序列；纠错单元，用于用上述预定频数的k-mer序列对上述二代序列进行纠错；其中，上述纠错包括：将上述预定频数的k-mer序列比对回上述二代序列，若比对上，则从该k-mer头尾两端对应在上述二代序列上的位置将上述二代序列截断；保留截断后的长度大于预定长度的序列，舍弃截断后的长度小于上述预定长度的序列。
进一步地，上述装置还包括：压缩单元，用于压缩分类过的k-mer序列。
本发明的方法和装置，主要对二代序列的k-mer进行分类，并统计每种k-mer的频数，并筛选出低频k-mer从而对二代序列进行纠错。能够大大降低二代序列基于德布鲁因组装所需要的内存峰值，从而能够降低基因组组装的难度。
附图说明
图1示出本发明一个实施例的二代序列k-mer纠错方法的流程图；
图2示出本发明一个实施例的二代序列组成和获取二代序列的k-mer序列示意图，每个二代序列包括四行，第一行是序列id，第二行是序列的碱基信息，第三行是“+”号，第四行是第二行每个碱基对应的测序质量值；
图3示出本发明一个实施例中用低频k-mer序列对二代序列进行纠错的原理示意图；
图4示出本发明一个实施例的二代序列k-mer纠错装置的结构框图；
图5示出本发明一个实施例的k-39序列频数分布图。
具体实施方式
下面通过具体实施方式结合附图对本发明作进一步详细说明。
在本发明的一个实施例，提供一种二代序列k-mer纠错的方法，旨在降低二代序列的错误率，同时降低基于德布鲁因组装所需要的内存峰值，从而降低基因组组装的难度。
本发明实施例中，二代序列是指第二代测序技术产生的测序序列，也称读长(reads)，目前二代序列的读长主要有100bp、150bp和250bp三种。
k-mer，即长度为k的短序列，其是从二代序列中截取出来的，前后相继的k-mer错位一个碱基。k-mer序列的长度范围一般是17至75个碱基长度，k-mer序列的长度最好是3的整数倍，例如39个碱基长度的k-mer，可以称为k-39。
图1示出了本发明一个实施例的二代序列k-mer纠错方法的流程图。
如图1所示，在步骤102中，获取二代序列的k-mer序列。
结合二代序列读长，选取一定长度的k-mer对二代序列进行分割并保存，k-mer值一般选取17到75的范围，另外为了方便后续的压缩，一般k-mer值选取为3的整数倍。
图2示出了本发明一个实施例中的二代序列组成和获取二代序列的k-mer序列，每个二代序列包括四行，第一行是序列id，第二行是序列的碱基信息，第三行是“+”号，第四行是第二行每个碱基对应的测序质量值。选取39个碱基长度的k-mer(k-39)对二代序列进行分割并保存，前后相继的k-mer错位一个碱基。
如图1所示，在步骤104中，对k-mer序列进行分类。
由于测序错误，会引入很多测序深度只有1层(1×)的k-mer序列，如果全部k-mer序列放在一起统计k-mer频率的话，特别是比较大的基因组，需要非常大的内存，故本申请设计了分类方法以降低内存。例如，选取k-mer序列中连续n位碱基，该连续n位碱基可以是k-mer序列的前n位碱基，也可以是距离k-mer序列首部或尾部一定距离的连续n位碱基，并根据这几位碱基的组合将k-mer序列分成多个文件保存。比如，截取k-mer序列的前3位碱基序列，其中含n的k-mer序列去掉，由于每位碱基有a、t、c、g四种可能，故3位碱基序列有64种可能，可将k-mer序列归类到64份文件中。如此类推，如果截取首位碱基的话可分为4份文件，如果截取首2位碱基的话可分为16份文件，即4的n次方份文件，其中n表示截取的碱基的位数。
如图1所示，作为可选步骤，在步骤106中，压缩k-mer序列。
由于每条二代序列划分为多个k-mer序列会大大增加存储，故对k-mer序列进行压缩可以大大降低存储空间。如表1所示，每三个碱基可用1位ascii字符替代，故k-mer序列的存储可降到原来的三分之一。而在步骤102中，已经说明了如果要压缩，k-mer的大小选择3的整数倍，故刚好可压缩k-mer序列。
表1
如图1所示，在步骤108中，统计k-mer序列的频数。
分别对分类好的文件进行k-mer序列频数统计，如果不需要压缩就对步骤104的文件进行统计，需要压缩就对步骤106的文件进行统计。
如图1所示，在步骤110中，筛选并获得低频k-mer序列。
根据步骤108得到的频数文件，统计频数为1至m次的k-mer序列个数，并可从频数统计文件中筛选并获得低频k-mer序列。所谓“低频”，是指低于预定频数。由于k-mer序列符合泊松分布，可以认为低于一定覆盖深度的k-mer序列是由于测序错误造成的。本发明实施例中，预定频数是根据具体应用项目确定的，根据不同项目要求，可以将预定频数确定在不同的数值或者数值范围，例如1-100，优选2-50，更优选3-20，特别优选4-10。在本发明一个实施例中，用k-39得到的频数统计图，频数小于4的可以认为是低频kmer序列。
如图1所示，在步骤112中，用低频k-mer序列对二代序列进行纠错。
将低频k-mer序列比对回二代序列，若比对上，则从该k-mer头尾两端对应在二代序列上的位置将二代序列截断；保留截断后的长度大于预定长度的序列，舍弃截断后的长度小于预定长度的序列。预定长度可以是80-120个碱基长度(bp)，例如，在本发明一个实施例中，预定长度是100个碱基长度，如果截断后的序列长度大于100bp则保留，否则舍弃。
图3示出了本发明一个实施例中用低频k-mer序列对二代序列进行纠错的原理。将低频k-mer比对到二代序列上；若比对上，从低频k-mer头尾比对上的位置截断二代序列；获得截断后的二代序列hc1和c2t，如果hc1长度大于100bp则保留hc1序列，否则舍弃；如果c2t长度大于100bp则保留，否则舍弃。
相应于上述实施例的二代序列k-mer纠错方法，本发明实施例还提供一种二代序列k-mer纠错装置，如图4所示，包括：获取单元402，用于获取二代序列的k-mer序列；分类单元404，用于对k-mer序列进行分类，得到多个分类文件；统计单元408，用于对每个分类文件，统计k-mer序列的频数；筛选单元410，用于筛选并获得低于预定频数的k-mer序列；纠错单元412，用于用预定频数的k-mer序列对二代序列进行纠错；其中，纠错包括：将预定频数的k-mer序列比对回二代序列，若比对上，则从该k-mer头尾两端对应在二代序列上的位置将二代序列截断；保留截断后的长度大于预定长度的序列，舍弃截断后的长度小于预定长度的序列。
进一步地，本发明实施例的装置还包括：压缩单元406，用于压缩分类过的k-mer序列。
本领域技术人员可以理解，上述实施方式中各种方法的全部或部分步骤可以通过程序来指令相关硬件完成，该程序可以存储于计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。
下面提供某个虫子基因组大小约为1.65g具体应用例，来说明本发明实施例的二代序列k-mer纠错方法和其技术效果。在该实施例中，实现二代序列用k-mer进行纠错，具体步骤如下：
(一)获取二代序列的k-mer序列
以k＝39获取二代序列的k-mer序列。去掉含n的k-mer序列，得到k＝39的k-mer序列文件。
(二)对k-mer序列进行分类
截取k-mer序列的头2位，并把头2位相同的k-mer序列归类到同一个文件，共产生16份文件。
(三)压缩k-mer序列
对归类好的16份文件中的k-mer序列进行压缩，其中每3个碱基按照表1中的对应关系替换成1个ascii码进行，完成压缩。
(四)统计k-mer序列频数
用哈希列表单独对每份压缩文件k-mer序列进行次数统计。得到16份k-mer序列频数表，格式为两列，第一列为压缩的k-mer序列，第二列是对应的频数。再对16份k-mer序列频数表中每个频数出现的次数进行统计。得到频数统计表，格式为第一列频数深度，第二列为此频数深度的k-mer序列总数。画成如图5所示的频数分布图，由于k-mer序列频数符合泊松分布，从图5中可看出深度小于等于4的为低频k-mer序列。
(五)筛选并获得低频k-mer序列
从步骤(四)中的频数统计表中可得出k-mer频数小于4的属于低频k-mer序列。将频数小于4的k-mer序列获得到同一个文件中。
(六)用低频k-mer序列对二代序列进行纠错
将低频k-mer序列比对回二代序列，不允许错配，如果比对上，认为此二代序列可能有测序错误，将此序列截断，然后根据截断后的序列长度决定保留还是舍弃序列。具体地，如果截断后的序列长度大于100bp则保留，否则舍弃。
(七)未纠错和纠错后的二代序列组装内存对比
我们用soapdenovo软件(此软件可以从网上免费获得，网址为http://soap.genomics.org.cn/soapdenovo.html)对未纠错的和已纠错的二代序列进行分别组装，如表2所示，已纠错的组装内存峰值是未纠错的内存峰值的49.92％，效果十分明显。
表2
以上内容是结合具体的实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。