R 语言作图（例如 ggplot2）时怎么添加物种图片（剪影）？又见之前的一篇推送，

👉 临摹 Ecology Letter 文章的系统发育树

（点击 👉 即可跳转。）

1 包和数据

加载包，

library(tidyverse)
library(ggplot2)
library(ggforce)
library(vegan)
library(indicspecies)
library(BiodiversityR)
library(rphylopic)

数据来自以下论文，

👉 Big rodents disperse small seeds and spores in Neotropical wetlands

All data supporting the results of this paper are available via figshare: https://doi.org/10.6084/m9.figshare.25742442 (Hoffmann et al., 2024).

读取 .csv 数据，

# 数据
dat <- read.csv("Data_Big Rodents.csv", header = TRUE)

# 若无法下载以上数据# 可使用此虚拟数据# 去掉注释运行即可# set.seed(1)
# n1 <- 39
# n2 <- 43
# plant <- as.data.frame(
#   sapply(
#     paste0("plant", 1:27),
#     function(i) {
#       c(
#         sample(x = c(rep(0, 70), rpois(30, 2)), size = n1),
#         sample(x = c(rep(0, 70), rpois(30, 6)), size = n2)
#       )
#     },
#     simplify = TRUE
#   )
# )
# samples <- paste0("site", 1:n)
# species <- c(rep("capybara", n1), rep("nutria", n2))
# dat <- cbind(samples, species, plant)

# 查看列名
colnames(dat)

 [1] "samples"                 "species"                
 [3] "alternanthera"           "amaranthus"             
 [5] "spirodela_intermedia"    "hydrocotyle"            
 [7] "eclipta_prostrata"       "stellaria_media"        
 [9] "cyperus_squarrosus"      "cyperus_odoratus"       
[11] "cyperus_cf_brevifolius"  "cyperus_cf_iria"        
[13] "cyperus"                 "eleocharis_cf_maculosa" 
[15] "fimbristylis_autumnalis" "juncus_cf_capitatus"    
[17] "juncus_cf_bufonius"      "juncus_cf_tenuis"       
[19] "marsilea_cf_ancylopoda"  "ludwigia_leptocarpa"    
[21] "ludwigia_peploides"      "eragrostis_hypnoides"   
[23] "eragrostis"              "cynodon"                
[25] "lolium"                  "bacopa_monnieri"        
[27] "azolla"                  "boehmeria_cylindrica"   
[29] "salvinia"

# 显示前 5 列，前 10 行
dat[1:10, 1:5]

   samples  species alternanthera amaranthus spirodela_intermedia
1    CA1_1 capybara             0          0                    0
2    CA1_3 capybara             0          0                    0
3    CA1_4 capybara             3          0                    0
4    CA1_5 capybara             0          0                    0
5    CA1_6 capybara             0          0                    0
6    CA1_7 capybara             0          1                    0
7    CA1_9 capybara             0          0                    0
8   CA1_10 capybara             0          0                    0
9    CA2_1 capybara             0          0                    0
10   CA2_2 capybara             0          0                    0

1.1 选择植物数据

去掉分类变量的前两列，

plant.dat <- dat[, -c(1:2)]

1.2 指示物种

# 函数
?indicspecies::multipatt()

# Runs the combination analysis
set.seed(1)
mpt <- indicspecies::multipatt(
  plant.dat,
  dat$species,
  control = how(nperm = 999)
)

# Species with significant association to one combination
summary(mpt)

 Multilevel pattern analysis
 ---------------------------

 Association function: IndVal.g
 Significance level (alpha): 0.05

 Total number of species: 27
 Selected number of species: 10 
 Number of species associated to 1 group: 10 

 List of species associated to each combination: 

 Group capybara  #sps.  3 
                        stat p.value   
juncus_cf_bufonius     0.506   0.002 **
juncus_cf_capitatus    0.480   0.002 **
cyperus_cf_brevifolius 0.358   0.014 * 

 Group nutria  #sps.  7 
                        stat p.value    
salvinia               0.953   0.001 ***
eclipta_prostrata      0.550   0.001 ***
marsilea_cf_ancylopoda 0.506   0.001 ***
lolium                 0.457   0.005 ** 
azolla                 0.431   0.005 ** 
cyperus_odoratus       0.425   0.026 *  
hydrocotyle            0.404   0.032 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# 筛选
mpt.sign2 <- filter(mpt$sign, p.value < 0.05)
(s.capybara2 <- rownames(mpt.sign2[mpt.sign2$s.capybara == 1, ]))

[1] "cyperus_cf_brevifolius" "juncus_cf_capitatus"    "juncus_cf_bufonius"

(s.nutria2 <- rownames(mpt.sign2[mpt.sign2$s.nutria == 1, ]))

[1] "hydrocotyle"            "eclipta_prostrata"      "cyperus_odoratus"      
[4] "marsilea_cf_ancylopoda" "lolium"                 "azolla"                
[7] "salvinia"

2 PCoA

需要说明的是，此处的 PCoA 分析结果与原论文不完全一致（但是结果的性质没改变）。

其实 MDS 就是 PCoA。

# 拟合 PCoA 的函数
?cmdscale()
# cmdscale(d, k = 2, ...)
# d - a distance structure
# k - the dimension of the space
# ...

plant.dat2 <- plant.dat
# 排除空行
# plant.dat2 <- plant.dat[-which(rowSums(plant.dat) == 0), ]

# Bray–Curtis 距离
# 函数 vegdist 来自 vegan 包
d <- vegdist(plant.dat2, method = "bray")

# 拟合 PCoA
pcoa <- stats::cmdscale(d = d, k = 2, eig = TRUE)

# 各轴解释率
(Dim1.pct <- pcoa$eig[1] / sum(pcoa$eig))

[1] 0.2630981

(Dim2.pct <- pcoa$eig[2] / sum(pcoa$eig))

[1] 0.1711436

# 添加植物物种坐标点
plant.scores <- as.data.frame(BiodiversityR::add.spec.scores(pcoa, plant.dat2)$cproj)
plant.scores$plant.nms <- rownames(plant.scores)

3 可视化

3.1 整理数据

# 提取 PCoA 坐标
Dim1 <- vegan::scores(pcoa)[, 1]
Dim2 <- vegan::scores(pcoa)[, 2]

# 整理到一个数据框中
species <- dat$species
species.clrs <- ifelse(species == "capybara", "#67b715", "#1b81e5")
PCoA.df <- data.frame(Dim1, Dim2, species, species.clrs)

3.2 ggplot2 可视化

# 绘图
FONT <- "Myriad Pro"
CEX <- 18

XLIM <- c(-0.6, 0.6)
YLIM <- c(-0.6, 0.3)

p0 <- ggplot(
  data = PCoA.df,
  aes(Dim1, Dim2, color = species)
) +
  geom_point() +
  ggforce::geom_mark_ellipse(
    expand = unit(0, "mm"),
    show.legend = FALSE
  ) +
  geom_text(
    data = filter(plant.scores, plant.scores$plant.nms %in% c(s.capybara2, s.nutria2)),
    aes(Dim1, Dim2, label = plant.nms),
    size = 4,
    family = FONT,
    inherit.aes = FALSE
  ) +
  scale_color_manual(
    name = NULL,
    values = unique(species.clrs)
  ) +
  scale_y_continuous(
    breaks = round(seq(-0.6, 0.3, 0.3), 1),
    expand = expansion(mult = c(0, 0.06))
  ) +
  labs(
    x = sprintf("PCoA 1 (%.1f%%)", Dim1.pct * 100),
    y = sprintf("PCoA 2 (%.1f%%)", Dim2.pct * 100)
  ) +
  coord_fixed(
    xlim = XLIM,
    ylim = YLIM,
    clip = "off"
  ) +
  theme_classic(
    # theme_bw(
    base_family = FONT,
    base_size = CEX,
    # base_rect_size = 0.4,
    base_line_size = 0.4
  ) +
  theme(
    panel.grid = element_blank(),
    axis.title = element_text(size = CEX),
    axis.text = element_text(size = CEX, color = "black"),
    axis.ticks.length = unit(8, "pt"),
    legend.position = c(0, 0),
    legend.justification = c(0, 0),
    legend.background = element_blank(),
    legend.text = element_text(size = CEX)
  )

Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.

p0

# ggview::ggview(p0, width = 6, height = 5, device = "png", dpi = 900)
ggsave(
  "poca-p0.png",
  p0,
  width = 6,
  height = 5,
  dpi = 900
)

3.3 添加动物形象

从 phylopic.org 网站上获取剪影，

# 可以到网站上获取 uuid
# 也可以使用 get_uuid 函数获取 uuid
# get_uuid(name = "Hydrochoerus hydrochaeris") # capybara
# get_uuid(name = "Myocastor coypus") # "nutria" also known as "coypu"

# 或者
# 直接使用拉丁名
p1 <- p0 +
  add_phylopic(
    # uuid = "097336ef-db5f-4edc-a2e4-2c29a10aa512",
    name = "Hydrochoerus hydrochaeris",
    x = 0.3,
    y = 0.1,
    ysize = 0.4,
    fill = unique(species.clrs)[1],
    alpha = 0.6,
    verbose = TRUE
  ) +
  add_phylopic(
    # uuid = "5ee3b11f-9cd8-4aaf-a59e-ac1ddc18d5e8",
    name = "Myocastor coypus",
    x = -0.3,
    y = 0.1,
    ysize = 0.3,
    fill = unique(species.clrs)[2],
    alpha = 0.6,
    verbose = TRUE
  ) +
  theme(
    legend.position = "none"
  )
p1

Organism silhouettes are from PhyloPic (https://www.phylopic.org/; T. Michael Keesey, 2023) and were added using the rphylopic R package ver. 1.4.0 (Gearty & Jones, 2023). Silhouette was made by Steven Traver, 2012 (CC0 1.0). Silhouette was contributed by Steven Traver.

Organism silhouettes are from PhyloPic (https://www.phylopic.org/; T. Michael Keesey, 2023) and were added using the rphylopic R package ver. 1.4.0 (Gearty & Jones, 2023). Silhouette was made by Mattia Menchetti, 2014 (CC0 1.0). Silhouette was contributed by Mattia Menchetti.

# ggview::ggview(p1, width = 5.5, height = 4.5, device = "png", dpi = 900)
ggsave(
  "poca-p1.png",
  p1,
  width = 6,
  height = 5,
  dpi = 900
)

Organism silhouettes are from PhyloPic (https://www.phylopic.org/; T. Michael Keesey, 2023) and were added using the rphylopic R package ver. 1.4.0 (Gearty & Jones, 2023). Silhouette was made by Steven Traver, 2012 (CC0 1.0). Silhouette was contributed by Steven Traver.
Organism silhouettes are from PhyloPic (https://www.phylopic.org/; T. Michael Keesey, 2023) and were added using the rphylopic R package ver. 1.4.0 (Gearty & Jones, 2023). Silhouette was made by Mattia Menchetti, 2014 (CC0 1.0). Silhouette was contributed by Mattia Menchetti.

4 基础绘图可视化 base plot

4.1 绘图 1

标记 species 名称，按照 species 着色，

# 绘图参数
par(
  mar = c(4, 4, 1, 1),
  family = FONT,
  cex = 1.8,
  tcl = -0.4,
  yaxs = "i",
  xpd = TRUE
)

# 绘图 1
plot(
  Dim1,
  Dim2,
  type = "n",
  xlab = sprintf("PCoA 1 (%.1f%%)", Dim1.pct * 100),
  ylab = sprintf("PCoA 2 (%.1f%%)", Dim2.pct * 100),
  xlim = XLIM,
  ylim = YLIM
)
text(
  Dim1,
  Dim2,
  species,
  cex = 0.6,
  col = species.clrs
)

4.2 绘图 2

标记 samples 名称，还是按照 species 着色，

# 绘图参数
par(
  mar = c(4, 4, 1, 1),
  family = FONT,
  cex = 1.8,
  tcl = -0.4,
  yaxs = "i",
  xpd = TRUE
)

# 绘图 2
samples <- dat[-73, "samples"]
plot(
  Dim1,
  Dim2,
  type = "n",
  xlab = sprintf("PCoA 1 (%.1f%%)", Dim1.pct * 100),
  ylab = sprintf("PCoA 2 (%.1f%%)", Dim2.pct * 100),
  xlim = XLIM,
  ylim = YLIM
)
text(
  Dim1,
  Dim2,
  samples,
  cex = 0.6,
  col = species.clrs
)

4.3 绘图 3

包围每种颜色（去掉 for loop 注释即可），

# 绘图参数
par(
  mar = c(4, 4, 1, 1),
  family = FONT,
  cex = 1.8,
  tcl = -0.4,
  yaxs = "i",
  xpd = TRUE
)

# 绘图 3
plot(
  Dim1,
  Dim2,
  type = "n",
  xlab = sprintf("PCoA 1 (%.1f%%)", Dim1.pct * 100),
  ylab = sprintf("PCoA 2 (%.1f%%)", Dim2.pct * 100),
  xlim = XLIM,
  ylim = YLIM
)
points(
  Dim1,
  Dim2,
  pch = 16,
  cex = 0.6,
  col = species.clrs
)
# for (i.species in unique(PCoA.df$species)) {
#   temp.df <- subset(PCoA.df, species == i.species)
#   polygon(
#     temp.df[chull(temp.df$Dim1, temp.df$Dim2), ],
#     border = temp.df$species.clrs
#   )
# }
add_phylopic_base(
  # uuid = "097336ef-db5f-4edc-a2e4-2c29a10aa512",
  name = "Hydrochoerus hydrochaeris",
  x = 0.3,
  y = 0.1,
  ysize = 0.4,
  fill = unique(species.clrs)[1],
  alpha = 0.6,
  verbose = TRUE
)

Organism silhouettes are from PhyloPic (https://www.phylopic.org/; T. Michael Keesey, 2023) and were added using the rphylopic R package ver. 1.4.0 (Gearty & Jones, 2023). Silhouette was made by Steven Traver, 2012 (CC0 1.0). Silhouette was contributed by Steven Traver.

add_phylopic_base(
  # uuid = "5ee3b11f-9cd8-4aaf-a59e-ac1ddc18d5e8",
  name = "Myocastor coypus",
  x = -0.3,
  y = 0.1,
  ysize = 0.3,
  fill = unique(species.clrs)[2],
  alpha = 0.6,
  verbose = TRUE
)

Organism silhouettes are from PhyloPic (https://www.phylopic.org/; T. Michael Keesey, 2023) and were added using the rphylopic R package ver. 1.4.0 (Gearty & Jones, 2023). Silhouette was made by Mattia Menchetti, 2014 (CC0 1.0). Silhouette was contributed by Mattia Menchetti.

legend(
  "bottomleft",
  legend = unique(species),
  pch = 19,
  col = unique(unique(species.clrs)),
  bty = "n"
)

ecologyR

🚀打赏可获取本公众号代码合集（见置顶文章）📌统计案例📌统计制图📌显著性标记📌结构方程模型可视化工具📌SEM教程与案例📌论文代码复现📌地图可视化

R 代码（修正）｜在 PCoA 图里添加物种图片｜添加水豚剪影｜卡皮巴拉