plot-acute-ph.Rmd

---
title: "Acute PH (SFig. 2.4)"
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
)
```

```{r wall-time-start, cache=FALSE, include=FALSE}
# Track time spent on performing this analysis
start_time <- Sys.time()
```

## Introduction
Here we generate publication ready plots of the analysis of the acute PH 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(janitor)

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-ph"
output_path <- "output/mouse-acute-ph"

# 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"))
```

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

### Proliferation
```{r proliferation}
prolif_df <- read_delim(here(data_path, "prolif.txt"), delim = "\t") %>%
  clean_names() %>%
  rename(y = prolif, time = time_days) %>%
  mutate(time = ordered(time))

prolif_summary <- prolif_df %>%
  group_by(time) %>%
  summarise(mean_se(y)) %>%
  ungroup()

prolif <- prolif_df %>%
  ggplot(aes(x = time, y = y)) +
  geom_boxplot() +
  # geom_errorbar(aes(ymin = ymin, ymax = ymax), width = 0.5) +
  # geom_col() +
  labs(x = "Time in days", y = "Proliferation in % (24h period)") +
  my_theme(grid = "y", fsize = fz) +
  stat_compare_means(
    data = prolif_df, label = "p.signif",
    ref.group = "0", hide.ns = T
  ) 

prolif
```

### Lobular area
```{r lobular-area}
lob_df <- read_delim(here(data_path, "lobarea.txt"), delim = "\t") %>%
  clean_names() %>%
  rename(y = lobarea, time = time_days) %>%
  mutate(time = ordered(time))

lob_summary <- lob_df %>%
  group_by(time) %>%
  summarise(mean_se(y)) %>%
  ungroup()

lob <- lob_df %>%
  ggplot(aes(x = time, y = y)) +
  geom_boxplot() +
  # geom_errorbar(aes(ymin = ymin, ymax = ymax), width = 0.5) +
  # geom_col() +
  labs(x = "Time in days", y = expression(paste("Lobule area in ", mm^2))) +
  my_theme(grid = "y", fsize = fz) +
  stat_compare_means(
    data = lob_df, label = "p.signif",
    ref.group = "0", hide.ns = TRUE
  ) 

lob
```

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

keys <- key_mm %>%
  filter(treatment == "Partial Hepatectomy" & class == "Acute") %>%
  distinct(time = value, label = time_label2) %>%
  drop_na() %>%
  add_row(time = 0, label = "Control") %>%
  mutate(
    time = ordered(time),
    label = fct_reorder(label, as.numeric(time))
  )

pca_plot <- pca_result$coords %>%
  mutate(time = ordered(round(as.numeric(as.character(time)) * 24, 0))) %>%
  inner_join(keys, by = "time") %>%
  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(0, 0, 0, 0)
  ) +
  scale_color_manual(values = aachen_color(c(
    "blue", "purple", "violet",
    "bordeaux", "red", "orange",
    "maygreen", "green", "turquoise",
    "petrol", "magenta", "black"
  ))) +
  guides(color = guide_legend(nrow = 2))

pca_plot
```

## Volcano plot
```{r volcano-plot}
df <- readRDS(here(output_path, "limma_result.rds")) %>%
  filter(contrast_reference == "hepatec") %>%
  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 == "hepatec") %>%
  inner_join(key_mm) %>%
  select(-class) %>%
  rename(class = time_label2)

top_genes_df <- df %>%
  group_by(class, sign(logFC)) %>%
  slice_max(order_by = abs(statistic), 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}
stem_res <- readRDS(here(output_path, "stem_result.rds")) %>%
  filter(key == "hepatec") %>%
  filter(p <= 0.05) %>%
  mutate(profile = fct_reorder(str_c(
    "STEM ID: ",
    as.character(profile)
  ), p)) %>%
  filter(!time %in% c(168, 336, 672, 2016))

# 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
### Supplementary Figure 2.4
```{r s-fig-2-2, fig.width=7.87, fig.height=11.69}
sfig2_4 <- (design + plot_spacer()) /
  (((prolif / lob) | 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", "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_4

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

### Supplementary Figure 2.5
```{r s-fig-2-3, fig.width=7.87, fig.height=11.69}
sfig2_5 <- (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_5

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

```{r wall-time-end, cache=FALSE, include=FALSE}
duration <- abs(as.numeric(difftime(Sys.time(), start_time, units = "secs")))
t = print(sprintf("%02d:%02d", duration %% 3600 %/% 60,  duration %% 60 %/% 1))
```
Time spend to execute this analysis: `r t` minutes.