R包vegan的MRPP分析判断群落结构差异
示例数据概要
示例数据共涉及了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检验群落结构差异
全局检验
首先在所有分组水平上，使用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值，可判断不同环境间群落物种组成差异水平的大小，以及评估同一环境中群落的变异程度等。