一种PacBio测序数据组装得到的基因组序列的纠错方法

技术领域

本发明涉及生物信息技术领域,更具体的说，它涉及一种PacBio测序数据组装得到的基因组序列的纠错方法。

背景技术

PacBio是一家测序仪公司，提供第三代测序技术测序平台，他们的测序仪产生的数据，在业内叫PacBio数据或PacBio测序数据；

Illumina是一家美国的测序仪公司，提供第二代测序技术测序平台，他们的测序仪产生的数据，在业内叫Illumina 数据或Illumina测序数据。

PacBio第三代测序技术具有超长读长、无PCR扩增、极小GC偏向等优势，越来越多的基因组是采用三代PacBio测序数据组装。但PacBio单次测序的错误率约为 15％，目前主要采用组装前对测序数据进行纠错，组装后序列不再纠错。然而，组装后序列还存在很多错误，包括单碱基错误和碱基插入缺失错误。单碱基错误和碱基插入缺失错误都对后续分析造成很大影响，比如，如果这种错误存在于基因区域，可能导致基因预测不出来或预测出错误基因；如果错误存在于重复序列区域，可能导致序列分化时间估算错误等。

发明内容

本发明的目的是解决以上提出的问题，提供一种PacBio测序数据组装后序列的纠错方法，最大程度的减少组装序列的错误。

本发明是通过以下技术方案实现的：

本发明为一种PacBio测序数据组装得到的基因组序列的纠错方法，包括以下步骤：

步骤一：使用比对软件将Illumina测序数据比对到PacBio测序数据组装得到的基因组序列上；

步骤二：根据步骤一的比对结果文件提取可能存在错误的位置和对应位置的碱基类型信息；

步骤三：根据步骤一的比对结果文件提取可能存在错误的位置的碱基类型的覆盖深度信息；

步骤四：根据可能存在错误的位置的原碱基类型的覆盖深度与对应位置其他类型碱基的覆盖深度的比值小于0.5，对PacBio测序数据组装得到的基因组序列该位置的碱基用该位置覆盖深度最大的其他类型碱基进行替换纠错，得到新的基因组序列，反之就不替换纠错。

作为优化，所述步骤一使用的Illumina测序数据样本DNA，与PacBio测序数据样本DNA来自同一样本的DNA。

作为优化，所述步骤二包含质控，所述质控是在提取出可能存在错误的位置和对应位置的碱基类型信息前去除reads比对错误数大于read长度的3％或者reads 不能完全比对上的比对信息。

作为优化，所述步骤三包含过滤，所述的过滤所述的过滤是在提取可能存在错误的位置的碱基类型的覆盖深度信息的同时去除覆盖深度低于3的错误位置信息。

作为优化，所述步骤二和步骤三中的错误的位置的碱基类型，是指单碱基错误和小于6bp的碱基插入缺失错误。

作为优化，所述步骤一中的Illumina测序数据，采用的是全基因组鸟枪法小片段构建的文库测序的数据。

作为优化，所述步骤一中的Illumina测序数据，由Hiseq2500测序仪测序而得，所述步骤一中的PacBio测序数据，由PacBio RSII测序仪测序而得。

作为优化，所述步骤一中采用的比对软件为BWA。

本发明的有益效果如下：

         本发明的方法实现了PacBio测序数据组装后序列的纠错，PacBio测序数据组装序列后主要的错误(包括单碱基错误和碱基插入缺失错误)被移除，有效的提高了组装序列的准确度；

因为组装序列是后续分析的基础，在后续分析中，有助于提高基因的结构预测准确度，重复序列预测的准确度，序列比较分析的准确性，明显降低了后续研究的错误风险。

附图说明

图1：本发明的主要流程示意图。

具体实施方式

下面结合附图和例子对本发明的实施例进行进一步详细说明：

本实施例为一种PacBio测序数据组装后序列的纠错方法，包括以下步骤：

步骤一：使用比对软件BWA将某一物种(比如白菜)Illumina测序数据比对到同一物种同一样品PacBio测序数据组装得到的基因组序列上。

步骤二：根据步骤一比对结果文件的第3列比对上序列名称信息，第4列的比对位置信息，第6列标记的插入缺失信息和第13列标记的比对不一致碱基信息，提取可能存在错误的位置和对应位置的碱基类型信息，比对结果文件信息格式为一般行业人员所熟知的；例如，比对结果文件第3列为Chr1，第4列为1120，第 6列为125M(完全比对上)，第13列为42C82，则提取可能存在错误的位置为Chr1 的第1162碱基位置，对应位置的碱基类型信息为“C”。

步骤三：根据步骤一比对结果文件的第3列比对上序列名称信息，第4列的比对位置信息，第6列标记的插入缺失信息和第13列标记的比对不一致碱基信息，在整个比对结果文件中统计可能存在错误的位置的碱基类型的覆盖深度信息，比对结果文件信息格式为一般行业人员所熟知的；例如，统计比对序列Chr1的第1162 碱基为C的共有20条reads，没有错误的比对到该位置的reads为0条。

步骤四：根据步骤三的统计，得到比对序列Chr1的第1162碱基为C的共有20条reads，没有错误的比对到该位置的reads为0条，0/20＝0，而0<0.5，则PacBio 测序数据组装得到的基因组序列的Chr1序列第1162碱基替换成“C”。

步骤一使用的Illumina测序数据样本DNA，与PacBio测序数据样本DNA来自同一样本的DNA。

步骤二包含质控，质控在步骤一之后，步骤二提取可能存在错误的位置和对应位置的碱基类型信息之前，质控是在提取出可能存在错误的位置和对应位置的碱基类型信息前去除reads比对错误数大于read长度的3％或者reads不能完全比对上的比对信息。

步骤三包含过滤，过滤与提取可能存在错误的位置的碱基类型的覆盖深度信息同时进行，过滤是在提取可能存在错误的位置的碱基类型的覆盖深度信息的同时去除覆盖深度低于3的错误位置信息。

步骤二和步骤三中的错误的位置的碱基类型，是指单碱基错误和小于6bp的碱基插入缺失错误。

步骤一中的Illumina测序数据，采用的是全基因组鸟枪法小片段构建的文库测序的数据。

步骤一中的Illumina测序数据，使用的是Hiseq2500测序仪测序而得，所述步骤一中的PacBio测序数据，使用的是PacBio RSII测序仪测序而得。

PacBio是一家测序仪公司，他们的测序仪产生的数据，称为PacBio测序数据。

Illumina是一家美国的测序仪公司，他们的测序仪产生的数据，称为Illumina测序数据。

BWA是对比软件的名称，无中文名称，在行业内直接用英文表达。

以上所述的仅是本发明的优选实施方式，应当指出，对于本技术领域中的普通技术人员来说，在不脱离本发明核心技术特征的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。