实例教程 | 功能多样性三元图的复现代码

文摘   科学   2024-10-11 09:00   江苏  
在生态学家工具箱中的众多多样性指数中,可划分为加性项的指数尤其有用,因为不同的组成部分可与塑造群落结构的不同生态过程相关联。Ricotta2023年在《Methods in Ecology and Evolution》期刊上提出了一种新的加性多样性分解法,将给定群落的多样性结构划分为三个互补部分:功能多样性、功能冗余和物种优势。这三个部分的总和为1。因此,它们可以用来在三元图中描绘群落结构。有关功能多样性三元图及其划分原理的更多细节可查看上一篇推文(Methods in Ecology and Evolution | 功能多样性的三元图),本篇教程将尝试复现该文的结果,以期为大家提供一个绘制功能多样性的三元图的工作流程示例。教程中的数据和完整代码可从github下载(https://github.com/biodiversity-monitoring/WeChat/tree/main/code/R/FDTD)。

  • 划分方法

长期以来,三元图一直被用于地球科学领域,以显示总和为1100%的三个成分的比例。在生态学中,Grime19741977)引入了三元图,根据植物的主要适应策略将其分为竞争物种、耐胁迫物种和杂草物种(CSR)。同样,Ricotta提出了一种基于三个互补成分的群落多样性三角表示法:Rao功能多样性(Q)、功能冗余(R)和辛普森优势度(D)。由于所有这三个多样性成分都可以用功能相似性/差异来解释,因此用DRQ的分数来表示特定群落,反映了群落中所有物种之间的相似性结构。

DRQ图中,三角形的顶点对应于三个分量DRQ,每个地点用一个点表示,其位置由三个加性分量的实际值决定。三角形的每个角都代表一种情况,即一个分量的值等于1,另外两个分量的值为0。每个分量的值随着与相应角的距离增加而线性减少。例如,当D值较大时,点会靠近D角,反映出研究地点物种优势度高、功能多样性和冗余度低的情况。同样,如果一个点靠近R角,则相应地点的功能冗余度非常高。最后,如果一个点靠近Q角,则该地点的功能多样性高,冗余度和优势度低。值得注意的是,三元图中的三个多样性组成部分所捕捉的变异类型各不相同。Q侧重于个体水平的功能变异,R仅考虑种间功能变异,而D则考虑物种多度的变异。
具体的推导公式可以查看上一篇译文,这里不再赘述。
  • 示例数据

示例数据来自于原文的附录,我对其他进行了简单整理,分为两个.csv文件(Species_abundances_and_functional_traits_in_grazed_plots.csvSpecies_abundances_and_functional_traits_in_ungrazed_plots.csv)。其包括了托斯卡纳(意大利中部)8块放牧和7块非放牧草地植物群落的实际多度数据,以及比叶面积、叶干物质含量、叶氮和碳含量4个功能性状数据。

  • 计算代码

在开始计算前,我们需要先设置好路径,并将示例数据放到路径下。

# 设置路径setwd("...")getwd()# 导入相应的R包library(adiv)library(ggtern)
接着导入示例数据。
# 导入数据grazed_data <- read.csv("Species_abundances_and_functional_traits_in_grazed_plots.csv")ungrazed_data <- read.csv("Species_abundances_and_functional_traits_in_ungrazed_plots.csv")# 查看数据结构head(grazed_data)##                      Species   G1   G2 G3 G4 G5 G6   G7   G8   SLA  LDMC     N     C## 1        Acer monspessulanum  0.0  0.5  0  0  1  0  0.0  0.0 0.134 0.906 0.471 0.808## 2            Acinos arvensis  0.5  2.0  1  0  0 16 28.0  0.5 0.074 0.511 0.082 0.745## 3      Arrhenatherum elatius  0.0  1.0  0  0  0  0  0.0  0.0 0.683 0.503 0.659 0.580## 4             Bromus erectus 10.0 11.0 17 32  9 11  7.0  9.0 0.299 0.587 0.403 0.547## 5            Bromus sterilis  0.0  0.5  0  0  0  0  0.0  1.0 0.959 0.237 0.666 0.644## 6       Bupleurum baldense**  1.0  1.0  0  0  0  0  0.0  0.0 0.394 0.535 0.109 0.463head(ungrazed_data)##                     Species    N1 N2   N3   N4 N5 N6    N7   SLA  LDMC     N     C## 1       Acer monspessulanum  60.0  0  0.0  1.0  0  0   0.0 0.131 0.915 0.269 0.868## 2           Acinos arvensis   0.0  0  0.0  0.0  2  4   0.0 0.040 0.475 0.000 0.722## 3     Arrhenatherum elatius   0.5  0  0.5  0.0  0  3   0.0 0.849 0.313 0.603 0.384## 4            Bromus erectus   0.5 18 86.0 61.0 58 65   2.0 0.039 0.757 0.297 0.540## 5           Bromus sterilis   0.5  4  0.0  0.0  0  0   0.0 0.385 0.317 0.495 0.391## 6      Bupleurum baldense**   0.0  0  0.5  1.0  2  2   0.0 0.394 0.535 0.109 0.463
从上面的数据结构中我们发现群落数据和性状数据都在同一个数据框中,因此我们需要进一步其拆分,以方便后面的计算。
# 去掉NA值grazed_data <- na.omit(grazed_data)ungrazed_data <- na.omit(ungrazed_data)# 提取群落数据com1 <- grazed_data[, c("G1", "G2", "G3", "G4", "G5", "G6", "G7", "G8")]rownames(com1) <- grazed_data$Speciescom1 <- t(com1)com1com2 <- ungrazed_data[, c("N1", "N2", "N3", "N4", "N5", "N6", "N7")]rownames(com2) <- ungrazed_data$Speciescom2 <- t(com2)com2# 提取性状数据traits1 <- grazed_data[, c("SLA", "LDMC", "N", "C")]rownames(traits1) <- grazed_data$Speciestraits1traits2 <- ungrazed_data[, c("SLA", "LDMC", "N", "C")]rownames(traits2) <- ungrazed_data$Speciestraits2

数据准备好后,就可以开始计算功能多样性的各个成分了

功能多样性的计算参照原文附录中的说明:所有性状数据均按其最小值和最大值在[0,1]范围内进行线性缩放。根据缩放值计算两种处理中所有物种对之间的欧氏功能距离。然后将功能距离除以两种处理中的最大值。所有计算均使用R软件包“adiv”的函数“speciesdiv”(用于计算辛普森多样性)和“QE”(用于计算Rao’s Q)进行。

根据这个流程,我封装了一个函数来完成功能多样性三个组分的计算:

FD_DRQ <- function(traits, com) {  # 线性缩放性状数据到[0,1]范围  traits_scaled <- apply(traits, 2, function(x) (x - min(x)) / (max(x) - min(x)))  # 计算物种间的功能欧几里得距离  dist_func <- dist(traits_scaled, method = "euclidean")  # 将距离除以最大值进行缩放  dist_func_scaled <- dist_func / max(dist_func)  # 使用adiv包计算Simpson优势度(D)、功能冗余(R)Rao's Q  S <- speciesdiv(com, "GiniSimpson")  colnames(S) <- "S"  Q <- QE(com, dist_func_scaled, "QE")  colnames(Q) <- "Q"  R <- S - Q   colnames(R) <- "R"  D <- 1 - S  colnames(D) <- "D"  # 创建DRQ数据框  DRQ <- data.frame(D = D, R = R, Q = Q)  return(DRQ)}
接下来只需要调用FD_DRQ函数即可,三元图的绘制这里我采用了ggtern包来实现。

FD_DRQ(traits1, com1)## D R Q## G1 0.1862207 0.4914946 0.3222847## G2 0.1305590 0.5057338 0.3637071## G3 0.2122165 0.3986428 0.3891407## G4 0.1564776 0.4990360 0.3444864## G5 0.1167779 0.5095429 0.3736792## G6 0.1355194 0.4800773 0.3844033## G7 0.1079387 0.4404172 0.4516440## G8 0.1598319 0.5786248 0.2615433FD_DRQ(traits2, com2)## D R Q## N1 0.5548090 0.2752250 0.1699660## N2 0.2117788 0.5008637 0.2873574## N3 0.2738727 0.4316556 0.2944718## N4 0.3571282 0.4149361 0.2279357## N5 0.2606493 0.4696835 0.2696672## N6 0.2490568 0.5032315 0.2477117## N7 0.5371157 0.2976231 0.1652613 DRQ <- FD_DRQ(traits1, com1)ggtern(data = DRQ, aes(x = D, y = R, z = Q)) +geom_point(shape=21,size=2,fill="red") +theme_custom() +labs(x = "D", y = "R", z = "Q")
最后,我们也可以把两个数据集绘制在一起,做一个比较。
# 合并两个数据集data1 <- FD_DRQ(traits1, com1)data1$group <- "Grazed"data2 <- FD_DRQ(traits2, com2)data2$group <- "Ungrazed"data <- rbind(data1, data2)data
# 创建三元图ggtern(data = data, aes(x = D, y = R, z = Q, fill = group, shape = group)) + geom_point(size = 2) + theme_custom() + labs(x = "D", y = "R", z = "Q") + scale_fill_manual(values = c("Grazed" = "red", "Ungrazed" = "blue")) + scale_shape_manual(values = c("Grazed" = 21, "Ungrazed" = 22)) + theme(legend.background = element_blank(), legend.key = element_blank())
如果还想进一步美化成图以更接近原图(如标签、图例、配色等),可以选择导出eps格式的矢量图在Adobe Illustrator中进行修改。
  • 参考来源

  1. Ricotta, C., Podani, J., Schmera, D., Bacaro, G., Maccherini, S., & Pavoine, S. (2023). The ternary diagram of functional diversity. Methods in Ecology and Evolution, 14, 1168–1174. https://doi.org/10.1111/2041-210X.14100

  2. Ricotta, C., Podani, J., Schmera, D., Bacaro, G., Maccherini, S., & Pavoine, S. (2023). Data from: The ternary diagram of functional diversity. Dryad Digital Repository. https://doi.org/10.5061/dryad.7pvmcvdzc

Biodiversity Monitoring
生物多样性;监测保护;群落生态;生态统计;R语言;python。 主要分享一些前沿的文献和方法实例,更新看心情和时间。
 最新文章