原创 阿童木 微生态

本文由阿童木根据实践经验而整理，希望对大家有帮助。

原创微文，欢迎转发转载。

导读

本系列的上一篇推文，即“开篇”中已经描述了宏基因组研究的基本思路和方法。先回顾一下，首先是收集样本和样本信息，然后是抽提样本中的总DNA，接着是对DNA进行文库构建，最后是将包含微生物群DNA信息的文库进行高通量测序。通过以上步骤能获得以G为单位的宏基因组测序数据。接下来非常重要的一个工作就是解读这些数据，即利用生物信息学软件分析测序数据中微生物群落的组成、基因、功能、通路等等。

为了提高分析的准确性和效率，在解读这些数据前还要完成一件非常重要的工作——测序数据过滤：1）除去建库过程中在样本DNA上添加测序接头等序列；2）除去会影响后续生物信息学分析准确性的测序过程中产生的低质量的碱基和序列；3）除去样本中可能混杂的宿主DNA序列。本篇内容分为以下三块：1）认识高通量测序数据；2）测序数据过滤的具体方法；3）测序数据的质量检查和统计。

高通量测序数据

微生太采用Illumina Novaseq平台对样本总DNA进行PE150（双端150bp长度）测序。测序仪器会把检测到的DNA序列信息以FASTQ格式的文件保存，格式如下：

FASTQ文件每四行对应一条测序序列的信息：

第一行以符号“@”起始，接着是序列编号等信息；

第二行为实际测得的一条碱基序列，也称一条Read；

第三行通常是连接符“+”，也可以有其他信息；

第四行是第二行中对应的每个碱基的测序质量。

过滤

测序获得的原始数据（Raw Data）中难免会存在一些低质量数据和非目的数据。为了保证后续分析结果的准确性和可靠性，需要根据碱基的质量信息对原始数据进行处理——过滤：1）首先使用Cutadapt彻底清除原始数据中的Illumina接头序列；2）再用PrinSEQ去除低质量的序列片段和可疑序列；3）最后通过Bowtie2将序列比对到宿主的基因组，保留没有比对到宿主基因组的微生物序列做后续分析。过滤的具体方法和参数如下：

1） 去除测序接头序列

建库时需要在待测序列的两端加上测序接头。测序接头能和测序试剂中的flowcell结合辅助测序顺利进行。因为测序接头是人为加上去的序列，所以分析前需要把接头去掉。Cutadapt是除接头最常用的工具，参数如下：

2)  用prinseq进一步过滤

去除长度小于50bp的序列；

去除平均质量分低于1/1000的序列；

去除序列复杂度低于70（用entropy的方式计算）的序列；

去除 N 碱基（未知碱基）达到20bp数量的序列；

去除每条序列5'端10bp的碱基，参数如下：

3)  去除宿主序列

从人体、动物体或植物等宿主中采集的微生物群标本难免会参杂宿主本身的细胞。、经核酸提取、建库和测序最终宿主DNA序列会和目的DNA序列混在一起。现在已知的物种基本都有基因组参考数据库，因此我们可以通过序列比对的方法轻松的识别和去除宿主DNA。例如，使用Bowtie2和人类基因组参考数据库hg38（GRCh38）可以识别和去除人体微生物群样本中的人类DNA，参数如下：

质量检查和统计

每份DNA样品经Illumina Novaseq PE150测序都会得到两份数据：上游序列Read1.fastq和下游序列Read2.fastq。序列过滤前和过滤后，都要用FastQC统计这些序列的碱基质量并将结果进行统计和可视化，以便于评价序列质量和分析序列过滤有无效果。高通量测序中常用Q30（1/1000错误率）评价碱基质量的好坏，不同Q值的意义如下：

1) FastQC碱基质量统计：

 上图是序列过滤前后FastQC质检的统计结果：1）Raw（原始）数据中的Read1和Read2的碱基的质量值都很高，上下游序列的所有碱基质量的中位数（蓝色线）均在Q30以上；2）Read2的末端个别碱基的质量波动较大，但是经过滤所有碱基都能达到很高的质量水平。

2) 过滤前后序列统计信息：

 上图是10个示例测序数据过滤前后的统计信息，由图可见：1）测序采用的PE150测序；2）所有样本的测序数据量均 > 5G，部分样本能接近10G；3）过滤后的Clean数据在95%以上；4）宿主序列占总序列的百分比不到0.2%。可见原始数据的质量还是很高的。

通过序列过滤，我们能获得更加准确、可靠的数据。拿到这样的数据，我们才能放心的进行下游的生物信息学分析和统计学分析。下期将为大家带来微生物分类鉴定的篇章。

为了研究样品物种组成及多样性信息，我们用Kraken2对所有样品的全部的有效序列进行注释分类。Kraken2是基于Kraken后研发的的最新版本（2018）。Kraken系列软件是一种使用精确k-mer匹配的分类系统，可实现高精度和快速分类速度。该分类器将查询序列中的每个k聚体与包含给定k聚体的所有基因组的最低共同祖先（LCA）匹配。Kraken2在Kraken的基础上提供了重大改进，达到了数据库构建时间更短，数据库大小更小和分类速度更快的效果。在使用完Kraken2对宏基因组测序数据进行分类。我们继续用Bracken对Kraken2得到的分类结果进行分类后贝叶斯重新估算丰度来估算宏基因组样本的物种水平或属级丰度。

一、为何选择Kraken2联合Bracken？

1. Kraken被引量已达到1300+，远超同类宏基因组分类注释软件。

很多引用Kraken和Bracken的文章发表在Nature，Nature Genetics等高水平杂志，Rob Knight发起的美国肠道计划相关文章也采用了Kraken联合Bracken的方法。

2. 基于k-mer精确匹配的算法，Kraken不仅速度快，而且准确性更高；

3. 最新版本的Kraken2继承并优化了Kraken；

数据库构建时间更短，数据库大小更小和分类速度更快的效果。

4. Bracken基于贝叶斯算法大大提高了本套方法的敏感性；

详情请见

Kraken2官网：https://ccb.jhu.edu/software/kraken2/ 

Bracken官网：https://ccb.jhu.edu/software/bracken/ 

二、序列注释程度分析

基于Bracken结果的绝对丰度及注释信息，对每个样品在7个分类水平（界、门、纲、目、科、属、种）（Kingdom, Phylum, Class, Order, Family, Genus, Species )上的序列数目占总序列数的比例进行统计，可以有效的评估样本的物种注释分辨率（注释到属/种的比例越高表示样本的注释效果越好）。使用R软件ggplot2包中的ggplot函数可对示例样本在各分类水平注释的相对程度进行堆叠图绘制，可视化结果如下。

1）结果展示

注释到各分类水平的序列在各样品中的数量分布柱形图

2）结果说明

横坐标 (Sample Name) 是样品名，纵坐标 (Sequence Number Percent) 表示注释到该水平的序列数目占总注释数据的比率，柱状图自上而下的颜色顺序对应于右侧的图例颜色顺序。每个分类水平最高值为1，代表100%的序列都得到了至少在这个级别的注释。

三、物种组成分析

基于Bracken结果的绝对丰度及注释信息，计算每个样品中的物种分别在界、门、纲、目、科、属、种7个分类水平的相对丰度。有了注释和相对丰度信息，再加上样品表型信息（metadata），就可以利用QIIME2中的qiime taxa barplot功能对所有的数据进行分类整理和可视化，结果如下。

1）结果展示

各分类水平的物种在各样品中的丰度分布柱形图

2）结果说明

横坐标 (Sample Name) 是样品名，纵坐标 (Relative Abundance) 表示相对丰度。微生太网页版结题报告中有交互式网页，客户可以在其中打开7种分类水平的物种相对丰度图，并对样本或者物种在图片中呈现的顺序根据元数据（包括分组信息）或者相对丰度的大小进行调节，也可以下载所有分类的丰度表。Level 1，2，3，4，5，6，7依次界、门、纲、目、科、属、种7个分类水平。

3）结果统计

对示例数据分析和统计，物种组成分析能检测到的物种有：1）物种丰度：Archaea (39537), Bacteria (20076282), Phages (26), Plant sand Fungi (5259), Viruses (66)；2）物种占比情况：Archaea (0.20%), Bacteria (99.78%), Phages (0.00%), Plants and Fungi (0.03%), Viruses (0.00%)。

四、高丰度物种分析

在课题研究过程中，人们经常喜欢探究自己的宏基因组样品中有哪些特殊物种，其中高丰度（优势）物种分析正是此类研究最常用的手段，因此我们设计了“丰度前20的菌门在各个样品中分布”的分析。使用R软件ggplot2包中的ggplot函数可进行数据可视化。

1）结果展示

丰度前20的菌门在各个样品中的分布柱形图

 2）结果说明

通常在种/属水平因为分类种类过多而在图例中无法全部展示所有分类，特增加本图是对图4-2进行补充。横坐标(Sample Name)是样品名，纵坐标(Sequence Number Percent)表示注释到该门水平的序列数目占总注释数据的比率，柱状图自上而下的颜色顺序对应于右侧的图例颜色顺序。在门水平没有注释的序列被归为unclassified一类。图例中最多显示最优势的20个种类,余下的相对丰度较低的物种被归类为Other在图中展示。

总结

物种鉴定的好坏直接关系到下游的统计分析和功能验证能否顺利的进行，因此宏基因组分析中的物种组成分析至关重要。Kraken2和Bracken这两个工具是微生太从众多物种分类方法中挑选而出的工具，这套物种鉴定方法在速度、精确度、准确度和敏感性都要高于同类方法。

利用Kraken2联合Bracken的方法能快速、准确的从宏基因组测序数据中获取物种的注释和丰度信息（详见第二篇），获得这些数据之后就可以根据观测（表型）数据对样本进行分组统计学分析，寻找与自己感兴趣的观测有关的微生物。本篇将介绍微生太宏基因组结题报告中的五种物种分组统计分析，包含：1）聚类分析；2）LEfSe差异分析；3）ANCOM差异分析；4）共有物种分析；5）多样性分析。

一、聚类分析

为了研究不同样品间的相似性，可以通过对样品进行聚类分析从而构建样品的聚类。通过选取感兴趣的物种（默认选取物种绝对丰度排名前20），根据样品的物种组成或物种在各样本的丰度分布情况进行聚类，以此考察不同样品或者分组间的相似或差异性（以门水平聚类分析为例）。使用R软件pheatmap包中的pheatmap函数可进行聚类分析。

R软件：https://www.r-project.org/ 

1）结果展示

【图1】

2）结果说明

纵轴为样品名称信息，同时也包括了分组信息。横轴为物种注释名称（本图为门水平）。图中上方的聚类树为物种在各样本中分布的相似度聚类，左侧的聚类树为样品聚类树，中间的热图是物种的相对丰度热图，颜色与相对丰度的关系见图上方的刻度尺。

二、LEfSe差异分析

LEfSe的统计方法可以寻找每一个分组的特征微生物（默认为LDA>2的微生物），也就是相对于其他分组，在这个组中丰度较高的微生物。LEfSe不仅有linux本地版，也有windows网页版，windows用户也能在网页中进行LEfSe分析。

LEfSe网页版：http://huttenhower.sph.harvard.edu/galaxy/ 

1） 结果展示

【图2】

2）结果说明

每一横向柱形体代表一个物种，柱形体的长度对应LDA值，LDA值越高则差异越大。柱形的颜色对应该物种是那个分组的特征微生物，特征微生物（在对应分组中的丰度相对较高）。

三、ANCOM差异分析

ANCOM(Analysis of composition of microbiomes)是另外一种比较微生物组学数据中物种在组间的显著性差异的分析方法。ANCOM分析不依赖于数据的分布假设，克服了相对丰度分析法的一些限制，能够有效降低结果的假阳性，更准确的分析各分类水平的物种在组间的差异情况。使用QIIME2可进行ANCOM分析。

QIIME2：https://qiime2.org/ 

1）结果展示

【图3】

2）结果说明

在ANCOM分析中，W值是一个衡量组间差异显著性的统计量（类似F值，t值），W值越高，代表该物种在组间的差异显著性越高。图中的每一个点都代表了一个比较的物种，纵坐标代表W值，横坐标clr值代表组间样品丰度的差异程度，数字绝对值越高代表相对丰度差异越大。因此在图中的点越靠近右上角（或左上角），则代表该物种与其他物种（靠近横轴0位置的物种）相比更具有显著性差异。

四、共有物种分析

对于分组较少（小于等于5）的数据，可通过绘制韦恩图（Venn diagram）统计分析不同样品组之间特有或共有物种的类型和数目，直观的展示样本分组间物种组成相似性及重叠情况。Venn图中，每个花瓣内的数字是对应分组特有的物种数目，每个圆内的数字总和是该分组的总物种数目。使用R软件VennDiagram包中的venn.diagrm函数可进行venn图绘制。

1）结果展示

【图4】

2）结果说明

韦恩图显示的是不同分组间共有或特有的物种数量，每个圆代表一个分组，面积由越大物种越多。

五、多样性分析

目前适用于生态学研究的降维分析主要是主成分分析 (PCA，Principal Component Analysis)和无度量多维标定法（NMDS，Non-Metric Multi-Dimensional Scaling）分析。其中，PCA是基于线型模型的一种降维分析，它应用方差分解的方法对多维数据进行降维，从而提取出数据中最主要的元素和结构。PCA 能够提取出最大程度反映样品间差异的两个坐标轴，从而将多维数据的差异反映在二维坐标图上，进而揭示复杂数据背景下的简单规律。而NMDS是非线性模型，其目的是为了克服线性模型的缺点，更好地反映生态学数据的非线性结构，应用NMDS分析，根据样本中包含的物种信息，以点的形式反映在多维空间上，而不同样本间的差异程度则是通过点与点间的距离体现，能够反映样本的组间或组内差异等。基于不同分类层级的物种丰度表，我们进行了PCA和NMDS分析，如果样品的物种组成越相似，则它们在PCA和NMDS图中的距离则越接近。

PCA是主坐标分析（PCoA）的一种特殊情况（基于欧式距离的PCoA就是PCA）。本报告利用PCoA分析方法，从多维数据中提取出最主要元素和能够最大程度反映样品间差异的三个坐标轴，从而将多维数据的差异反映在三维坐标图上，进而揭示复杂数据背景下的简单规律。我们基于Bray Curtis距离来进行PCoA分析，并选取贡献率最大的主坐标组合进行作图展示，图中样品的距离越接近，表示样品的物种组成结构越相似。微生太宏基因组分析报告中提供了PCoA的2D和3D图，方便选择使用。使用QIIME2可进行多样性分析。

1）结果展示

【图5】

2）结果说明

横坐标（Axis 1）表示第一主成分，百分比则表示第一主成分对样品差异的贡献值；纵坐标（Axis 2）表示第二主成分，百分比表示第二主成分对样品差异的贡献值；Axis 3坐标表示第三主成分，百分比表示第三主成分对样品差异的贡献值。客户可以根据元数据的分组信息可以个性化展示其他分类的距离关系。举例：可在右边设置区域的Select a Color Category下拉单选择分组信息，则同一个组的样品使用同一种颜色表示。

本篇内容分为以下三块：（1）NMDS和PCoA分析：上一篇的多样性分析中介绍了利用QIIME2进行PCoA分析研究菌群多样性，本篇将补充利用R进行NMDS和PCoA分析的内容；（2）RDA/CCA分析：这是以多元线性回归为基础的另一种PCA分析，常被用于宏基因组相关分析。例如，RDA/CCA分析可用来反映菌群与环境因子之间的关系，检测环境因子、样品、菌群（或功能）三者之间或者两两之间的关系；（3）相关分析：通过计算环境因子与物种间的Spearman相关系数，分析环境因子与宏组学数据（如，物种、KEGG功能、抗性基因等）之间的相关程度和显著性。

一、NMDS和PCoA

PCoA（主坐标分析）和NMDS（无度量多维标定法）是微生态研究中常用的降低数据维度的分析方法。基于欧式距离的PCoA就是我们熟悉的PCA。PCA是基于物种丰度矩阵的排序分析，PCoA是基于由物种丰度计算得到的距离矩阵。NMDS是非线性模型，它能克服线性模型的一些缺点，更好地反映生态学数据的非线性结构。PCoA和NMDS能通过从高维数据中提取出能够最大程度反映样品间差异的主要成分，揭示复杂数据背景下的简单规律。使用R软件VEGAN可进行微生物群PCoA和NMDS分析。计算距离矩阵的算法有很多，例如，Euclidean，Bray-Curtis，Jaccard，还有(un)weighted Unifrac，Bray-Curtis最为常用。有了矩阵数据就可以进行PCoA和NMDS分析。利用R软件基础包或ggplot2可进行结果可视化。

结果展示：

种水平物种NMDS和PCoA图

【图1】

种水平物种3D PCoA图

【图2】

结果说明：

横坐标表示第一主成分，百分比则表示第一主成分对样品差异的贡献值；纵坐标表示第二主成分，百分比表示第二主成分对样品差异的贡献值；第三坐标表示第三主成分，百分比表示第三主成分对样品差异的贡献值。一个点代表一个样本；颜色或者椭圆反应样本分组。

二、RDA/CCA分析

RDA/CCA分析是响应变量矩阵与解释变量之间多元多重线性回归的拟合值矩阵的PCA分析。RDA基于线性模型，CCA基于单峰模型。选择RDA还是CCA呢？首先要进行DCA分析，观察最大轴的值是否大于4，如果大于4.0，就选CCA，否则选RDA。RDA/CCA可分析微生物群与环境因子之间的相关性。其中，RDA/CCA分析给出的p值可以反应环境因子对微生物群落变异的影响是否显著。用R语言VEGAN包可进行RDA/CCA分析，用ggplot2可对结果进行可视化。

结果展示：

种水平物种RDA/CCA排序图1

【图3】

种水平物种CCA/RDA排序图2

【图4】

结果说明：

图三中的点代表细菌，图四中的点代表样本。环境因子用箭头表示，箭头的长度表示相关程度（箭头越长，相关性越大），箭头之间的夹角表示环境因子之间的相关关系（锐角为正相关，钝角为负相关）。图三中的每个点代表一个物种，点越大，物种丰度越高（红点为高丰度物种，灰点为低丰度物种）。图四中样本之间的距离大小表示菌群结构差异程度（距离越大，菌群差异程度越大）。

三、相关分析

进行相关分析一般首先要从所有物种中挑选出丰度较高的物种，接着计算物种丰度与环境因子相关矩阵，然后选出与环境因子相关最多的物种，最后用R语言pheatmap包可绘制相关分析热图。除了物种，也可以用GO、KEGG、抗性基因等宏基因组功能信息进行相关性分析，功能分析的内容将在后续篇章中讲解。

结果展示：

种水平物种与环境因子的相关性分析

【图5】

结果说明：

X轴上为环境因子，Y轴为物种。利用R语言pheatmap包中的函数可以计算衡量相关性大小的R值和校正错误发现率的p值。R值的大小在图中以不同颜色展示，右侧图例是不同R值的颜色区间。* 0.01≤ p 阈值，表示该功能在相应分组中丰度较高，而在其他分组中丰度较低，可作为候选Biomarker。下图是Pathway的LEfSe分析LDA柱形图。

结果说明：

每一横向柱形体代表一个Pathway，柱形体的长度对应LDA值，LDA值越高则差异越大。柱形的颜色对应代谢通路ID是哪个分组的特征Pathway，特征Pathway表示在对应分组中丰度相对较高的Pathway。为了方便客户了解不同KEGG功能层级上的组间差异，微生太对各个功能层级进行LEfSe分析。

三、Pathway可视化

利用KEGG数据库还可以绘制自己感兴趣或者组间有显著丰度差异的Pathway图（如下），在通路图中可以标注各个分组的特征基因（Gene biomarker）。

另外，点击报告中的网页文件可查看每个通路矩形节点的详细信息，鼠标悬在矩形框上可以看到对应通路节点有哪些KO在样本中检测到，以及这些KO是哪个分组的特征KO。网页形式如下： 

四、功能物种来源分析

LEfSe分析得到的差异Pathway/Module/KO结合HUMAnN2的分析结果，可进行功能物种来源分析，分析特征Pathway/Module/KO在每个样品的物种中的分布情况。下图是KO物种来源组成柱形图。其中横坐标对应样本，以及样本分组，不同分组用不同颜色标出。纵坐标对应各个样本该功能的相对丰度，不同物种来源用不同颜色标出。

五、聚类分析

从Pathway/Module/KO的相对丰度表中筛选出最大丰度排名前20的功能，接着可以通过对多个样品的功能进行丰度聚类来考察Pathway/Module/KO丰度在不同样品或者分组间的相似或差异性。下图是Pathway聚类热图。

六、 Circos图分析

Circos图可以用于展示每个样本中Pathway/Module/KO（丰度前10）的比例，以及各个功能在各个样本中的比例。下图是Pathway在各个样本中的分布Circos图。其中左半圈为丰度最高的十个Pathway，每个Pathway内，不同颜色代表不同样本来源的比例；右边半圈为样本，样本内不同颜色代表不同Pathway的比例。

七、 相关分析

相关分析可用于分析环境因子或其它组学数据与宏组学数据之间的相关性。先选出丰度较高的Pathway/Module/KO，再计算环境因子与抗性基因间的Spearman相关系数进行，最后选出与环境因子相关最多的20个Pathway，可以得到如下的Pathway与环境因子之间的相互关系热图。其中X轴上为环境因子，Y轴为Pathway。R值（相关系数）在图中以不同颜色展示，右侧图例是不同R值的颜色区间。* 0.01≤ P 阈值，表示该GO在相应分组中丰度较高，而在其他分组中丰度较低。下图是GO的LEfSe分析LDA柱形图。

三、功能物种来源分析

LEfSe分析得到的差异GO结合HUMAnN2的分析结果，可进行功能物种来源分析，分析GO在每个样品的物种中的分布情况。下图是GO物种来源组成柱形图,其中横坐标对应样本，以及样本分组，不同分组用不同颜色标出。纵坐标对应各个样本该功能的相对丰度，不同物种来源用不同颜色标出。

四、聚类分析

从GO的相对丰度表中筛选出丰度排名前20的GO，接着通过对多个样品的GO进行丰度聚类来考察GO丰度在不同样品或者分组间的相似或差异性。下图是GO聚类热图。

五、Circos图分析

Circos图可以用于展示每个样本中GO的丰度，以及各个GO分别在各个样本中的比例。下图是GO在各个样本中的分布Circos图。其中左半圈为丰度最高的十个GO，每个GO内，不同颜色代表不同样本来源的比例；右边半圈为样本，样本内不同颜色代表不同GO的比例。

六、 相关分析

相关分析可用于分析环境因子或其它组学数据与宏组学数据之间的相关性。选出丰度较高的GO，计算环境因子与抗性基因间的Spearman相关系数进行，再选出与环境因子相关最多的GO，可以得到如下的GO与环境因子之间的相互关系热图。其中X轴上为环境因子，Y轴为GO。R值（相关系数）在图中以不同颜色展示，右侧图例是不同R值的颜色区间。* 0.01≤ P 阈值，表示该Pathway在相应分组中丰度较高，而在其他分组中丰度较低。下图是Pathway的LEfSe分析LDA柱形图。

三、功能物种来源分析

LEfSe分析得到的差异Pathway结合HUMAnN2的分析结果，可进行功能物种来源分析，分析特征Pathway在每个样品的物种中的分布情况。下图是Pathway物种来源组成柱形图,其中横坐标对应样本，以及样本分组，不同分组用不同颜色标出。纵坐标对应各个样本该功能的相对丰度，不同物种来源用不同颜色标出。

四、聚类分析

从Pathway的相对丰度表中筛选出丰度排名前20的Pathway，接着通过对多个样品的Pathway进行丰度聚类来考察Pathway丰度在不同样品或者分组间的相似或差异性。下图是Pathway聚类热图。

五、Circos图分析

Circos图可以用于展示每个样本中Pathway的丰度，以及各个Pathway分别在各个样本中的比例。下图是Pathway在各个样本中的分布Circos图。其中左半圈为丰度最高的十个Pathway，每个Pathway内，不同颜色代表不同样本来源的比例；右边半圈为样本，样本内不同颜色代表不同Pathway的比例。

六、 相关分析

相关分析可用于分析环境因子或其它组学数据与宏组学数据之间的相关性。先选出丰度较高的Pathway，再计算环境因子与抗性基因间的Spearman相关系数进行，最后选出与环境因子相关最多的Pathway，可以得到如下的Pathway与环境因子之间的相互关系热图。其中X轴上为环境因子，Y轴为Pathway。R值（相关系数）在图中以不同颜色展示，右侧图例是不同R值的颜色区间。* 0.01≤ P 阈值，表示该COG在相应分组中丰度较高，而在其他分组中丰度较低。下图是COG的LEfSe分析LDA柱形图。

三、功能物种来源分析

LEfSe分析得到的差异COG结合HUMAnN2的分析结果，可进行功能物种来源分析，分析特征COG在每个样品的物种中的分布情况。下图是COG物种来源组成柱形图,其中横坐标对应样本，以及样本分组，不同分组用不同颜色标出。纵坐标对应各个样本该功能的相对丰度，不同物种来源用不同颜色标出。

四、聚类分析

从COG的相对丰度表中筛选出丰度排名前20的COG，接着通过对多个样品的COG进行丰度聚类来考察COG丰度在不同样品或者分组间的相似或差异性。下图是COG聚类热图。

五、Circos图分析

Circos图可以用于展示每个样本中COG的丰度，以及各个COG分别在各个样本中的比例。下图是COG在各个样本中的分布Circos图。其中左半圈为丰度最高的十个COG，每个COG内，不同颜色代表不同样本来源的比例；右边半圈为样本，样本内不同颜色代表不同COG的比例。

六、 相关分析

相关分析可用于分析环境因子或其它组学数据与宏组学数据之间的相关性。先选出丰度较高的COG，再计算环境因子与抗性基因间的Spearman相关系数进行，最后选出与环境因子相关最多的COG，可以得到如下的COG与环境因子之间的相互关系热图。其中X轴上为环境因子，Y轴为COG。R值（相关系数）在图中以不同颜色展示，右侧图例是不同R值的颜色区间。* 0.01≤ P 阈值，表示该碳水化合物活性酶在相应分组中丰度较高，而在其他分组中丰度较低。下图是碳水化合物活性酶的LEfSe分析LDA柱形图。

三、功能物种来源分析

LEfSe分析得到的差异碳水化合物活性酶结合HUMAnN2的分析结果，可进行功能物种来源分析，分析特征碳水化合物活性酶在每个样品的物种中的分布情况。下图是碳水化合物活性酶物种来源组成柱形图,其中横坐标对应样本，以及样本分组，不同分组用不同颜色标出。纵坐标对应各个样本该功能的相对丰度，不同物种来源用不同颜色标出。

四、聚类分析

从碳水化合物活性酶的相对丰度表中筛选出丰度排名前20的碳水化合物活性酶，接着通过对多个样品的碳水化合物活性酶进行丰度聚类来考察碳水化合物活性酶丰度在不同样品或者分组间的相似或差异性。下图是碳水化合物活性酶聚类热图。

五、Circos图分析

Circos图可以用于展示每个样本中碳水化合物活性酶的丰度，以及各个碳水化合物活性酶分别在各个样本中的比例。下图是碳水化合物活性酶在各个样本中的分布Circos图。其中左半圈为丰度最高的十个碳水化合物活性酶，每个碳水化合物活性酶内，不同颜色代表不同样本来源的比例；右边半圈为样本，样本内不同颜色代表不同碳水化合物活性酶的比例。 

六、相关分析

相关分析可用于分析环境因子或其它组学数据与宏组学数据之间的相关性。先选出丰度较高的碳水化合物活性酶，再计算环境因子与抗性基因间的Spearman相关系数进行，最后选出与环境因子相关最多的碳水化合物活性酶，可以得到如下的碳水化合物活性酶与环境因子之间的相互关系热图。其中X轴上为环境因子，Y轴为碳水化合物活性酶。R值（相关系数）在图中以不同颜色展示，右侧图例是不同R值的颜色区间。* 0.01≤ P 阈值，表示该EC酶在相应分组中丰度较高，而在其他分组中丰度较低。图2是EC酶的LEfSe分析LDA柱形图。

三、功能物种来源分析

LEfSe分析得到的差异EC酶结合HUMAnN2的分析结果，可进行功能物种来源分析，分析特征EC酶在每个样品的物种中的分布情况。图3是EC酶物种来源组成柱形图,其中横坐标对应样本，以及样本分组，不同分组用不同颜色标出。纵坐标对应各个样本该功能的相对丰度，不同物种来源用不同颜色标出。

四、聚类分析

从EC酶的相对丰度表中筛选出丰度排名前20的EC酶，接着通过对多个样品的EC酶进行丰度聚类来考察EC酶丰度在不同样品或者分组间的相似或差异性。下图图4是EC酶聚类热图。

图4 EC酶聚类热图

五、Circos图分析

Circos图可以用于展示每个样本中EC酶的丰度，以及各个EC酶分别在各个样本中的比例。下图是EC酶在各个样本中的分布Circos图。其中左半圈为丰度最高的十个EC酶，每个EC酶内，不同颜色代表不同样本来源的比例；右边半圈为样本，样本内不同颜色代表不同EC酶的比例。

六、相关分析

相关分析可用于分析环境因子或其它组学数据与宏组学数据之间的相关性。先选出丰度较高的EC酶，再计算环境因子与抗性基因间的Spearman相关系数进行，最后选出与环境因子相关最多的EC酶，可以得到如下的EC酶与环境因子之间的相互关系热图。其中X轴上为环境因子，Y轴为EC酶。R值（相关系数）在图中以不同颜色展示，右侧图例是不同R值的颜色区间。* 0.01≤ P 阈值，表示该抗性基因在相应分组中丰度较高，而在其他分组中丰度较低。下图是抗性基因的LEfSe分析LDA柱形图。

三、聚类分析

从抗性基因的相对丰度表中筛选出丰度排名前20的抗性基因，接着通过对多个样品的抗性基因进行丰度聚类来考察抗性基因丰度在不同样品或者分组间的相似或差异性。下图是抗性基因聚类热图。

四、Circos图分析

Circos图可以用于展示每个样本中抗性基因的丰度，以及各个抗性基因分别在各个样本中的比例。下图是抗性基因在各个样本中的分布Circos图。其中左半圈为丰度最高的十个抗性基因，每个抗性基因内，不同颜色代表不同样本来源的比例；右边半圈为样本，样本内不同颜色代表不同抗性基因的比例。 

 五、相关分析

相关分析可用于分析环境因子或其它组学数据与宏组学数据之间的相关性。先选出丰度较高的抗性基因，再计算环境因子与抗性基因间的Spearman相关系数进行，最后选出与环境因子相关最多的抗性基因，可以得到如下的抗性基因与环境因子之间的相互关系热图。其中X轴上为环境因子，Y轴为抗性基因。R值（相关系数）在图中以不同颜色展示，右侧图例是不同R值的颜色区间。* 0.01≤ P