|
"R實(shí)戰(zhàn)"專(zhuān)題·第24篇
今天與大家分享一種分析方法，與Mantel test非常相似，常用于評(píng)估物種與環(huán)境或功能的關(guān)聯(lián)度，被稱(chēng)為普魯克分析（Procrustes Analysis）,也稱(chēng)為普氏分析。為什么分享呢？因?yàn)楫?dāng)我們只聚焦于兩個(gè)數(shù)據(jù)集一致性時(shí)，普氏分析似乎更為直觀，略勝一籌。簡(jiǎn)要介紹該分析方法的原理和適用范圍、R語(yǔ)言實(shí)現(xiàn)過(guò)程、提取相應(yīng)結(jié)果進(jìn)行繪圖。
普魯克分析（Procrustes analysis）的簡(jiǎn)介 普魯克分析（Procrustes analysis）是一種用于形狀分布的分析方法。數(shù)學(xué)上：通過(guò)不斷迭代，尋找標(biāo)準(zhǔn)形狀（canonical shape），并利用最小二乘法尋找每個(gè)樣本形狀到這個(gè)標(biāo)準(zhǔn)形狀的仿射變化方式。普氏分析可基于不同多元數(shù)據(jù)集的排序構(gòu)型（≥2組），通過(guò)平移、旋轉(zhuǎn)、縮放等轉(zhuǎn)換方式，實(shí)現(xiàn)最大疊合（maximal superimposition），用于各數(shù)據(jù)集間的對(duì)比分析。排序方法可選擇PCA、PCoA、NMDS等。換句話說(shuō)，普氏分析的作用可以看作是一種對(duì)原始數(shù)據(jù)的預(yù)處理，目的是為了獲取更好的局部變化模型作為后續(xù)模型學(xué)習(xí)的基礎(chǔ)。如下圖所示，每一個(gè)人臉特征點(diǎn)可以用一種單獨(dú)的顏色表示；經(jīng)過(guò)歸一化變化，人臉的結(jié)構(gòu)越來(lái)越明顯，即臉部特征簇的位置越來(lái)越接近他們的平均位置；經(jīng)過(guò)一系列迭代，尺度和旋轉(zhuǎn)的歸一化操作，這些特征簇變得更加緊湊，它們的分布越來(lái)越能表達(dá)人臉表情的變化。（剔除剛性部分、保留柔性部分） 
對(duì)于微生態(tài)領(lǐng)域，該方法通常用于分析自相同樣本不同數(shù)據(jù)集之間是否存在相似性關(guān)系。例如物種與環(huán)境、物種與功能基因等等，可以分析兩個(gè)數(shù)據(jù)集之間的關(guān)聯(lián)度/是否存在潛在的一致性。
下面我們將利用R語(yǔ)言vegan包中procrustes和protest兩個(gè)函數(shù)來(lái)實(shí)現(xiàn)普氏分析，并利用ggplot2包完成對(duì)應(yīng)分析的繪圖。普氏分析的R語(yǔ)言實(shí)現(xiàn)過(guò)程 1) 加載R包和自帶數(shù)據(jù)集（varespec和varechem），并對(duì)數(shù)據(jù)進(jìn)行相應(yīng)的距離轉(zhuǎn)換 # Procrustes Anaylsis ----------
# 評(píng)估物種群落結(jié)構(gòu)與環(huán)境因子間是否具存在顯著的相關(guān)性（兩個(gè)數(shù)據(jù)集）
# 加載R包
library(vegan)
# 加載數(shù)據(jù)
data(varespec)
head(varespec)
data(varechem)
head(varechem)
# 需要先對(duì)數(shù)據(jù)計(jì)算其樣本間的距離
# 一般物種群落使用bray-curtis距離，而環(huán)境因子使用歐式距離
spe.dist <- vegdist(varespec) # 默認(rèn)Bray-Curtis
env.dist <- vegdist(scale(varechem), "euclid")查看數(shù)據(jù)格式，需要注意：無(wú)論分析何種數(shù)據(jù)，兩個(gè)數(shù)據(jù)表格的樣品均要一一對(duì)應(yīng)；數(shù)據(jù)的格式如下，就是普通的寬格式的表格，第一行為物種或環(huán)境因子名，第一列為樣本名。 
2）由于兩個(gè)數(shù)據(jù)集的屬性不同，需要分別對(duì)兩個(gè)數(shù)據(jù)集進(jìn)行降維分析，并提取前兩軸坐標(biāo)（具有數(shù)據(jù)集代表性的線性組合）進(jìn)行比較。 # 在進(jìn)行普魯克分析前，需要先對(duì)兩個(gè)數(shù)據(jù)進(jìn)行降為分析，這里使用的是NMDS
# 也可以使用PCA或者PCoA，然后進(jìn)行普魯克分析，計(jì)算顯著性
mds.s <- monoMDS(spe.dist)
mds.e <- monoMDS(env.dist)3) 進(jìn)行普氏分析，用到procrustes函數(shù)，我們先大概了解一下該函數(shù)。 
4) 通過(guò)對(duì)該函數(shù)中各參數(shù)的了解，可知X為目標(biāo)矩陣也就是降維后的環(huán)境（功能基因等）坐，Y為降維后的物種數(shù)據(jù)的坐標(biāo)，因?yàn)楹罄m(xù)普氏分析中旋轉(zhuǎn)和縮放操作是針對(duì)Y，將Y匹配給X。 另外，當(dāng)symmetric=FALSE時(shí)，處于"非對(duì)稱(chēng)"模式，X和Y的分配值調(diào)換后，普氏分析的偏差平方和（M2）也會(huì)隨之改變。而當(dāng)symmetric=TRUE時(shí)，從而給出更合適比例的對(duì)稱(chēng)統(tǒng)計(jì)?！皩?duì)稱(chēng)”模式下，X和Y的分配值調(diào)換后，普氏分析的偏差平方和（M2）不會(huì)發(fā)生改變，但注意旋轉(zhuǎn)仍將是非對(duì)稱(chēng)的。# 以對(duì)稱(chēng)模式為例進(jìn)行普氏分析（symmetric = TRUE）
pro.s.e <- procrustes(mds.s,mds.e, symmetric = TRUE)
summary(pro.s.e)可得知：通過(guò)最小二乘求得組坐標(biāo)之間的偏差平方和（M2統(tǒng)計(jì)量）為0.65
5) 進(jìn)一步通過(guò)各樣本之間的殘差來(lái)判斷物種群落與環(huán)境因子的一致性。 plot(pro.s.e, kind = 2)
residuals(pro.s.e)從圖1來(lái)看，如果樣本中物種與環(huán)境一致性（相似性）越近，則對(duì)應(yīng)的殘差越小，反之物種與環(huán)境的相似性越遠(yuǎn)，則殘差越大（三條輔助線對(duì)應(yīng)的位置分別為殘差25%、50%和75%）；圖2為對(duì)應(yīng)的樣本殘差值。 
6) 事后檢驗(yàn)，對(duì)偏差平方和M2統(tǒng)計(jì)量進(jìn)行置換999次的普氏檢驗(yàn)。由于置換999次檢驗(yàn)會(huì)存在細(xì)微的誤差，我們?cè)O(shè)定了種子數(shù)，避免重復(fù)存在差異。 # 普氏分析中M2統(tǒng)計(jì)量的顯著性檢驗(yàn)
set.seed(1)
pro.s.e_t <- protest(mds.s,mds.e, permutations = 999)
pro.s.e_t
提取對(duì)應(yīng)的結(jié)果：
# 偏差平方和（M2統(tǒng)計(jì)量）
pro.s.e_t$ss
# 對(duì)應(yīng)p值結(jié)果
pro.s.e_t$signif999次置換檢驗(yàn)后顯示p<0.001，結(jié)果是非常顯著的 
基于ggplot2包ggplot函數(shù)將其結(jié)果繪制成圖，并利用export包導(dǎo)出圖片。首先提取降維后的數(shù)據(jù)軸1和2的坐標(biāo)，并且提取轉(zhuǎn)換的坐標(biāo)；然后進(jìn)行繪制。library(ggplot2)
# 獲得x和y軸的坐標(biāo)及旋轉(zhuǎn)過(guò)的坐標(biāo)
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 繪制兩點(diǎn)間的直線
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的優(yōu)勢(shì)和不足
Forcino(2015)也證實(shí)了Jackson(1995)和Peres-Neto和Jackson(2001)的說(shuō)法，即Protest是一種強(qiáng)大而準(zhǔn)確用于度量數(shù)據(jù)集間相似性的分析方法，甚至優(yōu)于Mantel test。然而，在某些情況下，Mantel test產(chǎn)生的結(jié)果與Protest不同。例如，Manteltest更接近于三個(gè)更相似的數(shù)據(jù)集比較中的兩個(gè)數(shù)據(jù)集的評(píng)估評(píng)分的變化幅度。因此，我們?cè)诜治鰰r(shí)要根據(jù)我們的目標(biāo)和實(shí)際情況進(jìn)行選擇，特別是聚焦于兩個(gè)數(shù)據(jù)集的一致性時(shí)，Procrustes可能是更優(yōu)的方案。 參考鏈接： 1.https://www.cnblogs.com/nsnow/p/4745730.html 2.https://www./article/8301917564/ 3.https://www.jianshu.com/p/ec8faacac230 4.https://www./science/article/pii/S0031018215001546 (Evaluating the effectiveness of the Mantel test and Procrustes randomization test for exploratory ecological similarity among paleocommunitie）
|
