CN108573127B - 一种核酸第三代测序原始数据的处理方法及其应用 - Google Patents

Abstract

本申请公开了一种核酸第三代测序原始数据的处理方法及其应用。本申请的核酸第三代测序原始数据的处理方法，包括将第二代短序列数据比对到第三代自纠错数据上，统计比对结果中第三代自纠错数据的单碱基覆盖深度，将单碱基覆盖深度低于阈值的区域屏蔽为N，采用第二代测序的补洞软件对N屏蔽区域进行补洞，以获得单碱基错误率较低的核酸第三代测序数据。本申请的核酸第三代测序原始数据的处理方法，利用第二代短序列数据与第三代长序列数据进行比对，并利用第二代测序的补洞软件对比对结果中单碱基覆盖深度较低的N屏蔽区进行补齐，有效的降低了第三代测序数据中的单碱基错误率，提高了测序质量。

Description

一种核酸第三代测序原始数据的处理方法及其应用

技术领域

本申请涉及核酸测序数据处理领域，特别是涉及一种核酸第三代测序原始数据的处理方法及其应用。

背景技术

随着第二代测序技术(Next-generation sequencing，NGS)的成熟和普及，测序成本大大降低，其中二代测序仪Hiseq2500一次运行就可产出600Gb的数据量，相当于人类基因组的200倍。二代测序技术虽然可以快速产出大量数据，但是其有一个致命的缺点就是测序读长过短，众所周知，基因组组装最重要的指标就是N50以及基因组的完整度，但是传统二代测序技术由于其读长过短，其组装算法大多都是基于德布鲁因图论(de Bruijngraph)的思想，从而使得组装中遇到的最大的挑战就是解决高重复及高杂合基因组。

把组装出的重叠群(Contig)或骨架序列(Scaffold)从大到小排列，当其累计长度刚刚超过全部组装序列总长度50％时，最后一个重叠群(Contig)或骨架序列(Scaffold)的大小即为N50的大小，N50对评价基因测序的完整性有重要意义。把组装出的重叠群(Contig)或骨架序列(Scaffold)从大到小排列，当其累计长度刚刚超过全部组装序列总长度90％时，最后一个重叠群(Contig)或骨架序列(Scaffold)的大小即为N90的大小。

在一条骨架序列(Scaffold)中，重叠群(Contig)之间无序列信息且被表示为N的区域，被称为洞(gap)。

K-mer是指将一条长度为L的序列，按照长度K由L序列的5’到3’端挨个碱基进行分割，从而得到L-K+1个长度为K的核苷酸序列。例如，长度为90bp的短序列，按照17bp从头到尾挨个碱基进行分割，可以得到74个17bp的连续序列，即17-mer序列。

第三代测序技术也称为单分子实时测序技术(Single Molecule Real Time,SMRT)。PacBio第三代测序仪具有超长读长、无PCR扩增偏差的单分子测序、直接分析碱基修饰等技术优势，已经快速应用于基因组de novo组装、转录组学研究等领域，其平均读长10-15Kb，最长读长可超过40Kb。第三代测序技术的优势就在于单分子测序，对于高杂合、高重复、或者高GC，不存在任何偏好，所以三代测序技术可以显著的提升重叠群(Contig)组装水平。

但是，第三代数据最大的问题就在于其极高的单碱基错误率，高达15％。其错误类型主要是插入缺失，并且随机分布，可以通过生物信息的手段进行一定程度的修正，目前发布的三代组装流程中，如SMRT、Falcon、Pbcr、Canu，都具有第三代数据自纠错的功能，可以将20倍以上的三代数据的错误率从15％降低到3％左右，但是3％的错误率对于基于OLC(Overlap-Layout-Consensus)算法的组装软件干扰依然很大。

因此，亟需一种错误率更低的第三代测序数据的处理方法，以提高核酸第三代测序的质量。

发明内容

本申请的目的是提供一种新的核酸第三代测序原始数据的处理方法，及其应用。

为了实现上述目的，本申请采用了以下技术方案：

本申请公开了一种核酸第三代测序原始数据的处理方法，包括将第二代短序列数据比对到第三代自纠错数据上，统计比对结果中第三代自纠错数据的单碱基覆盖深度，将单碱基覆盖深度低于阈值的区域屏蔽为N，采用第二代测序的补洞软件对N屏蔽区域进行补洞，以获得单碱基错误率较低的核酸第三代测序数据。

需要说明的是，本申请的关键在于，利用精准度相对较高的第二代短序列数据与单碱基错误率较高的第三代长序列数据进行比对，将单碱基覆盖深度低于阈值的区域屏蔽为N，并将含有N的第三代长序列数据作为含洞的骨架序列，利用第二代测序的补洞软件进行补齐，以获得错误率更低的第三代测序数据，本申请的一种实现方式中，错误率仅为0.65％。

优选的，本申请的核酸第三代测序原始数据处理方法，具体包括以下步骤，

(a)根据接头序列信息、低质量碱基占整条序列的百分比、含N个数占整条序列的百分比以及复制序列情况，对相同核酸待测样品的第二代测序原始数据进行过滤处理，获得第二代短序列数据；可以理解，如果是采用事先存储的第二代短序列数据，则可以省略该步骤，例如本申请的实施例中直接调用网上下载的第二代短序列数据，则可以不用该步骤；

(b)根据接头序列、序列长度和序列质量值，对相同核酸待测样品的第三代序列原始数据进行过滤处理，获得第三代长序列数据；可以理解，如果是采用事先存储的第三代长序列数据，则可以省略该步骤，例如本申请的实施例中直接调用网上下载的第三代长序列数据，则可以不用该步骤；

(c)对步骤(b)获得的第三代长序列数据，通过多重序列比对的方法将20倍以上的第三代数据的错误率降低到3％左右，获得第三代自纠错数据；

其中，20倍以上的第三代数据，就是步骤(b)得到的第三代数据，多重序列比对的方法进行纠错，实际上就是这批第三代数据自我比对的结果，通过自我比对降低错误率。

(d)将步骤(a)获得的第二代短序列数据比对到步骤(c)的第三代自纠错数据上，统计第三代自纠错数据中单碱基覆盖深度以及长序列的覆盖率；

(e)将单碱基覆盖深度低于阈值的区域屏蔽为N，获得屏蔽后的第三代自纠错数据；

(f)将步骤(e)中获得的含有N的屏蔽后的第三代自纠错数据视为含洞的骨架序列，采用至少两种不同的第二代测序短序列补洞软件，分别对屏蔽后的第三代自纠错数据进行补洞，获得第三代二次自纠错数据，并提取获得所有补洞序列；

其中，补洞是基于第二代测序短序列数据对屏蔽后的第三代自纠错数据中N屏蔽区域进行补齐，此次采用的第二代测序短序列数据已经经过过滤及纠错处理，利用高精准的第二代测序短序列数据对第三代数据的屏蔽区域进行补齐，从而提高第三代数据的准确性。

(g)将所有补洞序列与其相应的被屏蔽处的原始序列进行一一比对，相似度低于相似性阈值的区域还原为原始序列，相似度高于或等于相似性阈值的区域替换为补洞序列，并且，将没有补洞或者补洞不完全的区域还原为原始序列，即获得最终的核酸第三代测序数据。

优选的，在将步骤(a)获得的第二代短序列数据比对到步骤(c)的第三代自纠错数据上之前，还包括步骤(a2)根据步骤(a)获得的第二代短序列数据构建K-mer频谱，通过K-mer频谱对所述第二代短序列数据进行纠错，获得高精准的第二代短序列数据，再将该高精准的第二代短序列数据用于步骤(d)，与步骤(c)的第三代自纠错数据比对。

需要说明的是，步骤(a2)的主要目的是将第二代短序列数据的精准度进一步提高，精准度越高，会进一步减小第二代短序列数据本身单碱基错误对本申请处理方法的干扰。

优选的，K-mer频谱的构建以及K-mer频谱对第二代短序列数据进行纠错采用的软件为Soapdenovo2软件包及其纠错模块SOAPec-2.0。

优选的，步骤(e)中，具体的将单碱基覆盖深度低于3的区域屏蔽为N。

需要说明的是，步骤(e)将单碱基覆盖深度低于阈值的区域屏蔽为N，其中阈值的选择直接影响屏蔽区域的大小，可以理解，如果阈值太大，屏蔽区域也会相应的增加，从而导致补洞纠错的准确性有所降低；而如果阈值太小，屏蔽区域相应减小，后续补洞纠错的量也相应减小，本申请的处理方法对第三代测序数据的错误率纠正效果也越差。综合考虑，本申请优选的方案中阈值为3，即将单碱基覆盖深度低于3的区域屏蔽为N，屏蔽区覆盖率约为32.54％至41.95％。

优选的，步骤(g)中，相似性阈值为95％。

优选的，在步骤(e)之前，还包括步骤(d2)将双末端覆盖深度为0的碱基切除，并删除覆盖率低于25％的长序列。

需要说明的是，步骤(d2)的目的主要是去除无法进行纠错的区域及序列，减少这部分数据的干扰，使得最终得到的第三代测序数据更加准确。可以理解，步骤(d2)可以根据所处理的核酸对象或测序准确度要求而选择。不过，需要注意的是，如果第三代自纠错数据未经过步骤(d2)处理，则两端覆盖深度为0的区域需保留一定长度不被屏蔽，例如保留约75bp不被屏蔽。当然，步骤(d2)处理在提高最终核酸第三代测序数据准确性的同时，也会造成数据损失。损失数据包括：第一，切除的双末端覆盖深度为0的碱基；第二，删除的覆盖率低于25％的长序列。对于所删除的长序列，本申请的优选方案选取覆盖率低于25％作为判断阈值，可以理解，如果阈值过低，达不到剔除数据提高准确性的效果，或者准确性提高效果相应减小；如果阈值过高，则删除的数据量增加，造成大量数据损失，同样影响测序质量。综合考虑，本申请优选的删除覆盖率低于25％的长序列，其数据损失约为16.99％至23.99％；而采用本申请的处理方法，最终得到的核酸第三代测序数据错误率仅为0.65％。

优选的，步骤(c)采用的自纠错软件为MHAP、Falcon、HGAP或Canu，优选为MHAP。

优选的，步骤(d)和(e)采用的软件为Soap2或Bwa。

需要说明的是，步骤(d)将第二代短序列数据比对到第三代自纠错数据上，通过比对结果统计第三代自纠错数据中单碱基覆盖深度以及长序列的覆盖率，相当于将第三代长序列视为基因组序列的重叠群(Contig)，该步骤可以使用短序列比对软件Soap2来实现，并利用软件soap.coverage2.27来统计单碱基覆盖深度，Soap2以及soap.coverage2.27均可以通过http://soap.genomics.org.cn/免费获得。该步骤也可以使用短序列比对软件bwamem结合bwa depth来实现。

优选的，步骤(f)中采用了两种补洞软件进行补洞，具体为KGF和Gapcloser。

优选的，步骤(g)中采用的比对软件为Blast+。

本申请的另一面公开了本申请的处理方法在核酸第三代测序平台或系统中的应用。

本申请的再一面公开了一种核酸测序平台或系统，该核酸测序平台或系统采用本申请的处理方法对测序获得的原始数据进行处理。

可以理解，本申请的第三代测序数据处理方法，能够有效的提高核酸测序质量，降低单碱基错误率，这对单碱基错误率较高的第三代测序来说，是很好的补充和纠正，大大提高了第三代测序的准确性；因此，完全可以将本申请的处理方法，或者按照本申请的处理方法步骤或思路，将其整合到核酸测序平台或系统中，进而提供一种高通量、高准确性的核酸测序平台、核酸测序系统或核酸测序仪。

由于采用以上技术方案，本申请的有益效果在于：

本申请的核酸第三代测序原始数据的处理方法，利用第二代短序列数据与第三代长序列数据进行比对，并利用第二代测序的补洞软件对比对结果中单碱基覆盖深度较低的N屏蔽区进行补齐，有效的降低了第三代测序数据中的单碱基错误率，提高了测序质量。

附图说明

图1是本申请实施例中核酸第三代测序原始数据处理方法的流程框图；

图2是本申请实施例核酸第三代测序原始数据处理方法中，获得第三代二次自纠错数据后的处理流程框图。

具体实施方式

本申请的核酸第三代测序原始数据处理方法，其关键在于利用精准度相对较高的第二代短序列数据与单碱基错误率较高的第三代长序列数据进行比对，并利用第二代测序的补洞软件对比对结果中单碱基覆盖深度较低的N屏蔽区进行补齐。其具体步骤，如图1所示，图1中，101为步骤(a)、102为步骤(a2)、103为步骤(b)、104为步骤(c)、105为步骤(d)、106为步骤(d2)、107为步骤(e)、108为步骤(f)、109为步骤(g)。如图1所示，本申请的处理方法可包含以下4种纠错方案：

方案一：包括101、102、103、104、105、106、107、108、109步骤；

方案二：包括101、103、104、105、106、107、108、109步骤；

方案三：包括101、102、103、104、105、107、108、109步骤；

方案四：包括101、103、104、105、107、108、109步骤；

有上述四种方案可知，步骤102以及步骤106为可选步骤，其中步骤102会损失部分第二代测序数据，但是可以使第二代序列更加准确，步骤106也会删除部分难以纠正的第三代数据，所以方案一纠错精度更高，而方案四能够保证数据量的完整性，如果后续用于组装且数据量足够的话，可以采用方案一。

下面参照图1和图2对本申请进行更全面的描述：

步骤101，使用原始的第二代测序数据，根据接头序列信息、低质量碱基占整条序列的百分比、含N个数占整条序列的百分比以及复制序列情况，对第二代测序原始数据进行过滤处理，从而获得高质量二代数据，即第二代短序列数据。该步骤采用常规的测序分析软件和平台即可完成。

步骤102，即步骤(a2)，将高质量二代数据进行纠错，此步骤可选，此步骤的主要目的是将二代数据的精准度进一步提高，二代数据的精准度越高，会进一步减小二代数据单碱基错误对本方法的干扰。具体的，根据第二代短序列数据构建K-mer频谱(K-merFrequency Spectrums,缩写KMF)，从而通过K-mer频谱对第二代短序列数据的短序列进行纠错，得到高精准二代数据，即高精准的第二代短序列数据。本申请的实施例中所用软件为Soapdenovo2软件包中纠错模块SOAPec-2.0，可从网站http://soap.genomics.org.cn/soapdenovo.html中免费获得。

步骤103，根据接头序列、序列长度和序列质量值，对原始三代序列数据进行过滤处理，最终获得质量值较高的三代长序列数据，即第三代长序列数据。同样的，可以采用常规的测序分析软件和平台完成。

步骤104，由于三代测序主要采用的是单分子实时测序技术，导致测序得到的三代数据存在高达15％的单碱基插入缺失错误，可以通过多重序列比对的方法将20倍以上的三代数据的错误率降低到3％左右，获得第三代自纠错数据。本申请的处理方法主要是针对自纠错后的第三代自纠错数据进行的二次纠错，从而得到更高准确性的第三代测序数据。目前，主流的三代数据组装软件都含有三代数据自纠错功能，如MHAP、Falcon、HGAP、Canu等，本申请实施例中所用软件为MHAP。

步骤105，将第二代短序列数据比对到第三代自纠错数据上，通过比对结果统计三代数据的单碱基覆盖深度，相当于将三代长序列视为基因组序列的重叠群(Contig)，该步骤可以使用短序列比对软件Soap2来实现，并利用软件soap.coverage2.27来统计单碱基覆盖深度，Soap2以及soap.coverage2.27均可以通过http://soap.genomics.org.cn/免费获得。该步骤也可以使用短序列比对软件bwamem结合bwa depth来实现。

步骤106，根据三代数据的单碱基覆盖情况，将两端覆盖度为0的碱基切除，并将覆盖度过低的长序列抛弃。增加此步骤的目的主要是去除本纠错方法无法进行纠错的区域及序列，减少这部分数据的干扰，可以使得最终得到的三代数据更加准确。

步骤107，根据三代数据的单碱基覆盖情况，将单碱基覆盖深度低于阈值的区域屏蔽为N，获得屏蔽后的第三代自纠错数据。如果三代数据未经过步骤106的处理，则两端覆盖度为0的区域保留一定长度不被屏蔽。

步骤108，采用第二代测序短序列补洞软件，对屏蔽后的第三代自纠错数据进行补洞。此步骤为本申请处理方法的核心步骤，该步骤主要是基于二代短序列的补洞算法对三代数据进行纠错，相当于将含有N序列的第三代自纠错数据视为基因组组装得到的含洞的骨架序列(Scaffold)。该步骤前后使用两种不同的二代短序列补洞软件，分别是KGF和Gapcloser，结合这两种补洞算法，可以将绝大多数的洞进行补齐，从而获得初级的第三代二次自纠错数据。下面详细介绍这两种不同的补洞算法：

KGF采用的版本是KGF1.19，KGF1.19可从Githup免费获得，其网址为，https://github.com/gigascience/paper-zhang2014/tree/master/Genome_assembly/GapCloser/krskgf。KGF补洞算法主要基于二代数据中短序列间的末端配对(paired-end)关系以及短序列与洞边缘的序列比对关系进行补洞，将二代短序列比对到包含洞的骨架序列(Scaffold)中，并将可以比对到洞相邻重叠群(Contig)的配对短序列收集起来，将收集起来的短序列根据重叠群(Overlap)关系进行组装，组装后序列再与洞两边相邻的重叠区(Contig)进行比对，根据比对结果，对洞进行填补。

Gapcloser的版本为Gapcloser_v1.12，为Soapdenovo2软件包中一个模块，可从SF.net(https://sourceforge.net/projects/soapdenovo2/files/GapCloser/)免费获得，Gapcloser补洞算法首先将二代短序列通过两个索引表储存到内存中，这两个索引表分别包含按照字典顺序排序的正向序列和反向互补序列信息，然后该算法进入迭代补洞阶段，每一迭代都会通过快速检索上述索引表提取出每个重叠群(Contig)比对上的末端配对短序列，沿着洞两端进行单碱基延伸，如果该位点可以得到延伸，说明提取出来的短序列超过80％支持该碱基，否则，该位点将被定义为分歧位点，从而导致本次迭代延伸停止。通过迭代处理，Gapcloser可以将较大洞进行补全。

步骤109，从第三代二次自纠错数据中提取补洞序列，将补洞序列与原位置屏蔽区域序列进行比对、还原、替换处理，得到终极的核酸第三代测序数据。经过步骤108，可以获得初级的第三代二次自纠错数据，该数据不成熟需要进行处理的主要原因有两方面，一方面第三代二次自纠错数据仍然存在部分N序列未被补全；另一方面，本申请的处理方法主要针对那些无法通过步骤104的三代自纠错方法修正的单碱基插入缺失，因为本申请的纠错过程是利用第二代数据进行补洞处理，难免会由于第二代数据结构以及补洞算法的原因，导致第三代数据发生结构上的改变。因此，需要对第三代二次自纠错数据进行处理，具体的，如图2所示，图2展示了初级的第三代二次自纠错数据的处理过程，该处理过程可以保证通过该方法纠正过来的第三代长序列没有发生结构上的变化。

步骤201，根据第三代长序列中原屏蔽区域的补洞情况，提取出第三代二次自纠错数据中相应屏蔽区域的补洞序列，该步骤可以追踪补洞软件的日志文件实现。

步骤202，将提取出的补洞序列，与对应区域的原屏蔽序列进行精确比对，本申请的实施例中采用的是精确比对软件blast-2.5.0+中blastn，比对结果以blast+的m6格式输出。

步骤203，根据比对结果，可以得知第三代数据屏蔽区域中补洞序列与原屏蔽区的原始序列的相似性，将相似性低于相似性阈值的屏蔽区还原为原始序列。

步骤204，根据比对结果，将相似度高于或等于相似性阈值的屏蔽区域替换为相应位置的补洞序列。

步骤205，经过步骤203和步骤204处理后，屏蔽后的第三代数据大部分屏蔽区域已经通过还原、替换处理得到有效序列信息，但是仍然存在一些屏蔽区域在补洞过程中无法补全，从而无法提取补洞序列与原屏蔽序列进行比对，这部分屏蔽区域直接还原为原屏蔽区域的原始序列，最终得到本申请的核酸第三代测序数据。

下面通过具体实施例和附图对本申请作进一步详细说明。以下实施例仅对本申请进行进一步说明，不应理解为对本申请的限制。

实施例

本例以有效参考序列大小为35Mb，实际大小为51Mb，的人22号染色体为分析对象，采用本申请的核酸第三代测序原始数据处理方法进行纠错处理。具体如下：

1数据处理

1)通过NCBI分别获得人样品的第二代数据和第三代数据进行测试，网址如下：ftp://ftp.ncbi.nih.gov/giab/ftp/data/AshkenazimTrio/HG002_NA24385_son/。

2)利用针对单分子实时测序长序列的全局比对软件blasr将第三代数据比对到人类基因组hg19上，提取第22号染色体上的第三代测序长序列，获得数据量3.14Gb，N50为11.46Kb。即本申请的处理方法中步骤(b)的第三代长序列数据。

需要说明的是，本例的第三代数据直接从第22号染色体数据上比对提取，因此，不需要步骤(b)的“根据接头序列、序列长度和序列质量值，对相同核酸待测样品的第三代序列原始数据进行过滤处理”。但是，常情况下，测序样品都是无参考序列的，不能够直接比对提取获得，所以需要步骤(b)。

3)利用短序列比对软件soap2将第二代数据比对到人类基因组hg19上，提取第22号染色体上的具有末端配对关系的成对短序列，最终获得5.58M对插入片段为350bp的末端配对序列，序列读长为250bp，总数据量为2.79Gb。

需要说明的是，本例的第二代数据是直接从第22号染色体数据上比对提取的，因此，不需要经过过滤纠错处理就可以保证二代数据的精准性，即不需要步骤(a)和步骤(a2)。可以理解，本例采用的是人22号染色体序列进行的试验，其本身有参考序列，因此无需过滤纠错处理；而通常情况下，测序样品都是无参考序列的，所以，还是需要经过过滤纠错处理，提高第二代数据准确性，即需要步骤(a)和步骤(a2)。

2Pacbio数据自纠错

将提取得到的人22号染色体的第三代长序列数据，通过MHAP软件进行自纠错，得到第三代自纠错数据，数据量为1.53Gb，N50为11.26Kb。即本申请的处理方法的步骤(c)。

3对Pacbio自纠错数据再次纠错

1)将自纠错的第三代自纠错数据按照20Mb大小进行分割，得到74份数据；需要说明的是，该步骤可以增加分割份数，即增加并行任务数，以加快程序运行；

2)利用短序列比对软件soap2将二代数据比对到每一份三代数据上并利用配套软件soap.coverage2.27统计三代数据单碱基覆盖深度，其中soap2参数设置为-m 0-x 1000-p 6-v 0-r 2，每份三代数据的覆盖率范围是55.24％至63.26％，例如，第一份三代数据大小L＝20916609bp，单碱基深度不为0的大小C＝12652235bp，数据覆盖率为C/L*100＝60.4889％；即本申请的处理方法的步骤(d)；

3)根据三代数据单碱基覆盖深度，将双末端覆盖深度为0的碱基切除，并删除覆盖率低于25％的三代长序列，每份三代数据的损失率范围是16.99％至23.99％，例如，第一份三代数据原大小L＝20916609bp，处理后数据大小T＝16389626bp，数据损失率为(1-T/L)*100＝21.643％；即本申请的处理方法的步骤(d2)；

4)处理后的三代长序列，根据覆盖深度，将覆盖深度低于3的碱基屏蔽为N，同时为了提高补洞效率，屏蔽为N的连续区域不得短于100bp，双末端要保留75bp的长度不被屏蔽，每份数据的屏蔽率范围是32.54％至41.95％，例如，第一份处理后数据大小T＝16389626bp，屏蔽区域大小M＝5424418bp，数据屏蔽率为M/T*100＝33.0967％；即本申请的处理方法的步骤(e)；

5)先后使用补洞软件KGF和Gapcloser对屏蔽后的每份三代数据进行补洞，获得初级三代二次纠错数据，即第三代二次自纠错数据；即本申请的处理方法的步骤(f)；

6)使用软件blast-2.5.0+中blastn将提取得到的初级三代二次纠错数据的补洞序列与对应区域的原屏蔽序列进行精确比对，以比对结果中相似性95％为阈值，将未比对上、相似性低于阈值以及未补全的屏蔽区域还原为原序列，相似度高于或等于阈值的屏蔽区域替换为补洞序列，每份三代数据的屏蔽区域的补洞序列替换率范围是79.96％至83.25％，替换率从一个侧面展示了本纠错方法的纠错效率；即本申请的处理方法的步骤(g)；

例如，第一份数据屏蔽区域大小M＝5424418bp，其中被替换为补洞序列的大小为F＝4445515bp，数据替换率F/M*100＝81.9537％。

8)将再次纠错后的74份三代数据进行合并，得到最终的核酸第三代测序数据，数据量为1.21Gb，N50为11.07Kb。可以理解，由于数据量比较大，为了提高运行速度和处理效率，本例将三代数据分割成了74份，74份数据分别按照本申请的第三代测序数据处理方法进行处理，而后合并得到最终的核酸第三代测序数据。

4Pacbio数据评估

1)错误率评估

利用bwamem算法分别将第三代长序列数据、第三代自纠错数据、本例处理后的最终的核酸第三代测序数据，比对到人参考基因组hg19第22号染色体序列上，计算比对结果中插入、缺失及错配，从而统计数据错误率，其中bwamem的-x参数设置为pacbio，结果如表1所示。

表1三代数据错误率统计表

标题	原数据	自纠错数据	AMF纠错数据
				总序列条数	353,418	76505	65995
总碱基数	3,137,854,583	1,000,078,976	829,495,184
				序列比对条数	352,955	76,505	65,995
碱基比对总数	2,288,816,181	911,760,814	818,977,833
				插入碱基数/％	173,856,628/7.60％	4,214,820/0.46％	965,558/0.12％
缺失碱基数/％	89,191,342/3.90％	8,823,386/0.97％	2,480,745/0.30％
				错配碱基数/％	95,549,815/4.17％	11,696,948/1.28％	1,864,127/0.23％
序列条数比对率	99.87％	100.00％	100.00％
				碱基个数比对率	72.94％	91.17％	98.73％
单碱基错误率	15.67％	2.71％	0.65％

表1中，原数据是指“1数据处理”中“2)”小节得到的第三代长序列数据，自纠错数据是指“2Pacbio数据自纠错”得到的第三代自纠错数据，AMF纠错数据是指本例处理后的最终的核酸第三代测序数据。需要说明的是，AMF是Align-Mask-Fillgap的缩写，由于本申请的处理方法核心思想就在于比对、屏蔽和补洞，因此，本申请的处理方法缩写为AMF纠错处理，所得到的数据缩写为AMF纠错数据。

表1的结果显示，原三代数据的单碱基错误率为15.67％，三代自纠错数据的单碱基错误率为2.71％，而本例的AMF纠错数据的单碱基错误率仅为0.65％，AMF纠错处理后三代数据的错误率得到显著的修正。

2)组装测试

使用基于OLC(Overlap-Layout-Consensus)算法的基因组组装软件CeleraAssembler(Eugene W.Myers,2000)对自纠错后的三代数据，即第三代自纠错数据，和AMF纠错处理后三代数据，即本例获得的最终的核酸第三代测序数据，分别进行组装测试。按照两组数据中序列的长度进行排序，并选取20倍，约700M，较长的序列作为CA的输入数据，在保证数据量以及参数相同的条件下分别进行组装，组装结果如表2所示。

表2三代数据组装结果统计表

表2的结果显示，经过AMF纠错处理后的三代数据要比仅仅自纠错后的三代数据噪音区域更少，对基于OLC算法的组装软件干扰更小，更容易组装得到组装指标较佳的基因组。

3)组装结果soap比对

利用短序列比对软件soap2将“1数据处理”第“3)”小节得到的第二代数据比对到上述两种组装结果上，通过统计二代短序列的比对率以及组装结果的覆盖率来评估其组装准确性，结果如表3所示。

表3组装结果soap比对统计表

标题	自纠错数据	AMF纠错数据
			二代数据总序列条数	11,152,260	11,152,260
序列比对条数/％	5421334/48.61％	5,893,582/52.85％
			数据总碱基数	36,418,229	34,077,730
比对覆盖碱基数	32,054,978	33,181,971
			比对覆盖率	88.019％	97.3714％

表3的结果显示，经过AMF纠错处理后的三代数据组装结果，与第二代数据比对的覆盖率97.3714％远高于仅仅自纠错后的三代数据组装结果的覆盖率88.019％，可见，AMF纠错处理后的三代数据组装结果准确率远高于仅仅自纠错后的三代数据组装结果。

以上内容是结合具体的实施方式对本申请所作的进一步详细说明，不能认定本申请的具体实施只局限于这些说明。对于本申请所属技术领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本申请的保护范围。

Claims (9)

1.一种核酸第三代测序原始数据的处理方法，其特征在于：包括将第二代短序列数据比对到第三代自纠错数据上，统计比对结果中第三代自纠错数据的单碱基覆盖深度，将单碱基覆盖深度低于阈值的区域屏蔽为N，采用第二代测序的补洞软件对N屏蔽区域进行补洞，以获得单碱基错误率较低的核酸第三代测序数据。

2.根据权利要求1所述的处理方法，其特征在于：具体包括以下步骤，

（a）根据接头序列信息、低质量碱基占整条序列的百分比、含N个数占整条序列的百分比以及复制序列情况，对相同核酸待测样品的第二代测序原始数据进行过滤处理，获得第二代短序列数据；

（b）根据接头序列、序列长度和序列质量值，对相同核酸待测样品的第三代序列原始数据进行过滤处理，获得第三代长序列数据；

（c）对步骤（b）获得的第三代长序列数据，通过多重序列比对的方法将20倍以上的第三代数据的错误率降低到3%左右，获得第三代自纠错数据；

（d）将步骤（a）获得的第二代短序列数据比对到步骤（c）的第三代自纠错数据上，统计第三代自纠错数据中单碱基覆盖深度以及长序列的覆盖率；

（e）将单碱基覆盖深度低于阈值的区域屏蔽为N，获得屏蔽后的第三代自纠错数据；

（f）将步骤（e）中获得的含有N的屏蔽后的第三代自纠错数据视为含洞的骨架序列，采用至少两种不同的第二代测序短序列补洞软件，分别对屏蔽后的第三代自纠错数据进行补洞，获得第三代二次自纠错数据，并提取获得所有补洞序列；

（g）将所有补洞序列与其相应的被屏蔽处的原始序列进行一一比对，相似度低于相似性阈值的区域还原为原始序列，相似度高于或等于相似性阈值的区域替换为补洞序列，并且，将没有补洞或者补洞不完全的区域还原为原始序列，即获得最终的核酸第三代测序数据。

3.根据权利要求2所述的处理方法，其特征在于：在将步骤（a）获得的第二代短序列数据比对到步骤（c）的第三代自纠错数据上之前，还包括步骤（a2）根据步骤（a）获得的第二代短序列数据构建K-mer频谱，通过K-mer频谱对所述第二代短序列数据进行纠错，获得高精准的第二代短序列数据，再将该高精准的第二代短序列数据用于步骤（d），与步骤（c）的第三代自纠错数据比对。

4.根据权利要求3所述的处理方法，其特征在于：所述K-mer频谱的构建以及K-mer频谱对所述第二代短序列数据进行纠错采用的软件为Soapdenovo2软件包及其纠错模块SOAPec-2.0。

5.根据权利要求2所述的处理方法，其特征在于：所述步骤（e）中，具体的将单碱基覆盖深度低于3的区域屏蔽为N；所述步骤（g）中，相似性阈值为95%。

6.根据权利要求2-5任一项所述的处理方法，其特征在于：在步骤（e）之前，还包括步骤（d2）将双末端覆盖深度为0的碱基切除，并删除覆盖率低于25%的长序列。

7.根据权利要求6所述的处理方法，其特征在于：所述步骤（c）采用的自纠错软件为MHAP、Falcon、HGAP或Canu。

8.根据权利要求6所述的处理方法，其特征在于：所述步骤（d）和（e）采用的软件为Soap2或Bwa；所述步骤（f）中采用了两种补洞软件进行补洞，具体为KGF和Gapcloser；所述步骤（g）中采用的比对软件为Blast+。

9.一种核酸测序平台或系统，其特征在于：所述核酸测序平台或系统采用权利要求1-8任一项所述的处理方法对测序获得的原始数据进行处理。

Images