R包vegan的相似性分析（ANOSIM）判断群落结构差异
ANOSIM示意图，通过将每个类别内的距离与类别之间的平均距离进行比较，当类内距离较小且类间距离较大时，表明分类明显。
例如在群落分析中，常使用ANOSIM结合排序分析一起，评估不同环境群落组成结构的相似/差异程度，组间群落差异是否显著不同于组内差异。
本篇以某16S扩增子测序所得的细菌群落数据为例，展示R包vegan进行ANOSIM检验群落结构组成差异的一般过程。
示例数据及R代码的百度盘链接：
https://pan.baidu.com/s/11U3kECzMkPriCKl43Wk60w
示例数据概要
示例数据共涉及了40个16S测序样本（细菌群落），均来自土壤。这40个土壤样本共来自5个不同的采样地点，每个地点各包含8个样本（生物学重复）。
此处希望通过ANOSIM，查看不同采样地点的土壤细菌群落组成结构是否具有显著的不同。
文件“otu_table.txt”为OTU水平的物种丰度表格。每一列为一个样本，每一行为一种OTU，交叉区域为每种OTU在各样本中的丰度。
文件“group.txt”为样本分组信息。第一列（names）为各样本名称；第二列（site）为各样本的分组信息，即这些样本所属的采样地点（s1，地点1；s2，地点2；s3……以此类推）。
文件“bray.txt”为提前计算得到的距离矩阵文件（此处为细菌群落间的Bray-curtis距离）。每一列为一个样本，每一行为一个样本，交叉区域为样本间的Bray-curtis距离（取值范围0-1，越接近于1表明样本间细菌群落组成差异越大）。
ANOSIM检验群落结构差异
全局检验
首先在所有分组水平上，使用ANOSIM检验整体差异，即查看来自5个不同采样地点所获得的土壤细菌群落结构组成在整体上是否不具备一致性，各组间群落差异是否显著不同于各组内的差异。
vegan中执行ANOSIM的函数为anosim()。输入anosim()的数据即可以为变量矩阵+分组信息，也可以为距离矩阵+分组信息。
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)
##ANOSIM 分析（所有分组间比较，即整体差异）
#根据 group$site 这一列分组进行 ANOSIM 分析检验组间差异，基于 999 次置换，详情 ?anosim
#（1）直接输入 OTU 丰度表，在参数中指定距离类型
#使用 Bray-Curtis 距离测度
anosim_result <- anosim(otu, group$site, distance = 'bray', permutations = 999)
#（2）或者首先根据丰度表计算群落间距离，仍然使用 Bray-Curtis 距离测度，详情 ?vegdist
dis1 <- vegdist(otu, method = 'bray')
#再将所的距离数据作为输入
anosim_result <- anosim(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)
anosim_result <- anosim(dis1, group$site, permutations = 999)
特别注意，为了避免识别错误，不要使用纯数字作为分组名称。对于置换次数，视数据量而定。
以上展示了3种处理过程，由于所使用距离类型相同（均为Bray-Curtis距离），分组信息、置换检验次数（均为999次）等也一致，若置换检验次数已达到饱和，则结果也是一致的。
#查看结果，上述 3 条命令所计算的内容一致，以其中一个为例summary(anosim_result)
#或者
names(anosim_result)
anosim_result$signif #p 值
anosim_result$statistic #R 值
summary()后屏幕输出了包含运行命令行、置换检验次数、差异计算结果等信息。
其中，可主要关注两个重要统计值，R值（ANOSIM statistic R）和p值（Significance）。R值可以得出组间与组内比较的差异程度，其取值范围（-1，1）；R>0，说明组间差异大于组内差异，即组间差异显著；R<0，说明组内差异大于组间差异；R值的绝对值越大表明相对差异越大。p值越低表明这种差异检验结果越显著，一般以0.05为显著性水平界限。
使用plot()函数能够将计算结果可视化。
#作图展示#pdf(paste('anosim.all.pdf', sep = ''), width = 10, height = 5)
png(paste('anosim.all.png', sep = ''), width = 800, height = 400)
plot(anosim_result, col = c('gray', 'red', 'green', 'blue', 'orange', 'purple'))
dev.off()
ANOSIM分析基于两两样本之间的距离值排序获得的秩，这样任一两两组的比较可以获得三个分类的数据。以箱线图的形式展示组间与组内的秩的分布，横坐标表示所有样品 （Between）以及各分组（本示例为s1、s2、s3、s4、s5），纵坐标表示距离（本示例使用Bray-Curtis距离）的秩。当Between组相对与其他每个分组的秩较高时，则表明组间差异大于组内差异。同时图的上方标注了R值与p值两个重要统计指标，便于我们直观地对组间差异是否显著不同于各组内的差异进行判断。
据此，我们可知5个采样地点的土壤细菌群落结构在整体水平上是不一致的，各组间群落差异是否显著不同于各组内的差异。
组间两两比较
可进一步使用ANOSIM评估两两环境（两两采样地）之间的土壤细菌群落结构组成的差异，包括是否具有显著性，以及组间和组内的变异程度等。
##ANOSIM 分析（使用循环处理，进行小分组间比较，如两组间）
#推荐使用 OTU 丰度表作为输入数据，每次筛选分组后重新计算样本距离，避免由于样本数减少可能导致的距离变动而造成误差
group_name <- unique(group$site)
dir.create('anosim_two', recursive = TRUE)
anosim_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, ]
anosim_result_otu_ij <- anosim(otu_ij, group_ij$site, permutations = 999, distance = 'bray') #Bray-Curtis 距离测度，基于 999 次置换
anosim_result_two <- rbind(anosim_result_two, c(paste(group_name[i], group_name[j], sep = '/'), 'Bray-Curtis', anosim_result_otu_ij$statistic, anosim_result_otu_ij$signif))
#每次循环输出图片
#pdf(paste('anosim_two/anosim.', group_name[i], '_', group_name[j], '.pdf', sep = ''), width = 7, height = 5)
png(paste('anosim_two/anosim.', group_name[i], '_', group_name[j], '.png', sep = ''), width = 600, height = 400)
plot(anosim_result_otu_ij, col = c('gray', 'red', 'blue'))
dev.off()
}
}
#带 R 值和 p 值的表格
anosim_result_two <- data.frame(anosim_result_two, stringsAsFactors = FALSE)
names(anosim_result_two) <- c('group', 'distance', 'R', 'P_value')
#可选添加 p 值校正过程，例如 Benjamini 校正
anosim_result_two$P_value <- as.numeric(anosim_result_two$P_value)
anosim_result_two$P_adj_BH <- p.adjust(anosim_result_two$P_value, method = 'BH')
write.table(anosim_result_two, 'anosim_two/ANOSIM.result_two.txt', row.names = FALSE, sep = '\t', quote = FALSE, na = '')
以上示例，在每一步循环中挑选出两个分组，并进行两两分组间的ANOSIM。输出的ANOSIM统计总表中，包含了计算所依据的距离类型（Distance）、ANOSIM分析的R值（R）、显著性p值（P_value）以及BH校正后p值（P_adj_BH）等多项内容。
每步循环的ANOSIM箱线图保存至“anosim_two”中。
其余的图片结果，如“anosim.s1_s2.png”，通过该图我们可知s1和s2采样地间的土壤细菌群落组成无较大差异，或者说组内差异大于了组间差异。具体表现在Between组相对与其他两个分组的秩较低，以及p值不显著（0.122，高于0.05水平）。
类似地，通过“anosim.s2_s4.png”，我们可知s2和s4采样地间的土壤细菌群落组成具有显著差异，即两组间群落差异大于了各组内的差异。具体表现在Between组的秩明显高于其他两个分组的秩，R值较高（0.586），且p值很显著（0.001水平，不过这仅为校正前的p值，可再结合统计表查看校正后的p值是否同样显著）。
即可获知来自5个采样地点的土壤细菌群落中，两两之间的相似/相异水平，以及组间和组内的变异程度，并帮助寻找在排序图中肉眼观察不明显的组间差异。