示例資料共涉及了40個16S測序樣本（細菌群落），均來自土壤。這40個土壤樣本共來自5個不同的取樣地點，每個地點各包含8個樣本（生物學重複）。
此處希望透過MRPP，檢視不同取樣地點的土壤細菌群落組成結構是否具有顯著的不同。
檔案“otu_table.txt”為OTU水平的物種丰度表格。每一列為一個樣本，每一行為一種OTU，交叉區域為每種OTU在各樣本中的丰度。
檔案“group.txt”為樣本分組資訊。第一列（names）為各樣本名稱；第二列（site）為各樣本的分組資訊，即這些樣本所屬的取樣地點（s1，地點1；s2，地點2；s3……以此類推）。
檔案“bray.txt”為提前計算得到的距離矩陣檔案（此處為細菌群落間的Bray-curtis距離）。每一列為一個樣本，每一行為一個樣本，交叉區域為樣本間的Bray-curtis距離（取值範圍0-1，越接近於1表明樣本間細菌群落組成差異越大）。
全域性檢驗
首先在所有分組水平上，使用MRPP檢驗整體差異，即檢視來自5個不同取樣地點所獲得的土壤細菌群落結構組成在整體上是否不具備一致性。
vegan中執行MRPP的函式為mrpp()。輸入mrpp()的資料即可以為變數矩陣+分組資訊，也可以為距離矩陣+分組資訊。
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)
##MRPP 分析（所有分組間比較，即整體差異）
#根據 group$site 這一列分組進行 MRPP 分析檢驗組間差異，基於 999 次置換，詳情 ?mrpp
#（1）直接輸入 OTU 丰度表，在引數中指定距離型別
#使用 Bray-Curtis 距離測度
mrpp_result <- mrpp(otu, group$site, distance = 'bray', permutations = 999)
#（2）或者首先根據丰度表計算群落間距離，仍然使用 Bray-Curtis 距離測度，詳情 ?vegdist
dis1 <- vegdist(otu, method = 'bray')
#再將所的距離資料作為輸入
mrpp_result <- mrpp(dis1, group$site, permutations = 999)
#（3）若是已經提供好了距離矩陣，則直接使用現有的距離矩陣進行分析即可
#現有的距離矩陣，這裡為 Bray-Curtis 距離測度
dis <- read.delim('bray.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
dis1 <- as.dist(dis)
mrpp_result <- mrpp(dis1, group$site, permutations = 999)
特別注意，為了避免識別錯誤，不要使用純數字作為分組名稱。對於置換次數，視資料量而定。
以上展示了3種處理過程，由於所使用距離型別相同（均為Bray-Curtis距離），分組資訊、置換檢驗次數（均為999次）等也一致，若置換檢驗次數已達到飽和，則結果也是一致的。
#檢視結果，上述 3 條命令所計算的內容一致，以其中一個為例
mrpp_result
螢幕輸出統計結果概要，包含了執行命令列、樣本數、置換檢驗次數、差異計算結果等資訊。其中，Significance of delta，即顯著性p值，小於0.05說明差異顯著；Chance corrected within-group agreement A，即A值，大於0說明組間差異大於組內差異，小於0說明組內差異大於組間差異；observe delta值越小說明組內差異小，expect delta值越大說明組間差異大。
對於各部分運算細節，可使用summary()檢視，並提取對應的結果。
#或者選擇提取需要的部分
summary(mrpp_result)
mrpp_result$Pvalue #如 p 值
可選提取多個重要的指標後，統一輸出至某檔案。
#可選輸出
write.table(data.frame(Group = 'all', Distance = 'Bray-Curtis', A = mrpp_result$A,
Observe_delta = mrpp_result$delta, Expect_delta = mrpp_result$E.delta, P_value = mrpp_result$Pvalue),
file = 'MRPP.result_all.txt', row.names = FALSE, sep = '\t', quote = FALSE)
示例輸出結果中，包含了計算所依據的距離測度（Distance）、A值（A）、observe delta值（Observe_delta）、expect delta值（Expect_delta）、顯著性p值（P_value）等多項內容。
據此，我們可知5個取樣地點的土壤細菌群落結構在整體水平上是不一致的。
組間兩兩比較
可進一步使用MRPP評估兩兩環境（兩兩取樣地）之間的土壤細菌群落結構組成的一致性等。
##MRPP 分析（使用迴圈處理，進行小分組間比較，如兩組間）
#推薦使用 OTU 丰度表作為輸入資料，每次篩選分組後重新計算樣本距離，避免由於樣本數減少可能導致的距離變動而造成誤差
group_name <- unique(group$site)
mrpp_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, ]
mrpp_result_otu_ij <- mrpp(otu_ij, group_ij$site, permutations = 999, distance = 'bray') #Bray-Curtis 距離測度，基於 999 次置換
mrpp_result_two <- rbind(mrpp_result_two, c(paste(group_name[i], group_name[j], sep = '/'), 'Bray-Curtis', mrpp_result_otu_ij$A, mrpp_result_otu_ij$delta, mrpp_result_otu_ij$E.delta, mrpp_result_otu_ij$Pvalue))
}
}
mrpp_result_two <- data.frame(mrpp_result_two, stringsAsFactors = FALSE)
names(mrpp_result_two) <- c('group', 'distance', 'A', 'Observe_delta', 'Expect_delta', 'P_value')
#可選新增 p 值校正過程，例如 Benjamini 校正
mrpp_result_two$P_value <- as.numeric(mrpp_result_two$P_value)
mrpp_result_two$P_adj_BH <- p.adjust(mrpp_result_two$P_value, method = 'BH')
#輸出
write.table(mrpp_result_two, 'MRPP.result_two.txt', row.names = FALSE, sep = '\t', quote = FALSE, na = '')
兩組間細菌群落的MRPP分析結果中，包含了計算所依據的距離測度（Distance）、A值（A）、observe delta值（Observe_delta）、expect delta值（Expect_delta）、顯著性p值（P_value）以及BH校正後p值（P_adj_BH）等多項重要內容。
透過顯著性p值，可檢視兩組間差異是否顯著；透過A值、observe delta值以及observe delta值，可判斷不同環境間群落物種組成差異水平的大小，以及評估同一環境中群落的變異程度等。
