一种用于16SrRNA基因的细菌群落组成和多样性分析的自动
化方法
技术领域：
[0001] 本
发明一般有关于分子
生物学技术领域，特别涉及高通量测序数据分析技术领域，并且更具体地说，涉及一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法。背景技术：
[0002] 新一代高通量测序技术大幅度降低了测序的时间和成本，使得大规模测序逐渐成为常规的研究和检测手段，测序产生的数据量急剧增加。如何高效地分析这些数据，已成为迫切需要解决的问题。
[0003] 目前高通量测序数据分析工具很多，进行分析序列信息的生物信息学工具纷繁复杂，对于分析菌群微生态的大规模测序数据，也已开发产生多种成熟的分析工具可供使用。其中，针对于菌群多样性和组成谱检测这一广泛应用的研究策略，绝大多数基于
微生物核糖体RNA(rRNA)基因的序列分析工具并不能一次性满足研究人员的大部分分析需求，往往需要多次调用不同的分析工具(甚至不同的计算机平台系统)以完成所有相关分析。
[0004] 当前用于rRNA基因序列分析的工具分为两大类：基于Web的工具和本地安装的工具。基于Web的工具托管在
服务器上，为用户提供Web界面。例如Ribosomal Database Project(RDP)的分类鉴定工具等；研究人员可以通过互联网上传测得的序列数据，并使用Web界面配置可选参数进行分析。但是对于一次上传的数据量存在限制，并受网络环境所约束，不适合大规模数据分析，并且远程用户无法根据自己的需求来自定义后端分析
软件，并不方便。本地化的16S rRNA基因分析工具包括Mothur和QIIME等。使用这些本地工具时，不需要将数据上传到远程服务器，只需在Linux服务器/工作站合理配置安装即能投入使用。这些工具大多是开源形式，允许研究人员根据需要自定义软件。然而，很多情况下，根据不同的分析需求，实现一个完整的分析流程仍需要对众多工具进行整合。如何能正确高效地选择并整合这些工具已成为迫切需求。
[0005] 现有的16S rRNA基因的细菌群落组成和多样性分析流程中，分析内容较为
基础和简单，未包含目前流行的多项主流分析项目：Specaccum物种累积曲线、LEfSe分析、Wilcoxon秩和检验等。在一些情况下，无法满足研究人员的分析需求。此外，原有分析流程在进行后续PCA、PCoA等比较分析时，未对数据进行拉平处理，由此将会引入测序深度不一导致的分析误差。
发明内容
[0006] 本发明的一个目的是为了克服现有工具中固有的弱点并结合其优势，提供一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法，该方法具有可靠的流程，已确保分析结果的准确。并满足多种分析需求。
[0007] 本发明的另一个目的是为了提供一种16S rRNA测序数据分析流程，实现各个环节的高效自动化管理和分析，从而节省时间成本，减轻研究人员的数据分析负担。
[0008] 为了实现上述目的，本发明所采用的技术方案如下：
[0009] 一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法，其提供的16S rRNA测序数据分析流程以测序原始序列数据作为输入,调用业界标准的分析工具(如：
Mothur、QIIME等)，最终对数据进行
可视化，并得到易于解读的分析结果，具体包括以下步骤：
[0010] 1)通过原始序列的测序
质量值、模糊
碱基数目、序列长度、引物序列和barcode序列的匹配度信息，对原始序列进行过滤和质量控制，并检查和剔除嵌合体，获得高质量序列；
[0011] 2)对步骤1)获得的高质量序列的长度分布进行统计；
[0012] 3)对步骤1)获得的高质量序列按97％的序列相似度进行归并和OTU划分，并选取每个OTU中丰度最高的序列作为该OTU的代表序列，随后，根据每个OTU在每个样本中所包含的序列数，构建OTU在各样本中丰度的矩阵文件；
[0013] 4)通过将OTU代表序列与对应
数据库的模板序列相比对，获取每个OTU所对应的分类学信息；
[0014] 5)将丰度值低于全体样本测序总量0.001％(十万分之一)的OTU去除，并将去除了稀有OTU的此丰度矩阵用于后续的一系列分析；
[0015] 6)根据获得的OTU丰度矩阵，计算各样本组共有OTU的数量，并通过Venn图直观地呈现各样本组所共有和独有OTU所占的比例；
[0016] 7)对OTU丰度矩阵中每个样本的序列总数在不同测序深度下依次随机抽样，以每个深度下
抽取到的序列数及其对应的OTU数绘制稀疏曲线；
[0017] 8)对OTU丰度矩阵中每个样本所对应的OTU总数绘制Specaccum物种累积曲线；
[0018] 9)对OTU及其对应的丰度值经Log2对数转换绘制各样本的丰度等级曲线；
[0019] 对OTU丰度矩阵中的全体样本根据最低测序深度统一进行随机重抽样(即序列拉平处理)，随后，分别对每个样本计算四种多样性指数；
[0020] 10)根据OTU划分和分类地位鉴定结果，可以获得每个样本在各分类
水平(界/
门/纲/目/科/属/种等)的具体组成；
[0021] 11)获取各样本在
指定分类水平上的组成和丰度分布表，并通过饼图、柱状图或面积图呈现分析结果，根据研究对象是单个或多个群落样本，绘图结果可能会以不同方式进行展示；
[0022] 12)获取各样本在指定分类水平上的组成和绝对丰度分布表，调用Metastats的统计学
算法，对指定分类水平的各个分类单元在样本组之间的序列量即绝对丰度差异进行两两比较检验；
[0023] 13)获取各样本在指定分类水平上的组成和相对丰度分布表，进行LEfSe分析，筛选关键的生物标记物；
[0024] 14)获取各样本在指定分类水平上的组成和相对丰度分布表，对各分类单元在两组样本中的丰度分布差异进行Wilcoxon秩和检验或Welch’s t检验，从而获得在两组中存在显著性差异的分类单元；
[0025] 15)获取各样本在指定分类水平上的组成和相对丰度分布表，对各分类单元在两个样本中的丰度分布差异进行Fisher’s检验，从而获得在两个样本中存在显著性差异的分类单元；
[0026] 16)获取各样本在指定分类水平上的组成和相对丰度分布表，对各分类单元在多组样本中的丰度分布差异进行ANOVA方差分析/Kruskal-Wallis H检验，从而获得在多组样本中存在显著性差异的分类单元；
[0027] 18)对前述OTU代表序列，通过PyNAST和MAFFT等工具进行多序列比对，之后通过FastTree工具构建OTU代表序列的系统发育树，该文件以Newick格式保存；
[0028] 19)根据前述OTU丰度矩阵和OTU划分和分类地位鉴定结果，将每个样本所含有的OTU的丰度信息和分类学组成数据映射至NCBI Taxonomy所提供的微生物分类等级树，统一呈现所有样本在各分类水平的具体组成；
[0029] 20)获取各样本在指定分类水平上的组成和相对丰度分布表，对样本总体在各分类水平的组成构建等级树，同时以不同
颜色区分各分类单元，并通过
节点大小反映它们的丰度分布；
[0030] 21)获取各样本在指定分类水平上的组成和丰度分布表，通过Krona软件进行群落分类学组成的交互展示；
[0031] 22)根据前述OTU丰度矩阵和OTU划分和分类地位鉴定结果，构建交互式OTU热图；
[0032] 23)获取各样本在指定分类水平上的组成和相对丰度分布表，对丰度前50位的分类单元进行
聚类分析并绘制热图；
[0033] 24)获取各样本在指定分类水平上的组成和相对丰度分布表，对指定分类水平的群落组成结构进行PCA主成分分析，并且以二维和三维图像描述样本间的自然分布特征；
[0034] 25)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行PcoA主坐标分析，并且以二维和三维图像描述样本间基于微生物系统发育关系的群落空间分布特征；
[0035] 26)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行NMDS非度量多维尺度分析，通过二维或三维排序图描述群落样本的结构分布；
[0036] 27)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，对加权及非加权距离矩阵分别进行UPGMA聚类分析；
[0037] 28)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，对加权及非加权距离矩阵分别进行组内和组间距离均值进行t检验，并通过1000次蒙特卡罗置换检验判断统计显著性，从而衡量组内和组间距离差异；
[0038] 30)获取各样本在指定分类水平上的组成和相对丰度分布表，通过三元
相图分析显示出不同物种在各样本中的相对丰度差异；
[0039] 31)获取各样本在指定分类水平上的组成和相对丰度分布表，对指定分类水平的相对丰度矩阵进行RDA冗余分析，通过1000次置换检验确定统计显著性，并生成包含“样本—分类单元—影响因素”三种元素排序图；
[0040] 32)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行CAP主坐标典型相关分析分析，并作1000次置换检验，确定组间差异是否具有统计学显著性，使用CAP分析得到的前两维数据生成CAP约束排序图；
[0041] 33)获取各样本在指定分类水平上的组成和相对丰度分布表和样本分组数据构建PLS-DA判别模型，并根据分析得到的前两维数据生成PLS-DA约束排序图；
[0042] 34)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行Adonis/PERMANOVA分析，并作999次置换检验确定组间差异是否具有统计学显著性；
[0043] 35)根据前述OTU丰度矩阵和OTU代表序列的系统发育树，基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行ANOSIM分析，通过对样本距离等级排序来判断样本组内和组间差异的大小，并通过置换检验评价原始样本组间差异的统计学显著性；
[0044] 36)根据前述对OTU丰度矩阵中的全体样本根据最低测序深度统一进行随机重抽样后，对重抽样OTU丰度矩阵使用
随机森林算法挑取丰度分布在不同组间存在显著差异的OTU。挑取OTU时使用1000棵随机森林
决策树进行建模，并以10倍交叉验证估计“基线”误差的大小随机森林分析；
[0045] 37)获取各样本在指定分类水平上的组成和相对丰度分布表，计算丰度位于前50位的分类单元之间的Spearman等级相关系数，对其中|rho|>0.8且P值<0.01的相关优势属构建关联网络，并进行可视化。
[0046] 38)根据前述对OTU丰度矩阵中的全体样本根据最低测序深度统一进行随机重抽样后，去除绝对丰度<2的OTU，进行二分网络分析，并进行可视化，样本和OTU根据
弹簧镶嵌模型排布。
[0047] 39)对前述获得的高质量序列，使用PICRUSt软件，根据KEGG数据库中微生物代谢功能的类别对群落样本进行预测，并通过条形图展示各样本中编码的功能基因类别及其丰度。
[0048] 本发明的有益效果是：
[0049] 改进后的流程包含了目前流行的主流分析项目，同时分析内容实现模
块化，
数据挖掘分析的方法更多样、更深入，可以根据不同的需要结合不同的分析模块内容，先后顺序的流程安排也更合理；此外，消除了测序深度不一导致的分析误差，使分析结果更全面、准确、可靠。
[0050] 改进后的流程通过自动化脚本控制运行，一站式输出可视化图形和分析结果,简化了数据分析过程中的人工操作,提高了运行效率。研究人员只需导入数据及根据分析需求简单调整脚本即可完成操作，从而节省时间，减轻分析的工作负担。
附图说明
[0051] 图1是根据本发明构造的的一个细菌16S rRNA测序数据分析
流程图。
具体实施方式
[0052] 在具体实施方式中，该方法如图1所示包括以下步骤：
[0053] 步骤1：对于Illumina Miseq平台原始双端测序数据，以原始下机数据作为输入数据执行MiSeqQuality16S.pl脚本，窗口大小为10bp，步长为1bp，从5’端第一个碱基
位置开始移动，要求窗口中碱基平均质量≥Q20(即碱基平均测序准确率≥99％)，从第一个平均质量值低于Q20的窗口处截断序列，并要求截断后的序列长度≥150bp，且不允许存在模糊碱基N。随后，利用FLASH软件，对通过质量初筛的双端序列根据重叠碱基进行
配对连接：要求Read 1和Read 2两条序列的重叠碱基长度≥10bp，且不允许碱基错配。最后，以FLASH拼接输出序列和barcode序列作为输入文件，将执行MiSeqRawDataSplit.pl脚本，把连接后的序列识别分配入对应样本(要求Index序列完全匹配)，从而获得每个样本的有效序列。
[0054] 对于Roche 454FLX+平台原始测序数据，以原始下机数据和barcode序列作为输入文件，以原始下机数据和barcode序列作为输入数据通过sfffile工具把序列识别分配入对应样本(要求Index序列完全匹配)，以sfffile输出序列作为输入文件通过sffinfo工具输出FASTA格式文件和QUAL格式文件，随后通过QIIME进行过滤，要求碱基平均质量≥Q25，不允许存在模糊碱基，从而获得每个样本的有效序列。
[0055] 以上述有效序列和样本分组信息作为输入文件，执行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代表序列的系统发育树。
[0056] 步骤2：通过BIOM文件处理工具和modify_otu_table.R脚本，以步骤1输出OTU丰度矩阵作为输入文件将丰度值低于全体样本测序总量0.001％(十万分之一)的OTU去除。
[0057] 步骤:3：以步骤2输出OTU丰度矩阵和分组信息作为输入文件，通过creat_venn_map.py、plot_venn_graph.R脚本计算各样本(组)共有OTU的数量，并进行可视化。
[0058] 步骤4：以步骤2输出OTU丰度矩阵作为输入文件，通过create_table2_4_in_report.R、bar.R脚本，计算各样本中能分类至门、纲、目、科、属、种各分类水平的OTU数及可视化。
[0059] 步骤5：以步骤2输出OTU丰度矩阵作为输入文件，使用BIOM工具和QIIME软件绘制稀疏曲线。
[0060] 步骤6：以步骤2输出OTU丰度矩阵作为输入文件，通过species_curve.R脚本绘制Specaccum物种累积曲线。
[0061] 步骤7：以步骤2输出OTU丰度矩阵作为输入文件，通过plot_rank_abundance_graph.R脚本绘制各样本的丰度等级曲线。
[0062] 步骤8：以步骤2输出OTU丰度矩阵作为输入文件，通过normalize_rarefied_otu_table.sh脚本对全体样本根据最低测序深度统一进行随机重抽样100次，并计算平均值，剔除平均值<1的OTU，从而最大限度消除测序深度不一导致的分析误差并保证结果的准确性。
[0063] 步骤9：以步骤8输出OTU丰度矩阵作为输入文件，通过QIIME软件分别对每个样本计算四种多样性指数。
[0064] 步骤10：以步骤8输出OTU丰度矩阵作为输入文件，通过bar.R脚本对门、纲、目、科、属、种六个分类水平各自含有的微生物类群数可视化。
[0065] 步骤11：以步骤2输出OTU丰度矩阵作为输入文件，通过QIIME软件和bar-phylum.R脚本获取各样本在门、纲、目、科、属五个分类水平上的组成和丰度分布表及可视化。
[0066] 步骤12：以步骤11输出OTU丰度矩阵作为输入文件，通过metats.sh、diff_group2.R脚本，调用Metastats的统计学算法，对指定分类水平的各个分类单元在样本(组)之间的序列量(即绝对丰度)差异进行两两比较检验。
[0067] 步骤13：以步骤11输出OTU丰度矩阵作为输入文件，通过LEfSe本地工具进行LEfSe分析，筛选关键的生物标记物。
[0068] 步骤14：以步骤11输出OTU丰度矩阵和分组信息作为输入文件，通过wilcox_box.R脚本对各分类单元在两组样本中的丰度分布差异进行Wilcoxon秩和检验，从而获得在两组中存在显著性差异的分类单元并输出可视化结果。
[0069] 步骤15：以步骤11输出OTU丰度矩阵和分组信息作为输入文件，通过STAMP软件对各分类单元在两组样本中的丰度分布差异进行Welch’s t检验，从而获得在两组中存在显著性差异的分类单元并输出可视化结果。
[0070] 步骤16：以步骤11输出OTU丰度矩阵和分组信息作为输入文件，通过STAMP软件对各分类单元在两个样本中的丰度分布差异进行Fisher’s检验，从而获得在两个样本中存在显著性差异的分类单元并输出可视化结果。
[0071] 步骤18：以步骤11输出OTU丰度矩阵和分组信息作为输入文件，通过STAMP软件对各分类单元在多组样本中的丰度分布差异进行ANOVA方差分析/Kruskal-Wallis H检验，从而获得在多组样本中存在显著性差异的分类单元并输出可视化结果。
[0072] 步骤19：以步骤2输出OTU丰度矩阵作为输入文件，通过MEGAN软件将每个样本所含有的OTU的丰度信息和分类学组成数据映射至NCBI Taxonomy所提供的微生物分类等级树，统一呈现所有样本在各分类水平的具体组成。
[0073] 步骤20：以步骤11输出OTU丰度矩阵作为输入文件，通过graphlan.sh脚本对样本总体在各分类水平的组成构建等级树，同时以不同颜色区分各分类单元，并通过节点大小反映它们的丰度分布。
[0074] 步骤21：以步骤11输出OTU丰度矩阵作为输入文件，通过divide_sample_from_taxa.py脚本进行群落分类学组成的交互展示。
[0075] 步骤22：以步骤2和步骤11输出OTU丰度矩阵和分组信息作为输入文件，通过QIIME软件和Heatmap.R绘制热图。
[0076] 步骤23：以步骤11输出OTU丰度矩阵和分组信息作为输入文件，通过plot_pca_graph.R、plot_3d_pca_graph.R脚本对指定分类水平的群落组成结构进行PCA主成分分析，并且可视化。
[0077] 步骤24：以步骤8输出OTU丰度矩阵、步骤2输出系统发育数和分组信息作为输入文件，通过QIIME软件进行基于Unifrac距离计算样本间的距离矩阵，由加权及非加权距离矩阵分别进行PcoA主坐标分析，并且以二维和三维图像描述样本间基于微生物系统发育关系的群落空间分布特征；进行NMDS非度量多维尺度分析，通过二维或三维排序图描述群落样本的结构分布；进行UPGMA聚类分析；对加权及非加权距离矩阵分别进行组内和组间距离均值进行t检验，并通过1000次蒙特卡罗置换检验判断统计显著性，从而衡量组内和组间距离差异。
[0078] 步骤25：以步骤23输出nmds坐标和分组信息作为输入文件，通过plot_nmds.R脚本进行可视化。
[0079] 步骤26：以步骤11输出OTU丰度矩阵作为输入文件，通过make_ternary_plot.R脚本进行
三元相图分析显示出不同物种在各样本中的相对丰度差异。
[0080] 步骤27：以步骤11输出OTU丰度矩阵、环境因子和分组信息作为输入文件，通过rda.R脚本对指定分类水平的相对丰度矩阵进行RDA冗余分析，通过1000次置换检验确定统计显著性，并生成包含“样本—分类单元—影响因素”三种元素排序图。
[0081] 步骤28：以步骤23输出距离矩阵和分组信息作为输入文件，通过R软件BiodiversityR
软件包和ggplot2软件包进行CAP主坐标典型相关分析分析，并作1000次置换检验，确定组间差异是否具有统计学显著性。使用CAP分析得到的前两维数据生成CAP约束排序图。
[0082] 步骤29：以步骤11输出OTU丰度矩阵、环境因子和分组信息作为输入文件，通过PLSDA.R脚本构建PLS-DA判别模型并作图。
[0083] 步骤30：以步骤23输出距离矩阵和分组信息作为输入文件，通过QIIME软件进行Adonis/PERMANOVA分析，并作999次置换检验确定组间差异是否具有统计学显著性。
[0084] 步骤31：以步骤23输出距离矩阵和分组信息作为输入文件，通过QIIME软件进行ANOSIM分析，通过对样本距离等级排序来判断样本组内和组间差异的大小，并通过置换检验评价原始样本组间差异的统计学显著性。
[0085] 步骤32：以步骤8输出OTU丰度矩阵和分组信息作为输入文件，通过QIIME软件使用随机森林算法挑取丰度分布在不同组间存在显著差异的OTU。挑取OTU时使用1000棵随机森林决策树进行建模，并以10倍交叉验证估计误差的大小并完成随机森林分析。
[0086] 步骤33：以步骤11输出OTU丰度矩阵作为输入文件，通过network_spearman.R脚本和Cytoscape软件计算丰度位于前50位的分类单元之间的Spearman等级相关系数，对其中|rho|>0.8且P值<0.01的相关优势属构建关联网络及可视化。
[0087] 步骤34：以步骤8输出OTU丰度矩阵作为输入文件，通过QIIME软件和Cytoscape软件进行二分网络分析，并进行可视化，样本和OTU根据弹簧镶嵌模型排布。
[0088] 步骤35：以步骤1获得的高质量序列作为输入文件，使用PICRUSt软件，根据KEGG数据库中微生物代谢功能的类别对群落样本进行预测，并通过条形图展示各样本中编码的功能基因类别及其丰度。
[0089] 上述步骤中除少数步骤外，都可以通过linux
操作系统，以命令行的形式整合到shell脚本中，从而方便一次性执行，实现整个分析流程的自动化，提高分析效率。也可以根据需要单独执行。
