R包vegan的主坐标分析（PCoA）
物种数据的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中没有），仅观测排序图中样方的相对位置，二者结果是一致的（如果比较未经标尺缩放前的最原始坐标，刻度也应当是一致的）。
R包ade4的模糊主成分分析（FPCA）及模糊对应分析（FCA）
对应分析（CA）和去趋势对应分析（DCA）在群落分析中的应用