科研星球

R语言统计与绘图:通过相似性或相异指数的数值分布比较群落Beta多样性高低

在基于高通量测序的微生物群落分析中,若提到如何描述不同群落β多样性是否存在差异或者评估组内或组间的差异程度,我们通常可以想到很多方法,这些都是描述群落β多样性特征的常见分析。例如,一般我们首先会基于微生物分类群丰度数据计算不同群落样本间的成对相似或相异矩阵,然后通过PCAPCoA等观察不同组的群落样本是否沿特征轴具有明显的分界,或者通过聚类分析观察不同组的群落样本是否划分到不同的聚类簇中,以及通过一些统计检验方法如PERMANOVA等确定不同群落在物种组成方面的差异显著性,等等……


那么,如果想通过上述提到的相似或相异矩阵,去比较几组群落β多样性的高低,该如何实现呢?因为通常无论相似性或相异指数,都是在反映两两群落样本之间的差异程度,又该如何通过它们描述整体群落的β多样性的高低水平呢?

  

通过相似性或相异指数的数值分布确定Beta多样性高低


  

既然这种样本间的成对相似或相异矩阵是反映群落β多样性的一种度量,那么它的数值差异肯定也可以用于描述β多样性的高低,例如以下两篇文献中这样。


下载.jpeg


Zhang等(2018)比较了中国51个大豆农田的土壤细菌群落、根际细菌群落以及根内细菌群落β多样性。基于Bray-Curtis相异指数的β多样性分析表明,根际细菌群落的β多样性最高,因为该组Bray-Curtis相异指数的数值分布相对其它两组更高。

下载 (1).jpeg

Jiao等(2018)根据300个土壤样品的Bray-Curtis相异矩阵,估计了细菌、古细菌和真菌群落之间的β多样性差异。根据各组Bray-Curtis相异指数的数值分布高低,可判断真菌群落具有最高的β多样性,古细菌群落次之,细菌群落的β多样性最低。


当然这种方法有个前提,就是每组中包含多个群落样本,这样就获得了各组内群落之间的相似性或相异指数分布。之后再通过不同组内的相似性或相异指数分布直接比较不同组间数值的高低差异,就可以确定β多样性的高低水平了。

 

那么,为什么可以这样比较呢?鱼老弟开始也是很迷惑,就去问师兄,师兄给解释了一番后才想明白。简单来说,在定义上β多样性反映了在一个梯度上从一个生境到另一个生境所发生的种的多样性变化的速率和范围,它是研究群落之间的种多度关系。所以可知,β多样性越高就代表了物种多样性改变程度越大。而Bray-Curtis相异指数越大,则意味者群落之间物种组成的相似性越低,因此也就可以近似理解,一组更高的Bray-Curtis相异指数分布就意味着该组群落的β多样性更高。

类似地,更换为其它类型的与β多样性有关的相似性或相异指数也可以用来描述类似的情形。

   


一个通过相异指数的数值分布比较Beta多样性高低的分析和作图示例


  

接下来,就模仿上述两篇文献,展示这种通过相似性或相异指数确定不同组群落β多样性高低的方法在R语言中的计算和作图过程。

下文中所使用的示例数据和R代码的百度盘链接(提取码,iw17):

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

若百度盘失效,也可在GitHub的备份中获取:

https://github.com/lyao222lll/sheng-xin-xiao-bai-yu-2021

  


示例数据



示例数据共涉及了48个由16S测序所得的细菌群落数据,这48个样本分别来自3个不同的环境(Env1Env2Env3),每个环境下各包含16个样本。

文件“otu_table.txt”记录了所有环境样本的细菌群落丰度(OTU水平)。每一列为一个样本,每一行为一种OTU,交叉区域为每种OTU在各样本中的丰度。

下载 (2).jpeg

文件“group.txt”记录了48个样本分组信息,其中group1列为样本来源的3个不同的环境分组(Env1Env2Env3),group2列不重要请无视吧……


下载 (3).jpeg


  


计算与Beta多样性有关的相异指数



我们期望观察和比较三种不同的环境中(Env1Env2Env3),微生物群落的β多样性的高低。为此,我们首先可以基于物种丰度数据计算群落相似性或相异指数,因为它们是常用于反映群落β多样性的度量值。

对于常见的群落相似性或相异指数的定义以及计算方法,可参考前文“β多样性、群落相似性与距离”。这里以微生物数据中比较常用的Bray-Curtis距离为例,获取距离矩阵。

otu <- t(otu)
 
#计算与 Beta 多样性有关的群落相异指数,例如使用 vegan 包计算 Bray-curtis 距离,详情加载 vegan 包后 ?vegdist
dis <- vegan::vegdist(otu, method = 'bray')
 
#以矩阵形式输出
dis <- as.matrix(dis)
write.table(dis, 'Bray-curtis.txt', sep = '\t', col.names = NA, quote = FALSE)

下载 (4).jpeg

输出文件“Bray-Curtis.txt”记录了所有环境样本之间成对的Bray-Curtis相异矩阵。

备注:这里直接通过上述提供的物种丰度表计算获得Bray-Curtis相异矩阵。如果您希望计算其它类型的相似性或相异指数,如微生物群落中另一常见的Unifrac相异指数,也是可以的,具体使用哪种相似性或相异指数根据实际情况选择。此外,当然如果您在一开始就准备好了这个样本间成对的相似性或相异矩阵文件(而不再通过分类群的丰度表从头计算),也可以直接使用它。

  


绘制箱线图初步观察群落Beta多样性高低水平



总之,我们将已经准备好的,代表群落β多样性的相异矩阵文件读入到R中,开始作图并进行统计计算。

由于我们期望观察和比较三种不同的环境中(Env1Env2Env3),微生物群落的β多样性的高低,因此在差异分析之前,不妨首先绘制一个箱线图进行探索性的分析,大致观察一下Bray-Curtis相异指数的整体分布。

 
#读取样本分组信息
group <- read.delim('group.txt', stringsAsFactors = FALSE)
 
##例如,比较 Env1、Env2、Env3 三组之间,群落的 Beta 多样性差异
#根据分组获得组内距离矩阵
env1 <- subset(group, group1 == 'Env1')$samples
dis_env1 <- dis[env1,env1]
 
env2 <- subset(group, group1 == 'Env2')$samples
dis_env2 <- dis[env2,env2]
 
env3 <- subset(group, group1 == 'Env3')$samples
dis_env3 <- dis[env3,env3]
 
#将矩阵转化为向量,以便用于作图和统计
dis_env1 <- as.vector(as.dist(dis_env1))
dis_env2 <- as.vector(as.dist(dis_env2))
dis_env3 <- as.vector(as.dist(dis_env3))
 
#构建作图数据集
dat <- data.frame(
    dis = c(dis_env1, dis_env2, dis_env3),
    group = factor(c(
        rep('Env1', length(dis_env1)), 
        rep('Env2', length(dis_env2)), 
        rep('Env3', length(dis_env3))
    ), levels = c('Env1', 'Env2', 'Env3'))
)
 
#使用 ggplot2 绘制各组内 Bray-curtis 距离指数分布的箱线图
library(ggplot2)
 
p <- ggplot(dat, aes(group, dis)) +
geom_boxplot(aes(fill = group), width = 0.6) +
scale_fill_manual(values = c('#CD5B45', '#228B22', '#00688B')) +
theme(panel.grid = element_blank(), panel.background = element_blank(), 
    axis.line = element_line(colour = 'black'), legend.position = 'none') +
labs(x = NULL, y = 'Bray-Curtis dissimilarity\n')
 
p

下载 (5).jpeg

根据箱线图中展示的Bray-Curtis相异指数的数值分布,我们可以初步评估,似乎环境Env2中的微生物群落的β多样性更高,而环境Env1中的微生物群落的β多样性更低。

  


通过统计检验确定群落Beta多样性高低差异



当然仅使用眼睛判断肯定不靠谱,需通过某些统计手段确定β多样性的高低水平。包括比较三组环境中微生物群落β多样性的整体差异,再分别比较两两环境的差异。我本来想直接使用方差分析,但是在评估时发现数据违反了等方差假设,不满足方差分析的条件,因此采用了非参数方法:Kruskal-Wallis检验进行三组的整体比较,并进而通过Wilcoxon秩和检验继续比较两两差异。


#三组的整体差异分析,使用 Kruskal-Wallis Test 执行,详情 ?kruskal.test
kruskal.test(dis~group, data = dat)
 
#如果整体显著再进行两两分组的比较,使用 Wilcoxon 秩和检验执行双侧检验,详情 ?wilcox.test
wilcox.test(dis_env1, dis_env2, alternative = 'two.sided')
wilcox.test(dis_env2, dis_env3, alternative = 'two.sided')
wilcox.test(dis_env1, dis_env3, alternative = 'two.sided')

下载 (6).jpeg

根据不同组Bray-Curtis相异指数的数值比较,我们可判断三种不同的环境中(Env1Env2Env3)微生物群落β多样性的高低水平存在显著不同,包括两两环境之间亦是如此。

那么,组间差异确定了,怎样再确定哪种环境中微生物群落β多样性更高呢?回顾t检验方差分析中,我们可以通过均值的高低去判断;而对于非参数的Wilcoxon秩和检验方法而言,则不是比较的均值,而是中位数(《环境与生态统计:R语言的应用》第一版中文译本,第68页提到“Wilcoxon秩和检验检验的是两组样本的中位数是不是相等”)。

#考虑到 Wilcoxon 秩和检验体现了中位数的差异,因此计算三组数据的中位数以评估 Beta 多样性的高低水平
median(dis_env1)
median(dis_env2)
median(dis_env3)

下载 (7).jpeg

综合上述统计结果,可知三种不同的环境中(Env1Env2Env3),两两环境之间的微生物群落β多样性的高低水平均存在显著不同,且反映β多样性的Bray-Curtis相异指数的中位数Env2>Env3>Env1。最后确定了结论,环境Env2中的微生物群落的β多样性最高,而环境Env1中的微生物群落的β多样性最低。

#基于上述统计结果,判断好组间差异后,将差异分析结果添加到箱线图中
p +
annotate('text', label = 'Kruskal-Wallis Test', x = 1, y = 0.56, size = 3) +
annotate('text', label = sprintf('italic(P) < %.3f', 0.001), x = 1, y = 0.53, size = 3, parse = TRUE) +
annotate('text', label = 'c', x = 1, y = max(dis_env1)+0.05, size = 3) +
annotate('text', label = 'a', x = 2, y = max(dis_env2)+0.05, size = 3) +
annotate('text', label = 'b', x = 3, y = max(dis_env3)+0.05, size = 3)

下载 (8).jpeg

如此,三种不同的环境中(Env1Env2Env3),微生物群落β多样性的高低水平一目了然。

   


* 关于其它的Beta多样性指数


  

除此之外,β多样性也可以计算一个具体的数值(β多样性指数),类似α多样性指数那样,可以直接用于比较数值大小。师兄给的推荐,LegendreDe Caceres2013)总结了多种度量β多样性指数的方法,大家有兴趣可以阅读大佬的文章。


   


参考文献


  


钱松环境与生态统计:R语言的应用(曾思育 译)高等教育出版社, 2011.
Baselga A. Multiple site dissimilarity quantifies compositional heterogeneity among several sites, while average pairwise dissimilarity may be misleading. Ecography, 2013, 36(2): 124-128.
Jiao S, Chen W, Wang J et al. Soil microbiomes with distinct assemblies through vertical soil profiles drive the cycling of multiple nutrients in reforested ecosystems. Microbiome, 2018, 6: 146.
Legendre P, De Cáceres M, Beta diversity as the variance of community data: dissimilarity coefficients and partitioning. Ecol Lett, 2013, 16: 951-63.
Zhang B, Zhang J, Liu Y et al. Biogeography and ecological processes affecting root-associated bacterial communities in soybean fields across China. Sci Total Environ, 2018, 627: 20-27.



没有账号?