R包vegan的置換多元方差分析（PERMANOVA）判斷群落結構差異

R包vegan的置換多元方差分析（PERMANOVA）判斷群落結構差異

置換多元方差分析（Permutational multivariate analysis of variance，PERMANOVA），又稱非引數多因素方差分析（nonparametric multivariate analysis of variance）、或者ADONIS分析，其本質是基於F統計的方差分析，依據距離矩陣對總方差進行分解的非引數多元方差分析方法。使用PERMANOVA可分析不同分組因素對樣品差異的解釋度，並使用置換檢驗進行顯著性統計。

前文已經對R語言執行PERMANOVA的方法作了初步簡介。

PERMANOVA也常應用在物種多度資料的β多樣性分析中。例如在PCA、PCoA、NMDS等排序分析中，透過排序圖顯示了不同樣方之間群落結構區分明顯。但由於這些非約束排序僅為探索性分析，不涉及統計檢驗，無法提供一個明確的指示值，告訴我們所觀察到的差異是否是有效的，或者差異的程度。這時，即可透過PERMANOVA確定群落結構差異的顯著性。同樣地，對於無法在排序圖中區分的分組，也可透過PERMANOVA評估潛在的差異。

本篇以某16S擴增子測序所得的細菌群落資料為例，展示R包vegan進行PERMANOVA檢驗群落結構組成差異的一般過程。

示例資料及R程式碼的百度盤連結：

https://pan.baidu.com/s/1L7T71u1iSiiZpo9y6SeQzQ

示例資料概要

示例資料共涉及了40個16S測序樣本（細菌群落），均來自土壤。這40個土壤樣本共來自5個不同的取樣地點，每個地點各包含8個樣本（生物學重複）。

此處希望透過PERMANOVA，檢視不同取樣地點的土壤細菌群落組成結構是否具有顯著的不同。

檔案“otu_table.txt”為OTU水平的物種丰度表格。每一列為一個樣本，每一行為一種OTU，交叉區域為每種OTU在各樣本中的丰度。

檔案“group.txt”為樣本分組資訊。第一列（names）為各樣本名稱；第二列（site）為各樣本的分組資訊，即這些樣本所屬的取樣地點（s1，地點1；s2，地點2；s3……以此類推）。

檔案“bray.txt”為提前計算得到的距離矩陣檔案（此處為細菌群落間的Bray-curtis距離）。每一列為一個樣本，每一行為一個樣本，交叉區域為樣本間的Bray-curtis距離（取值範圍0-1，越接近於1表明樣本間細菌群落組成差異越大）。

PERMANOVA檢驗群落結構差異

全域性檢驗

首先在所有分組水平上，使用PERMANOVA檢驗整體差異，即檢視來自5個不同取樣地點所獲得的土壤細菌群落結構組成在整體上是否不具備一致性。

vegan中執行PERMANOVA的函式為adonis()。輸入adonis()的資料即可以為變數矩陣+分組資訊，也可以為距離矩陣+分組資訊。

library(vegan)

##讀入檔案

#OTU 丰度表

otu <- read.delim('otu_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)

otu <- data.frame(t(otu))

#樣本分組檔案

group <- read.delim('group.txt', sep = '\t', stringsAsFactors = FALSE)

##PERMANOVA 分析（所有分組間比較，即整體差異）

#根據 group$site 這一列分組進行 PERMANOVA 分析檢驗組間差異，基於 999 次置換，詳情 ?adonis

#（1）直接輸入 OTU 丰度表，在引數中指定距離型別

#使用 Bray-Curtis 距離測度

adonis_result <- adonis(otu~site, group, distance = 'bray', permutations = 999)

#（2）或者首先根據丰度表計算群落間距離，仍然使用 Bray-Curtis 距離測度，詳情 ?vegdist

dis1 <- vegdist(otu, method = 'bray')

#再將所的距離資料作為輸入

adonis_result <- adonis(dis1~site, group, permutations = 999)

#（3）若是已經提供好了距離矩陣，則直接使用現有的距離矩陣進行分析即可

#現有的距離矩陣，這裡為 Bray-Curtis 距離測度

dis <- read.delim('bray.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)

dis1 <- as.dist(dis)

adonis_result <- adonis(dis1~site, group, permutations = 999)

特別注意，為了避免識別錯誤，不要使用純數字作為分組名稱。對於置換次數，視資料量而定。

以上展示了3種處理過程，由於所使用距離型別相同（均為Bray-Curtis距離），分組資訊、置換檢驗次數（均為999次）等也一致，若置換次數足夠大，則結果也是一致的。

#檢視結果，上述 3 條命令所計算的內容一致，以其中一個為例

adonis_result

螢幕輸出了執行命令列、置換檢驗次數、差異計算結果等資訊。其中可主要關注中間的列表。site行即為本示例site分組的組間差異檢驗結果，Residuals為殘差，Total為site和Residuals的加和。

對於各項統計內容：Df，自由度，其值=所比較的分組數量-1；SumsOfSqs，即Sums of squares，總方差，又稱離差平方和；MeanSqs，即Mean squares，均方（差）；F.Model，F檢驗值；R2，即Variation (R²)，方差貢獻，表示不同分組對樣品差異的解釋度，即分組方差與總方差的比值，R²越大表示分組對差異的解釋度越高；Pr (>F)，顯著性p值，預設p<0.05即存在顯著差異。

對於各部分運算細節，可使用summary()命令檢視結果“adonis_result_dis”中所含內容，並提取相關結果檢視。

#或者例如

summary(adonis_result)

adonis_result$aov.tab

可將主要的統計結果轉化為資料框的型別，並輸出至檔案。

#可選輸出

otuput <- data.frame(adonis_result$aov.tab, check.names = FALSE, stringsAsFactors = FALSE)

otuput <- cbind(rownames(otuput), otuput)

names(otuput) <- c('', 'Df', 'Sums of squares', 'Mean squares', 'F.Model', 'Variation (R2)', 'Pr (>F)')

write.table(otuput, file = 'PERMANOVA.result_all.txt', row.names = FALSE, sep = '\t', quote = FALSE, na = '')

結果中包含了自由度（Df）、平方和（Sums of squares）、均方差（Mean squares）、F檢驗值（F.Model）、方差貢獻（Variation (R²)）、顯著性p值（Pr (>F)）等多項重要內容。

據此，我們可知5個取樣地點的土壤細菌群落結構在整體水平上是不一致的。

組間兩兩比較

可進一步使用PERMANOVA評估兩兩環境（兩兩取樣地）之間的土壤細菌群落結構組成的差異，包括是否具有顯著性，以及組間和組內的變異程度等。

##PERMANOVA 分析（使用迴圈處理，進行小分組間比較，如兩組間）

#推薦使用 OTU 丰度表作為輸入資料，每次篩選分組後重新計算樣本距離，避免由於樣本數減少可能導致的距離變動而造成誤差

group_name <- unique(group$site)

adonis_result_two <- NULL

for (i in 1:(length(group_name) – 1)) {

for (j in (i + 1):length(group_name)) {

group_ij <- subset(group, site %in% c(group_name[i], group_name[j]))

otu_ij <- otu[group_ij$names, ]

adonis_result_otu_ij <- adonis(otu_ij~site, group_ij, permutations = 999, distance = 'bray') #Bray-Curtis 距離測度，基於 999 次置換

adonis_result_two <- rbind(adonis_result_two, c(paste(group_name[i], group_name[j], sep = '/'), 'Bray-Curtis', unlist(data.frame(adonis_result_otu_ij$aov.tab, check.names = FALSE)[1, ])))

}

adonis_result_two <- data.frame(adonis_result_two, stringsAsFactors = FALSE)

names(adonis_result_two) <- c('group', 'distance', 'Df', 'Sums of squares', 'Mean squares', 'F.Model', 'Variation (R2)', 'Pr (>F)')

#可選新增 p 值校正，例如 Benjamini 校正

adonis_result_two$'Pr (>F)' <- as.numeric(adonis_result_two$'Pr (>F)')

adonis_result_two$P_adj_BH <- p.adjust(adonis_result_two$'Pr (>F)', method = 'BH')

#輸出

write.table(adonis_result_two, 'PERMANOVA.result_two.txt', row.names = FALSE, sep = '\t', quote = FALSE, na = '')

兩組間細菌群落PERMANOVA分析的結果中，包含了計算所依據的距離測度（Distance）、自由度（Df）、平方和（Sums of squares）、均方差（Mean squares）、F檢驗值（F.Model）、方差貢獻（Variation (R²)）、顯著性p值（Pr (>F)）以及BH校正後p值（P_adj_BH）等多項重要內容。