CN105279391A - 一种宏基因组16S rRNA的高通量测序数据处理及分析流程控制方法 - Google Patents

Abstract

本发明公开了一种宏基因组16S？rRNA的高通量测序数据处理及分析流程控制方法，其首先由系统生成自定义参数配置文件，再根据用户设定参数后的自定义参数文件和高通量数据处理流程模块生成与数据流程对应的批处理可执行文件；由系统执行批处理可执行文件，实现数据流程自动化，最终生成结果报告文件。从而能高效的帮助生物信息分析人员完成一套标准化的高通量数据分析流程，甚至可以让不懂高通量数据分析的科研人员自己完成高通量数据分析。从而可以达到优化科研人员的工作效率，降低科研成本的目的。本发明不仅仅可以用于宏基因组16S？rRNA高通量数据分析流程，也可用于18S和ITS等高通量测序分析流程，甚至可以在高通量测序领域通用，其实现方法简单，应用范围较为广泛。

Description

一种宏基因组16S rRNA的高通量测序数据处理及分析流程控制方法

技术领域

本发明涉及药物基因组学和计算生物学领域，具体涉及一种宏基因组16SrRNA的高通量测序数据处理及分析流程控制方法。

背景技术

宏基因组学在微生物研究中占据了非常重要的地位，宏基因组是以环境中微生物的基因组的总和为研究对象。16SrRNA(smallsubunitribosomalRNA)基因是对原核微生物进行系统化分类研究时最常用的分子标志物，广泛用于微生物生态学研究中。近年来随着高通量测序技术及数据分析方法等不断进步，大量基于16SrRNA基因的研究促进了微生物生态学的快速发展，例如：气候变化、水处理工程系统、大气污染、极端环境、人体肠道、石油污染修复和生物冶金，甚至和人体健康也密切关联。然而使用16SrRNA作为分子标志物时也存在诸多问题，例如水平基因转移、多拷贝的异质性、基因扩增效率的差异、数据分析方法的选择等，这些问题影响了微生物群落组成和多样性分析时的准确性，尤其是与高通量测序技术相关的大数据处理及分析流程控制，给相关科研工作者带来了挑战和困难，成为该领域目前急需解决的问题。

发明内容

本发明要解决的技术问题是克服现有技术宏基因组16SrRNA高通量测序数据处理中不准确性、以及分析流程中步骤繁琐、费时费力等缺陷，提供一种16SrRNA的高通量测序数据处理及分析流程控制方法。

为解决上述技术方案，本发明提供一种宏基因组16SrRNA的高通量测序数据处理及分析流程控制方法，其特征在于，其包括如下步骤：

(1)自定义参数配置文件的生成步骤；导入宏基因组16SrRNA高通量测序原始序列数据，经过筛选和拼接得到理论上有效的16SrRNA可变区全长序列，在此基础上进行生物信息学参数分析；

(2)输入步骤：用户根据需要，输入设定的各参数配置文件；

(3)分析步骤：根据参数配置文件，宏基因组高通量数据处理流程模块生成对应的自动化分析流程；

(4)执行及输出步骤：执行所描述的自动化分析流程，获得并输出宏基因组16SrRNA分析结果报告。

本发明的优选技术方案中，所述的步骤(1)中，具体包括如下步骤：

(A)导入宏基因组16SrRNA高通量测序原始序列文件，

(B)对所述的宏基因组16SrRNA高通量测序原始序列文件进行质量控制与统计，并剔除低质量序列数据，获得经过筛选的序列数据；

(C)将所述的经过筛选的数据进行拼接，组装成全长的16SrRNA可变区序列；

(D)将拼接结果进行质量控制，并去除嵌合体,得到理论上有效16SrRNA的全长序列。

本发明的优选技术方案中，所述的步骤(C)中，使用PANDseq拼接软件，对重叠区域进行比对打分，比对打分值低于0.6时将被去除，重叠区域小于5bp或者重叠区域大于2个mismatch也就去除，根据拼接结果选择有效序列在400～480bp之间的序列用于下一步分析

本发明的优选技术方案中，所述的步骤(D)中，先UCHIME软件在de-novo模式下去除嵌合体序列，然后USEARCH软件在有参模式进一步去除嵌合体序列，最终得到理论上有效的16SrRNA可变区全长序列。

本发明的优选技术方案中，所述的步骤(1)中，生物信息学参数分析包括对于获得的16SrRNA可变区全长序列进行聚类；包括输入指令采用使用UCLUST方法进行OTU聚类，OTU中序列相似性设为97％，得到OTU列表及OTU代表性序列。

本发明的优选技术方案中，所述的步骤(1)中，包括进一步对OTU代表性序列进行物种分类分析。所述的物种分类分析包括，物种进化分析，物种丰富度分析，物种鉴定分析和α多样性指数分析。

本发明的优选技术方案中，系统将多样品OTU代表性序列进行聚类与差异性分析，包括β多样性分析和多样品聚类分析。

对每个OTU选择一条代表性序列，使用RDPclassifier对代表性序列进行物种分类注释，从而得到每个样本的群落组成。

在本发明的一个实施方案中，使用RDPclassifier贝叶斯算法对97％相似水平的OTU代表序列进行分类学分析，并在各个水平统计每个样本的群落组成，比对数据库为Silva_11116SrRNAdatabase(http://www.arb-silva.de/)。

本发明的方法还可以对多个样品进行样品聚类分析，如采用Qiime平台，使用UPGMA(Unweightedpairgroupmethodwitharithmeticmean)聚类方法，基于weightedunifrac和unweightedunifrac距离矩阵，将样品进行聚类。

β多样性值为两个样本间的相异系数，反映不同样本间的多样性的差异，利用各样品序列间的进化和丰度信息计算样品间的距离，反映样品间是否有显著地微生物群落差异。在本发明的一个实施方案中，采用Qiime平台，首先利用来自不同环境样品的OTU代表序列构建一个进化树，Unifrac度量标准根据构建的进化树枝的长度计量两个不同环境样品之间的差异。UniFrac分析分为weightedunifrac和unweightedunifrac两种度量方法，两者之间差异在于是否计入不同环境样品的序列相对丰度。weightedunifrac算法在计算树枝长度时将序列的丰度信息进行加权计算，因此unweightedunifrac可以检测样品间变化的存在，而weightedunifrac可以更进一步定量的检测样品间不同谱系上发生的变异。

在本发明的方法中，使用Qiime平台，采用对序列进行随机抽样的方法，以抽到的有效序列数进行OTU的分析，并分别分别使用ACE算法、Chao算法、Shannon算法、Simpson算法、Good'sCoverage计算各α多样性指数。

Ace：用来估计群落中含有OTU数目的指数，由Chao提出，是生态学中估计物种总数常用指数之一。(http://www.mothur.org/wiki/Ace)

$S_{ACE}=\left\{\begin{array}{ccc}{S}_{abund}+\frac{S_{rare}}{C_{ACE}}+\frac{n_1}{C_{ACE}}{\widehat{\mathrm{\&gamma;}}}_{ACE}^2,& for& {\widehat{\mathrm{\&gamma;}}}_{ACE}<0.80\\ S_{abund}+\frac{S_{rare}}{C_{ACE}}+\frac{n_1}{C_{ACE}}{\widetilde{\mathrm{\&gamma;}}}_{ACE}^2,& rof& {\widehat{\mathrm{\&gamma;}}}_{ACE}\&GreaterEqual;0.80\end{array}\right.$

n_i:表示含有i条序列的OTU数目；

abunf:设定的一个OTU丰度阈值；

S_rare：低于或等于该丰度阈值的OTU数目；

S_abund:高于该丰度阈值的OTU数目；

Chao：是用Chao1算法估计样品中所含OTU数目的指数，Chao在生态学中常用来评估物种总数。(http://www.mothur.org/wiki/Chao)

$S_{chao1}=S_{obs}+\frac{n_1(n_1-1)}{2(n_2+1)}$

S_chao1:最终评估的OTU数目；

S_obs:实际测出的OTU数目；

n₁:表示含有1条序列的OTU数目；

n₂:表示含有2条序列的OTU数目；

Shannon:常用于反映α多样性指数，用来估算样品中微生物多样性。Shannon值越大，说明群落多样性越高。(http://www.mothur.org/wiki/Shannon)

$H_{shannon}=-\underset{i=1}{\overset{S_{obs}}{\&Sigma;}}\frac{n_i}{N}\ln \frac{n_i}{N}$

S_obs:实际测出的OTU数目；

n_i:表示含有i条序列的OTU数目；

N:所有测得序列数。

Simpson：辛普森多样性指数，由EdwardHughSimpson(1949)提出，在生态学中常用来定量的描述一个区域的生物多样性。Simpson指数越大，说明群落多样性越低。(http://www.mothur.org/wiki/Simpson)

$D_{\mathrm{si}mpson}=\frac{{\&Sigma;}_{i=1}^{S_{obs}}{n}_i(n_i-1)}{N(N-1)}$

S_obs:实际测出的OTU数目；

n_i:表示含有i条序列的OTU数目；

N:所有测得序列数。

Good'sCoverage:是指各样本文库的覆盖率，其数值越高，则样本中序列没有被测出的概率越低。(http://www.mothur.org/wiki/Coverage)

$C=1-\frac{n_1}{N}$

n₁:表示含有1条序列的OTU数目；

N:所有测得序列数。

Rank-abundance曲线可反映物种丰度和物种均匀度两个方面，物种丰度由曲线在横轴上的长度来反映，曲线在横轴上的范围越大，物种的丰度越高；物种均匀度由曲线的形状(平滑度)来反映，曲线越平坦，表示物种的均匀度越高。

在本发明的一个实施方案中，使用qiime平台，采用随机抽样法，以抽到的序列数与它们所能代表OTU的数目构建稀释曲线。

利用本发明，将宏基因组16SrRNA各分析步骤模块分和流程分，能够单独运行一个模块或流程中的局部分析模块，并进行模块内规定数据分析流程的快速执行。从而通过不同模块的选取，帮助科研人员和检测人员迅速完成一套高通量数据的前期数据质控、功能分析和结果报告。该工具能够优化生物信息分析人员和科研人员的工作时间，显著提高工作效率，降低科研成本，本发明的分析流程思路清晰，其实现方法简单，可广泛应用于微生物学研究工作中，也可用于人类肠道菌分析以及临床相关应用。

本发明的方法首先由系统生成自定义参数配置文件，再根据用户设定参数后的自定义参数文件和高通量数据处理流程模块生成与数据流程对应的批处理可执行文件；由系统执行批处理可执行文件，实现数据流程自动化，最终生成结果报告文件。从而能高效的帮助生物信息分析人员完成一套标准化的高通量数据分析流程，甚至可以让不懂高通量数据分析的科研人员自己完成高通量数据分析。从而可以达到优化科研人员的工作效率，降低科研成本的目的。本发明不仅仅可以用于宏基因组16SrRNA高通量数据分析流程，也可用于18S和ITS等高通量测序分析流程，甚至可以在高通量测序领域通用，其实现方法简单，应用范围较为广泛。

附图说明

图1是宏基因组16S自动化分析流程；

图2是宏基因组16S生物信息学分析步骤。

具体实施方式

以下结合具体实施例对上述方案做进一步说明。应理解，这些实施例是用于说明本发明而不是限制本发明的范围。实施例中采用的实施条件可以根据具体应用要求的条件做进一步调整，未注明的实施条件通常为常规实验中的条件。

首先对原始数据进行过滤处理，然后去除嵌合体序列，得到有效序列后进行聚类分析，每一个聚类称为一个物种操作单元(OperationalTaxonomicUnits，OTU)，对OTU的代表序列作分类学分析。基于OTU分析结果，可以对各个样本进行多种α多样性指数分析，以及对测序深度的检测；基于分类学信息，可以在各个分类水平上进行群落结构的统计分析；基于多个样本可以进行Unifrac等分析，在上述分析的基础上，可以进行一系列统计学和可视化分析。

1.测序数据(PFdata)进行数据量和测序质量的统计。

Table1.1StatisticsofPFdata

列名解释：

2.测序数据质量优化

高通量测序中通常会出现一些点突变等测序错误，而且序列末端的质量比较低，为了得到更高质量及更准确的生物信息分析结果，需要对测序原始数据进行优化处理。

分析软件：Pandaseq(v2.7)、Trimmomatic(v0.30)、Usearch(v8.0)

优化步骤及参数：

使用pandaseq(v2.7)将两条序列进行比对，根据比对的末端重叠区进行拼接，拼接时保证至少有20bp的重叠区，去除拼接结果中含有N的序列；

使用Trimmomatic(v0.30)去除引物和接头序列，去除两端质量值低于20的碱基，去除长度小于400bp的序列；

使用usearch(v8.0)将上面拼接过滤后的序列与数据库进行比对，去除其中的嵌合体序列(chimerasequence)，得到最终的有效数据。

结果展示：

effective_stat.txt各样本有效数据统计表

Table2.1Statisticsofeffectivedata

列名解释：

3.OTU分析

OTU是在群体遗传学研究中，为了便于分析，人为给某一个分类单元(属、种、分组等)设置的统一标志。在生物信息分析中，测序得到的每一条序列来自于一个菌种，要了解一个样本测序结果中的菌种、属等数目信息，就需要对序列进行归类操作。通过归类操作，将序列按照彼此的相似性归类为许多小组，一个小组就是一个OTU。通常在97％的相似水平下对所有序列进行OTU划分并进行生物信息统计分析。

软件平台：Qiime平台(v1.7)

分析方法：使用UCLUST方法进行OTU聚类，OTU中序列相似性设为97％，得到OTU列表及OTU代表性序列；

●otu_table.xls各样本OTU中序列数统计表

Table3.1OTUtable

列名解释：

4.物种分类分析

为了得到OTU对应的物种分类信息，对每个OTU选择一条代表性序列，使用RDPclassifier对代表性序列进行物种分类注释，从而得到每个样本的群落组成。

软件平台：使用Qiime平台(v1.7)

分析方法：使用RDPclassifier贝叶斯算法对97％相似水平的OTU代表序列进行分类学分析，并在各个水平统计每个样本的群落组成，比对数据库为Silva_11116SrRNAdatabase(http://www.arb-silva.de/)。

Table4.1Taxonomytreefile

列名解释：

●otu_table_mc2_w_tax_sorted_L2每个样本在门水平下个物种的百分比统计表

Table4.2TaxaStatisticsatPhylumlevel

列名解释：

●03_Taxonomy/taxa_summary_by_sample/Sample_tax_stat.xls每个样本在不同分类水平下(门、纲、目、科、属)物种种类数目统计表

Table4.3StatisticsofTaxonomicComposition

注：上表中统计的物种种类数目为物种百分比统计表中不为0的分类数目，包含other的类目。

列名解释：

5.稀释曲线

基于ObservedOTUs数，构建的稀释曲线(rarefactioncurve)，是用于描述随着样品量的加大，可能检测到的物种种类随之增加的状况，是调查样品的物种组成和预测样品中物种丰度的有效工具，在生物多样性和群落调查中，被广泛用于判断样品量是否充分以及估计物种丰富度。因此，通过稀释曲线不仅可以判断样品量是否充分，在样品量充分的前提下，运用稀释曲线还可以对物种丰富度进行预测。

分析软件：Qiime平台(v1.7)

分析方法：采用对序列进行随机抽样的方法，以抽到的序列数与它们所能代表OTU的数目构建稀释曲线。

6.α多样性指数分析

群落生态学中，α多样性主要关注单样本的多样性分析，可以反映微生物群落中物种的数目，通过一系列统计学指数的分析来估计环境群落的物种丰度和多样性。

计算菌群丰度(Communityrichness)的指数有：

ACE：用来估计群落中含有OTU数目的指数，由Chao提出，是生态学中估计物种总数常用指数之一。(http://www.mothur.org/wiki/Ace)

Chao：是用Chao1算法估计样品中所含OTU数目的指数，Chao在生态学中常用来评估物种总数。(http://www.mothur.org/wiki/Chao)

计算菌群多样性(Communitydiversity)的指数有：

Shannon:常用于反映α多样性指数，用来估算样品中微生物多样性。Shannon值越大，说明群落多样性越高。(http://www.mothur.org/wiki/Shannon)

Simpson：辛普森多样性指数，由EdwardHughSimpson(1949)提出，在生态学中常用来定量的描述一个区域的生物多样性。Simpson指数越大，说明群落多样性越低。(http://www.mothur.org/wiki/Simpson)

测序深度(Coverage)指数有：

Good'sCoverage:是指各样本文库的覆盖率，其数值越高，则样本中序列没有被测出的概率越低。(http://www.mothur.org/wiki/Coverage)

各指数计算公式可查阅报告的“4.2结果说明”部分。

分析软件：Qiime平台(v1.7)

分析方法：采用对序列进行随机抽样的方法，以抽到的有效序列数进行OTU的分析，并分别计算各α多样性指数。

●alpha_rarefaction.xlsα多样性指数

Table6.1Collationofalphadiversityresults

列名解释：

7.β多样性分析

β多样性值为两个样本间的相异系数，反映不同样本间的多样性的差异，利用各样品序列间的进化和丰度信息计算样品间的距离，反映样品间是否有显著地微生物群落差异，可通过UniFrac分析实现。

分析软件：Qiime平台(v1.7)

分析方法：首先利用来自不同环境样品的OTU代表序列构建一个进化树，Unifrac度量标准根据构建的进化树枝的长度计量两个不同环境样品之间的差异。

UniFrac分析分为weightedunifrac和unweightedunifrac两种度量方法，两者之间差异在于是否计入不同环境样品的序列相对丰度。weightedunifrac算法在计算树枝长度时将序列的丰度信息进行加权计算，因此unweightedunifrac可以检测样品间变化的存在，而weightedunifrac可以更进一步定量的检测样品间不同谱系上发生的变异。

结果展示：

weighted_unifrac.txtweightedunifrac距离矩阵

Table7.1Weightedunifracdistance

●unweighted_unifrac.txtunweightedunifrac距离矩阵

Table7.2Unweightedunifracdistance

注：表中行列表示不同样本，表中的数字表示样本间的相异系数，值越小，表示两个样本间的物种多样性差异越小。

8.PCoA分析

PCoA(PrincipalCo-ordinatesAnalysis)分析即主坐标分析，是一种研究数据相似性或差异性的可视化方法，它与PCA类似，通过一系列的特征值和特征向量进行排序后，选择主要排在前几位的特征值，找到距离矩阵中最主要的坐标，结果是数据矩阵的一个旋转，它没有改变样品点之间的相互位置关系，只是改变了坐标系统。两者的区别为PCA是基于样本的相似系数矩阵来寻找主坐标，而PCoA是基于距离矩阵来寻找主坐标。

分析软件：Qiime平台(v1.7)

分析方法：基于weightedunifrac和unweightedunifrac距离矩阵进行PCoA作图分析。

9.UPGMATree

样本聚类分析利用各样品序列间的进化信息来比较环境样品在特定的进化谱系中是否有显著的微生物群落差异。

分析软件：Qiime平台(v1.7)

分析方法：使用UPGMA(Unweightedpairgroupmethodwitharithmeticmean)聚类方法，基于weightedunifrac和unweightedunifrac距离矩阵，将样品进行聚类。

以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实例的限制，上述实例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等同物界定。

Claims (7)

1.一种宏基因组16SrRNA的高通量测序数据处理及分析流程控制方法，其特征在于，其包括如下步骤：

(1)自定义参数配置文件的生成步骤；导入宏基因组16SrRNA高通量测序原始序列数据，经过筛选和拼接得到理论上有效的16SrRNA可变区全长序列，在此基础上进行生物信息学参数分析；

(2)输入步骤：用户根据需要，输入设定的各参数配置文件；

(3)分析步骤：根据参数配置文件，宏基因组高通量数据处理流程模块生成对应的自动化分析流程；

(4)执行及输出步骤：执行所描述的自动化分析流程，获得并输出宏基因组16SrRNA分析结果报告。

2.根据权利要求1所述的宏基因组16SrRNA的高通量测序数据处理及分析流程控制方法，其特征在于，所述的步骤(1)中，具体包括如下步骤：

(A)导入宏基因组16SrRNA高通量测序原始序列文件，

(B)对所述的宏基因组16SrRNA高通量测序原始序列文件进行质量控制与统计，并剔除低质量序列数据，获得经过筛选的序列数据；

(C)将所述的经过筛选的数据进行拼接，组装成全长的16SrRNA可变区序列；

(D)将拼接结果进行质量控制，并去除嵌合体,得到理论上有效16SrRNA的全长序列。

3.根据权利要求2所述的宏基因组16SrRNA的高通量测序数据处理及分析流程控制方法，其特征在于，所述的步骤(C)中，使用PANDseq拼接软件，对重叠区域进行比对打分，比对打分值低于0.6时将被去除，重叠区域小于5bp或者重叠区域大于2个mismatch也就去除，根据拼接结果选择有效序列在400～480bp之间的序列用于下一步分析。

4.根据权利要求2所述的宏基因组16SrRNA的高通量测序数据处理及分析流程控制方法，其特征在于，所述的步骤(D)中，先UCHIME软件在de-novo模式下去除嵌合体序列，然后USEARCH软件在有参模式进一步去除嵌合体序列，最终得到理论上有效的16SrRNA可变区全长序列。

5.根据权利要求1所述的宏基因组16SrRNA的高通量测序数据处理及分析流程控制方法，其特征在于，所述的步骤(1)中，生物信息学参数分析包括对于获得的16SrRNA可变区全长序列进行聚类；包括输入指令采用使用UCLUST方法进行OTU聚类，OTU中序列相似性设为97％，得到OTU列表及OTU代表性序列。

6.根据权利要求5所述的宏基因组16SrRNA的高通量测序数据处理及分析流程控制方法，其特征在于，所述的步骤(1)中，包括进一步对OTU代表性序列进行物种分类分析。所述的物种分类分析包括，物种进化分析，物种丰富度分析，物种鉴定分析和α多样性指数分析。

7.根据权利要求5所述的宏基因组16SrRNA的高通量测序数据处理及分析流程控制方法，其特征在于，系统将多样品OTU代表性序列进行聚类与差异性分析，包括β多样性分析和多样品聚类分析。