技术贴 | 微生太宏基因组报告解读 | 第四篇：PCoA、NMDS、RDA/CCA、相关分析
本篇内容分为以下三块：（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 <0.05，** 0.001≤ p < 0.01，*** p < 0.001。
参考：
1. 微生太宏基因组结题报告模板
2. 数量生态学笔记||冗余分析（RDA）概述
https://www.jianshu.com/p/00f69e8bd5ef
3. PCA、PCoA、NMDS、CCA、RDA傻傻分不清楚
http://www.dxy.cn/bbs/topic/32534684?sf=2&dn=10