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）等多项重要内容。

即可获知来自5个采样地点的土壤细菌群落中，两两之间的差异水平。据此评估各组群落间的相似/相异程度，以及帮助寻找在排序图中肉眼观察不明显的组间差异。

友情链接