在微生物群落研究的過程中，我們經常需要評估微生物群落結構與環境因子整體之間是否具有顯著的相關性，此時，通常使用的方式是Mantel test和普氏分析。
當然除了分析群落結構與環境因子的相關性之外，這兩個分析還可以用於分析同一樣品不同類型微生物群落之間的相關性，比如同一樣品的稀有和豐富物種或者同一樣品細菌和真菌群落結構的相關性。
不同微生物群落之間，該分析更多的還是用於分析微生物群落組成結構與其它功能基因組之間的關係，比如細菌組成結構與抗生素抗性組的相關性。
最後還有一種不太常用的用法，就是分析配對的兩種不同類型樣品微生物群落的相關性，比如河流或海洋同一位置水體和沉積物細菌群落組成結構的相關性，從而分析這兩種關聯樣品類型中微生物群落的轉移。
這兩種方法的具體原理我這裡就不講了，網上有很多，感興趣的朋友可以自行搜索一下，本推文主要介紹分析的實現方法以及對結果可視化的美化。
輸入數據既然是評估兩組數據的相關性，那麼輸入數據就是兩個數據表，分別對應要分析的兩組數據，本推文已一套樣本的細菌群落豐度數據及其對應的環境因子數據為例。
⚠️無論分析何種數據，兩個數據表格的樣品均要一一對應。
數據的格式如下，就是普通的寬格式的表格，第一行和第一列分別對應樣本和物種/環境因子，文件為制表符分割的txt格式。
注意行為樣本還是列為樣本並不重要，最終分析的時候是行為樣本，至少輸入數據後注意是否進行轉置即可。
這裡需要使用「vegan」包，首先載入該分析包，並將待分析數據導入。
在進行Mantel test和普氏分析之前，需要先對輸入的數據計算其樣本間的距離，一般物種組成數據使用bray-curtis距離，而環境因子數據使用歐式距離，具體代碼如下：
library(vegan)
data <- read.csv("class.txt", head=TRUE,sep="\t",row.names = 1)
data <- data[which(rowSums(data) > 0),]
data <- t(data)
s.dist <- vegdist(data,method = "bray")
data <- read.csv("environment.txt", head=TRUE,sep="\t",row.names = 1)
r.dist <- vegdist(data)
隨後進行Mantel test
mantel(s.dist,r.dist)
這裡注意記錄結果，主要是Mantel statistic r和Significance兩個數值。
默認是使用pearson相關性，也可以根據情況使用spearman相關性。
mantel(s.dist,r.dist,method = "spearman")
在進行普氏分析之前，需要先對兩個數據進行降為分析，這裡使用的是NMDS，也可以使用PCoA。
mds.s <- monoMDS(s.dist)
mds.r <- monoMDS(r.dist)
隨後進行普氏分析，並計算結果顯著性。
pro.s.r <- procrustes(mds.s,mds.r)
protest(mds.s,mds.r)
注意記錄結果，主要是Procrustes Sum of Squares (m12 squared)和Significance兩個數值。
結果可視化接下來使用「ggplot2」包進行結果的可視化，首先載入繪圖包並建立繪圖數據。
library(ggplot2)
Y <- cbind(data.frame(pro.s.r$Yrot), data.frame(pro.s.r$X))
X <- data.frame(pro.s.r$rotation)
Y$ID <- rownames(Y)
之後開始繪圖。
p <- ggplot(Y) +
geom_segment(aes(x = X1, y = X2, xend = (X1 + MDS1)/2, yend = (X2 + MDS2)/2),
arrow = arrow(length = unit(0, 'cm')),
color = "#B2182B", size = 1) +
geom_segment(aes(x = (X1 + MDS1)/2, y = (X2 + MDS2)/2, xend = MDS1, yend = MDS2),
arrow = arrow(length = unit(0, 'cm')),
color = "#56B4E9", size = 1) +
geom_point(aes(X1, X2), fill = "#B2182B", size = 4, shape = 21) +
geom_point(aes(MDS1, MDS2), fill = "#56B4E9", size = 4, shape = 21) +
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 = X[1,2]/X[1,1], size = 0.3) +
geom_abline(intercept = 0, slope = X[2,2]/X[2,1], size = 0.3) +
annotate('text', label = 'Procrustes analysis:\n M2 = 0.8358, p-value = 0.035\nMantel test:\n r = 0.1703, p-value = 0.04',
x = -1.5, y = 1.2, size = 4,hjust = 0) +
theme(plot.title = element_text(size=16,colour = "black",hjust = 0,face = "bold"))
然後保存圖片。
png("proc.png",width = 3600,height = 3200,res = 600)
p
dev.off()
pdf(file = "prco.pdf",width = 6,height = 5.4)
p
dev.off()
⚠️因為我懶得做一個公式連結自動標註統計學檢驗結果，所以是手動記錄結果只手自己輸入的，因此更換數據之後，在annotate那個代碼裡，你們要自己修改內容以及文本插入的坐標軸位置。
關注公眾號「紅皇后學術」，後臺回復「普氏分析」獲取示例數據和繪圖代碼！
另感謝杰哥提供了一個基本的繪圖代碼🙏！！