一种基于三代测序技术的宏基因组测序数据的自动化分析方法[发明专利]

(12)发明专利申请

(10)申请公布号 CN 112133368 A(43)申请公布日 2020.12.25

(21)申请号 202011090802.8(22)申请日 2020.10.13

(71)申请人 南开大学

地址 300110 天津市南开区卫津路94号(72)发明人 刘健　田妹　陈娇　

(74)专利代理机构 北京科家知识产权代理事务

所(普通合伙) 11427

代理人 徐思波(51)Int.Cl.

G16B 20/30(2019.01)G16B 30/10(2019.01)G16B 40/00(2019.01)G16B 45/00(2019.01)

权利要求书2页 说明书5页 附图1页

(54)发明名称

一种基于三代测序技术的宏基因组测序数据的自动化分析方法(57)摘要

本发明公开了一种基于三代测序技术的宏基因组测序数据自动化分析方法，包括以下步骤：1）原始三代测序数据进行质量控制，得到clean reads；2）对clean reads拼接组装得到contigs序列；3）对contigs进行纠错优化；4）对优化后的contigs进行物种注释；5）基于物种丰度矩阵对样本多样性进行统计分析；6）基于物种丰度矩阵对样本组间差异显著物种进行统计分析；7）对优化后的contigs进行分箱；8）对分箱得到的bins进行基因注释；9）基于基因丰度矩阵对样本组间差异显著基因进行统计分析；10）基于基因注释结果，完成功能注释以及物种注释；本发明提供了从三代宏基因组测序数据处理到物种组成分析、基因组成分析与功能注释的分析方法，解决了宏基因组数据自动化精准分析难题。

CN 112133368 ACN 112133368 A

权　利　要　求　书

1/2页

1.一种基于三代测序技术的宏基因组测序数据的自动化分析方法，其特征在于，包括以下步骤：

1) 对原始测序数据进行质量控制得到clean reads；2) 对经过质量控制的clean reads进行拼接组装得到consensus序列；3) 对拼接组装得到的consensus序列进行纠错优化；4) 对经过纠错优化的consensus进行物种注释；5) 基于物种丰度矩阵对样本多样性进行统计分析；6) 基于物种丰度矩阵对样本组间差异显著物种进行统计分析；7) 对经过纠错优化的consensus进行分箱；8) 对分箱得到的bins进行基因注释；9) 基于基因丰度矩阵对样本组间差异显著基因进行统计分析；10) 基于基因注释的结果，对序列进行功能注释以及物种注释。

2.根据权利要求1所述的一种基于三代测序技术的宏基因组测序数据的自动化分析方法，其特征在于，所述步骤1）具体过程如下：

a. 根据测序实验信息、序列长度、碱基质量以及序列的GC百分比对原始测序数据进行过滤筛选，通过设置阈值的方式过滤掉超出阈值范围的序列；

b. 对上述经过质量过滤的序列，以环境中可能存在的污染源的序列作为参考序列，通过长序列比对，对非微生物DNA进行过滤筛除，得到clean reads。

3.根据权利要求1所述的一种基于三代测序技术的宏基因组测序数据的自动化分析方法，其特征在于，所述步骤2）具体过程如下：

a、Overlap：对所有clean reads进行两两比对，找到片段间的重叠信息；b、Layout：根据得到的重叠信息将存在的重叠片段建立一种组合关系，形成重叠群，即contig；

c、在重叠群中寻找一条最优的序列路径，并获得与路径对应的序列，即Consensus。4.根据权利要求1所述的一种基于三代测序技术的宏基因组测序数据的自动化分析方法，其特征在于，所述步骤3）具体过程如下：

a、将clean reads与拼接组装产生的contigs进行比对对齐；b、根据每个位点的一致性信息进行修正；c、将上述步骤b得到的结果作为更新的contigs，进行多轮纠正。

5.根据权利要求1所述的一种基于三代测序技术的宏基因组测序数据的自动化分析方法，其特征在于，所述步骤4）具体过程如下：

a、将经过纠错优化的contigs与基因组数据库进行比对；b、根据比对结果对contigs进行物种注释，并统计各物种的丰度；c、对物种丰度进行热图可视化；d、基于物种丰度矩阵，可视化物种相对丰度柱状图；e、将物种进化分支树进行可视化。

6.根据权利要求1所述的一种基于三代测序技术的宏基因组测序数据的自动化分析方法，其特征在于，所述步骤5）具体过程如下：

a、根据各样本的物种丰度矩阵，对指定分类水平的群落组成结构进行CCA主成分分析，

2

CN 112133368 A

权　利　要　求　书

2/2页

并且以二维和三维图像描述样本间的自然分布特征；

b、根据各样本的物种丰度与群落物种的群落发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行PcoA主坐标分析，并且以二维和三维图像描述样本间基于微生物系统发育关系的群落空间分布特征；

c、根据各样本的物种丰度与群落物种的群落发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行NMDS主坐标分析，并且以二维和三维图像描述样本间基于微生物系统发育关系的群落空间分布特征;

d、根据各样本的丰度矩阵绘制Specaccum物种累积曲线。

7.根据权利要求1所述的一种基于三代测序技术的宏基因组测序数据的自动化分析方法，其特征在于，所述步骤6）具体过程如下：

a、在多组样本间采用非参数因子Kruskal-Wallis秩和检验检测不同分组建丰度差异显著的物种；

b、利用Wilcoxon秩和检验检查对显著差异物种类中的所有亚种比较；c、使用线性判别分析（LDA）对数据进行降维和评估差异显著的物种并将每个物种的影响程度进行量化，并提供可视化结果。

8.根据权利要求1所述的一种基于三代测序技术的宏基因组测序数据的自动化分析方法，其特征在于，所述步骤7）具体过程如下：

a、将clean reads比对到组装得到的contigs上，统计contigs的丰度，根据contigs的丰度及其GC含量进行分箱；

b、对上述步骤a分箱结果根据完成度与污染度进行过滤；c、将上述步骤b的过滤进过进行重组装，得到最终分箱结果bins。

9.根据权利要求1所述的一种基于三代测序技术的宏基因组测序数据的自动化分析方法，其特征在于，所述步骤8）具体过程如下：

a、基于动态规划算法对bins进行基因结构预测；b、计算基因丰度矩阵。

10.根据权利要求1所述的一种基于三代测序技术的宏基因组测序数据的自动化分析方法，其特征在于，所述步骤9）具体过程如下：

a、在多组样本间采用非参数因子Kruskal-Wallis秩和检验检测不同分组建丰度差异显著的物种；

b、利用Wilcoxon秩和检验检查对显著差异物种类中的所有亚种比较；c、使用线性判别分析（LDA）对数据进行降维和评估差异显著的物种并将每个物种的影响程度进行量化，并提供可视化结果；

d、对上述基因结构预测得到的编码基因序列进行去冗余；e、将得到的非冗余的编码基因序列与功能蛋白数据库进行比对，对编码基因序列进行功能注释以及物种注释；

f、根据功能注释的结果，统计针对每个功能的丰度，并进行柱状图可视化。

3

CN 112133368 A

说　明　书

1/5页

一种基于三代测序技术的宏基因组测序数据的自动化分析

方法

技术领域

[0001]本发明一般有关三代测序技术领域，具体为一种基于三代测序技术的宏基因组测序数据的自动化分析方法。

背景技术

[0002]基于基因组序列的生物信息学已经应用在生物界的方方面面，从动物、植物再到微生物，人们投入大量资金组织实施着各种基因组计划。测序技术也因此不断更新发展，虽然二代测序在当下使用的更为频繁普遍，但由于其读长太短的特点，使得最终得到的基因组序列缺失大量片段，基因组信息不完整，对于后续的分析造成了很大的困扰。以PacBio公司的SMRT技术和Oxford Nanopore Technologies公司的纳米孔单分子技术为代表的新一代测序技术被称为第三代测序技术，与前两代测序技术相比，其最大的特点就是单分子测序，测序过程无需进行PCR扩增，并且理论上可以测定无限长度的核酸序列。相比于二代测序技术200-500bp的测序读长，三代测序技术的读长可达几十kb，甚至100kb，并且三代测序技术的速度也远远快于二代测序技术。[0003]宏基因组是基因组学一个新兴的科学研究方向，宏基因组学的研究对象是整个微生物群落，宏基因组学是研究直接从环境样本中提取的基因组遗传物质的学科。传统的微生物研究依赖于实验室培养，宏基因组学的兴起填补了无法在传统实验室中培养的微生物研究的空白。DNA测序技术的不断进步以及测序通量和分析方法的改进使得人们得以一窥未知的基因组科学领域。此外，微生物通常是以群落方式共生与某一环境中，它们的很多特性是基于整个群落环境及个体间的相互影响的，因此宏基因组学相比于做单个个体的研究更能发现其特性。

[0004]虽然三代测序技术解决了二代测序技术读长短的问题，但其错误率较高，并且有关三代测序技术的分析工具也较少。在宏基因组领域，虽然已经发明了一些基于三代测序技术的分析工具，但是一个相对完善的分析流程还未被设计与开发。发明内容

[0005]本发明的目的在于提供一种基于三代测序技术的宏基因组测序数据的自动化分析方法，以解决上述背景技术中提出的需求。[0006]为实现上述目的，本发明提供如下技术方案：一种基于三代测序技术的宏基因组测序数据的自动化分析方法,包括以下步骤：

1）对原始测序数据进行质量控制得到clean reads；2）对经过质量控制的clean reads进行拼接组装得到contigs序列；3）对拼接组装得到的contigs序列进行纠错优化；4）对经过纠错优化的contigs进行物种注释；5）基于物种丰度矩阵对样本多样性进行统计分析；

4

CN 112133368 A

说　明　书

2/5页

6）基于物种丰度矩阵对样本组间差异显著物种进行统计分析；7）对经过纠错优化的contigs进行分箱；8）对分箱得到的bins进行基因注释；9）基于基因丰度矩阵对样本组间差异显著基因进行统计分析；10）基于基因注释的结果，对序列进行功能注释以及物种注释。[0007]优选的，所述步骤1）具体过程如下：

a、根据测序实验信息、序列长度、碱基质量以及序列的GC百分比对原始测序数据进行过滤筛选，通过设置阈值的方式过滤掉超出阈值范围的序列；

b、对上述经过质量过滤的序列，以环境中可能存在的污染源的序列作为参考序列，通过序列比对，对非微生物DNA进行过滤筛除，得到clean reads。[0008]优选的，所述步骤2）具体过程如下：

a、对所有clean reads进行两两比对，找到片段间的重叠信息；b、根据得到的重叠信息将存在的重叠片段建立一种组合关系，形成重叠群；c、在重叠群众找一条最优的序列路径，并获得与路径对应的序列，从而得到最终的contigs。

[0009]优选的，所述步骤3）具体过程如下：

a、将clean reads与拼接组装产生的contigs进行比对对齐；b、根据每个位点的一致性信息进行修正；c、将上述步骤b得到的结果作为更新的contigs，进行多轮纠正。[0010]优选的，所述步骤4）具体过程如下：

a、将经过纠错优化的contigs与基因组数据库进行比对；b、根据比对结果对contigs进行物种注释，并统计各物种的 丰度；c、对物种丰度进行热图可视化；d、基于物种丰度矩阵，可视化物种相对丰度柱状图；e、将物种进化分支树进行可视化。[0011]优选的，所述步骤5）具体过程如下：

a、根据各样本的物种丰度矩阵，对指定分类水平的群落组成结构进行CCA主成分分析，并且以二维和三维图像描述样本间的自然分布特征；

b、根据各样本的物种丰度与群落物种的群落发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行PcoA主坐标分析，并且以二维和三维图像描述样本间基于微生物系统发育关系的群落空间分布特征；

c、根据各样本的物种丰度与群落物种的群落发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行NMDS主坐标分析，并且以二维和三维图像描述样本间基于微生物系统发育关系的群落空间分布特征；

d、根据各样本的丰度矩阵绘制Specaccum物种累积曲线。[0012]优选的，所述步骤6）具体过程如下：

a、在多组样本间采用非参数因子Kruskal-Wallis秩和检验检测不同分组建丰度差异显著的物种；

b、利用Wilcoxon秩和检验检查对显著差异物种类中的所有亚种比较；

5

CN 112133368 A

说　明　书

3/5页

c、使用线性判别分析（LDA）对数据进行降维和评估差异显著的物种并将每个物种的影响程度进行量化，并提供可视化结果。[0013]优选的，所述步骤7）具体过程如下：

a、将clean reads比对到组装得到的contigs上，统计contigs的丰度，根据contigs的丰度及其GC含量进行分箱；

b、对上述步骤a分箱结果根据完成度与污染度进行过滤；c、将上述步骤b的过滤进过进行重组装，得到最终分箱结果bins。[0014]优选的，所述步骤8）具体过程如下：

a、基于动态规划算法对bins进行基因结构预测；b、计算基因丰度矩阵。[0015]优选的，所述步骤9）具体过程如下：

a、在多组样本间采用非参数因子Kruskal-Wallis秩和检验检测不同分组建丰度差异显著的物种；

b、利用Wilcoxon秩和检验检查对显著差异物种类中的所有亚种比较；c、使用线性判别分析（LDA）对数据进行降维和评估差异显著的物种并将每个物种的影响程度进行量化，并提供可视化结果。[0016]优选的，所述步骤10）具体过程如下：

a、对上述基因结构预测得到的编码基因序列进行去冗余；b、将得到的非冗余的编码基因序列与功能蛋白数据库进行比对，对编码基因序列进行功能注释以及物种注释；

c、根据功能注释的结果，统计针对每个功能的丰度，并进行柱状图可视化。[0017]与现有技术相比，本发明的有益效果是：解决了当前宏基因组领域对于基于三代测序技术的自动化分析流程的需求，为研究人员提供便利。附图说明

[0018]图1为本发明流程图。

具体实施方式

[0019]下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。[0020]请参阅图1，本发明提供一种技术方案：本发明提供如下技术方案：一种基于三代测序技术的宏基因组分析的自动化分析方法,包括以下步骤：

1）对原始测序数据进行质量控制得到clean reads；2）对经过质量控制的clean reads进行拼接组装得到contigs序列；3）对拼接组装得到的contigs序列进行纠错优化；4）对经过纠错优化的contigs进行物种注释；5）基于物种丰度矩阵对样本多样性进行统计分析；

6

CN 112133368 A

说　明　书

4/5页

6）基于物种丰度矩阵对样本组间差异显著物种进行统计分析；7）对经过纠错优化的contigs进行分箱；8）对分箱得到的bins进行基因注释；9）基于基因丰度矩阵对样本组间差异显著基因进行统计分析；10）基于基因注释的结果，对序列进行功能注释以及物种注释。[0021]本发明中，所述步骤1）具体过程如下：

a、根据测序实验信息、序列长度、碱基质量以及序列的GC百分比对原始测序数据进行过滤筛选，通过设置阈值的方式过滤掉超出阈值范围的序列；

b、对上述经过质量过滤的序列，以环境中可能存在的污染源的序列作为参考序列，通过序列比对，对非微生物DNA进行过滤筛除，得到clean reads。[0022]具体地可以使用NanoFilt软件进行过滤。[0023]本发明中，所述步骤2）具体过程如下：

a、对所有clean reads进行两两比对，找到片段间的重叠信息；b、根据得到的重叠信息将存在的重叠片段建立一种组合关系，形成重叠群；c、在重叠群众找一条最优的序列路径，并获得与路径对应的序列，从而得到最终的contigs。

[0024]具体地地可以使用Flye软件进行实现。[0025]本发明中，所述步骤3）具体过程如下：

a、将clean reads与拼接组装产生的contigs进行比对对齐；b、根据每个位点的一致性信息进行修正；c、将上述步骤b得到的结果作为更新的contigs，进行多轮纠正。[0026]具体地可以使用medaka软件实现。[0027]本发明中，所述步骤4）具体过程如下：

a、将经过纠错优化的contigs与基因组数据库进行比对；b、根据比对结果对contigs进行物种注释，并统计各物种的丰度；c、对物种丰度进行热图可视化；d、基于物种丰度矩阵，可视化物种相对丰度柱状图；e、将物种进化分支树进行可视化。

[0028]具体地可以使用DIAMOND比对软件以及R语言实现可视化。[0029]本发明中，所述步骤5）具体过程如下：

a、根据各样本的物种丰度矩阵，对指定分类水平的群落组成结构进行CCA主成分分析，并且以二维和三维图像描述样本间的自然分布特征；

b、根据各样本的物种丰度与群落物种的群落发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行PcoA主坐标分析，并且以二维和三维图像描述样本间基于微生物系统发育关系的群落空间分布特征；

c、根据各样本的物种丰度与群落物种的群落发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行NMDS主坐标分析，并且以二维和三维图像描述样本间基于微生物系统发育关系的群落空间分布特征；

d、根据各样本的丰度矩阵绘制Specaccum物种累积曲线。

7

CN 112133368 A[0030]

说　明　书

5/5页

具体地，上述步骤可以使用R语言实现。

[0031]本发明中，所述步骤6）具体过程如下：

a、在多组样本间采用非参数因子Kruskal-Wallis秩和检验检测不同分组建丰度差异显著的物种；

b、利用Wilcoxon秩和检验检查对显著差异物种类中的所有亚种比较；c、使用线性判别分析（LDA）对数据进行降维和评估差异显著的物种并将每个物种的影响程度进行量化，并提供可视化结果。

[0032]具体地上述步骤可以使用LEfSe软件实现。[0033]本发明中，所述步骤7）具体过程如下：

a、将clean reads比对到组装得到的contigs上，统计contigs的丰度，根据contigs的丰度及其GC含量进行分箱；

b、对上述步骤a分箱结果根据完成度与污染度进行过滤；c、将上述步骤b的过滤进过进行重组装，得到最终分箱结果bins。[0034]具体地，上述步骤可以使用metaWrap软件实现。[0035]本发明中，所述步骤8）具体过程如下：

a、基于动态规划算法对bins进行基因结构预测；b、计算基因丰度矩阵。[0036]具体地，可以使用metaProdigal软件实现。[0037]本发明中，所述步骤9）具体过程如下：

a、在多组样本间采用非参数因子Kruskal-Wallis秩和检验检测不同分组建丰度差异显著的物种；

b、利用Wilcoxon秩和检验检查对显著差异物种类中的所有亚种比较；c、使用线性判别分析（LDA）对数据进行降维和评估差异显著的物种并将每个物种的影响程度进行量化，并提供可视化结果。[0038]具体地，上述步骤可以使用LEfSe软件实现。[0039]本发明中，所述步骤10）具体过程如下：

a、对上述基因结构预测得到的编码基因序列进行去冗余；b、将得到的非冗余的编码基因序列与功能蛋白数据库进行比对，对编码基因序列进行功能注释以及物种注释；

c、根据功能注释的结果，统计针对每个功能的丰度，并进行柱状图可视化。[0040]具体地，可以使用minimap2软件进行比对，使用R语言进行可视化。[0041]综上所述，本发明通过整合分析结果更为准确的分析方法以及较为全面的分析模块的方法，开发了基于三代测序的宏基因组测序数据的自动化分析方法，从而解决了当前宏基因组领域对于基于三代测序技术的自动化分析流程的需求，为研究人员提供便利。[0042]对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为所涉及的权利要求。

8

CN 112133368 A

说　明　书　附　图

1/1页

图1

9