三代测序数据纠错的方法、装置和计算机可读存储介质与流程

三代测序数据纠错的方法、装置和计算机可读存储介质与流程

文档序号：15616049发布日期：2018-10-09 21:24

导航： X技术> 最新专利>计算;推算;计数设备的制造及其应用技术

本发明涉及生物信息技术领域，具体涉及三代测序数据纠错的方法、装置和计算机可读存储介质。

背景技术：

以Pacbio为代表的第三代测序平台，其测序读长(reads)长(平均10～15k)且无GC偏向性的优势，使其在基因组组装等方面得到了广泛的应用，但其过高的错误率(15％～20％)，使得组装算法的复杂性大大提高。对于组装策略而言，输入的读长的错误率、长度和测序深度是影响组装效果的主要指标。因此，在利用第三代测序数据进行组装时，充分利用数据特征，对第三代测序读长进行纠错，降低其数据错误率，是数据前处理的一个重要步骤。

以Hiseq为代表的第二代测序，相对于第三代测序来说是一个成本更低、准确率更高的测序方式。Hiseq数据的错误率要比PacBio数据的错误率低1～1.5个数量级。因而，利用第二代测序数据对第三代测序数据进行纠错，在可以将第三代数据的错误率降低到与第二代测序相当的水平。从而更有利用基因组组装及其他相关的技术应用。

目前的纠错技术包括：三代测序数据自纠错技术和二代测序数据对三代测序数据纠错的技术。

三代测序数据自纠错技术，以Quiver和PBcR MHAP为代表。以Quiver为例，其技术路线为：将长的三代PacBio读长作为参考序列(reference)，将其他读长比对到参考序列上。然后利用序列间的一致性去推断比对区域的一致性序列，用得到的一致性序列替换原序列从而得到纠错后的读长(Chen-Shan Chin,et al(2013),Nonhybrid,finished microbial genome assemblies from long-read SMRT sequencing data-nature methods,Supplementary Note 1,p13～p16)。PBcR MHAP与Quiver类似，其优势在于只用到了三代测序数据，缺点在于需要较高的数据深度。

二代测序数据对三代测序数据纠错的技术，以PacbioToCA和ECtools为代表。以PacbioToCA为例，其技术路线为：将二代短读长比对到三代读长上，然后将比对到一起的二代和三代读长合并，生成一致序列。然后截断并分离二代比对出现间隙(gap)的位点，作为纠错后的读长(Sergey Koren et al(2012)Hybrid error correction and de novo assembly of single-molecule sequencing reads-nature biotechnology)。该方案综合利用了二代测序数据和三代测序数据，但没有考虑同一基因组多个相似片段间的差异。

总而言之，现有技术方案存在如下缺陷：在典型应用场景下通常会导致大量的数据损失；会导致读长长度的缩短，不利于充分利用三代数据的读长优势；纠错结果为纯序列格式，无质量值系统，无法评估纠错结果中各碱基的错误率；并且自纠错技术需要一定深度的三代数据才能完成纠错，对低深度数据不适用。

技术实现要素：

本发明提供一种三代测序数据纠错的方法、装置和计算机可读存储介质。

根据第一方面，一种实施例中提供一种三代测序数据纠错的方法，包括：

利用二代测序数据和/或三代测序数据，组装出一个初步的参考基因组；

将上述二代测序数据和上述三代测序数据比对到上述参考基因组上；

对于上述三代测序数据比对结果中每个比对片段上的每个碱基位置，推断并赋予该碱基位置一个最大可能性的碱基型和质量值；

对于上述三代测序数据的读长中有多个比对片段和/或未比对上的片段，根据上述最大可能性的碱基型和质量值，将上述多个比对片段和/或未比对上的片段整合为一条读长。

进一步地，上述推断并赋予该碱基位置一个最大可能性的碱基型和质量值，通过最大后验模型、最大似然模型或隐马尔可夫模型实现；

优选地，上述最大后验模型包括：

对于每个碱基位置，在给定该碱基位置的基因组拷贝数、同一基因组位置上其它二代和三代比对片段的基因型以及二代和三代测序错误的先验概率的条件下，计算该碱基位置上各种可能的碱基型出现的后验概率；

将该碱基位置的碱基型推断为拥有最大后验概率的碱基型；将该碱基位置的质量值赋值为各种可能的碱基型的后验概率中的最大值除以第二大值然后取常用对数并乘以-10所得的结果；

更优选地，对于碱基位置L，假定其基因组拷贝数为n，在该位置上二代测序的基因型为R＝(r1，r2，r3，…，rk)，三代测序的基因型为S＝(s1，s2，s3，…，sl)，某三代测序读长的测序结果为obs,obs∈S；该位置所有可能的等位基因型为所有可能的基因型为则上述三代测序读长在该位置的碱基型为：

其中S′＝S-obs。

进一步地，对于上述三代测序数据的读长中有多个比对片段和/或未比对上的片段，上述最大可能性的碱基型和质量值，通过如下步骤确定：

对于未比对上的片段，其每个碱基位置的碱基型赋值为该位置观察的碱基型，质量值赋值为其测序的质量值；

对于同一碱基位置被多个比对片段覆盖，该位置的碱基型赋值为拥有最高质量值的碱基型，其质量值赋值为各碱基型质量值中的最大值减去次大值。

进一步地，在上述组装出一个初步的参考基因组之后且在上述将上述二代测序数据和上述三代测序数据比对到上述参考基因组上之前，还包括：处理上述参考基因组的组装结果，使其片段长度和基因组复杂度适于上述二代测序数据和上述三代测序数据比对。

根据第二方面，一种实施例中提供一种三代测序数据纠错的装置，包括：

组装装置，用于利用二代测序数据和/或三代测序数据，组装出一个初步的参考基因组；

比对装置，用于将上述二代测序数据和上述三代测序数据比对到上述参考基因组上；

推断装置，用于对于上述三代测序数据比对结果中每个比对片段上的每个碱基位置，推断并赋予该碱基位置一个最大可能性的碱基型和质量值；

整合装置，用于对于上述三代测序数据的读长中有多个比对片段和/或未比对上的片段，根据上述最大可能性的碱基型和质量值，将上述多个比对片段和/或未比对上的片段整合为一条读长。

进一步地，上述推断并赋予该碱基位置一个最大可能性的碱基型和质量值，通过最大后验模型、最大似然模型或隐马尔可夫模型实现；

优选地，上述最大后验模型包括：

对于每个碱基位置，在给定该碱基位置的基因组拷贝数、同一基因组位置上其它二代和三代比对片段的基因型以及二代和三代测序错误的先验概率的条件下，计算该碱基位置上各种可能的碱基型出现的后验概率；

将该碱基位置的碱基型推断为拥有最大后验概率的碱基型；将该碱基位置的质量值赋值为各种可能的碱基型的后验概率中的最大值除以第二大值然后取常用对数并乘以-10所得的结果；

更优选地，对于碱基位置L，假定其基因组拷贝数为n，在该位置上二代测序的基因型为R＝(r1,r2,r3,…,rk)，三代测序的基因型为S＝(s1,s2,s3,…,sl)，某三代测序读长的测序结果为obs,obs∈S；该位置所有可能的等位基因型为所有可能的基因型为则上述三代测序读长在该位置的碱基型为：

其中S′＝S-obs。

进一步地，对于上述三代测序数据的读长中有多个比对片段和/或未比对上的片段，上述最大可能性的碱基型和质量值，通过如下步骤确定：

对于未比对上的片段，其每个碱基位置的碱基型赋值为该位置观察的碱基型，质量值赋值为其测序的质量值；

对于同一碱基位置被多个比对片段覆盖，该位置的碱基型赋值为拥有最高质量值的碱基型，其质量值赋值为各碱基型质量值中的最大值减去次大值。

进一步地，在上述组装装置和上述比对装置之间，还包括：

处理装置，用于处理上述参考基因组的组装结果，使其片段长度和基因组复杂度适于上述二代测序数据和上述三代测序数据比对。

根据第三方面，一种实施例中提供一种三代测序数据纠错的装置，该装置包括：

存储器，用于存储程序；

处理器，用于通过执行上述存储器存储的程序以实现如第一方面的方法。

根据第四方面，一种实施例中提供一种计算机可读存储介质，包括程序，上述程序能够被处理器执行实现如第一方面的方法。

本发明对三代测序数据深度没有限制，可以在任意数据深度下完成数据纠错，从而实现对低平均深度(小于20×)的三代测序数据的纠错；除数据本身的插入缺失纠正引起的读长长度变化外，不引入额外的数据损失和读长长度损失，更多的数据量和更长的读长有利于下游的数据分析；引入纠错结果的质量值体系，使得纠错结果的单碱基质量可以评价，现有技术方案未能实现此功能。

附图说明

图1为本发明一个实施例中提供的三代测序数据纠错的方法流程图；

图2为本发明一个实施例中三代测序数据的读长中有多个比对片段和/或未比对上的片段的示意图；

图3为本发明一个实施例中提供的三代测序数据纠错的装置的结构框图。

具体实施方式

下面通过具体实施方式结合附图对本发明作进一步详细说明。在以下的实施方式中，很多细节描述是为了使得本发明能被更好的理解。然而，本领域技术人员可以毫不费力的认识到，其中部分特征在不同情况下是可以省略的，或者可以由其他元件、材料、方法所替代。在某些情况下，本发明相关的一些操作并没有在说明书中显示或者描述，这是为了避免本发明的核心部分被过多的描述所淹没，而对于本领域技术人员而言，详细描述这些相关操作并不是必要的，他们根据说明书中的描述以及本领域的一般技术知识即可完整了解相关操作。

另外，说明书中所描述的特点、操作或者特征可以以任意适当的方式结合形成各种实施方式。同时，方法描述中的各步骤或者动作也可以按照本领域技术人员所能显而易见的方式进行顺序调换或调整。因此，说明书和附图中的各种顺序只是为了清楚描述某一个实施例，并不意味着是必须的顺序，除非另有说明其中某个顺序是必须遵循的。

本发明实施例涉及三代测序数据纠错的方法和装置，使用二代测序数据和三代测序数据，其中在测序平台选择上，二代测序数据可以来源于Hiseq以及其他平台，例如454、BGISEQ；三代测序数据可以来源于Pacbio以及其他平台，例如Nanopore。因此，在本发明实施例中，对数据来源对测序平台的选择没有限制。

如图1所示，根据本发明的一个实施例，提供一种三代测序数据纠错的方法，包括：

步骤110：利用二代测序数据和/或三代测序数据，组装出一个初步的参考基因组。

该步骤中，参考基因组的组装可以使用单纯的二代测序数据，或者单纯的三代测序数据，也可以使用二代测序数据和三代测序数据进行混合组装。这些数据的组装方法，可以按照本领域通常用的技术进行。例如，使用DBG2OLC软件(参考Chengxi Ye等人,DBG2OLC:Efficient Assembly of Large Genomes Using Long Erroneous Reads of the Third Generation Sequencing Technologies,Scientific Reports 6:31900(2016)；来源http://www.nature.com/articles/srep31900)进行二代测序数据和三代测序数据混合组装。

二代测序数据可以来源于Hiseq或BGISEQ等测序平台；三代测序数据可以来源于Pacbio或Nanopore等测序平台。在本发明一个实施例中，二代测序数据来源于Hiseq测序平台，三代测序数据来源于Pacbio测序平台。

测序数据的深度没有特别的限制，尤其是对于三代测序数据而言，可以是小于20×的低平均深度数据，例如，在本发明一个实施例中，使用68×Hiseq数据作为二代测序数据，使用15×Pacbio数据作为三代测序数据。

经过该步骤的组装以后，组装结果体现为多个重叠群，一般情况下，不作任何处理即适合用于作为比对的参考基因组。例如，在本发明一个实施例中，组装结果数据统计显示组装重叠群(contig)N50达到289kb，自比对显示无长片段高重复序列，因此可以不作任何处理即用作比对的参考基因组。然而，在某些情况下，例如，基因组特别复杂和/或组装效果特别差的情况下，需要对组装结果进行优化。具体而言，基因组特别复杂(例如重复大于70％)时，需要对组装结果做序列层面的去冗余；组装效果特别差(N50小于1k)时，通常需要进行优化组装参数的处理。这些处理方法，可以参考现有去冗余或优化组装参数的技术。

步骤120：将二代测序数据和三代测序数据比对到上述参考基因组上。

通过将二代测序数据和三代测序数据分别比对到参考基因组上，能够间接实现二代测序数据和三代测序数据的比对，便于进行后续步骤的执行。

该步骤的比对方法，可以按照现有技术实现，例如，在本发明一个实施例中，利用bwa软件分别将二代测序数据(Hiseq数据)和三代测序数据(Pacbio数据)比对到参考基因组上。

步骤130：对于三代测序数据比对结果中每个比对片段上的每个碱基位置，推断并赋予该碱基位置一个最大可能性的碱基型和质量值。

该步骤可以通过几种模型实现，例如最大后验模型、最大似然模型或隐马尔可夫模型等。

对于最大后验模型而言，具体可以包括：对于每个碱基位置，在给定该碱基位置的基因组拷贝数、同一基因组位置上其它二代和三代比对片段的基因型以及二代和三代测序错误的先验概率的条件下，计算该碱基位置上各种可能的碱基型出现的后验概率；然后将该碱基位置的碱基型推断为拥有最大后验概率的碱基型，也就是最大可能性的碱基型；将该碱基位置的质量值赋值为各种可能的碱基型的后验概率中的最大值除以第二大值然后取常用对数并乘以-10所得的结果，也就是最大可能性的质量值。

在一个更优选的实施例中，对于碱基位置L，假定其基因组拷贝数为n，在该位置上二代测序的基因型为R＝(r1,r2,r3,…,rk)，其中，r1,r2,r3,…,rk中的每个表示一条二代测序读长在这个位置上的基因型，三代测序的基因型为S＝(s1,s2,s3,…,sl)，其中，s1,s2,s3,…,sl中的每个表示一条三代测序读长在这个位置上的基因型，某三代测序读长的测序结果为obs,obs∈S；该位置所有可能的等位基因型为所有可能的基因型为则上述三代测序读长在该位置的碱基型为：

其中S′＝S-obs。

如图2所示，某三代测序数据的读长中可能有多个比对片段(R1、R2和R3)，分别比对到在参考基因组的物理位置上彼此不相邻的多个重叠群(C1、C2和C3)上，并且可能有未比对上的片段(R4)，没有比对到参考基因组的任何重叠群上。在这种情况下，最大可能性的碱基型和质量值，可以通过如下步骤确定：

(1)对于未比对上的片段(R4)，其每个碱基位置的碱基型赋值为该位置观察的碱基型(即测序结果)，质量值赋值为其测序的质量值。

(2)对于同一碱基位置被多个比对片段覆盖，如图2所示，比对片段(R1)的X至Y之间的序列与比对片段(R2)的P至Q之间的序列是重叠的，也就是说，这条三代测序读长有一部分同时比对到重叠群(C1)和重叠群(C2)，那么X至Y(即P至Q)之间的多个碱基被两个比对片段(R1和R2)覆盖。在这种情况下，被多个比对片段覆盖的同一碱基位置的碱基型赋值为拥有最高质量值的碱基型，其质量值赋值为各碱基型质量值中的最大值减去次大值。

步骤140：对于三代测序数据的读长中有多个比对片段和/或未比对上的片段，根据最大可能性的碱基型和质量值，将多个比对片段和/或未比对上的片段整合为一条读长。

如图2所示，多个比对片段(R1、R2、R3和R4)经过步骤130，各碱基位置将具有唯一的纠错后碱基型和质量值，串联各碱基位置即可得纠错后的读长(图2中R)，即将多个比对片段和/或未比对上的片段整合为一条读长。

本领域技术人员可以理解，上述实施方式中各种方法的全部或部分功能可以通过硬件的方式实现，也可以通过计算机程序的方式实现。当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘、光盘、硬盘等，通过计算机执行该程序以实现上述功能。例如，将程序存储在设备的存储器中，当通过处理器执行存储器中程序，即可实现上述全部或部分功能。另外，当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序也可以存储在服务器、另一计算机、磁盘、光盘、闪存盘或移动硬盘等存储介质中，通过下载或复制保存到本地设备的存储器中，或对本地设备的系统进行版本更新，当通过处理器执行存储器中的程序时，即可实现上述实施方式中全部或部分功能。

如图3所示，根据本发明的一个实施例，还提供一种三代测序数据纠错的装置，包括：

组装装置310，用于利用二代测序数据和/或三代测序数据，组装出一个初步的参考基因组；比对装置320，用于将二代测序数据和三代测序数据比对到上述参考基因组上；推断装置330，用于对于三代测序数据比对结果中每个比对片段上的每个碱基位置，推断并赋予该碱基位置一个最大可能性的碱基型和质量值；整合装置340，用于对于三代测序数据的读长中有多个比对片段和/或未比对上的片段，根据最大可能性的碱基型和质量值，将多个比对片段和/或未比对上的片段整合为一条读长。

与本发明实施例的方法相对应，本发明实施例的装置中，推断并赋予该碱基位置一个最大可能性的碱基型和质量值，可以通过最大后验模型、最大似然模型或隐马尔可夫模型实现。

其中，在本发明的一个实施例中，最大后验模型包括：对于每个碱基位置，在给定该碱基位置的基因组拷贝数、同一基因组位置上其它二代和三代比对片段的基因型以及二代和三代测序错误的先验概率的条件下，计算该碱基位置上各种可能的碱基型出现的后验概率；将该碱基位置的碱基型推断为拥有最大后验概率的碱基型；将该碱基位置的质量值赋值为各种可能的碱基型的后验概率中的最大值除以第二大值然后取常用对数并乘以-10所得的结果。

在一个更优选的实施例中，对于碱基位置L，假定其基因组拷贝数为n，在该位置上二代测序的基因型为R＝(r1,r2,r3,…,rk)，三代测序的基因型为S＝(s1,s2,s3,…,sl)，某三代测序读长的测序结果为obs,obs∈S；该位置所有可能的等位基因型为所有可能的基因型为则上述三代测序读长在该位置的碱基型为：

其中S′＝S-obs。

在一个更优选的实施例中，对于三代测序数据的读长中有多个比对片段和/或未比对上的片段，最大可能性的碱基型和质量值，通过如下步骤确定：对于未比对上的片段，其每个碱基位置的碱基型赋值为该位置观察的碱基型，质量值赋值为其测序的质量值；对于同一碱基位置被多个比对片段覆盖，该位置的碱基型赋值为拥有最高质量值的碱基型，其质量值赋值为各碱基型质量值中的最大值减去次大值。

在一个更优选的实施例中，在组装装置310和比对装置320之间，还包括：处理装置，用于处理参考基因组的组装结果，使其片段长度和基因组复杂度适于二代测序数据和三代测序数据比对。

在本发明一个实施例中，还提供一种三代测序数据纠错的装置，该装置包括：存储器，用于存储程序；处理器，用于通过执行存储器存储的程序以实现如本发明实施例的方法。

在本发明一个实施例中，还提供一种计算机可读存储介质，包括程序，该程序能够被处理器执行实现如本发明实施例的方法。

以下通过实施例详细说明本发明的技术方案和效果，应当理解，实施例仅是示例性的，不能理解为对本发明保护范围的限制。

实施例1：人22号染色体PacBio测序数据混合纠错

本实施例的原始数据为68×Hiseq PE250测序数据和90×Pacbio测序数据。在本实施例中，我们抽取了15×Pacbio测序数据，作为示例数据来测试本发明在低深度数据下的表现。

1、利用68×Hiseq PE250测序数据和15×Pacbio测序数据，用DBG2OLC软件做混合组装，得到一个初步的参考基因组。

2、数据统计结果显示，组装重叠群(contig)N50为289kb，自比对显示无长片段高重复序列，组装结果不作任何处理即适合用于作为比对的参考序列。

3、利用bwa软件分别将Hiseq测序数据和Pacbio测序数据比对到参考序列上。

4、利用最大后验模型推断每个碱基位置的碱基型，并评估其质量值。

对于碱基位置L，假定其基因组拷贝数为n，在该位置上二代测序的基因型为R＝(r1,r2,r3,…,rk)，其中，r1,r2,r3,…,rk中的每个表示一条二代测序读长在这个位置上的基因型，三代测序的基因型为S＝(s1,s2,s3,…,sl)，其中，s1,s2,s3,…,sl中的每个表示一条三代测序读长在这个位置上的基因型，某三代测序读长的测序结果为obs,obs∈S；该位置所有可能的等位基因型为所有可能的基因型为则上述三代测序读长在该位置的碱基型为：

其中S′＝S-obs。

通过贝叶斯(Bayesian)公式展开及合并化简，可得：

 

式中，P(obs|g)为Pacbio测序的先验错误率，P(g|G)服从几何分布，P(G)服从多项分布(multinomial distribution)，P(R)P(S)为常数，L(G|R)和L(G|S)可通过基因分型算法或软件求得。

实际应用中，可以利用GATK软件分别对Hiseq测序数据和Pacbio测序数据识别基因型(call genotype)，然后利用各基因型的PL值计算得到L(G|R)和L(G|S)。然后利用公式计算得到各等位基因型的后验概率。该碱基位置的碱基型推断为拥有最大后验概率的等位基因的碱基型。该碱基位置的质量值赋值为各等位基因型的后验概率中的最大值除以第二大值然后取常用对数并乘以-10。

5、根据比对结果，查找确定每条三代测序读长是否有多个比对片段和/或未比对上的片段。对于未比对上的片段，其每个碱基位置的碱基型赋值为该位点观察的碱基型，质量值赋值为其测序的质量值。若同一个碱基位置被多个比对片段覆盖，则该位置碱基型赋值为拥有最高质量值的碱基型，其质量值赋值为各碱基型质量值中的最大值减去次大值。通过该步骤，每条三代测序读长的各个碱基位置将有唯一的纠错后碱基型和质量值。最后串联各碱基位置即可得纠错后的读长。

6、将纠错结果比对到人参考基因组hg19的22号染色体上。在忽略样品基因组与参考基因组差异的前提下，用比对错误率来估计纠错后的错误率，以评价纠错效果。

比较例1

1、利用PBcR MHAP对30×Pacbio测序数据做自纠错。

2、将纠错结果比对到人参考基因组hg19的22号染色体上。在忽略样品基因组与参考基因组差异的前提下，用比对错误率来估计纠错后的错误率，以评价纠错效果。

实施例1和比较例1的结果如表1所示。结果显示，实施例1在Pacbio测序数据深度低一半的情况下，数据的留存率、纠正后错误率以及读长N50相比比较例1的自纠错，有明显提升。

表1

 

 

以上应用了具体个例对本发明进行阐述，只是用于帮助理解本发明，并不用以限制本发明。对于本发明所属技术领域的技术人员，依据本发明的思想，还可以做出若干简单推演、变形或替换。

完整全部详细技术资料下载