Composing two box plots with significance lines in R's ggplot2 and patchwork (CC310)

October 28, 2024 • PD Schloss • 9 min read

Pat recreates a two paneled figure with box plots that have lines to indicate significance using ggplot2 and patchwork in R. He customizes the appearance using geom_boxplot, stat_summary, factor, annotate, scale_y_continuous, coord_cartesian, labs, and theme functions. You can find the code he developed in this episode at https://www.riffomonas.org/code_club/2024-10-28-boxplots. The original figure was published in the journal mBio (https://journals.asm.org/doi/10.1128/mbio.01966-24); Pat recreated Figure 1. The newsletter describing how I would go about generating the figure can be found at https://shop.riffomonas.org/posts/fun-with-boxplots.

Code

library(tidyverse)
library(ggtext)
library(glue)
library(patchwork)

set.seed(19760620)

coral_physiology <- tibble(
  fragment = 1:72,
  temperature = rep(c(25, 32), each = 36),
  density = c(rnorm(n = 36, mean = 32, sd = 10),
              rnorm(n = 36, mean = 5, sd = 1)),
  photosynthesis = c(rnorm(n = 36, mean = 4.5, sd = 1),
                     rnorm(n = 36, mean = 2, sd = 1)),
  respiration = c(rnorm(n = 36, mean = -2, sd = 0.75),
                  rnorm(n = 36, mean = -1, sd = 0.75))
)


panel_a <- coral_physiology %>%
  ggplot(aes(x = as.character(temperature), y = density)) +
  geom_boxplot(fill = "gray", color = "darkgray",
               staplewidth = 0.3,
               width = 0.5) +
  stat_summary(geom = "crossbar",
               fun = median, fun.min = median, fun.max = median,
               width = 0.5,
               linewidth = 0.15) +
  annotate(geom = "segment",
           x = 1.1, xend = 1.9, y = 53, yend = 53) +
  annotate(geom = "text",
           x = 1.5, y = 54, label = "*", size = 8) +
  labs(x = "Treatments", y = "Symbionts cm<sup>-2</sup>",
       subtitle = "x10<sup>5</sup>") +
  scale_y_continuous(breaks = seq(0, 55, 5), 
                     expand = c(0, 0)) +
  coord_cartesian(clip = "off", ylim = c(0, 55), xlim = c(1, 2)) +
  annotate("segment", x = 0.0, xend = 3, y = -14, yend = -14) +
  theme_classic() +
  theme(axis.title.y = element_markdown(size = 20),
        axis.ticks.y = element_line(color = "gray"),
        axis.text.y = element_text(size = 20),
        axis.ticks.x = element_blank(),
        axis.text.x = element_text(size = 22),
        axis.title.x = element_text(size = 22),
        plot.subtitle = element_markdown(size = 17))


panel_b <- coral_physiology %>%
  select(-density) %>%
  pivot_longer(cols = c(photosynthesis, respiration)) %>%
  mutate(pretty_name = if_else(name == "photosynthesis",
                               glue("{temperature} light"),
                               glue("{temperature} dark")),
         pretty_name = factor(pretty_name, 
                              levels = c("25 light", "32 light",
                                         "25 dark", "32 dark"))) %>%
  ggplot(aes(x = pretty_name, y = value)) +
  geom_boxplot(fill = "gray", color = "darkgray", staplewidth = 0.3) +
  stat_summary(geom = "crossbar",
               fun = median, fun.min = median, fun.max = median,
               width = 0.75,
               linewidth = 0.15) +
  annotate(geom = "segment",
           x = 1.1, xend = 1.9, y = 6.9, yend = 6.9) +
  annotate(geom = "text",
           x = 1.5, y = 7, label = "*", size = 8) +
  annotate(geom = "segment",
           x = 3.1, xend = 3.9, y = 1.5, yend = 1.5) +
  annotate(geom = "text",
           x = 3.5, y = 1.6, label = "*", size = 8) +
  
  labs(x = "Treatments",
       y = "&mu;mol O<sub>2</sub>cm<sup>-2</sup>h<sup>-1</sup>") +
  scale_y_continuous(limits = c(-3, 7.5), expand = c(0, 0),
                     breaks = seq(-4, 7, 1)) +
  coord_cartesian(clip = "off") +
  theme_classic() +
  theme(axis.title.y = element_markdown(size = 17),
        axis.ticks.y = element_line(color = "gray"),
        axis.text.y = element_text(size = 17),
        axis.ticks.x = element_blank(),
        axis.text.x = element_text(size = 22),
        axis.title.x = element_text(size = 22))


panel_a / panel_b + plot_annotation(tag_levels = "A") &
  theme(plot.tag = element_text(size = 24, face = "bold"))
ggsave("coral_physiology.png", width = 7.25, height = 9)