普鲁克分析(Procrustes analysis)的简介
普氏分析的R语言实现过程
1) 加载R包和自带数据集(varespec和varechem),并对数据进行相应的距离转换
# Procrustes Anaylsis ----------
# 评估物种群落结构与环境因子间是否具存在显著的相关性(两个数据集)
# 加载R包
library(vegan)
# 加载数据
data(varespec)
head(varespec)
data(varechem)
head(varechem)
# 需要先对数据计算其样本间的距离
# 一般物种群落使用bray-curtis距离,而环境因子使用欧式距离
spe.dist <- vegdist(varespec) # 默认Bray-Curtis
env.dist <- vegdist(scale(varechem), "euclid")
查看数据格式,需要注意:无论分析何种数据,两个数据表格的样品均要一一对应;数据的格式如下,就是普通的宽格式的表格,第一行为物种或环境因子名,第一列为样本名。
2)由于两个数据集的属性不同,需要分别对两个数据集进行降维分析,并提取前两轴坐标(具有数据集代表性的线性组合)进行比较。
# 在进行普鲁克分析前,需要先对两个数据进行降为分析,这里使用的是NMDS
# 也可以使用PCA或者PCoA,然后进行普鲁克分析,计算显著性
mds.s <- monoMDS(spe.dist)
mds.e <- monoMDS(env.dist)
3) 进行普氏分析,用到procrustes函数,我们先大概了解一下该函数。
?procrustes
4) 通过对该函数中各参数的了解,可知X为目标矩阵也就是降维后的环境(功能基因等)坐,Y为降维后的物种数据的坐标,因为后续普氏分析中旋转和缩放操作是针对Y,将Y匹配给X。
# 以对称模式为例进行普氏分析(symmetric = TRUE)
pro.s.e <- procrustes(mds.s,mds.e, symmetric = TRUE)
summary(pro.s.e)
可得知:通过最小二乘求得组坐标之间的偏差平方和(M2统计量)为0.65
5) 进一步通过各样本之间的残差来判断物种群落与环境因子的一致性。
plot(pro.s.e, kind = 2)
residuals(pro.s.e)
从图1来看,如果样本中物种与环境一致性(相似性)越近,则对应的残差越小,反之物种与环境的相似性越远,则残差越大(三条辅助线对应的位置分别为残差25%、50%和75%);图2为对应的样本残差值。
6) 事后检验,对偏差平方和M2统计量进行置换999次的普氏检验。由于置换999次检验会存在细微的误差,我们设定了种子数,避免重复存在差异。
# 普氏分析中M2统计量的显著性检验
set.seed(1)
pro.s.e_t <- protest(mds.s,mds.e, permutations = 999)
pro.s.e_t
提取对应的结果:
# 偏差平方和(M2统计量)
pro.s.e_t$ss
# 对应p值结果
pro.s.e_t$signif
999次置换检验后显示p<0.001,结果是非常显著的
普氏分析结果的绘制
library(ggplot2)
# 获得x和y轴的坐标及旋转过的坐标
Pro_Y <- cbind(data.frame(pro.s.e$Yrot), data.frame(pro.s.e$X))
Pro_X <- data.frame(pro.s.e$rotation)
# 绘图
ggplot(data=Pro_Y) +
geom_segment(aes(x = X1, y = X2,
xend = (X1 + MDS1)/2, yend = (X2 + MDS2)/2),
# geom_segment 绘制两点间的直线
arrow = arrow(length = unit(0, 'cm')),
color = "#9BBB59", size = 1) +
geom_segment(aes(x = (X1 + MDS1)/2, y = (X2 + MDS2)/2,
xend = MDS1, yend = MDS2),
arrow = arrow(length = unit(0.2, 'cm')),
color = "#957DB1", size = 1) +
geom_point(aes(X1, X2), color = "#9BBB59", size = 3, shape = 16) +
geom_point(aes(MDS1, MDS2), color = "#957DB1", size = 3, shape = 16) +
theme(panel.grid = element_blank(), # 绘制背景
panel.background = element_rect(color = 'black',
fill = 'transparent'),
legend.key = element_rect(fill = 'transparent'),
axis.ticks.length = unit(0.4,"lines"),
axis.ticks = element_line(color='black'),
axis.line = element_line(colour = "black"),
axis.title.x=element_text(colour='black', size=14),
axis.title.y=element_text(colour='black', size=14),
axis.text=element_text(colour='black',size=12)) +
labs(x = 'Dimension 1', y = 'Dimension 2', color = '') +
labs(title="Correlation between community and environment") +
geom_vline(xintercept = 0, color = 'gray', linetype = 2, size = 0.3) +
geom_hline(yintercept = 0, color = 'gray', linetype = 2, size = 0.3) +
geom_abline(intercept = 0, slope = Pro_X[1,2]/Pro_X[1,1], size = 0.3) +
geom_abline(intercept = 0, slope = Pro_X[2,2]/Pro_X[2,1], size = 0.3) +
annotate('text', label = 'Procrustes analysis:\n
M2 = 0.6297, p-value = 0.001',
x = -0.3, y = 0.3, size = 4,hjust = 0) +
theme(plot.title = element_text(size=14,colour = "black",
hjust = 0.5,face = "bold"))
library(export)
graph2ppt(file="Procrustes.ppt", append=T, height=5, width=5)
Pocrustes的优势和不足
Forcino(2015)也证实了Jackson(1995)和Peres-Neto和Jackson(2001)的说法,即Protest是一种强大而准确用于度量数据集间相似性的分析方法,甚至优于Mantel test。然而,在某些情况下,Mantel test产生的结果与Protest不同。例如,Manteltest更接近于三个更相似的数据集比较中的两个数据集的评估评分的变化幅度。因此,我们在分析时要根据我们的目标和实际情况进行选择,特别是聚焦于两个数据集的一致性时,Procrustes可能是更优的方案。
参考链接:
1.https://www.cnblogs.com/nsnow/p/4745730.html
2.https://www.freesion.com/article/8301917564/
3.https://www.jianshu.com/p/ec8faacac230
4.https://www.sciencedirect.com/science/article/pii/S0031018215001546
(Evaluating the effectiveness of the Mantel test and Procrustes randomization test for exploratory ecological similarity among paleocommunitie)
联系客服