R包vegan的主坐标分析（PCoA）

主坐标分析（Principal Coordinates Analysis，PCoA）是一种基于距离的排序方法（distance-based ordination methods）。与主成分分析（PCA）或对应分析（CA）相比，它以对象间相异矩阵开始（而非原始的对象-变量组成矩阵），为每个对象分配低维空间中的位置，可以灵活地选择任意关联测度类型。

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

本篇以某16S扩增子测序所得细菌群落数据为例，简介R包vegan的PCoA过程。

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

https://pan.baidu.com/s/1uNRUUyuip-D6hoDorDvEog

物种数据的PCoA分析

网盘文件“otu_table.txt”为通过16S测序所得的细菌群落丰度表格，OTU水平。

接下来期望通过PCoA，评估3组群落的物种组成差异水平。

PCoA执行

由于PCoA识别的是距离测度数据，而非原始的样方-物种丰度对应关系，所以首先需要根据OTU的丰度组成，计算样方间距离。并用计算好的距离测度作为PCoA的输入数据。

如果已经提供了现有的距离矩阵文件，则可以直接读取该距离矩阵，直接作为PCoA的输入。

此处展示两种过程，若两种过程中所使用的距离测度是一致的，则最后所得结果也是相同的。vegan中，cmdscale()函数用于执行PCoA。

library(vegan)

##读取 OTU 丰度表，排序
otu <- read.delim('otu_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- data.frame(t(otu))

#根据物种组成计算样方距离，如 Bray-curtis 距离，详情 ?vegdist
bray_dis <- vegdist(otu, method = 'bray') #结果以 dist 数据类型存储

#输出距离矩阵
#write.table(as.matrix(bray_dis), 'bray_distance.txt', sep = '\t', col.names = NA, quote = FALSE)

#PCoA 排序，详情 ?cmdscale
pcoa <- cmdscale(bray_dis, k = (nrow(otu) - 1), eig = TRUE)

##或者读取已经准备好的距离矩阵文件，如 Bray-curtis 距离，排序
dis <- read.delim('bray_distance.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
bray_dis <- as.dist(dis) #转为 dist 数据类型

#PCoA 排序，详情 ?cmdscale
pcoa <- cmdscale(bray_dis, k = (nrow(dis) - 1), eig = TRUE)

如果使用非欧氏距离（比方说这里的Bray-curtis距离），PCoA一般会产生负特征根。负特征根的处理详见下文。 

主要内容提取

在排序结果中提取主要结果，如特征值，样方得分（排序坐标）等。

#各 PCoA 轴的特征值
pcoa_eig <- pcoa$eig
#先评估下负特征值（末尾几个轴）
barplot(pcoa_eig)

#如果负特征值影响甚微，则继续
#如果负特征值非常明显，则应当先校正

#这里我们先继续，负特征值的校正后面再提

#各 PCoA 轴的解释量
pcoa_exp <- pcoa$eig/sum(pcoa$eig)

#样方排序坐标
site <- pcoa$point
#或
site <- scores(pcoa)
#输出
write.table(site, 'pcoa_site.txt', sep = '\t', col.names = NA, quote = FALSE)

对于负特征值，不妨先观察下它们的数值。如果绝对值特别小，比方说本示例，忽略它们倒也无关紧要。

断棍模型等评估特征轴的方法，PCoA中仍然适用。

#断棍模型评估各轴特征值
n <- length(pcoa_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

barplot(t(cbind(100 * pcoa_eig/sum(pcoa_eig), bsm$p[n:1])), beside = TRUE, main = '% 变差', col = c('orange', 'bisque'), las = 2)
abline(h = mean(100 * pcoa_eig/sum(pcoa_eig)), col = 'red')
legend('topright', c('% 特征值', '断棍模型', '平均特征值'), pch = c(15, 15, NA), col = c('orange', 'bisque', 'red'), lwd = c(NA, NA, 1), bty = 'n')

被动添加物种变量

如前文所述，由于物种信息在相异矩阵的计算过程中丢失，原始PCoA结果中只有样方得分。

如若期望补充物种得分，物种变量可以直接与PCoA轴进行相关分析，或者通过与它们所在样方得分的多度加权平均与PCoA轴建立关联而投影到PCoA空间中。

#例如使用 wascores() 被动添加物种得分（坐标），丰度加权平均方法，详情 ?wascores
#断棍模型认为前 4 轴特征值具有代表性，这里就暂且先映射前 4 轴
species <- wascores(pcoa$points[,1:4], otu)

#输出
write.table(species, 'pcoa_species.txt', sep = '\t', col.names = NA, quote = FALSE)

PCoA排序图

最后根据PCoA所得的样方得分（排序坐标）绘图就可以了。可以同时将被动添加的物种得分也一并展示出。

#ordiplot() 简洁版排序图
par(mfrow = c(1, 2))
ordiplot(pcoa, type = 'text', main = 'PCoA仅样方')
ordiplot(pcoa, type = 'text', main = 'PCoA（带物种投影）')
points(species[ ,1:2], pch = 3, cex = 0.7, col = 'blue')
#ordiplot() 作图，可视化调整
#断棍模型认为前 4 轴特征值具有代表性，这里就展示前 4 轴
#前 4 轴解释量
pcoa1 <- paste('PCoA axis1 :', round(100*pcoa_exp[1], 2), '%')
pcoa2 <- paste('PCoA axis2 :', round(100*pcoa_exp[2], 2), '%')
pcoa3 <- paste('PCoA axis3 :', round(100*pcoa_exp[3], 2), '%')
pcoa4 <- paste('PCoA axis4 :', round(100*pcoa_exp[4], 2), '%')

#所有物种太多了，挑主要的展示，如 top10 丰度物种
#计算 top10 丰度物种
abundance <- apply(otu, 2, sum)
abundance_top10 <- names(abundance[order(abundance, decreasing = TRUE)][1:10])

#作图，前 4 轴两两组合展示下
par(mfrow = c(2, 2))

ordiplot(site[ ,1:2], type = 'none', main = 'PCoA仅样方', xlab = pcoa1, ylab = pcoa2)
points(site[ ,1:2], pch = 19, col = c(rep('red', 12), rep('orange', 12), rep('green3', 12)), cex = 1)
#text(site[ ,1:2], rownames(site), cex = 0.7)

ordiplot(site[ ,c(3,4)], type = 'none', main = 'PCoA仅样方', xlab = pcoa3, ylab = pcoa4)
points(site[ ,c(3,4)], pch = 19, col = c(rep('red', 12), rep('orange', 12), rep('green3', 12)), cex = 1)
#text(site[ ,c(3,4)], rownames(site), cex = 0.7)

ordiplot(site[ ,1:2], type = 'none', main = 'PCoA（带top10丰度物种投影）', xlab = pcoa1, ylab = pcoa2)
points(site[ ,1:2], pch = 19, col = c(rep('red', 12), rep('orange', 12), rep('green3', 12)), cex = 1)
#text(site[ ,1:2], rownames(site), cex = 0.7)
text(species[abundance_top10,1:2], abundance_top10, cex = 0.7, col = 'blue')
#ggplot2 作图，已知第3、4 轴无法区分样方差异，就只观测前 2 轴
#添加分组信息
site <- data.frame(pcoa$point)[1:2]
site$name <- rownames(site)
site$group <- c(rep('A', 12), rep('B', 12), rep('C', 12))

species_top10 <- data.frame(species[abundance_top10,1:2])
species_top10$name <- rownames(species_top10)

#ggplot2 作图
library(ggplot2)

p <- ggplot(data = site, aes(X1, X2)) +
geom_point(aes(color = group)) +
stat_ellipse(aes(fill = group), geom = 'polygon', level = 0.95, alpha = 0.1, show.legend = FALSE) + #添加置信椭圆，注意不是聚类
scale_color_manual(values = c('red3', 'orange3', 'green3')) +
scale_fill_manual(values = c('red', 'orange', 'green3')) +
theme(panel.grid.major = element_line(color = 'gray', size = 0.2), panel.background = element_rect(color = 'black', fill = 'transparent'),
plot.title = element_text(hjust = 0.5), legend.position = 'none') +
geom_vline(xintercept = 0, color = 'gray', size = 0.5) +
geom_hline(yintercept = 0, color = 'gray', size = 0.5) +
geom_text(data = species_top10, aes(label = name), color = 'blue', size = 3) + #添 top10 丰度物种标签
labs(x = pcoa1, y = pcoa2, title = 'PCoA（带top10丰度物种投影）') +
annotate('text', label = 'A', x = -0.25, y = 0.04, size = 5, colour = 'red3') +
annotate('text', label = 'B', x = -0.05, y = -0.15, size = 5, colour = 'orange3') +
annotate('text', label = 'C', x = 0.25, y = 0.17, size = 5, colour = 'green3')

p

#ggsave('pcoa.pdf', p, width = 5.5, height = 5.5)
ggsave('pcoa.png', p, width = 5.5, height = 5.5)

其它相关内容

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

校正负特征根

如前文所述，由于绝大多数情况下，负特征根产生在末端的几个轴中，一般不会带来很大影响，此时大可忽略它们。例如上述示例。

但少数情形中，负特征根的绝对值大小与前几轴的正特征根相当，此时必须重视。以下是一些常见解决方案。

#由于仅欧式距离不产生负特征值
#可将非欧式距离的 Bray-curtis 距离平方根转化，使它获得欧式几何特征
#将转化后的距离用于PCoA，从而避免负特征值的产生
library(ade4)
is.euclid(bray_dis) #判断是否为欧式距离，显示为否

bray_dis_sqrt <- sqrt(bray_dis) #平方根转化
is.euclid(bray_dis_sqrt) #判断是否为欧式距离，显示为是

#PCoA 函数中本身也提供了校正参数（add），避免负特征值的产生
#Cailliez 校正
pcoa <- cmdscale(bray_dis, k = (nrow(otu) - 1), eig = TRUE, add = TRUE)

PCoA与PCA和CA

如前文所述，对于PCoA，若输入的距离测度为欧几里得距离，则PCoA将产生与PCA相似的排序结果。类似地，输入卡方距离测度，获得与CA相似的结果。

#物种多度数据的 PCA，推荐先作个 Hellinger 转化
otu_hell <- decostand(otu, method = 'hellinger')
tb_pca <- rda(otu_hell)

#物种多度数据的 CA
ca <- cca(otu)

#计算欧几里得距离，并 PCoA
euclid <- vegdist(otu_hell, method = 'euclidean')
pcoa_euclid <- cmdscale(euclid, k = (nrow(otu_hell) - 1), eig = TRUE)

#先做卡方标准化，再计算欧几里得距离，即可得卡方距离，再 PCoA
otu_chi <- decostand(otu, method = 'chi.square')
chi <- vegdist(otu_chi, method = 'euclidean')
pcoa_chi <- cmdscale(chi, k = (nrow(otu_chi) - 1), eig = TRUE)

#比较
par(mfrow = c(2, 2))
ordiplot(tb_pca, type = 'text', display = 'site', scaling = 1, main = 'PCA')
ordiplot(ca, type = 'text', display = 'site', scaling = 1, main = 'CA')
ordiplot(pcoa_euclid, type = 'text', main = 'PCoA-euclidean')
ordiplot(pcoa_chi, type = 'text', main = 'PCoA-chi.square')
#PCoA-chi.square 图有点别扭？先不慌
#由于轴的方向没有意义，因此可允许我们反转轴的方向展示
par(mfrow = c(2, 2))
ordiplot(tb_pca, type = 'text', display = 'site', scaling = 1, main = 'PCA')
ordiplot(ca, type = 'text', display = 'site', scaling = 1, main = 'CA')
ordiplot(pcoa_euclid, type = 'text', main = 'PCoA-euclidean')

pcoa_chi_site <- data.frame(pcoa_chi$point)
pcoa_chi_site[2] <- -pcoa_chi_site[2]
ordiplot(pcoa_chi_site, type = 'text', main = 'PCoA-chi.square')

结果中，忽略坐标轴刻度（PCA、CA中存在标尺缩放，即scaling参数；PCoA中没有），仅观测排序图中样方的相对位置，二者结果是一致的（如果比较未经标尺缩放前的最原始坐标，刻度也应当是一致的）。

友情链接