R包vegan的群落PCA及tb-PCA分析

主成分分析（Principal Component Analysis，PCA）是基于特征向量的线性无约束排序方法。在群落分析中，分析对象是原始的定量数据，排序图展示样方之间的欧几里得距离。

PCA的基本概念及其在群落分析中的应用，可参考前文。

本篇以某16S扩增子测序所得细菌群落数据，以及测量的环境因子数据为例，在生态学分析的角度，简介在R中运行PCA的过程。

示例文件、R脚本等，已上传至百度盘：

https://pan.baidu.com/s/1RKOSji5gPHqSey_AlJSwEA

（1）文件“phylum_table.txt”为通过16S测序所得的细菌类群丰度表格，此处统计至细菌“门水平”（即界门纲目科属种的“门”这一水平）。（为了方便演示，我没使用OTU水平的）

（2）文件“env_table.txt”记录了多种土壤理化指标信息，即环境数据。

接下来期望通过PCA分析这些数据。

环境数据的PCA分析

数据集中存在9种定量环境变量，它们有什么相关性？各样方的环境特征是否存在相似性？

PCA执行

vegan包中，可通过rda()执行PCA，该函数的输入格式大致如下。

library(vegan)
 
#详情 ?rda
#带协变量的 RDA
rda(Y, X, W)
#不带协变量的 RDA，即常规 RDA
rda(Y, X)
#PCA
rda(Y)

其中，Y为响应变量矩阵，X为解释变量矩阵，W为偏RDA分析需要输入的协变量矩阵。X或W是在RDA分析中才需要提供的，对于PCA而言，只需提供响应变量矩阵Y即可，即由对象（行）和变量（列）组成的矩阵。

因此，对环境数据的PCA分析，方法如下。

#读取环境数据
env <- read.delim('env_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
 
##PCA 排序
#环境变量需要标准化，详情 ?rda
pca_env <- rda(env, scale = TRUE)
 
#I 型标尺
summary(pca_env, scaling = 1)
#II 型标尺
#summary(pca_env, scaling = 2)

rda()中scale是对数据进行标准化。“scale=F”，不对变量标准化，排序对象是离差矩阵（包含方差和协方差的变量之间的关联矩阵）；“scale=T”，对变量标准化，排序对象是相关矩阵（不同量纲的变量之间的相关系数矩阵）。由于各环境变量具有不同的量纲，为了使不同物理单位的变量趋向一致，且降低极端值的影响，因此需要标准化。

本篇末尾处，也举例比较了标准化前后PCA的差异，并简述了环境数据标准化的必要性。

PCA结果解读

summary()展示结果时，根据所重要要关注的内容，可选根据I型标尺或II型标尺缩放排序坐标后展示出。如果对排序样方之间的距离更感兴趣，则考虑I型标尺；如果对变量之间的相关关系更感兴趣，则考虑II型标尺。关于I型标尺和II型标尺，可参考前文。

这里展示了按I型标尺缩放后排序结果（参数scaling = 1），部分内容如下所示。

主要内容提取

可以将重要的信息提取出来，例如：

#各主成分（PCA轴）特征根，特征根是固定的，和标尺选择无关
pca_eig <- pca_env$CA$eig
#除以特征根总和，即可得各主成分（PCA轴）解释量
pca_exp <- pca_env$CA$eig / sum(pca_env$CA$eig)
 
#提取对象排序坐标，以 I 型标尺为例，以前两轴为例
#scores() 提取样方坐标
site.scaling1 <- scores(pca_env, choices = 1:2, scaling = 1, display = 'site')
#或者在 summary 中提取
site.scaling1 <- summary(pca_env, scaling = 1)$sites[ ,1:2]
site.scaling1
#write.table(site.scaling1, 'site.scaling1.txt', col.names = NA, sep = '\t', quote = FALSE)
 
#提取变量排序坐标，以 II 型标尺为例，以前两轴为例
#scores() 提取环境变量坐标
env.scaling2 <- scores(pca_env, choices = 1:2, scaling = 2, display = 'sp')
#或者 summary 中提取
env.scaling2 <- summary(pca_env, scaling = 2)$species[ ,1:2]
env.scaling2
#write.table(env.scaling2, 'env.scaling2.txt', col.names = NA, sep = '\t', quote = FALSE)

PCA双序图

最后根据PCA所得的样方排坐标以及环境变量的排序坐标，绘图就可以了。

通常，排序对象直接在对应坐标处绘制为点；变量则呈现为向量，由原点(0,0)起始，指向变量得分的对应坐标处，向量的方向表示了该变量数值增加的方向。虽然也可以将变量在顶点坐标处直接展示为点（而非箭头形式），但考虑到它们可能会与对象坐标点混淆，因此一般不推荐这样表示。

#简洁版双序图，详情 ?biplot
pc1 <- paste('PC1:', round(pca_exp[1]*100, 2), '%')
pc2 <- paste('PC2:',round(pca_exp[2]*100, 2), '%')
 
par(mfrow = c(1, 2))
biplot(pca_env, choices = c(1, 2), scaling = 1, main = 'I 型标尺', col = c('red', 'blue'), xlab = pc1, ylab = pc2)
biplot(pca_env, choices = c(1, 2), scaling = 2, main = 'II 型标尺', col = c('red', 'blue'), xlab = pc1, ylab = pc2)
env.scaling1 <- summary(pca_env, scaling = 1)$species[ ,1:2]

#ordiplot() 作图，I 型标尺为例，详情 ?ordiplot
png('pca_env.png', width = 600, height = 600, res = 100, units = 'px')

ordiplot(pca_env, type = 'n', choices = c(1, 2), scaling = 1, main = 'I 型标尺', xlab = pc1, ylab = pc2)
points(pca_env, dis = 'site', choices = c(1, 2), scaling = 1, pch = 21, bg = c(rep('red', 4), rep('orange', 4), rep('green3', 4)), col = NA, cex = 1.2)
arrows(0, 0, env.scaling1[ ,1], env.scaling1[ ,2], length = 0.1, lty = 1, col = 'blue')
text(pca_env, dis = 'sp', choices = c(1, 2), scaling = 1, col = 'blue', cex = 0.8)

dev.off()
#ggplot2 作图，I 型标尺为例
#提取样方和环境变量排序坐标，前两轴
site.scaling1 <- data.frame(summary(pca_env, scaling = 1)$sites[ ,1:2])
env.scaling1 <- data.frame(summary(pca_env, scaling = 1)$species[ ,1:2])
 
#添加分组信息
site.scaling1$sample <- rownames(site.scaling1)
site.scaling1$group <- c(rep('A', 4), rep('B', 4), rep('C', 4))
env.scaling1$group <- rownames(env.scaling1)
 
#ggplot2 作图
library(ggplot2)
 
p <- ggplot(site.scaling1, aes(PC1, PC2)) +
geom_point(aes(color = group)) +
scale_color_manual(values = c('red', 'orange', 'green3')) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.title = element_blank(), legend.key = element_rect(fill = 'transparent')) + 
labs(x = pc1, y = pc2) +
geom_vline(xintercept = 0, color = 'gray', size = 0.5) + 
geom_hline(yintercept = 0, color = 'gray', size = 0.5) +
geom_segment(data = env.scaling1, aes(x = 0, y = 0, xend = PC1, yend = PC2), arrow = arrow(length = unit(0.1, 'cm')), size = 0.3, color = 'blue') +
geom_text(data = env.scaling1, aes(PC1 * 1.1, PC2 * 1.1, label = group), color = 'blue', size = 3)
 
#ggsave('pca_env.pdf', p, width = 4.5, height = 4)
ggsave('pca_env.png', p, width = 4.5, height = 4)

物种组成数据的PCA分析

物种组成数据比较特殊。由于PCA本质上是一种展示对象之间欧几里得距离的线性方法，所以可能不适合原始的物种多度数据分析（特别是含有许多0值的情形，欧几里得距离对双零问题很敏感）。但是，可以事先对物种多度数据执行Hellinger转化，将转化后的物种多度数据用于PCA。原因可参考前文。

这种通过首先对原始物种组成数据作转化后，再执行PCA的方法，也称为基于转化的PCA（Transformation-based Principal Component Analysis，tb-PCA）。

更推荐的方法是使用主坐标分析（Principal Coordinate Analysis，PCoA）代替PCA，通过计算样方间的定量非对称距离测度（如Bray-curtis距离），实现多维物种空间中的样方排序。本篇不讨论PCoA，继续看tb-PCA。

我们通过tb-PCA分析物种组成数据，细菌门水平丰度表。

此时，rda()函数中的scale参数（标准化）未再使用“scale=TRUE”。由于物种数据的量纲单位都一致，且Hellinger转化本身即可降低高丰度物种的权重，所以我使用了“scale=FALSE”。

#读取物种数据
species <- read.delim('phylum_table.txt', sep = '\t', row.names = 1, stringsAsFactors = FALSE, check.names = FALSE)
 
#物种数据 Hellinger 预转化，详情 ?decostand
species_hel <- decostand(species, method = 'hellinger')
 
#使用转化后的物种数据执行 PCA，无需标准化
pca_sp <- rda(species_hel, scale = FALSE)
 
#summary(pca_sp, scaling = 1)
#summary(pca_sp, scaling = 2)
 
#这里的排序图就只简单展示样方，不展示物种了
par(mfrow = c(1, 2))
ordiplot(pca_sp, scaling = 1, display = 'site', type = 'text')
ordiplot(pca_sp, scaling = 1, display = 'site', type = 'points')
 
#其它细节可参考上文

其它细节信息，如结果解读、关键信息提取、排序图可视化等，参考上文内容即可，方法一致，就不再演示了。

有一点需要注意，由于物种数据在这里作了转化，所以排序图中的物种相关性并不完全等于原始数据的物种相关性。

其它相关内容

以下是一些常见的问题，在此列出来帮助理解PCA。

评估有解读价值的PCA轴

PCA只是一种探索性分析方法，主要目的是在低维空间尽可能多地展示数据的主要趋势特征。最终选择多少个观测轴、选择哪些轴并没有统一的标准，自己看着来。

通常的做法就是选取前2-4个排序轴，然后提取对象或变量在这些轴中的坐标，绘制排序图查看对象或变量沿这些轴的分布是否出现了某种规律。借此评估数据的趋势特征。

排位偏后的轴（通常就是第5轴之后）一般就不考虑了，因为它们的特征值太低。

可以通过柱形图由高到低展示所有的特征值，然后根据特征值做个评判，多少个排序轴值得保留和解读。（这也只是一种参考，并不是强制的）

#Kaiser-Guttman 准则
#计算特征值的平均值，可据此保留高于平均值的主成分轴
pca_eig[pca_eig > mean(pca_eig)]
 
#断棍模型
#将单位长度的“棍子”随机分成与 PCA 轴数一样多的几段，由高往低排序
#可据此选取特征根大于所对应断棍长度的轴，或者总和大于所对应断棍长度总和的前几轴
n <- length(pca_eig)
bsm <- data.frame(j=seq(1:n), p = 0)
bsm$p[1] <- 1/n
for (i in 2:n) bsm$p[i] <- bsm$p[i-1] + (1/(n + 1 - i))
bsm$p <- 100*bsm$p/n
bsm
 
#绘制每轴的特征根和方差百分比 
par(mfrow = c(2, 1))
barplot(pca_eig, main = '特征根', col = 'bisque', las = 2)
abline(h = mean(pca_eig), col = 'red')
legend('topright', '平均特征根', lwd = 1, col = 2, bty = 'n')
barplot(t(cbind(100 * pca_eig/sum(pca_eig), bsm$p[n:1])), beside = TRUE, main = '% 变差', col = c('bisque', 2), las = 2)
legend('topright', c('% 特征根', '断棍模型'), pch = 15, col = c('bisque', 2), bty = 'n')

评估高贡献的变量

哪些变量对数据特征具有较高的贡献度？

可以将变量向量垂直投影至两侧的坐标轴上，对于主要的PCA轴，根据变量投影到其上的相对长度，选择高贡献的变量。

还有就是通过平衡贡献圈判断，仅可在I型标尺中适用，平衡贡献圈概念可参考前文。向量长度超过平衡贡献圈半径的变量，其重要性相对较高。

II型标尺可帮助评估高度共线性、高度相关的变量，据此可以手动对变量去冗余。

#某大佬的函数
source('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/scripts/NumEcolR2/cleanplot.pca.R')
env <- read.delim('env_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
pca_env <- rda(env, scale = TRUE)
 
par(mfrow = c(1, 2))
cleanplot.pca(pca_env, scaling = 1)        #I 型标尺中加入平衡贡献圈，评估重要的变量
cleanplot.pca(pca_env, scaling = 2)

关于环境变量的标准化

PCA分析必须从相同量纲的变量矩阵开始，原因是需要将变量总方差分配给特征根。因此，变量必须有相同的物理单位，方差和才有意义（方差的单位是变量单位的平方）。或者，变量是无量纲的数据，例如标准化后的数据。

然后我们不妨比较一下，标准化前后的环境数据，执行PCA后的差异。

#对比环境变量标准化前后的差异
env <- read.delim('env_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
 
pca_env1 <- rda(env, scale = FALSE)    #不执行标准化
pca_env2 <- rda(env, scale = TRUE)      #执行标准化
 
#I 型标尺
par(mfrow = c(1, 2))
cleanplot.pca(pca_env1, scaling = 1)
cleanplot.pca(pca_env2, scaling = 1)
 
round(pca_env1$CA$eig/sum(pca_env1$CA$eig)*100, 2)
round(pca_env2$CA$eig/sum(pca_env2$CA$eig)*100, 2)

结果显而易见，标准化前（左）相对于标准化后（右）的环境数据，方差很大，受高权重数值（由于各变量物理单位不同，变量间数值大小差异巨大）的影响甚重，特别是第一主成分轴的特征值极高。所以环境数据标准化是很有必要的。

关于物种组成数据Hellinger的转化

一般而言，对于物种组成数据较少使用PCA这种线性模型去做，PCoA或NMDS等是比较推荐的。

如果仍要使用PCA的话，事实上，并不是说物种多度数据一定要作转化才行。如果物种多度数据比较均匀，且含有很少的0值（这是最关键的，避开了欧几里得距离对“双零”敏感这一问题），则不转化也是可以的。

比方说示例数据“phylum_table.txt”，由于是门水平（一个很大的类别）的类群统计，本身0值较少，其实也可以不执行Hellinger转化的。尽管如此，可能高丰度和低丰度种群的丰度差别还是非常大，数据并不“均匀”，此时可以通过在运行rda()时，通过“scale=TRUE”消除高丰度物种与稀有物种在重要性方面的差异，降低欧几里得距离对高数值的敏感性。

但如果0值过多，执行PCA时还是推荐Hellinger转化，我们换个数据说明下吧。

网盘文件“otu_table.txt”是OTU水平的物种丰度表，其中含有非常多的0值，然后我们比较下Hellinger转化前后的PCA差异。

#对比物种数据 Hellinger 转化前后的差异
otu <- read.delim('otu_table.txt', sep = '\t', row.names = 1, stringsAsFactors = FALSE, check.names = FALSE)
otu_hel <- decostand(otu, method = 'hellinger')
 
pca_sp1 <- rda(otu, scale = FALSE)      #使用 Hellinger 转化前的数据
pca_sp2 <- rda(otu_hel, scale = FALSE)        #使用 Hellinger 转化后的数据
pca_sp3 <- rda(otu, scale = TRUE)        #使用 Hellinger 转化前的数据，但执行标准化
 
#特征值提取
pca_exp1 <- pca_sp1$CA$eig / sum(pca_sp1$CA$eig)
pc1_sp1 <- paste('PC1:', round(pca_exp1[1]*100, 2), '%')
pc2_sp1 <- paste('PC2:',round(pca_exp1[2]*100, 2), '%')
 
pca_exp2 <- pca_sp2$CA$eig / sum(pca_sp2$CA$eig)
pc1_sp2 <- paste('PC1:', round(pca_exp2[1]*100, 2), '%')
pc2_sp2 <- paste('PC2:',round(pca_exp2[2]*100, 2), '%')
 
pca_exp3 <- pca_sp3$CA$eig / sum(pca_sp3$CA$eig)
pc1_sp3 <- paste('PC1:', round(pca_exp3[1]*100, 2), '%')
pc2_sp3 <- paste('PC2:',round(pca_exp3[2]*100, 2), '%')
 
#I 型标尺
par(mfrow = c(2, 2))
#ordiplot(pca_sp1, scaling = 1, display = 'site', type = 'text')
#ordiplot(pca_sp2, scaling = 1, display = 'site', type = 'text')
#ordiplot(pca_sp3, scaling = 1, display = 'site', type = 'text')
ordiplot(pca_sp1, dis = 'site', type = 'n', choices = c(1, 2), scaling = 1, main = 'Hellinger 前，不标准化', xlab = pc1_sp1, ylab = pc2_sp1)
points(pca_sp1, dis = 'site', choices = c(1, 2), scaling = 1, pch = 21, bg = c(rep('red', 12), rep('orange', 12), rep('green3', 12)), col = NA, cex = 1.2)
ordiplot(pca_sp2, dis = 'site', type = 'n', choices = c(1, 2), scaling = 1, main = 'Hellinger 后，不标准化', xlab = pc1_sp2, ylab = pc2_sp2)
points(pca_sp2, dis = 'site', choices = c(1, 2), scaling = 1, pch = 21, bg = c(rep('red', 12), rep('orange', 12), rep('green3', 12)), col = NA, cex = 1.2)
ordiplot(pca_sp3, dis = 'site', type = 'n', choices = c(1, 2), scaling = 1, main = 'Hellinger 前，标准化', xlab = pc1_sp3, ylab = pc2_sp3)
points(pca_sp3, dis = 'site', choices = c(1, 2), scaling = 1, pch = 21, bg = c(rep('red', 12), rep('orange', 12), rep('green3', 12)), col = NA, cex = 1.2)

Hellinger转化前，且未标准化的数据，存在很大的方差、组内离散度高、趋势不明显且略微呈现了“马蹄形效应”；Hellinger转化后或者未Hellinger转化但标准化后的数据，样方分布明显，但相较之下，Hellinger转化后的区分程度更高，主成分轴的解释量具代表性。

PCA图中添加聚类树

还可以对样方执行聚类，将聚类树与PCA排序图相结合展示，也不失为一种可视化好方法。

#排序图中添加聚类树
env <- read.delim('env_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
pca_env <- rda(env, scale = TRUE)
 
#I 型标尺
p <- ordiplot(pca_env, dis = 'site', type = 'n', choices = c(1, 2), scaling = 1)
points(pca_env, dis = 'site', choices = c(1, 2), scaling = 1, pch = 21, bg = c(rep('red', 4), rep('orange', 4), rep('green3', 4)), col = NA, cex = 1.2)
 
#例如 UPGMA 聚类
env_upgma <- hclust(dist(scale(env)), method = 'average')
ordicluster(p, env_upgma, col = 'gray')

排序图中的对象点或变量向量可以同比例缩放

有时发现vegan输出的排序图中，所展示的对象或变量坐标，与计算结果中输出的数值对应不起来？但是仔细看的话，尽管横纵坐标的数值不对应，但各对象之间或各变量之间的相对位置并没有变，似乎是同比例缩放了。

在绘制排序图时，有时会考虑将排序对象或变量的得分（即坐标）乘以一个常量，即对排序图中的点或向量坐标同比例放大或缩小一定的数值后展示，以产生易于解读的生态排序结果（可参见Legendre和Legendre，1998，404页）。要么同时缩放所有对象，要么同时缩放所有变量，不能只改变部分对象或变量的位置。这种做法对于解答生态学问题通常是没有影响的。

比方说A图使用原始样方及物种得分（排序坐标）绘图展示，但物种向量的长度过短影响观测。因此考虑将物种得分同比例放大4倍，即对物种向量顶点的原始横纵坐标均乘以4后表示，即得到B图的样式更加易于辨认。

那么，变量向量缩放了，怎么再通过平衡贡献圈评判重要的变量呢？很简单，平衡贡献圈的半径也同比例缩放一下就可以了。

比方说上述使用的函数“cleanplot.pca”（细心的同学可能也注意到了），就考虑到使用一个内部常数对结果重新标定，这时变量向量的箭头长度及平衡贡献圈的半径不完全等于原始值，而是等比例原始值。

参考资料

DanielBorcard, FranoisGillet, PierreLegendre, et al. 数量生态学:R语言的应用（赖江山 译）. 高等教育出版社, 2014.

Legendre P, Legendre L. Numerical Ecology. Second English edition. Developments in Environmental Modelling, 1998, 20, Elsevier

友情链接