plot-acute-apap.Rmd

---
title: "Acute APAP (Fig. 2)"
output: 
  workflowr::wflow_html:
    code_folding: hide
editor_options:
  chunk_output_type: console
---

```{r chunk-setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  autodep = TRUE,
  cache = TRUE,
  message = FALSE,
  warning = FALSE
)
```

## Introduction
Here we generate publication ready plots of the analysis of the acute APAP mouse model.

## Libraries and sources
These libraries and sources are used for this analysis.
```{r libs-and-src, message=FALSE, warning=FALSE, cache=FALSE}
library(tidyverse)
library(tidylog)
library(here)

library(AachenColorPalette)
library(cowplot)
library(lemon)
library(ggpubr)
library(patchwork)

source(here("code/utils-plots.R"))
```

Definition of global variables that are used throughout this analysis.
```{r analysis-specific-params, cache=FALSE}
# i/o
data_path <- "data/mouse-acute-apap"
output_path <- "output/mouse-acute-apap"

# graphical parameters
# fontsize
fz <- 7

# keys to annotate contrasts
key_mm <- readRDS(here("data/meta-chronic-vs-acute/contrast_annotation.rds"))
```

## Design
```{r design}
design <- ggdraw() +
  draw_image(here(data_path, "exp-design.pdf"))

design
```

## Histology
```{r histology}
histology <- ggdraw() +
  draw_image(here(data_path, "histology.png"))

histology
```

## Liver enyzmes
```{r liver-enzymes}
df <- read_csv2(here(data_path, "liver_enzymes.csv")) %>%
  mutate(time = fct_inorder(time)) %>%
  pivot_longer(col = -c(time, enzyme), names_to = "mouse", values_to = "y") %>%
  mutate(enzyme = factor(str_to_upper(enzyme), levels = c("ALT", "AST")))

df_summary <- df %>%
  group_by(time, enzyme) %>%
  summarise(mean_se(y)) %>%
  ungroup()

liver_enzymes <- df_summary %>%
  ggplot(aes(x = time, y = y)) +
  geom_errorbar(aes(ymin = ymin, ymax = ymax), width = 0.5) +
  geom_col() +
  facet_rep_wrap(~enzyme, scales = "free", ncol = 1) +
  labs(x = NULL, y = "U/L") +
  my_theme(grid = "y", fsize = fz) +
  stat_compare_means(
    data = df, label = "p.signif",
    ref.group = "Control", hide.ns = T
  ) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))

liver_enzymes
```

## PCA
```{r pca}
pca_result <- readRDS(here(output_path, "pca_result.rds"))

keys <- key_mm %>%
  filter(treatment == "APAP" & class == "Acute") %>%
  distinct(time = value, label = time_label2) %>%
  drop_na() %>%
  add_row(time = 0, label = "Control") %>%
  mutate(time = ordered(time))

pca_plot <- pca_result$coords %>%
  arrange(time) %>%
  inner_join(keys, by = "time") %>%
  mutate(label = fct_inorder(label)) %>%
  ggplot(aes(x = PC1, y = PC2, color = label, label = label)) +
  geom_point() +
  labs(
    x = paste0("PC1", " (", pca_result$var[1], "%)"),
    y = paste0("PC2", " (", pca_result$var[2], "%)"),
    color = "Time"
  ) +
  my_theme(fsize = fz) +
  theme(legend.position = "top",
        legend.box.margin = margin(10, 0, -20, 0)) +
  scale_color_manual(values = aachen_color(c(
    "blue", "purple", "violet",
    "bordeaux", "red", "orange",
    "maygreen", "green", "turquoise",
    "petrol"
  )))

pca_plot
```

## Volcano plot
```{r volcano-plot}
df <- readRDS(here(output_path, "limma_result.rds")) %>%
  filter(contrast_reference == "apap") %>%
  inner_join(key_mm, by = "contrast") %>%
  select(-contrast) %>%
  rename(contrast = time_label2) %>%
  mutate(contrast = fct_drop(contrast)) %>%
  mutate(regulation = fct_recode(regulation,
    Up = "up", Down = "down",
    n.s. = "ns"
  ))

deg_count <- df %>%
  add_count(contrast, regulation) %>%
  filter(regulation != "n.s.") %>%
  mutate(regulation = fct_drop(regulation)) %>%
  mutate(
    logFC = case_when(
      regulation == "Up" ~ 0.75 * max(logFC),
      regulation == "Down" ~ 0.75 * min(logFC)
    ),
    pval = 0.4
  ) %>%
  distinct(n, contrast, logFC, pval, regulation, value) %>%
  complete(contrast, nesting(regulation, logFC, pval), fill = list(n = 0))

main_time <- c(12, 24, 48)
# for main panel
volcano_main <- df %>%
  filter(value %in% main_time) %>%
  plot_volcano(nrow = 1) +
  geom_text(
    data = filter(deg_count, value %in% main_time),
    aes(y = pval, label = n), size = fz / (14 / 5),
    hjust = "inward", vjust = "inward",
    show.legend = F
  ) +
  theme(
    legend.position = "top",
    legend.box.margin = margin(-10, 0, -10, 0)
  ) +
  labs(color = "Regulation") +
  my_theme(grid = "y", fsize = fz)

# for supp panel
volcano_supp <- df %>%
  filter(!value %in% main_time) %>%
  plot_volcano(ncol = 1) +
  geom_text(
    data = filter(deg_count, !value %in% main_time),
    aes(y = pval, label = n), size = fz / (14 / 5),
    show.legend = F
  ) +
  theme(
    legend.position = "top",
    legend.box.margin = margin(-10, 0, -10, 0)
  ) +
  labs(color = "Regulation") +
  my_theme(grid = "y", fsize = fz)

volcano_main
volcano_supp
```

## Top DEGs
```{r top-degs}
df <- readRDS(here(output_path, "limma_result.rds")) %>%
  filter(contrast_reference == "apap") %>%
  inner_join(key_mm) %>%
  select(-class) %>%
  rename(class = time_label2)

top_genes_df <- df %>% 
  filter(regulation != "ns") %>%
  group_by(class, sign(logFC)) %>%
  slice_max(order_by = abs(logFC), n = 10, with_ties = F) %>%
  ungroup() %>% 
  nest(data = -c(class, value))

plots <- top_genes_df %>%
  mutate(p = pmap(., .f = plot_top_genes, fontsize = fz))

main_time <- c(12, 24, 48)

# for main panel
top_genes_main <- plots %>%
  filter(value %in% main_time) %>%
  pull(p) %>%
  wrap_plots() +
  plot_layout(nrow = 1)

# for supp panel
top_genes_supp <- plots %>%
  filter(!value %in% main_time) %>%
  pull(p) %>%
  wrap_plots() +
  plot_layout(ncol = 1)

top_genes_main
top_genes_supp
```

## Time series cluster
```{r ts-cluster}
profile_label <- tribble(
  ~profile, ~process,
  "STEM ID: 14", "Stress response",
  "STEM ID: 9", "Migration",
  "STEM ID: 12", "Proliferation",
  "STEM ID: 6", "Metabolism (1)",
  "STEM ID: 11", "Metabolism (2)",
  "STEM ID: 18", "Metabolism (3)",
  "STEM ID: 8", "Metabolism (4)"
) %>%
  mutate(
    profile = fct_inorder(profile),
    process = coalesce(process, profile),
    process = fct_inorder(process)
  )

stem_res <- readRDS(here(output_path, "stem_result.rds")) %>%
  filter(key == "apap") %>%
  filter(p <= 0.05) %>%
  mutate(profile = fct_reorder(str_c(
    "STEM ID: ",
    as.character(profile)
  ), p)) %>%
  inner_join(profile_label, by = "profile") %>%
  select(-profile) %>%
  rename(profile = process) %>%
  filter(!time %in% c(384, 192, 144))

# extract meta data of profiles
profile_anno <- stem_res %>%
  group_by(key, profile, p) %>%
  mutate(y = 1.2 * abs(max(value))) %>%
  ungroup() %>%
  mutate(max_time = max(time)) %>%
  distinct(key, profile, p, size, y, max_time) %>%
  mutate(label = str_c(size, " ", "genes"))

ts_cluster <- stem_res %>%
  plot_stem_profiles(model_profile = F, nrow = 1) +
  labs(x = "Time in Hours") +
  geom_text(
    data = profile_anno, aes(x = 0, y = y, label = label),
    inherit.aes = F, size = fz / (14 / 5), hjust = "inward"
  ) +
  my_theme(grid = "no", fsize = fz) +
  scale_x_continuous(
    breaks = unique(stem_res$time),
    guide = guide_axis(n.dodge = 3)
  )

ts_cluster
```

## Collage
### Figure 2
```{r fig-2, fig.width=7.87, fig.height=11.69}
fig2 <- (design + histology) /
  ((liver_enzymes | pca_plot) + plot_layout(widths = c(1,3))) /
  ((volcano_main / top_genes_main)) /
  ts_cluster +
  plot_layout(height = c(1.5, 2, 1.75, 0.75)) +
  plot_annotation(tag_levels = list(c("A", "B", "C", "D", "E", "F", "", "", "",
                                      "", "", "G"))) &
  theme(
    plot.tag = element_text(size = fz + 3, face = "bold"),
    legend.key.height = unit(11.5, "pt"),
    legend.key.width = unit(12.5, "pt")
  )

fig2

ggsave(here("figures/Figure 2.pdf"), fig2,
  width = 21, height = 29.7, units = c("cm")
)
ggsave(here("figures/Figure 2.png"), fig2,
  width = 21, height = 29.7, units = c("cm")
)
```

### Supplementary Figure 2.1
```{r s-fig-2-1, fig.width=7.87, fig.height=11.69}
sfig2_1 <- (volcano_supp | top_genes_supp) +
  plot_annotation(tag_levels = list(c("A", "B"))) &
  theme(
    plot.tag = element_text(size = fz + 3, face = "bold"),
    legend.key.height = unit(11.5, "pt"),
    legend.key.width = unit(12.5, "pt")
  )

sfig2_1

ggsave(here("figures/Supplementary Figure 2.1.pdf"), sfig2_1,
  width = 21, height = 29.7, units = c("cm")
)
ggsave(here("figures/Supplementary Figure 2.1.png"), sfig2_1,
  width = 21, height = 29.7, units = c("cm")
)
```