R语言做置换多元方差分析(PERMANOVA)

Rvegan执行群落结构差异检验之置换多元方差分析(PERMANOVA

 

置换多元方差分析(Permutational multivariate analysis of variancePERMANOVA),又称非参数多因素方差分析(nonparametric multivariate analysis of variance)、或者ADONIS分析,其本质是基于F统计的方差分析,依据距离矩阵对总方差进行分解的非参数多元方差分析方法。使用PERMANOVA可分析不同分组因素对样品差异的解释度,并使用置换检验进行显著性统计。在生态统计中,可以使用PERMANOVA(往往同时配合排序分析来使用,对于PERMANOVA,更常与PCoA排序分析放一起说明问题),查看不同环境的群落组成结构差异是否显著。

Rvegan包是生态统计分析的常用包之一,在本文,以16S扩增子测序所得的细菌群落数据为例,展示Rvegan进行PERMANOVA检验群落结构组成差异的一般过程。

 

首先介绍示例数据。我们此处共有4016S测序样本(细菌群落),均来自土壤。这40个土壤样本共来自5个不同的采样地点,每个地点各包含8个样本。

此处我们希望通过PERMANOVA,查看不同采样地点的土壤细菌群落组成结构是否具有显著的不同。

 

测试数据(data)、示例结果(result)、R脚本等,已上传至百度盘(提取码d1rv):

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

 


示例文件简要


文件“otu_table.txt”为OTU丰度表格,其内容展示如下。

每一列为一个样本,每一行为一种OTU,交叉区域为每种OTU在各样本中的丰度。

1.png

 

文件“bray.txt”为提前计算得到的样本距离矩阵文件(此处展示的是样本间Bray-curtis距离),其内容展示如下。

每一列为一个样本,每一行为一个样本,交叉区域为样本间的Bray-curtis距离(取值范围0-1,越接近于1表明样本间细菌群落组成差异越大)。

3.png


文件“group.txt”为样本分组信息,其内容展示如下。

第一列(names)为各样本名称;第二列(site)为各样本的分组信息,即这些样本所属的采样地点(s1,地点1s2,地点2s3……以此类推)。

注意:“names”列中样本名称的顺序一定要与OTU表或者距离矩阵中的样本名称顺序一致。

2.png

  


使用vegan包进行PERMANOVA检验群落结构差异


首先导入数据至R中。我们可选导入已经计算好的样本距离矩阵文件,也可使用OTU丰度表格文件,同时导入样本分组文件。

备注:当读入数据为距离矩阵时,需要将读取的数据框转化为dist类型,便于后续的函数识别;当读入本示例的OTU表格时,则首先要进行转置操作,即要求所得数据框格式:每一行为一个样本,每一列为物种信息。

##读入文件
#现有的距离矩阵
dis <- read.delim('bray.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
dis <- as.dist(dis)    #将导入的样本间距离转化为 dist 类型
#或者直接使用 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检验整体差异,即查看来自5个不同采样地点所获得的土壤细菌群落结构组成在整体上是否不具备一致性。

在该部分对以下步骤进行较为详细的说明展示。

 

导入vegan包,并使用vegan包中的adonis()命令执行PERMANOVAadonis()的详细说明可使用?adonis()查看。

library(vegan)
 
##PERMANOVA 分析(所有分组间比较,即整体差异)
#(1)若是已经提供好了距离矩阵,则直接使用现有的距离矩阵进行分析即可
adonis_result_dis <- adonis(dis~site, group, permutations = 999) #根据 group$site 这一列样本分组信息进行 PERMANOVA 分析,随机置换检验 999 次
 
#(2)若是使用 OTU 丰度表,则需要在计算时指定所依据的距离类型,这里依然使用 Bray-Curtis 距离
adonis_result_otu <- adonis(otu~site, group, permutations = 999, distance = 'bray')   #这条命令的详情和上述命令所表示的信息一致
 
#(3)或者首先根据丰度表计算样本距离,在将所的距离数据作为输入
dis1 <- vegdist(otu, method = 'bray')
adonis_result_dis1 <- adonis(dis1~site, group, permutations = 999)     #同上所述

使用adonis()命令进行PERMANOVA分析时,指定样本间距离数据(“dis”),以及各样本所属的分组信息(“group ”数据框的“site”列,为了避免识别错误,不要使用纯数字作为分组名称),且要注意分组信息中的分组名称与距离矩阵中样本名称要按顺序对应(需预先排列好二者的顺序)permutations参数用于指定随机置换检验的次数,这里设置置换999次(视数据量而定,数据量越大,置换次数也应当增大)。

若导入数据为OTU丰度表时,则也可以直接将OTU丰度表作为输入数据,此时需要使用distance参数指定所依据计算的距离类型;或者也可以首先使用vegan包中的vegdist()计算样本间距离,然后再将所得距离数据作为adonis()的输入。之后指定各样本所属的分组信息且要注意分组信息中的分组名称与OTU丰度表中样本名称要按顺序对应(需预先排列好二者的顺序),以及置换检验的次数等,进行PERMANOVA

 

以上展示了3种处理过程,由于所使用距离类型相同(均为Bray-Curtis距离),分组信息、置换检验次数(均为999次)等也一致,若置换检验次数已达到饱和,则结果也是一致的。

然后我们查看结果,评估组间群落组成差异。

#查看结果,上述 3 条命令所计算的内容一致,以其中一个为例
adonis_result_dis
#或者
summary(adonis_result_dis)

直接键入统计结果所赋值的对象名称后(例如本示例为“adonis_result_dis”),屏幕输出统计结果简要。包含了运行命令行、置换检验次数、差异计算结果等信息。

4.png

其中,可主要关注中间的列表。site行即为本示例site分组的组间差异检验结果,Residuals为残差,TotalsiteResiduals的加和。

对于各项统计内容:Df,自由度,其值=所比较的分组数量-1SumsOfSqs,即Sums of squares,总方差,又称离差平方和;MeanSqs,即Mean squares,均方(差);F.ModelF检验值;R2,即Variation (R2),方差贡献,表示不同分组对样品差异的解释度,即分组方差与总方差的比值,R2越大表示分组对差异的解释度越高;Pr (>F),显著性p值,小于0.05说明差异显著。

对于各部分运算细节,可使用summary()命令查看结果“adonis_result_dis”中所含内容,并提取相关结果查看。

5.png

将主要的统计结果转化为数据框的类型,并输出至文件。

#可选输出
otuput <- data.frame(adonis_result_dis$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 = '')

示例输出结果文件“PERMANOVA.result_all.txt”的内容如下,包含了自由度(Df)、平方和(Sums of squares)、均方差(Mean squares)、F检验值(F.Model)、方差贡献(Variation (R2))、显著性p值(Pr (>F))等多项重要内容。

6.png

 

据此,回顾本文开始的问题,我们可知5个采样地点的土壤细菌群落结构在整体水平上是不一致的。

 

两组间比较展示一个批量处理循环流程)


我们已知在整体水平上,组间群落组成是存在差异的。那么,究竟是哪些之间的差异所导致的呢?在这些物种组成具有显著差异的不同群落中,哪些群落之间的差异较小,而哪些群落之间具有更大的差异呢?

此外,有时候我们也并不关注整体水平,而更多是想关注某几个特定分组之间的群落差异是否显著,这时候也需要我们将这些特定分组的样本挑选出,单独进行比较。

 

接下来我们继续使用示例数据,使用PERMANOVA探求两两环境(两两采样地)之间的土壤细菌群落结构组成的差异,包括是否具有显著性,以及组间和组内的变异程度等。

在所比较的分组较多的情况下,可使用循环语句的方式逐一选择数据并处理(在该部分提供一个批处理循环示例作为参考),以节约手动一个个运行的时间。

在这里推荐使用OTU丰度表作为输入数据,每次筛选分组后重新计算样本距离;而不是读入现有的距离矩阵后,在矩阵中筛选样本行列。这样可避免由于样本数减少可能导致的距离变动而造成的误差(有时候因样本数的减少,导致了实际存在的物种数也相应地减少,这样重新计算的样本间距离与在之前的相比,可能会有不同)。

##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')     #随机置换检验 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)')

以上示例,在每一步循环中挑选出两个分组,并进行两两分组间的比较。读入数据OTU表,筛选出对应的样本数据并作为输入,距离类型使用Bray-Curtis距离,999次置换检验。在每次的统计结果中提取出带有方差贡献、F检验值、显著性p值等主要统计信息的列表数据,并将其中的组间差异检验结果那一行与前一个结果按行合并。

最后所得PERMANOVA结果数据框“adonis_result_two”的内容如下。各列的统计学意义与上文一致。

7.png

 

我们可选依据显著性p值水平,在结果中添加“*”标记后,输出统计结果。如下所示。

#可选标记 “*” 显著性
for (i in 1:nrow(adonis_result_two)) {
    if (adonis_result_two[i, 'Pr (>F)'] <= 0.001) adonis_result_two[i, 'Sig'] <- '***'
    else if (adonis_result_two[i, 'Pr (>F)'] <= 0.01) adonis_result_two[i, 'Sig'] <- '**'
    else if (adonis_result_two[i, 'Pr (>F)'] <= 0.05) adonis_result_two[i, 'Sig'] <- '*'
}
 
#输出
write.table(adonis_result_two, 'PERMANOVA.result_two.txt', row.names = FALSE, sep = '\t', quote = FALSE, na = '')

两组间细菌群落PERMANOVA分析输出结果“PERMANOVA.result_two.txt”的内容如下,与上文结果的内容一致,包含了计算所依据的距离类型(Distance)、自由度(Df)、平方和(Sums of squares)、均方差(Mean squares)、F检验值(F.Model)、方差贡献(Variation (R2))、显著性p值(Pr (>F))等多项重要内容。

8.png

 

我们即可得到来自5个采样地点的土壤细菌群落中,两两之间的差异水平,并可从中找出具有显著变异的细菌群落样本。

通过显著性p值,可查看两组间差异是否显著;并通过R2,评估各分组或环境因素对土壤细菌群落差异的解释度等。

 

补充内容


在进行PERMANOVA检验群落结构差异的同时,常结合排序分析(如PCAPCoANMDS,对于PERMANOVA更推荐配合PCoA)的可视化展示,一同说明问题,使结果更具说服力。

通过排序图,通过观察群落样本的坐标分布,可判断不同分组之间是否具有明显区分(组间差异),同一分组样本之间的离散程度(组内差异),哪些分组在排序图中所处的区域较为一致(组间相似),哪些分组与其它分组相比在图中单独位于一个较远的区域(差异较大的分组),哪些样本明显偏离了平均位置(离群点),各排序轴主要与那些环境因素存在潜在的相关性(群落差异主要受到那些环境的影响)及各排序轴的解释量高低(模型解释度)等。

再结合PERMANOVA的检验结果,验证以上观测,做出总结。

 


参考文档


vegan包中对PERMANOVA的说明

Permutational Multivariate Analysis of Variance Using Distance Matrices

 



链接地址:http://blog.sciencenet.cn/blog-3406804-1169075.html
发布了3 篇原创文章 · 获赞 28 · 访问量 4万+
展开阅读全文

没有更多推荐了,返回首页

©️2019 CSDN 皮肤主题: 大白 设计师: CSDN官方博客

分享到微信朋友圈

×

扫一扫,手机浏览