plot-acute-bdl.Rmd

---
title: "Acute BDL (SFig. 2.6)"
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 BDL 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-bdl"
output_path <- "output/mouse-acute-bdl"

# 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), names_to = "enzyme", values_to = "y") %>%
  mutate(enzyme = factor(str_to_upper(enzyme), levels = c("ALT", "AST", "ALP")))

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 == "Bile Duct Ligation" & class == "Acute") %>%
  distinct(time = value, label = time_label2) %>%
  drop_na() %>%
  mutate(time = ordered(time))

pca_plot <- pca_result$coords %>%
  inner_join(keys, by="time") %>%
  mutate(treatment = case_when(treatment == "bdl" ~ "BDL",
                               treatment == "sham" ~ "Sham Surgery")) %>%
  ggplot(aes(x=PC1, y=PC2, color=label, shape = treatment, label = label)) +
  geom_point() +
  labs(x = paste0("PC1", " (", pca_result$var[1], "%)"),
       y = paste0("PC2", " (", pca_result$var[2], "%)"),
       color = "Time", shape = "Group") +
  my_theme(fsize  = fz) +
  theme(legend.position = "top",
        legend.box.margin = margin(10, 0, -20, 10)) +
  scale_color_manual(values = aachen_color(c("blue", "bordeaux", "orange", 
                                             "green"))) +
  scale_shape_manual(values = c(19,15))

pca_plot
```

## Volcano plot
```{r volcano-plot}
df <- readRDS(here(output_path, "limma_result.rds")) %>%
  filter(contrast_reference == "bdl") %>%
  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(1,3,7)

# 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 == "bdl") %>%
  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(1,3,7)

# 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}
stem_res <- readRDS(here(output_path, "stem_result.rds")) %>%
  filter(key == "bdl") %>%
  filter(p <= 0.05) %>%
  mutate(profile = fct_reorder(str_c(
    "STEM ID: ",
    as.character(profile)
  ), p))

# 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 Days") +
  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 = 1)
  )

ts_cluster
```

## Collage
### Supplementary Figure 2.6
```{r s-fig-2-6, fig.width=7.87, fig.height=11.69}
sfig2_6 <- (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")
  )

sfig2_6

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

### Supplementary Figure 2.7
```{r s-fig-2-7}
sfig2_7 <- (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_7

ggsave(here("figures/Supplementary Figure 2.7.pdf"), sfig2_7,
  width = 21, height = 10, units = c("cm")
)
ggsave(here("figures/Supplementary Figure 2.7.png"), sfig2_7,
  width = 21, height = 10, units = c("cm")
)
```