CN106815492A - 一种用于16SrRNA基因的细菌群落组成和多样性分析的自动化方法 - Google Patents

一种用于16SrRNA基因的细菌群落组成和多样性分析的自动化方法 Download PDF

Info

Publication number
CN106815492A
CN106815492A CN201611187576.9A CN201611187576A CN106815492A CN 106815492 A CN106815492 A CN 106815492A CN 201611187576 A CN201611187576 A CN 201611187576A CN 106815492 A CN106815492 A CN 106815492A
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.)
Granted
Application number
CN201611187576.9A
Other languages
English (en)
Other versions
CN106815492B (zh
Inventor
薛正晟
寇文伯
王慧娟
姜丽荣
孙子奎
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
SHANGHAI PERSONAL BIOTECHNOLOGY CO Ltd
Original Assignee
SHANGHAI PERSONAL BIOTECHNOLOGY CO Ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by SHANGHAI PERSONAL BIOTECHNOLOGY CO Ltd filed Critical SHANGHAI PERSONAL BIOTECHNOLOGY CO Ltd
Priority to CN201611187576.9A priority Critical patent/CN106815492B/zh
Publication of CN106815492A publication Critical patent/CN106815492A/zh
Application granted granted Critical
Publication of CN106815492B publication Critical patent/CN106815492B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16BBIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00ICT 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)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
  • Information Retrieval, Db Structures And Fs Structures Therefor (AREA)

Abstract

本发明公开的一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法,其提供的16S rRNA测序数据分析流程以测序原始序列数据作为输入,调用业界标准的分析工具(如:Mothur、QIIME等),最终对数据进行可视化,并得到易于解读的分析结果。本发明包含了目前流行的主流分析项目,同时分析内容实现模块化,数据挖掘分析的方法更多样、更深入,可以根据不同的需要结合不同的分析模块内容,先后顺序的流程安排也更合理;此外,消除了测序深度不一导致的分析误差,使分析结果更全面、准确、可靠。

Description

一种用于16SrRNA基因的细菌群落组成和多样性分析的自动 化方法
技术领域:
本发明一般有关于分子生物学技术领域,特别涉及高通量测序数据分析技术领域,并且更具体地说,涉及一种用于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数据库中微生物代谢功能的类别对群落样本进行预测,并通过条形图展示各样本中编码的功能基因类别及其丰度。
CN201611187576.9A 2016-12-20 2016-12-20 一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法 Active CN106815492B (zh)

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 true CN106815492A (zh) 2017-06-09
CN106815492B 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)

Cited By (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107463800A (zh) * 2017-07-19 2017-12-12 东莞博奥木华基因科技有限公司 一种肠道微生物信息分析方法及系统
CN107475385A (zh) * 2017-08-21 2017-12-15 上海派森诺生物科技股份有限公司 一种基于smrt高通量测序技术的菌群多样性组成谱数据分析方法
CN109686406A (zh) * 2018-11-12 2019-04-26 山东省医学科学院基础医学研究所 一种系统发生树图制作方法及系统
CN109727644A (zh) * 2018-11-12 2019-05-07 山东省医学科学院基础医学研究所 基于微生物基因组二代测序数据的Venn图制作方法及系统
CN110111843A (zh) * 2018-01-05 2019-08-09 深圳华大基因科技服务有限公司 对核酸序列进行聚类的方法、设备及存储介质
CN110504006A (zh) * 2019-07-15 2019-11-26 广州奇辉生物科技有限公司 一种处理扩增子数据的方法、系统、平台及存储介质
CN111709219A (zh) * 2020-04-28 2020-09-25 上海欧易生物医学科技有限公司 单组学及多组学KEGG PATHWAY map表达热图个性化展示的方法及应用
CN112365928A (zh) * 2020-11-16 2021-02-12 赛福解码(北京)基因科技有限公司 生物信息数据分析和结果质控自动化方法和系统
CN112435713A (zh) * 2020-12-31 2021-03-02 申友基因组研究院(南京)有限公司 一种肠道微生物16SrRNA的NGS数据分析方法
CN112466391A (zh) * 2020-12-24 2021-03-09 广州基迪奥生物科技有限公司 一种16s测序和its测序数据关联分析方法及系统
CN112489726A (zh) * 2020-11-10 2021-03-12 哈尔滨因极科技有限公司 基于16s微生物扩增测序数据的分析方法、装置及设备
CN112908407A (zh) * 2021-02-02 2021-06-04 北京大学 一种用tRNA组学来质控蛋白生物合成体系的方法
CN113257362A (zh) * 2021-05-24 2021-08-13 自然资源部第三海洋研究所 一种生物环境样品的筛选方法
CN113793640A (zh) * 2021-09-17 2021-12-14 艾德范思(北京)医学检验实验室有限公司 基于二代测序的微生物16s扩增子数据分析的装置及方法
CN114703265A (zh) * 2022-03-28 2022-07-05 南京农业大学 一种基于16SrRNA扩增子测序检测土壤致病细菌生物污染的方法
CN114938947A (zh) * 2022-07-25 2022-08-26 天津医科大学总医院 一种多层级的阴道微生态评估体系构建方法和装置

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20050074431A1 (en) * 2003-10-01 2005-04-07 Martin Phyllis A.W. Chromobacterium suttsuga sp. nov. and use for control of insect pests
CN102382877A (zh) * 2010-08-30 2012-03-21 中国食品发酵工业研究院 一种研究大曲细菌群落结构多样性的方法
US20120129706A1 (en) * 2010-11-22 2012-05-24 Ashvini Chauhan Method of Assessing Soil Quality and Health

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20050074431A1 (en) * 2003-10-01 2005-04-07 Martin Phyllis A.W. Chromobacterium suttsuga sp. nov. and use for control of insect pests
CN102382877A (zh) * 2010-08-30 2012-03-21 中国食品发酵工业研究院 一种研究大曲细菌群落结构多样性的方法
US20120129706A1 (en) * 2010-11-22 2012-05-24 Ashvini Chauhan Method of Assessing Soil Quality and Health

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
A. MURAT EREN 等: "《Oligotyping: differentiating between closely related microbial taxa using 16S rRNA gene data》", 《METHODS IN ECOLOGY AND EVOLUTION》 *
郑艳玲 等: "《崇明东滩夏冬季表层沉积物细菌多样性研究》", 《中国环境科学》 *

Cited By (20)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107463800B (zh) * 2017-07-19 2018-05-11 东莞博奥木华基因科技有限公司 一种肠道微生物信息分析方法及系统
CN107463800A (zh) * 2017-07-19 2017-12-12 东莞博奥木华基因科技有限公司 一种肠道微生物信息分析方法及系统
CN107475385A (zh) * 2017-08-21 2017-12-15 上海派森诺生物科技股份有限公司 一种基于smrt高通量测序技术的菌群多样性组成谱数据分析方法
CN110111843B (zh) * 2018-01-05 2021-07-06 深圳华大基因科技服务有限公司 对核酸序列进行聚类的方法、设备及存储介质
CN110111843A (zh) * 2018-01-05 2019-08-09 深圳华大基因科技服务有限公司 对核酸序列进行聚类的方法、设备及存储介质
CN109686406A (zh) * 2018-11-12 2019-04-26 山东省医学科学院基础医学研究所 一种系统发生树图制作方法及系统
CN109727644A (zh) * 2018-11-12 2019-05-07 山东省医学科学院基础医学研究所 基于微生物基因组二代测序数据的Venn图制作方法及系统
CN110504006A (zh) * 2019-07-15 2019-11-26 广州奇辉生物科技有限公司 一种处理扩增子数据的方法、系统、平台及存储介质
CN111709219A (zh) * 2020-04-28 2020-09-25 上海欧易生物医学科技有限公司 单组学及多组学KEGG PATHWAY map表达热图个性化展示的方法及应用
CN112489726A (zh) * 2020-11-10 2021-03-12 哈尔滨因极科技有限公司 基于16s微生物扩增测序数据的分析方法、装置及设备
CN112365928A (zh) * 2020-11-16 2021-02-12 赛福解码(北京)基因科技有限公司 生物信息数据分析和结果质控自动化方法和系统
CN112365928B (zh) * 2020-11-16 2021-07-06 赛福解码(北京)基因科技有限公司 生物信息数据分析和结果质控自动化方法和系统
CN112466391A (zh) * 2020-12-24 2021-03-09 广州基迪奥生物科技有限公司 一种16s测序和its测序数据关联分析方法及系统
CN112435713A (zh) * 2020-12-31 2021-03-02 申友基因组研究院(南京)有限公司 一种肠道微生物16SrRNA的NGS数据分析方法
CN112908407A (zh) * 2021-02-02 2021-06-04 北京大学 一种用tRNA组学来质控蛋白生物合成体系的方法
CN113257362A (zh) * 2021-05-24 2021-08-13 自然资源部第三海洋研究所 一种生物环境样品的筛选方法
CN113793640A (zh) * 2021-09-17 2021-12-14 艾德范思(北京)医学检验实验室有限公司 基于二代测序的微生物16s扩增子数据分析的装置及方法
CN113793640B (zh) * 2021-09-17 2024-03-08 艾德范思(北京)医学检验实验室有限公司 基于二代测序的微生物16s扩增子数据分析的装置及方法
CN114703265A (zh) * 2022-03-28 2022-07-05 南京农业大学 一种基于16SrRNA扩增子测序检测土壤致病细菌生物污染的方法
CN114938947A (zh) * 2022-07-25 2022-08-26 天津医科大学总医院 一种多层级的阴道微生态评估体系构建方法和装置

Also Published As

Publication number Publication date
CN106815492B (zh) 2019-02-12

Similar Documents

Publication Publication Date Title
CN106815492B (zh) 一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法
Kassambara Practical guide to cluster analysis in R: Unsupervised machine learning
Jachner et al. Statistical methods for the qualitative assessment of dynamic models with time delay (R Package qualV)
CN107368700A (zh) 基于计算云平台的微生物多样性交互分析系统及其方法
CN108985380B (zh) 一种基于聚类集成的转辙机故障识别方法
CN112669899B (zh) 一种16s和宏基因组测序数据关联分析方法、系统及设备
CN108681742B (zh) 用于分析司机驾驶行为对车辆能耗敏感性的分析方法
CN114864004A (zh) 基于滑窗稀疏卷积去噪自编码器的缺失标记填充方法
CN116072227A (zh) 海洋营养成分生物合成途径挖掘方法、装置、设备和介质
CN111210085B (zh) 一种基于多视图集成学习的煤矿瓦斯浓度预警方法
CN117036781A (zh) 一种基于树综合多样性深度森林的图像分类方法
CN115034812B (zh) 基于大数据的钢铁行业销售量预测方法及装置
CN113793640B (zh) 基于二代测序的微生物16s扩增子数据分析的装置及方法
Peltonen et al. Parallel Coordinate Plots for Neighbor Retrieval.
US20220336057A1 (en) Efficient voxelization for deep learning
US11538555B1 (en) Protein structure-based protein language models
CN113870950B (zh) 一种稻瘟菌侵染水稻关键sRNA识别系统及识别方法
US20220336054A1 (en) Deep Convolutional Neural Networks to Predict Variant Pathogenicity using Three-Dimensional (3D) Protein Structures
EP4323989A1 (en) Efficient voxelization for deep learning
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
Gao et al. GBDT4CTRVis: visual analytics of gradient boosting decision tree for advertisement click-through rate prediction
CN117171676B (zh) 基于决策树的土壤微生物识别分析方法、系统及存储介质
CN117875726B (zh) 基于深度学习的价值链优化管控方法
Lee et al. Finding the optimal gene order in displaying microarray data
Langfelder et al. Package ‘WGCNA’

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

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

PE01 Entry into force 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

PC01 Cancellation of the registration of the contract for pledge of patent right
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

EC01 Cancellation of recordation of patent licensing contract