CN106815492B - 一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法 - Google Patents
一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法 Download PDFInfo
- Publication number
- CN106815492B CN106815492B CN201611187576.9A CN201611187576A CN106815492B CN 106815492 B CN106815492 B CN 106815492B CN 201611187576 A CN201611187576 A CN 201611187576A CN 106815492 B CN106815492 B CN 106815492B
- Authority
- CN
- China
- Prior art keywords
- otu
- sample
- abundance
- sequence
- matrix
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000004458 analytical method Methods 0.000 title claims abstract description 67
- 239000000203 mixture Substances 0.000 title claims abstract description 54
- 238000000034 method Methods 0.000 title claims abstract description 27
- 108020004465 16S ribosomal RNA Proteins 0.000 title claims abstract description 13
- 230000001580 bacterial effect Effects 0.000 title claims abstract description 11
- 108090000623 proteins and genes Proteins 0.000 title claims abstract description 11
- 238000012163 sequencing technique Methods 0.000 claims abstract description 29
- 230000008569 process Effects 0.000 claims abstract description 13
- 238000007405 data analysis Methods 0.000 claims abstract description 10
- 239000011159 matrix material Substances 0.000 claims description 100
- 238000001558 permutation test Methods 0.000 claims description 11
- 238000005070 sampling Methods 0.000 claims description 9
- 241000894007 species Species 0.000 claims description 9
- 238000007637 random forest analysis Methods 0.000 claims description 8
- 238000007689 inspection Methods 0.000 claims description 7
- 238000012797 qualification Methods 0.000 claims description 6
- 238000004422 calculation algorithm Methods 0.000 claims description 5
- 230000002452 interceptive effect Effects 0.000 claims description 5
- 238000000585 Mann–Whitney U test Methods 0.000 claims description 4
- 238000009825 accumulation Methods 0.000 claims description 4
- 230000008901 benefit Effects 0.000 claims description 4
- 244000005700 microbiome Species 0.000 claims description 4
- 238000011160 research Methods 0.000 claims description 4
- 241000196324 Embryophyta Species 0.000 claims description 3
- 238000000540 analysis of variance Methods 0.000 claims description 3
- 239000000090 biomarker Substances 0.000 claims description 3
- 239000003086 colorant Substances 0.000 claims description 3
- 238000010219 correlation analysis Methods 0.000 claims description 3
- 238000002790 cross-validation Methods 0.000 claims description 3
- 238000003066 decision tree Methods 0.000 claims description 3
- 238000005259 measurement Methods 0.000 claims description 3
- 230000007269 microbial metabolism Effects 0.000 claims description 3
- 244000005706 microflora Species 0.000 claims description 3
- 238000002887 multiple sequence alignment Methods 0.000 claims description 3
- 238000010587 phase diagram Methods 0.000 claims description 3
- 238000000513 principal component analysis Methods 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims description 2
- 238000011161 development Methods 0.000 claims description 2
- 238000006073 displacement reaction Methods 0.000 claims description 2
- 239000012530 fluid Substances 0.000 claims description 2
- 238000013507 mapping Methods 0.000 claims description 2
- 238000004321 preservation Methods 0.000 claims description 2
- 238000003908 quality control method Methods 0.000 claims description 2
- 238000012360 testing method Methods 0.000 claims description 2
- 239000004744 fabric Substances 0.000 claims 1
- 238000012216 screening Methods 0.000 claims 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims 1
- 238000007418 data mining Methods 0.000 abstract 1
- 239000002585 base Substances 0.000 description 11
- 238000012800 visualization Methods 0.000 description 8
- 230000008859 change Effects 0.000 description 2
- 230000007613 environmental effect Effects 0.000 description 2
- 238000012177 large-scale sequencing Methods 0.000 description 2
- 108091032973 (ribonucleotides)n+m Proteins 0.000 description 1
- 235000017284 Pometia pinnata Nutrition 0.000 description 1
- 238000012300 Sequence Analysis Methods 0.000 description 1
- 238000012896 Statistical algorithm Methods 0.000 description 1
- 239000003513 alkali Substances 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 101150039352 can gene Proteins 0.000 description 1
- 238000010835 comparative analysis Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000002068 genetic effect Effects 0.000 description 1
- 238000012165 high-throughput sequencing Methods 0.000 description 1
- 230000008676 import Effects 0.000 description 1
- 238000009434 installation Methods 0.000 description 1
- 230000004807 localization Effects 0.000 description 1
- 230000035800 maturation Effects 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 108700022487 rRNA Genes Proteins 0.000 description 1
- 210000003705 ribosome Anatomy 0.000 description 1
- 238000002864 sequence alignment Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 230000033772 system development Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
Landscapes
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Chemical & Material Sciences (AREA)
- Analytical Chemistry (AREA)
- Biophysics (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Health & Medical Sciences (AREA)
- Engineering & Computer Science (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Biotechnology (AREA)
- Evolutionary Biology (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Theoretical Computer Science (AREA)
- Information Retrieval, Db Structures And Fs Structures Therefor (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本发明公开的一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法，其提供的16S rRNA测序数据分析流程以测序原始序列数据作为输入,调用业界标准的分析工具(如：Mothur、QIIME等)，最终对数据进行可视化，并得到易于解读的分析结果。本发明包含了目前流行的主流分析项目，同时分析内容实现模块化，数据挖掘分析的方法更多样、更深入，可以根据不同的需要结合不同的分析模块内容，先后顺序的流程安排也更合理；此外，消除了测序深度不一导致的分析误差，使分析结果更全面、准确、可靠。
Description
技术领域：
本发明一般有关于分子生物学技术领域，特别涉及高通量测序数据分析技术领域，并且更具体地说，涉及一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法。
背景技术：
新一代高通量测序技术大幅度降低了测序的时间和成本，使得大规模测序逐渐成为常规的研究和检测手段，测序产生的数据量急剧增加。如何高效地分析这些数据，已成为迫切需要解决的问题。
目前高通量测序数据分析工具很多，进行分析序列信息的生物信息学工具纷繁复杂，对于分析菌群微生态的大规模测序数据，也已开发产生多种成熟的分析工具可供使用。其中，针对于菌群多样性和组成谱检测这一广泛应用的研究策略，绝大多数基于微生物核糖体RNA(rRNA)基因的序列分析工具并不能一次性满足研究人员的大部分分析需求，往往需要多次调用不同的分析工具(甚至不同的计算机平台系统)以完成所有相关分析。
当前用于rRNA基因序列分析的工具分为两大类：基于Web的工具和本地安装的工具。基于Web的工具托管在服务器上，为用户提供Web界面。例如Ribosomal DatabaseProject(RDP)的分类鉴定工具等；研究人员可以通过互联网上传测得的序列数据，并使用Web界面配置可选参数进行分析。但是对于一次上传的数据量存在限制，并受网络环境所约束，不适合大规模数据分析，并且远程用户无法根据自己的需求来自定义后端分析软件，并不方便。本地化的16S rRNA基因分析工具包括Mothur和QIIME等。使用这些本地工具时，不需要将数据上传到远程服务器，只需在Linux服务器/工作站合理配置安装即能投入使用。这些工具大多是开源形式，允许研究人员根据需要自定义软件。然而，很多情况下，根据不同的分析需求，实现一个完整的分析流程仍需要对众多工具进行整合。如何能正确高效地选择并整合这些工具已成为迫切需求。
现有的16S rRNA基因的细菌群落组成和多样性分析流程中，分析内容较为基础和简单，未包含目前流行的多项主流分析项目：Specaccum物种累积曲线、LEfSe分析、Wilcoxon秩和检验等。在一些情况下，无法满足研究人员的分析需求。此外，原有分析流程在进行后续PCA、PCoA等比较分析时，未对数据进行拉平处理，由此将会引入测序深度不一导致的分析误差。
发明内容
本发明的一个目的是为了克服现有工具中固有的弱点并结合其优势，提供一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法，该方法具有可靠的流程，已确保分析结果的准确。并满足多种分析需求。
本发明的另一个目的是为了提供一种16S rRNA测序数据分析流程，实现各个环节的高效自动化管理和分析，从而节省时间成本，减轻研究人员的数据分析负担。
为了实现上述目的，本发明所采用的技术方案如下：
一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法，其提供的16S rRNA测序数据分析流程以测序原始序列数据作为输入,调用业界标准的分析工具(如：Mothur、QIIME等)，最终对数据进行可视化，并得到易于解读的分析结果，具体包括以下步骤：
1)通过原始序列的测序质量值、模糊碱基数目、序列长度、引物序列和barcode序列的匹配度信息，对原始序列进行过滤和质量控制，并检查和剔除嵌合体，获得高质量序列；
2)对步骤1)获得的高质量序列的长度分布进行统计；
3)对步骤1)获得的高质量序列按97％的序列相似度进行归并和OTU划分，并选取每个OTU中丰度最高的序列作为该OTU的代表序列，随后，根据每个OTU在每个样本中所包含的序列数，构建OTU在各样本中丰度的矩阵文件；
4)通过将OTU代表序列与对应数据库的模板序列相比对，获取每个OTU所对应的分类学信息；
5)将丰度值低于全体样本测序总量0.001％(十万分之一)的OTU去除，并将去除了稀有OTU的此丰度矩阵用于后续的一系列分析；
6)根据获得的OTU丰度矩阵，计算各样本组共有OTU的数量，并通过Venn图直观地呈现各样本组所共有和独有OTU所占的比例；
7)对OTU丰度矩阵中每个样本的序列总数在不同测序深度下依次随机抽样，以每个深度下抽取到的序列数及其对应的OTU数绘制稀疏曲线；
8)对OTU丰度矩阵中每个样本所对应的OTU总数绘制Specaccum物种累积曲线；
9)对OTU及其对应的丰度值经Log2对数转换绘制各样本的丰度等级曲线；
对OTU丰度矩阵中的全体样本根据最低测序深度统一进行随机重抽样(即序列拉平处理)，随后，分别对每个样本计算四种多样性指数；
10)根据OTU划分和分类地位鉴定结果，可以获得每个样本在各分类水平(界/门/纲/目/科/属/种等)的具体组成；
11)获取各样本在指定分类水平上的组成和丰度分布表，并通过饼图、柱状图或面积图呈现分析结果，根据研究对象是单个或多个群落样本，绘图结果可能会以不同方式进行展示；
12)获取各样本在指定分类水平上的组成和绝对丰度分布表，调用Metastats的统计学算法，对指定分类水平的各个分类单元在样本组之间的序列量即绝对丰度差异进行两两比较检验；
13)获取各样本在指定分类水平上的组成和相对丰度分布表，进行LEfSe分析，筛选关键的生物标记物；
14)获取各样本在指定分类水平上的组成和相对丰度分布表，对各分类单元在两组样本中的丰度分布差异进行Wilcoxon秩和检验或Welch’s t检验，从而获得在两组中存在显著性差异的分类单元；
15)获取各样本在指定分类水平上的组成和相对丰度分布表，对各分类单元在两个样本中的丰度分布差异进行Fisher’s检验，从而获得在两个样本中存在显著性差异的分类单元；
16)获取各样本在指定分类水平上的组成和相对丰度分布表，对各分类单元在多组样本中的丰度分布差异进行ANOVA方差分析/Kruskal-Wallis H检验，从而获得在多组样本中存在显著性差异的分类单元；
18)对前述OTU代表序列，通过PyNAST和MAFFT等工具进行多序列比对，之后通过FastTree工具构建OTU代表序列的系统发育树，该文件以Newick格式保存；
19)根据前述OTU丰度矩阵和OTU划分和分类地位鉴定结果，将每个样本所含有的OTU的丰度信息和分类学组成数据映射至NCBI Taxonomy所提供的微生物分类等级树，统一呈现所有样本在各分类水平的具体组成；
20)获取各样本在指定分类水平上的组成和相对丰度分布表，对样本总体在各分类水平的组成构建等级树，同时以不同颜色区分各分类单元，并通过节点大小反映它们的丰度分布；
21)获取各样本在指定分类水平上的组成和丰度分布表，通过Krona软件进行群落分类学组成的交互展示；
22)根据前述OTU丰度矩阵和OTU划分和分类地位鉴定结果，构建交互式OTU热图；
23)获取各样本在指定分类水平上的组成和相对丰度分布表，对丰度前50位的分类单元进行聚类分析并绘制热图；
24)获取各样本在指定分类水平上的组成和相对丰度分布表，对指定分类水平的群落组成结构进行PCA主成分分析，并且以二维和三维图像描述样本间的自然分布特征；
25)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行PcoA主坐标分析，并且以二维和三维图像描述样本间基于微生物系统发育关系的群落空间分布特征；
26)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行NMDS非度量多维尺度分析，通过二维或三维排序图描述群落样本的结构分布；
27)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，对加权及非加权距离矩阵分别进行UPGMA聚类分析；
28)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，对加权及非加权距离矩阵分别进行组内和组间距离均值进行t检验，并通过1000次蒙特卡罗置换检验判断统计显著性，从而衡量组内和组间距离差异；
30)获取各样本在指定分类水平上的组成和相对丰度分布表，通过三元相图分析显示出不同物种在各样本中的相对丰度差异；
31)获取各样本在指定分类水平上的组成和相对丰度分布表，对指定分类水平的相对丰度矩阵进行RDA冗余分析，通过1000次置换检验确定统计显著性，并生成包含“样本—分类单元—影响因素”三种元素排序图；
32)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行CAP主坐标典型相关分析分析，并作1000次置换检验，确定组间差异是否具有统计学显著性，使用CAP分析得到的前两维数据生成CAP约束排序图；
33)获取各样本在指定分类水平上的组成和相对丰度分布表和样本分组数据构建PLS-DA判别模型，并根据分析得到的前两维数据生成PLS-DA约束排序图；
34)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行Adonis/PERMANOVA分析，并作999次置换检验确定组间差异是否具有统计学显著性；
35)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行ANOSIM分析，通过对样本距离等级排序来判断样本组内和组间差异的大小，并通过置换检验评价原始样本组间差异的统计学显著性；
36)根据前述对OTU丰度矩阵中的全体样本根据最低测序深度统一进行随机重抽样后，对重抽样OTU丰度矩阵使用随机森林算法挑取丰度分布在不同组间存在显著差异的OTU。挑取OTU时使用1000棵随机森林决策树进行建模，并以10倍交叉验证估计“基线”误差的大小随机森林分析；
37)获取各样本在指定分类水平上的组成和相对丰度分布表，计算丰度位于前50位的分类单元之间的Spearman等级相关系数，对其中|rho|>0.8且P值<0.01的相关优势属构建关联网络，并进行可视化。
38)根据前述对OTU丰度矩阵中的全体样本根据最低测序深度统一进行随机重抽样后，去除绝对丰度<2的OTU，进行二分网络分析，并进行可视化，样本和OTU根据弹簧镶嵌模型排布。
39)对前述获得的高质量序列，使用PICRUSt软件，根据KEGG数据库中微生物代谢功能的类别对群落样本进行预测，并通过条形图展示各样本中编码的功能基因类别及其丰度。
本发明的有益效果是：
改进后的流程包含了目前流行的主流分析项目，同时分析内容实现模块化，数据挖掘分析的方法更多样、更深入，可以根据不同的需要结合不同的分析模块内容，先后顺序的流程安排也更合理；此外，消除了测序深度不一导致的分析误差，使分析结果更全面、准确、可靠。
改进后的流程通过自动化脚本控制运行，一站式输出可视化图形和分析结果,简化了数据分析过程中的人工操作,提高了运行效率。研究人员只需导入数据及根据分析需求简单调整脚本即可完成操作，从而节省时间，减轻分析的工作负担。
附图说明
图1是根据本发明构造的的一个细菌16S rRNA测序数据分析流程图。
具体实施方式
在具体实施方式中，该方法如图1所示包括以下步骤：
步骤1：对于Illumina Miseq平台原始双端测序数据，以原始下机数据作为输入数据执行MiSeqQuality16S.pl脚本，窗口大小为10bp，步长为1bp，从5’端第一个碱基位置开始移动，要求窗口中碱基平均质量≥Q20(即碱基平均测序准确率≥99％)，从第一个平均质量值低于Q20的窗口处截断序列，并要求截断后的序列长度≥150bp，且不允许存在模糊碱基N。随后，利用FLASH软件，对通过质量初筛的双端序列根据重叠碱基进行配对连接：要求Read 1和Read 2两条序列的重叠碱基长度≥10bp，且不允许碱基错配。最后，以FLASH拼接输出序列和barcode序列作为输入文件，将执行MiSeqRawDataSplit.pl脚本，把连接后的序列识别分配入对应样本(要求Index序列完全匹配)，从而获得每个样本的有效序列。
对于Roche 454FLX+平台原始测序数据，以原始下机数据和barcode序列作为输入文件，以原始下机数据和barcode序列作为输入数据通过sfffile工具把序列识别分配入对应样本(要求Index序列完全匹配)，以sfffile输出序列作为输入文件通过sffinfo工具输出FASTA格式文件和QUAL格式文件，随后通过QIIME进行过滤，要求碱基平均质量≥Q25，不允许存在模糊碱基，从而获得每个样本的有效序列。
以上述有效序列和样本分组信息作为输入文件，执行seq_process.py脚本，调用QIIME首先识别疑问序列。要求序列长度≥150bp，剔除：1)5’端引物错配碱基数>1的序列；2)含有连续相同碱基数>8的序列。调用USEARCH检查和剔除嵌合体，获得高质量序列；调用R脚本计算高质量序列的长度分布，并输出可视化结果；对前述获得的高质量序列使用QIIME软件，调用UCLUST这一序列比对工具，按97％的序列相似度进行归并和OTU划分，并选取每个OTU中丰度最高的序列作为该OTU的代表序列。随后，根据每个OTU在每个样本中所包含的序列数，构建OTU在各样本中丰度的矩阵文件(即OTU table)；同时，在QIIME软件中使用默认参数，通过将OTU代表序列与对应数据库的模板序列相比对，获取每个OTU所对应的分类学信息；此外，在QIIME软件中对OTU代表序列，通过PyNAST或者MAFFT工具进行多序列比对，之后调用FastTree工具构建OTU代表序列的系统发育树。
步骤2：通过BIOM文件处理工具和modify_otu_table.R脚本，以步骤1输出OTU丰度矩阵作为输入文件将丰度值低于全体样本测序总量0.001％(十万分之一)的OTU去除。
步骤:3：以步骤2输出OTU丰度矩阵和分组信息作为输入文件，通过creat_venn_map.py、plot_venn_graph.R脚本计算各样本(组)共有OTU的数量，并进行可视化。
步骤4：以步骤2输出OTU丰度矩阵作为输入文件，通过create_table2_4_in_report.R、bar.R脚本，计算各样本中能分类至门、纲、目、科、属、种各分类水平的OTU数及可视化。
步骤5：以步骤2输出OTU丰度矩阵作为输入文件，使用BIOM工具和QIIME软件绘制稀疏曲线。
步骤6：以步骤2输出OTU丰度矩阵作为输入文件，通过species_curve.R脚本绘制Specaccum物种累积曲线。
步骤7：以步骤2输出OTU丰度矩阵作为输入文件，通过plot_rank_abundance_graph.R脚本绘制各样本的丰度等级曲线。
步骤8：以步骤2输出OTU丰度矩阵作为输入文件，通过normalize_rarefied_otu_table.sh脚本对全体样本根据最低测序深度统一进行随机重抽样100次，并计算平均值，剔除平均值<1的OTU，从而最大限度消除测序深度不一导致的分析误差并保证结果的准确性。
步骤9：以步骤8输出OTU丰度矩阵作为输入文件，通过QIIME软件分别对每个样本计算四种多样性指数。
步骤10：以步骤8输出OTU丰度矩阵作为输入文件，通过bar.R脚本对门、纲、目、科、属、种六个分类水平各自含有的微生物类群数可视化。
步骤11：以步骤2输出OTU丰度矩阵作为输入文件，通过QIIME软件和bar-phylum.R脚本获取各样本在门、纲、目、科、属五个分类水平上的组成和丰度分布表及可视化。
步骤12：以步骤11输出OTU丰度矩阵作为输入文件，通过metats.sh、diff_group2.R脚本，调用Metastats的统计学算法，对指定分类水平的各个分类单元在样本(组)之间的序列量(即绝对丰度)差异进行两两比较检验。
步骤13：以步骤11输出OTU丰度矩阵作为输入文件，通过LEfSe本地工具进行LEfSe分析，筛选关键的生物标记物。
步骤14：以步骤11输出OTU丰度矩阵和分组信息作为输入文件，通过wilcox_box.R脚本对各分类单元在两组样本中的丰度分布差异进行Wilcoxon秩和检验，从而获得在两组中存在显著性差异的分类单元并输出可视化结果。
步骤15：以步骤11输出OTU丰度矩阵和分组信息作为输入文件，通过STAMP软件对各分类单元在两组样本中的丰度分布差异进行Welch’s t检验，从而获得在两组中存在显著性差异的分类单元并输出可视化结果。
步骤16：以步骤11输出OTU丰度矩阵和分组信息作为输入文件，通过STAMP软件对各分类单元在两个样本中的丰度分布差异进行Fisher’s检验，从而获得在两个样本中存在显著性差异的分类单元并输出可视化结果。
步骤18：以步骤11输出OTU丰度矩阵和分组信息作为输入文件，通过STAMP软件对各分类单元在多组样本中的丰度分布差异进行ANOVA方差分析/Kruskal-Wallis H检验，从而获得在多组样本中存在显著性差异的分类单元并输出可视化结果。
步骤19：以步骤2输出OTU丰度矩阵作为输入文件，通过MEGAN软件将每个样本所含有的OTU的丰度信息和分类学组成数据映射至NCBI Taxonomy所提供的微生物分类等级树，统一呈现所有样本在各分类水平的具体组成。
步骤20：以步骤11输出OTU丰度矩阵作为输入文件，通过graphlan.sh脚本对样本总体在各分类水平的组成构建等级树，同时以不同颜色区分各分类单元，并通过节点大小反映它们的丰度分布。
步骤21：以步骤11输出OTU丰度矩阵作为输入文件，通过divide_sample_from_taxa.py脚本进行群落分类学组成的交互展示。
步骤22：以步骤2和步骤11输出OTU丰度矩阵和分组信息作为输入文件，通过QIIME软件和Heatmap.R绘制热图。
步骤23：以步骤11输出OTU丰度矩阵和分组信息作为输入文件，通过plot_pca_graph.R、plot_3d_pca_graph.R脚本对指定分类水平的群落组成结构进行PCA主成分分析，并且可视化。
步骤24：以步骤8输出OTU丰度矩阵、步骤2输出系统发育数和分组信息作为输入文件，通过QIIME软件进行基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行PcoA主坐标分析，并且以二维和三维图像描述样本间基于微生物系统发育关系的群落空间分布特征；进行NMDS非度量多维尺度分析，通过二维或三维排序图描述群落样本的结构分布；进行UPGMA聚类分析；对加权及非加权距离矩阵分别进行组内和组间距离均值进行t检验，并通过1000次蒙特卡罗置换检验判断统计显著性，从而衡量组内和组间距离差异。
步骤25：以步骤23输出nmds坐标和分组信息作为输入文件，通过plot_nmds.R脚本进行可视化。
步骤26：以步骤11输出OTU丰度矩阵作为输入文件，通过make_ternary_plot.R脚本进行三元相图分析显示出不同物种在各样本中的相对丰度差异。
步骤27：以步骤11输出OTU丰度矩阵、环境因子和分组信息作为输入文件，通过rda.R脚本对指定分类水平的相对丰度矩阵进行RDA冗余分析，通过1000次置换检验确定统计显著性，并生成包含“样本—分类单元—影响因素”三种元素排序图。
步骤28：以步骤23输出距离矩阵和分组信息作为输入文件，通过R软件BiodiversityR软件包和ggplot2软件包进行CAP主坐标典型相关分析分析，并作1000次置换检验，确定组间差异是否具有统计学显著性。使用CAP分析得到的前两维数据生成CAP约束排序图。
步骤29：以步骤11输出OTU丰度矩阵、环境因子和分组信息作为输入文件，通过PLSDA.R脚本构建PLS-DA判别模型并作图。
步骤30：以步骤23输出距离矩阵和分组信息作为输入文件，通过QIIME软件进行Adonis/PERMANOVA分析，并作999次置换检验确定组间差异是否具有统计学显著性。
步骤31：以步骤23输出距离矩阵和分组信息作为输入文件，通过QIIME软件进行ANOSIM分析，通过对样本距离等级排序来判断样本组内和组间差异的大小，并通过置换检验评价原始样本组间差异的统计学显著性。
步骤32：以步骤8输出OTU丰度矩阵和分组信息作为输入文件，通过QIIME软件使用随机森林算法挑取丰度分布在不同组间存在显著差异的OTU。挑取OTU时使用1000棵随机森林决策树进行建模，并以10倍交叉验证估计误差的大小并完成随机森林分析。
步骤33：以步骤11输出OTU丰度矩阵作为输入文件，通过network_spearman.R脚本和Cytoscape软件计算丰度位于前50位的分类单元之间的Spearman等级相关系数，对其中|rho|>0.8且P值<0.01的相关优势属构建关联网络及可视化。
步骤34：以步骤8输出OTU丰度矩阵作为输入文件，通过QIIME软件和Cytoscape软件进行二分网络分析，并进行可视化，样本和OTU根据弹簧镶嵌模型排布。
步骤35：以步骤1获得的高质量序列作为输入文件，使用PICRUSt软件，根据KEGG数据库中微生物代谢功能的类别对群落样本进行预测，并通过条形图展示各样本中编码的功能基因类别及其丰度。
上述步骤中除少数步骤外，都可以通过linux操作系统，以命令行的形式整合到shell脚本中，从而方便一次性执行，实现整个分析流程的自动化，提高分析效率。也可以根据需要单独执行。
Claims (1)
1.一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法，其提供的16SrRNA测序数据分析流程以测序原始序列数据作为输入,调用业界标准的分析工具，最终对数据进行可视化，并得到易于解读的分析结果，其特征在于，具体包括以下步骤：
1)通过原始序列的测序质量值、模糊碱基数目、序列长度、引物序列和barcode序列的匹配度信息，对原始序列进行过滤和质量控制，并检查和剔除嵌合体，获得高质量序列；
2)对步骤1)获得的高质量序列的长度分布进行统计；
3)对步骤1)获得的高质量序列按97％的序列相似度进行归并和OTU划分，并选取每个OTU中丰度最高的序列作为该OTU的代表序列，随后，根据每个OTU在每个样本中所包含的序列数，构建OTU在各样本中丰度的矩阵文件；
4)通过将OTU代表序列与对应数据库的模板序列相比对，获取每个OTU所对应的分类学信息；
5)将丰度值低于全体样本测序总量0.001％的OTU去除，并将去除了稀有OTU的此丰度矩阵用于后续的一系列分析；
6)根据获得的OTU丰度矩阵，计算各样本组共有OTU的数量，并通过Venn图直观地呈现各样本组所共有和独有OTU所占的比例；
7)对OTU丰度矩阵中每个样本的序列总数在不同测序深度下依次随机抽样，以每个深度下抽取到的序列数及其对应的OTU数绘制稀疏曲线；
8)对OTU丰度矩阵中每个样本所对应的OTU总数绘制Specaccum物种累积曲线；
9)对OTU及其对应的丰度值经Log2对数转换绘制各样本的丰度等级曲线；对OTU丰度矩阵中的全体样本根据最低测序深度统一进行随机重抽样(即序列拉平处理)，随后，分别对每个样本计算四种多样性指数；
10)根据OTU划分和分类地位鉴定结果，可以获得每个样本在各分类水平的具体组成；
11)获取各样本在指定分类水平上的组成和丰度分布表，并通过饼图、柱状图或面积图呈现分析结果，根据研究对象是单个或多个群落样本，绘图结果以不同方式进行展示；
12)获取各样本在指定分类水平上的组成和绝对丰度分布表，调用Metastats的统计学算法，对指定分类水平的各个分类单元在样本组之间的序列量即绝对丰度差异进行两两比较检验；
13)获取各样本在指定分类水平上的组成和相对丰度分布表，进行LEfSe分析，筛选关键的生物标记物；
14)获取各样本在指定分类水平上的组成和相对丰度分布表，对各分类单元在两组样本中的丰度分布差异进行Wilcoxon秩和检验或Welch’s t检验，从而获得在两组中存在显著性差异的分类单元；
15)获取各样本在指定分类水平上的组成和相对丰度分布表，对各分类单元在两个样本中的丰度分布差异进行Fisher’s检验，从而获得在两个样本中存在显著性差异的分类单元；
16)获取各样本在指定分类水平上的组成和相对丰度分布表，对各分类单元在多组样本中的丰度分布差异进行ANOVA方差分析/Kruskal-Wallis H检验，从而获得在多组样本中存在显著性差异的分类单元；
18)对前述OTU代表序列，通过PyNAST和MAFFT工具进行多序列比对，之后通过FastTree工具构建OTU代表序列的系统发育树，该文件以Newick格式保存；
19)根据前述OTU丰度矩阵和OTU划分和分类地位鉴定结果，将每个样本所含有的OTU的丰度信息和分类学组成数据映射至NCBI Taxonomy所提供的微生物分类等级树，统一呈现所有样本在各分类水平的具体组成；
20)获取各样本在指定分类水平上的组成和相对丰度分布表，对样本总体在各分类水平的组成构建等级树，同时以不同颜色区分各分类单元，并通过节点大小反映它们的丰度分布；
21)获取各样本在指定分类水平上的组成和丰度分布表，通过Krona软件进行群落分类学组成的交互展示；
22)根据前述OTU丰度矩阵和OTU划分和分类地位鉴定结果，构建交互式OTU热图；
23)获取各样本在指定分类水平上的组成和相对丰度分布表，对丰度前50位的分类单元进行聚类分析并绘制热图；
24)获取各样本在指定分类水平上的组成和相对丰度分布表，对指定分类水平的群落组成结构进行PCA主成分分析，并且以二维和三维图像描述样本间的自然分布特征；
25)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行PcoA主坐标分析，并且以二维和三维图像描述样本间基于微生物系统发育关系的群落空间分布特征；
26)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行NMDS非度量多维尺度分析，通过二维或三维排序图描述群落样本的结构分布；
27)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，对加权及非加权距离矩阵分别进行UPGMA聚类分析；
28)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，对加权及非加权距离矩阵分别进行组内和组间距离均值进行t检验，并通过1000次蒙特卡罗置换检验判断统计显著性，从而衡量组内和组间距离差异；
30)获取各样本在指定分类水平上的组成和相对丰度分布表，通过三元相图分析显示出不同物种在各样本中的相对丰度差异；
31)获取各样本在指定分类水平上的组成和相对丰度分布表，对指定分类水平的相对丰度矩阵进行RDA冗余分析，通过1000次置换检验确定统计显著性，并生成包含“样本—分类单元—影响因素”三种元素排序图；
32)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行CAP主坐标典型相关分析分析，并作1000次置换检验，确定组间差异是否具有统计学显著性，使用CAP分析得到的前两维数据生成CAP约束排序图；
33)获取各样本在指定分类水平上的组成和相对丰度分布表和样本分组数据构建PLS-DA判别模型，并根据分析得到的前两维数据生成PLS-DA约束排序图；
34)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行Adonis/PERMANOVA分析，并作999次置换检验确定组间差异是否具有统计学显著性；
35)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行ANOSIM分析，通过对样本距离等级排序来判断样本组内和组间差异的大小，并通过置换检验评价原始样本组间差异的统计学显著性；
36)根据前述对OTU丰度矩阵中的全体样本根据最低测序深度统一进行随机重抽样后，对重抽样OTU丰度矩阵使用随机森林算法挑取丰度分布在不同组间存在显著差异的OTU；挑取OTU时使用1000棵随机森林决策树进行建模，并以10倍交叉验证估计“基线”误差的大小随机森林分析；
37)获取各样本在指定分类水平上的组成和相对丰度分布表，计算丰度位于前50位的分类单元之间的Spearman等级相关系数，对其中|rho|>0.8且P值<0.01的相关优势属构建关联网络，并进行可视化；
38)根据前述对OTU丰度矩阵中的全体样本根据最低测序深度统一进行随机重抽样后，去除绝对丰度<2的OTU，进行二分网络分析，并进行可视化，样本和OTU根据弹簧镶嵌模型排布；
39)对前述获得的高质量序列，使用PICRUSt软件，根据KEGG数据库中微生物代谢功能的类别对群落样本进行预测，并通过条形图展示各样本中编码的功能基因类别及其丰度。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201611187576.9A CN106815492B (zh) | 2016-12-20 | 2016-12-20 | 一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201611187576.9A CN106815492B (zh) | 2016-12-20 | 2016-12-20 | 一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106815492A CN106815492A (zh) | 2017-06-09 |
| CN106815492B true CN106815492B (zh) | 2019-02-12 |
Family
ID=59109080
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201611187576.9A Active CN106815492B (zh) | 2016-12-20 | 2016-12-20 | 一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106815492B (zh) |
Families Citing this family (15)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107463800B (zh) * | 2017-07-19 | 2018-05-11 | 东莞博奥木华基因科技有限公司 | 一种肠道微生物信息分析方法及系统 |
| CN107475385A (zh) * | 2017-08-21 | 2017-12-15 | 上海派森诺生物科技股份有限公司 | 一种基于smrt高通量测序技术的菌群多样性组成谱数据分析方法 |
| CN110111843B (zh) * | 2018-01-05 | 2021-07-06 | 深圳华大基因科技服务有限公司 | 对核酸序列进行聚类的方法、设备及存储介质 |
| CN109686406A (zh) * | 2018-11-12 | 2019-04-26 | 山东省医学科学院基础医学研究所 | 一种系统发生树图制作方法及系统 |
| CN109727644B (zh) * | 2018-11-12 | 2021-09-07 | 山东省医学科学院基础医学研究所 | 基于微生物基因组二代测序数据的Venn图制作方法及系统 |
| CN110504006B (zh) * | 2019-07-15 | 2023-03-28 | 广州奇辉生物科技有限公司 | 一种处理扩增子数据的方法、系统、平台及存储介质 |
| CN111709219B (zh) * | 2020-04-28 | 2021-06-18 | 上海欧易生物医学科技有限公司 | 单组学及多组学KEGG PATHWAY map表达热图个性化展示的方法及应用 |
| CN112489726A (zh) * | 2020-11-10 | 2021-03-12 | 哈尔滨因极科技有限公司 | 基于16s微生物扩增测序数据的分析方法、装置及设备 |
| CN112365928B (zh) * | 2020-11-16 | 2021-07-06 | 赛福解码(北京)基因科技有限公司 | 生物信息数据分析和结果质控自动化方法和系统 |
| CN112466391A (zh) * | 2020-12-24 | 2021-03-09 | 广州基迪奥生物科技有限公司 | 一种16s测序和its测序数据关联分析方法及系统 |
| CN112908407B (zh) * | 2021-02-02 | 2021-10-01 | 北京大学 | 一种用tRNA组学来质控蛋白生物合成体系的方法 |
| CN113257362A (zh) * | 2021-05-24 | 2021-08-13 | 自然资源部第三海洋研究所 | 一种生物环境样品的筛选方法 |
| CN113793640B (zh) * | 2021-09-17 | 2024-03-08 | 艾德范思(北京)医学检验实验室有限公司 | 基于二代测序的微生物16s扩增子数据分析的装置及方法 |
| CN114703265A (zh) * | 2022-03-28 | 2022-07-05 | 南京农业大学 | 一种基于16SrRNA扩增子测序检测土壤致病细菌生物污染的方法 |
| CN114938947B (zh) * | 2022-07-25 | 2022-10-21 | 天津医科大学总医院 | 一种多层级的阴道微生态评估体系构建方法和装置 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102382877A (zh) * | 2010-08-30 | 2012-03-21 | 中国食品发酵工业研究院 | 一种研究大曲细菌群落结构多样性的方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7244607B2 (en) * | 2003-10-01 | 2007-07-17 | The United States Of America, As Represented By The Secretary Of Agriculture | Chromobacterium subtsugae sp. nov. for control of insect pests |
| US20120129706A1 (en) * | 2010-11-22 | 2012-05-24 | Ashvini Chauhan | Method of Assessing Soil Quality and Health |
-
2016
- 2016-12-20 CN CN201611187576.9A patent/CN106815492B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102382877A (zh) * | 2010-08-30 | 2012-03-21 | 中国食品发酵工业研究院 | 一种研究大曲细菌群落结构多样性的方法 |
Non-Patent Citations (2)
| Title |
|---|
| 《Oligotyping: differentiating between closely related microbial taxa using 16S rRNA gene data》;A. Murat Eren 等;《Methods in Ecology and Evolution》;20131231;第4卷(第12期);全文 |
| 《崇明东滩夏冬季表层沉积物细菌多样性研究》;郑艳玲 等;《中国环境科学》;20120131;第32卷(第2期);全文 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106815492A (zh) | 2017-06-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106815492B (zh) | 一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法 | |
| Czech et al. | Scalable methods for analyzing and visualizing phylogenetic placement of metagenomic samples | |
| Jachner et al. | Statistical methods for the qualitative assessment of dynamic models with time delay (R Package qualV) | |
| Birks | Overview of numerical methods in palaeolimnology | |
| CN103336916B (zh) | 一种测序序列映射方法及系统 | |
| Abdelaal et al. | Predicting cell populations in single cell mass cytometry data | |
| CN107368700A (zh) | 基于计算云平台的微生物多样性交互分析系统及其方法 | |
| Gustafsson et al. | Constructing and analyzing a large-scale gene-to-gene regulatory network Lasso-constrained inference and biological validation | |
| CN108985380B (zh) | 一种基于聚类集成的转辙机故障识别方法 | |
| Lancashire et al. | Classification of bacterial species from proteomic data using combinatorial approaches incorporating artificial neural networks, cluster analysis and principal components analysis | |
| Bhargava et al. | DNA barcoding in plants: evolution and applications of in silico approaches and resources | |
| CN105279391A (zh) | 一种宏基因组16S rRNA的高通量测序数据处理及分析流程控制方法 | |
| CN112669899B (zh) | 一种16s和宏基因组测序数据关联分析方法、系统及设备 | |
| CN108681742B (zh) | 用于分析司机驾驶行为对车辆能耗敏感性的分析方法 | |
| CN111710364B (zh) | 一种菌群标记物的获取方法、装置、终端及存储介质 | |
| CN102073882A (zh) | 高光谱遥感影像的dna计算光谱匹配分类方法 | |
| Liang et al. | MetaVelvet-DL: a MetaVelvet deep learning extension for de novo metagenome assembly | |
| CN114864004A (zh) | 基于滑窗稀疏卷积去噪自编码器的缺失标记填充方法 | |
| CN109448842A (zh) | 人体肠道微生态失衡的确定方法、装置及电子设备 | |
| CN113793640B (zh) | 基于二代测序的微生物16s扩增子数据分析的装置及方法 | |
| Peltonen et al. | Parallel Coordinate Plots for Neighbor Retrieval. | |
| Gorman et al. | Sparsification of large ultrametric matrices: insights into the microbial Tree of Life | |
| Abouabdallah et al. | Does clustering of DNA barcodes agree with botanical classification directly at high taxonomic levels? Trees in French Guiana as a case study | |
| Wagner et al. | On the effect of response transformations in sequential parameter optimization | |
| Boctor et al. | Comprehensive Guideline for Microbiome Analysis Using R |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| EE01 | Entry into force of recordation of patent licensing contract |
Application publication date: 20170609 Assignee: Shanghai Puchuang Longke Finance Leasing Co., Ltd Assignor: SHANGHAI PERSONAL BIOTECHNOLOGY CO., LTD. Contract record no.: X2019310000027 Denomination of invention: Automatic method for analyzing constitution and diversity of bacterial community of 16SrRNA gene Granted publication date: 20190212 License type: Exclusive License Record date: 20191225 |
|
| EE01 | Entry into force of recordation of patent licensing contract | ||
| PE01 | Entry into force of the registration of the contract for pledge of patent right | ||
| PE01 | Entry into force of the registration of the contract for pledge of patent right |
Denomination of invention: Automatic method for analyzing constitution and diversity of bacterial community of 16SrRNA gene Effective date of registration: 20191227 Granted publication date: 20190212 Pledgee: Shanghai Puchuang Longke Finance Leasing Co., Ltd Pledgor: SHANGHAI PERSONAL BIOTECHNOLOGY CO., LTD. Registration number: Y2019310000041 |
|
| PC01 | Cancellation of the registration of the contract for pledge of patent right | ||
| PC01 | Cancellation of the registration of the contract for pledge of patent right |
Date of cancellation: 20201228 Granted publication date: 20190212 Pledgee: Shanghai Puchuang Longke Finance Leasing Co.,Ltd. Pledgor: SHANGHAI PERSONAL BIOTECHNOLOGY Co.,Ltd. Registration number: Y2019310000041 |
|
| EC01 | Cancellation of recordation of patent licensing contract | ||
| EC01 | Cancellation of recordation of patent licensing contract |
Assignee: Shanghai Puchuang Longke Finance Leasing Co.,Ltd. Assignor: SHANGHAI PERSONAL BIOTECHNOLOGY Co.,Ltd. Contract record no.: X2019310000027 Date of cancellation: 20210108 |