R/plot_sleep.R

In anniliu1/communications: Automate publication-ready plotting for circadian team members

#' Sleep time boxplots and cumulative relative frequency plots
#'
#' This function group boxplot_sleep() and ecdfplot_sleep() can be used to generate publication-ready boxplots on total sleep, daytime sleep, ZT0-4 sleep time of \emph{Drosophila} and cumulative relative frequency plots on daytime sleep bout duration of \emph{Drosophila} under group enrichment or social isolation.
#'
#' @param data a dataframe
#' @param group1 a string name for group enrichment after either 1, 3, 5, or 7 days
#' (e.g. "Grp_1D")
#' @param group2 a string name for social isolation after either 1, 3, 5, or 7 days
#' (e.g. "Iso_1D")
#' @param sleep_time a string name for the sleep time type: "total_sleep", "daytime_sleep",
#' or "ZT04_sleep"
#' @param y_axis_name a string name for the y-axis title: "Total sleep (min)",
#' "Daytime sleep (min)", "ZT0-4 sleep (min)
#' @param y_axis_limit a vector of length 2 that specifies the range of y-axis:
#' c(0, 1300), c(0, 700), or c(0, 350)
#' @param mark_line_position a number that specifies the position of significance line
#'
#' @return An array of boxplots and cumulative relative frequency plots that characterize
#' sleep time of \emph{Drosophila} under group enrichment or social isolation
#'
#'
#' @describeIn Plot_sleep generate 12 boxplots on total sleep, daytime sleep, ZT0-4 sleep time of \emph{Drosophila} after 1, 3, 5, and 7 days of group enrichment or social isolation
#'
#'
#' @examples
#'# Install and load libraries
#'install.packages(c("tidyverse", "ggsignif", "ggpubr"))
#'x <- c("tidyverse", "ggsignif", "ggpubr")
#'lapply(x, require, character.only = TRUE)
#'
#'# Simulate a new dataset that contains total sleep, daytime sleep, and ZT0-4 sleep
#'set.seed(589)
#'sleep_data <- data.frame(
#'   group = c( # Replicate each group name 50 times
#'     rep(x = "Grp_1D", times = 50), rep(x = "Iso_1D", times = 50),
#'     rep(x = "Grp_3D", times = 50), rep(x = "Iso_3D", times = 50),
#'     rep(x = "Grp_5D", times = 50), rep(x = "Iso_5D", times = 50),
#'     rep(x = "Grp_7D", times = 50), rep(x = "Iso_7D", times = 50)
#'),
#'   total_sleep = c( # Generate normal distributions with different means and standard deviations
#'     abs(rnorm(n = 50, mean = 1000, sd = 50)), abs(rnorm(n = 50, mean = 1100, sd = 48)),
#'     abs(rnorm(n = 50, mean = 900, sd = 50)), abs(rnorm(n = 50, mean = 900, sd = 46)),
#'     abs(rnorm(n = 50, mean = 900, sd = 50)), abs(rnorm(n = 50, mean = 820, sd = 39)),
#'     abs(rnorm(n = 50, mean = 970, sd = 40)), abs(rnorm(n = 50, mean = 820, sd = 50))
#'),
#'  daytime_sleep = c(
#'    abs(rnorm(n = 50, mean = 380, sd = 70)), abs(rnorm(n = 50, mean = 400, sd = 65)),
#'    abs(rnorm(n = 50, mean = 380, sd = 70)), abs(rnorm(n = 50, mean = 380, sd = 73)),
#'    abs(rnorm(n = 50, mean = 320, sd = 70)), abs(rnorm(n = 50, mean = 280, sd = 80)),
#'    abs(rnorm(n = 50, mean = 390, sd = 70)), abs(rnorm(n = 50, mean = 230, sd = 79))
#'),
#'  ZT04_sleep = c(
#'    abs(rnorm(n = 50, mean = 60, sd = 68)), abs(rnorm(n = 50, mean = 70, sd = 72)),
#'    abs(rnorm(n = 50, mean = 60, sd = 68)), abs(rnorm(n = 50, mean = 52, sd = 71)),
#'    abs(rnorm(n = 50, mean = 70, sd = 80)), abs(rnorm(n = 50, mean = 50, sd = 82)),
#'    abs(rnorm(n = 50, mean = 70, sd = 78)), abs(rnorm(n = 50, mean = 48, sd = 82))
#'),
#'  sleep_bout = as.vector(sapply(1:8, function(i) { # Generate exponential distributions
#'    rexp(n = 50, r = 0.001)
#'}))
#')
#'
#'# Create data and comparison groups list for boxplots
#'comparison_group_list_boxplot <- list(
#'   list(data = sleep_data, data = sleep_data, data = sleep_data,
#'        data = sleep_data, data = sleep_data, data = sleep_data,
#'        data = sleep_data, data = sleep_data, data = sleep_data,
#'        data = sleep_data, data = sleep_data, data = sleep_data),
#'   list(group1 = "Grp_1D", group1 = "Grp_1D", group1 = "Grp_1D",
#'        group1 = "Grp_3D", group1 = "Grp_3D", group1 = "Grp_3D",
#'        group1 = "Grp_5D", group1 = "Grp_5D", group1 = "Grp_5D",
#'        group1 = "Grp_7D", group1 = "Grp_7D", group1 = "Grp_7D"),
#'   list(group2 = "Iso_1D", group2 = "Iso_1D", group2 = "Iso_1D",
#'        group2 = "Iso_3D", group2 = "Iso_3D", group2 = "Iso_3D",
#'        group2 = "Iso_5D", group2 = "Iso_5D", group2 = "Iso_5D",
#'        group2 = "Iso_7D", group2 = "Iso_7D", group2 = "Iso_7D"),
#'   list(sleep_time = "total_sleep", sleep_time = "daytime_sleep", sleep_time = "ZT04_sleep",
#'        sleep_time = "total_sleep", sleep_time = "daytime_sleep", sleep_time = "ZT04_sleep",
#'        sleep_time = "total_sleep", sleep_time = "daytime_sleep", sleep_time = "ZT04_sleep",
#'        sleep_time = "total_sleep", sleep_time = "daytime_sleep", sleep_time = "ZT04_sleep"),
#'   list(y_axis_name = "Total sleep (min)", y_axis_name = "Daytime sleep (min)", y_axis_name = "ZT0-4 sleep (min)",
#'        y_axis_name = "Total sleep (min)", y_axis_name = "Daytime sleep (min)", y_axis_name = "ZT0-4 sleep (min)",
#'        y_axis_name = "Total sleep (min)", y_axis_name = "Daytime sleep (min)", y_axis_name = "ZT0-4 sleep (min)",
#'        y_axis_name = "Total sleep (min)", y_axis_name = "Daytime sleep (min)", y_axis_name = "ZT0-4 sleep (min)"),
#'   list(y_axis_limit = c(0, 1300), y_axis_limit = c(0, 700), y_axis_limit = c(0, 350),
#'        y_axis_limit = c(0, 1300), y_axis_limit = c(0, 700), y_axis_limit = c(0, 350),
#'        y_axis_limit = c(0, 1300), y_axis_limit = c(0, 700), y_axis_limit = c(0, 350),
#'        y_axis_limit = c(0, 1300), y_axis_limit = c(0, 700), y_axis_limit = c(0, 350)),
#'   list(mark_line_position = 1250, mark_line_position = 670, mark_line_position = 330,
#'        mark_line_position = 1250, mark_line_position = 670, mark_line_position = 330,
#'        mark_line_position = 1250, mark_line_position = 670, mark_line_position = 330,
#'        mark_line_position = 1250, mark_line_position = 670, mark_line_position = 330)
#')
#'
#'# Iterate over data and comparison groups list to plot an array of boxplots
#'boxplot_array <- comparison_group_list_boxplot %>%
#'   pmap(.f = boxplot_sleep)
#' @export
#'
#' @md
boxplot_sleep <- function(data,
                       group1,
                       group2,
                       sleep_time,
                       y_axis_name,
                       y_axis_limit,
                       mark_line_position) {

  # Prepare data for specific sleep time
  data <- data %>%
    select(group, all_of(sleep_time)) %>% # Select group and specified sleep time columns
    rename(value = all_of(sleep_time)) # Rename specified sleep time into value

  # Prepare summary statistics (median, q1, and q3) for Boxplot plotting
  summary_data <- data %>%
    filter(group %in% c(group1, group2)) %>%
    dplyr::group_by(group) %>%
    summarise(
      median = median(value),
      q1 = quantile(value, probs = 0.25),
      q3 = quantile(value, probs = 0.75)
    )
  median_group1 <- summary_data %>%
    filter(group == group1) %>%
    pull(median)
  q1_group1 <- summary_data %>%
    filter(group == group1) %>%
    pull(q1)
  q3_group1 <- summary_data %>%
    filter(group == group1) %>%
    pull(q3)
  median_group2 <- summary_data %>%
    filter(group == group2) %>%
    pull(median)
  q1_group2 <- summary_data %>%
    filter(group == group2) %>%
    pull(q1)
  q3_group2 <- summary_data %>%
    filter(group == group2) %>%
    pull(q3)

  # Plot the boxplot
  plot <- ggplot(
    data = data %>%
      filter(group %in% c(group1, group2)),
    mapping = aes(x = group, y = value)
  ) +
    geom_boxplot(
      color = "white"
    ) +
    geom_jitter(
      aes(colour = group),
      size = 2,
      alpha = 0.2,
      width = 0.05
    ) +
    scale_color_manual(values = c("royalblue1", "indianred1")) +
    geom_segment(
      aes(
        x = 0.95, xend = 1.05,
        y = median_group1,
        yend = median_group1
      ),
      colour = "black"
    ) +
    geom_segment(
      aes(
        x = 0.95, xend = 1.05,
        y = q1_group1,
        yend = q1_group1
      ),
      colour = "black"
    ) +
    geom_segment(
      aes(
        x = 0.95, xend = 1.05,
        y = q3_group1,
        yend = q3_group1
      ),
      colour = "black"
    ) +
    geom_segment(
      aes(
        x = 1,
        xend = 1,
        y = q1_group1,
        yend = q3_group1
      ),
      colour = "black"
    ) +
    geom_segment(
      aes(
        x = 1.95, xend = 2.05,
        y = median_group2,
        yend = median_group2
      ),
      colour = "black"
    ) +
    geom_segment(
      aes(
        x = 1.95, xend = 2.05,
        y = q1_group2,
        yend = q1_group2
      ),
      colour = "black"
    ) +
    geom_segment(
      aes(
        x = 1.95, xend = 2.05,
        y = q3_group2,
        yend = q3_group2
      ),
      colour = "black"
    ) +
    geom_segment(
      aes(
        x = 2,
        xend = 2,
        y = q1_group2,
        yend = q3_group2
      ),
      colour = "black"
    ) +
    theme_classic() +
    theme(
      axis.text.x = element_text(
        size = 12,
        angle = 45,
        vjust = 1,
        hjust = 1
      ),
      legend.position = "none"
    ) +
    labs(
      title = "",
      x = "",
      y = y_axis_name
    ) +
    scale_y_continuous(
      n.breaks = 7,
      limits = y_axis_limit
    )

  # Add significance mark onto the boxplot
  boxplot <- plot +
    geom_signif(
      comparisons = list(c(group1, group2)),
      test = t.test,
      y_position = mark_line_position,
      tip_length = 0,
      vjust = 0.1,
      map_signif_level = c(
        "****" = 0.0001,
        "***" = 0.001,
        "**" = 0.01,
        "*" = 0.05
      )
    )
}

#' @describeIn Plot_sleep generate 4 cumulative relative frequency plots on daytime sleep bout duration of \emph{Drosophila} under group enrichment or social isolation
#'
#' @examples
#'
#'# Create data and comparison groups list for cumulative relative frequency plots
#'comparison_group_list_ecdfplot <- list( # Create data and comparison groups list
#'   list(data = sleep_data, data = sleep_data, data = sleep_data, data = sleep_data),
#'   list(group1 = "Grp_1D", group1 = "Grp_3D", group1 = "Grp_5D", group1 = "Grp_7D"),
#'   list(group2 = "Iso_1D", group2 = "Iso_3D", group2 = "Iso_5D", group2 = "Iso_7D")
#')
#'
#'# Iterate over data and comparison groups list to plot an array of cumulative relative frequency plots
#'ecdfplot_array <- comparison_group_list_ecdfplot %>%
#'   pmap(.f = ecdfplot_sleep)
#'
#'# Arrange boxplots and cumulative relative frequency plots together
#'plot_array <- ggarrange(
#'   ggarrange(boxplot_array[[1]], boxplot_array[[2]], boxplot_array[[3]], ecdfplot_array[[1]],
#'             ncol = 4, labels = "A", widths = c(1, 1, 1, 2)
#'),
#'   ggarrange(boxplot_array[[4]], boxplot_array[[5]], boxplot_array[[6]], ecdfplot_array[[2]],
#'             ncol = 4, labels = "B", widths = c(1, 1, 1, 2)
#'),
#'   ggarrange(boxplot_array[[7]], boxplot_array[[8]], boxplot_array[[9]], ecdfplot_array[[3]],
#'             ncol = 4, labels = "C", widths = c(1, 1, 1, 2)
#'),
#'   ggarrange(boxplot_array[[10]], boxplot_array[[11]], boxplot_array[[12]], ecdfplot_array[[4]],
#'             ncol = 4, labels = "D", widths = c(1, 1, 1, 2)
#'),
#'   nrow = 4
#') %>%
#'   annotate_figure(bottom = text_grob("Figure 3. Sleep time in Drosophila.",
#'                                      face = "italic"))
#'
#'# Save the plot
#'ggsave("plot_array.png", dpi = 600, height = 15, width = 8)

#' @export
#'
#' @md
ecdfplot_sleep <- function(data, group1, group2) {

  # Prepare data for plotting
  ecdf_data <- data %>%
    filter(group %in% c(group1, group2))

  ggplot(
    data = ecdf_data,
    mapping = aes(x = sleep_bout, color = group)
  ) +
    geom_point(stat = "ecdf", shape = 1) +
    labs(
      title = "",
      x = "Daytime sleep bout \nduration (min)",
      y = "Cumulative relative \nfrequency",
      color = ""
    ) +
    theme_classic() +
    theme(
      axis.text.x = element_text(
        size = 12,
        vjust = 1,
        hjust = 1
      ),
      legend.position = c(0.2, 0.4), # Adjust the legend position
      legend.background = element_blank()
    ) +
    scale_color_manual(values = c("royalblue1", "indianred1")) +
    scale_x_continuous(
      trans = "log10", # Log transform the x-axis
      limits = c(1, 1000)
    ) # Unify the location and number of ticks
}