【R分享|實戰】普魯克分析 (Procrustes Analysis) : 評估物種與環境/功能的關係

「 不求做的最好，但求做的更好。」   --科白君今天與大家分享一種分析方法，與Mantel test非常相似，常用於評估物種與環境或功能的關聯度，被稱為普魯克分析（Procrustes Analysis）,也稱為普氏分析。為什麼分享呢？因為當我們只聚焦於兩個數據集一致性時，普氏分析似乎更為直觀，略勝一籌。簡要介紹該分析方法的原理和適用範圍、R語言實現過程、提取相應結果進行繪圖。

普魯克分析（Procrustes analysis）的簡介

普魯克分析（Procrustes analysis）是一種用於形狀分布的分析方法。數學上：通過不斷迭代，尋找標準形狀（canonical shape），並利用最小二乘法尋找每個樣本形狀到這個標準形狀的仿射變化方式。普氏分析可基於不同多元數據集的排序構型（≥2組），通過平移、旋轉、縮放等轉換方式，實現最大疊合（maximal superimposition），用於各數據集間的對比分析。排序方法可選擇PCA、PCoA、NMDS等。換句話說，普氏分析的作用可以看作是一種對原始數據的預處理，目的是為了獲取更好的局部變化模型作為後續模型學習的基礎。如下圖所示，每一個人臉特徵點可以用一種單獨的顏色表示；經過歸一化變化，人臉的結構越來越明顯，即臉部特徵簇的位置越來越接近他們的平均位置；經過一系列迭代，尺度和旋轉的歸一化操作，這些特徵簇變得更加緊湊，它們的分布越來越能表達人臉表情的變化。（剔除剛性部分、保留柔性部分）

對於微生態領域，該方法通常用於分析自相同樣本不同數據集之間是否存在相似性關係。例如物種與環境、物種與功能基因等等，可以分析兩個數據集之間的關聯度/是否存在潛在的一致性。
下面我們將利用R語言vegan包中procrustes和protest兩個函數來實現普氏分析，並利用ggplot2包完成對應分析的繪圖。

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函數，我們先大概了解一下該函數。

4) 通過對該函數中各參數的了解，可知X為目標矩陣也就是降維後的環境（功能基因等）坐，Y為降維後的物種數據的坐標，因為後續普氏分析中旋轉和縮放操作是針對Y，將Y匹配給X。

另外，當symmetric=FALSE時，處於"非對稱"模式，X和Y的分配值調換後，普氏分析的偏差平方和（M2）也會隨之改變。而當symmetric=TRUE時，從而給出更合適比例的對稱統計。「對稱」模式下，X和Y的分配值調換後，普氏分析的偏差平方和（M2）不會發生改變，但注意旋轉仍將是非對稱的。

# 以對稱模式為例進行普氏分析（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，結果是非常顯著的

基於ggplot2包ggplot函數將其結果繪製成圖，並利用export包導出圖片。首先提取降維後的數據軸1和2的坐標，並且提取轉換的坐標；然後進行繪製。

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)

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）

分析中的數據可以添加我們的微信群獲得，

獲得途徑1，關注本公眾號，後臺回復 「客服微信」，小編將邀請您進群和我們一起交流和學習~

獲得途徑2，添加小編個人微信，小編將拉你進群。          

期待您的"分享"點讚"在看"