[发明专利]一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法

摘要

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

主权项

1.一种用于16S rRNA基因的细菌群落组成和多样性分析的自动化方法，其提供的16S rRNA测序数据分析流程以测序原始序列数据作为输入,调用业界标准的分析工具，最终对数据进行可视化，并得到易于解读的分析结果，其特征在于，具体包括以下步骤：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数据库中微生物代谢功能的类别对群落样本进行预测，并通过条形图展示各样本中编码的功能基因类别及其丰度。