plot-chronic-vs-acute.Rmd

---
title: "Chronic vs. acute (Fig. 3)"
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 for the comparison of chronic and acute mouse models.

## 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(VennDiagram)
library(scales)
library(lemon)
library(ComplexHeatmap)
library(ggwordcloud)
library(circlize)
library(patchwork)

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

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

# graphical parameters
# fontsize
fz <- 7
# color function for heatmaps
col_fun <- colorRamp2(
  c(-4, 0, 4),
  c(aachen_color("blue"), "white", aachen_color("red"))
)

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

## Time point of maximal liver damage
```{r time-point-max-liver-damage}
keys <- key_mm %>%
  filter(class == "Acute" & treatment_abbr %in% c("CCl4", "APAP", "PH")) %>%
  distinct(time = value, label = time_label2, treatment_abbr) %>%
  drop_na() %>%
  add_row(time = 0, label = "Control", treatment_abbr = "CCl4") %>%
  add_row(time = 0, label = "Control", treatment_abbr = "APAP") %>%
  add_row(time = 0, label = "Control", treatment_abbr = "PH") %>%
  mutate(time = ordered(time))

pca_dist <- readRDS(here(output_path, "pca_dist.rds")) %>%
  inner_join(keys, by = c("time", "treatment_abbr")) %>%
  arrange(time) %>%
  mutate(label = fct_inorder(label))

max_liver_damage <- pca_dist %>%
  ggplot(aes(x = label, y = dist, fill = max)) +
  geom_col() +
  my_theme(grid = "y", fsize = fz) +
  scale_fill_manual(values = aachen_color(c("black50", "green"))) +
  labs(
    x = NULL,
    y = "Absolute mean distance to control along PC1"
  ) +
  theme(
    legend.position = "none",
    axis.text.x = element_text(angle = 45, hjust = 1)
  ) +
  facet_rep_wrap(~treatment_abbr, scales = "free")

max_liver_damage
```

## Similarity of acute gene sets
```{r sim-acute-gene-sets}
keys <- key_mm %>%
  distinct(contrast, label = label2)

j <- readRDS(here(output_path, "gene_set_similarity.rds")) %>%
  separate(set1, into = c("tmp1", "tmp2", "tmp3", "contrast"), sep = "-") %>%
  select(-starts_with("tmp")) %>%
  inner_join(keys, by = "contrast") %>%
  rename(set1 = label) %>%
  select(-contrast) %>%
  separate(set2, into = c("tmp1", "tmp2", "tmp3", "contrast"), sep = "-") %>%
  select(-starts_with("tmp")) %>%
  inner_join(keys, by = "contrast") %>%
  rename(set2 = label) %>%
  select(-contrast)

acute_gs_sim <- j %>%
  mutate(set1 = fct_rev(set1)) %>%
  ggplot(aes(
    x = set1, y = set2, fill = similarity,
    label = round(similarity, 3)
  )) +
  geom_tile(color = "black", size = 0.2) +
  scale_fill_gradient(low = "white", high = aachen_color("green")) +
  labs(x = NULL, y = NULL, fill = "Jaccard\nIndex") +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    axis.line = element_blank(),
    axis.ticks = element_blank()
  ) +
  geom_text(size = (fz - 2) / (14 / 5)) +
  my_theme(grid = "no", fsize = fz)

acute_gs_sim
```

## Interstudy enrichment of acute studies
```{r interstudy-enrichment-acute}
keys <- key_mm %>%
  distinct(contrast, label = label2)

gsea_res <- readRDS(here(output_path, "interstudy_enrichment.rds")) %>%
  separate(signature,
    into = c("tmp1", "tmp2", "tmp3", "contrast"),
    sep = "-"
  ) %>%
  select(-starts_with("tmp")) %>%
  inner_join(keys, by = "contrast") %>%
  rename(signature = label) %>%
  select(-contrast) %>%
  separate(geneset, into = c("tmp1", "tmp2", "tmp3", "contrast"), sep = "-") %>%
  select(-starts_with("tmp")) %>%
  inner_join(keys, by = "contrast") %>%
  rename(geneset = label) %>%
  select(-contrast)

acute_interstudy_enrichment <- gsea_res %>%
  mutate(direction = fct_rev(str_to_title(direction))) %>%
  mutate(label = gtools::stars.pval(padj)) %>%
  ggplot(aes(x = signature, y = geneset, fill = ES)) +
  geom_tile() +
  geom_text(aes(label = label), size = fz / (14 / 5), vjust = 1) +
  facet_wrap(~direction, ncol = 1) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    axis.line = element_blank(),
    axis.ticks = element_blank()
  ) +
  scale_fill_gradient2(
    low = aachen_color("blue"), mid = "white",
    high = aachen_color("red")
  ) +
  my_theme(grid = "no", fsize = fz) +
  labs(x = "Signature", y = "Gene Set", fill = "ES") +
  guides(fill = guide_colorbar(title = "ES"))

acute_interstudy_enrichment
```

## Union of acute and chronic genes
```{r union-chronic-acute-gene-pools}
keys <- key_mm %>%
  distinct(contrast, label = label2)

c <- readRDS(here(output_path, "chronic_gene_pool.rds"))
a <- readRDS(here(output_path, "acute_gene_pool.rds"))

df <- bind_rows(c, a) %>%
  inner_join(keys, by = "contrast") %>%
  mutate(
    statistic = case_when(
      statistic >= 25 ~ 25,
      TRUE ~ statistic
    ),

    class = str_to_title(class)
  )

union_a <- df %>%
  filter(class == "Acute") %>%
  ggplot(aes(x = label, y = fct_reorder(gene, statistic, mean))) +
  geom_tile(aes(fill = statistic)) +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    legend.position = "none"
  ) +
  labs(y = "Pool of acute genes", x = NULL, fill = "t-statistic") +
  facet_rep_wrap(~class) +
  scale_fill_gradient2(
    low = aachen_color("blue"), mid = "white",
    high = aachen_color("red")
  ) +
  my_theme(grid = "no", fsize = fz)


union_c <- df %>%
  filter(class == "Chronic") %>%
  ggplot(aes(
    x = label, y = fct_reorder(gene, statistic, mean),
    fill = statistic
  )) +
  geom_tile() +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) +
  labs(y = "Pool of chronic genes", x = NULL, fill = "t-statistic") +
  facet_rep_wrap(~class) +
  scale_fill_gradient2(
    low = aachen_color("blue"), mid = "white",
    high = aachen_color("red")
  ) +
  my_theme(grid = "no", fsize = fz)

union_a +
  union_c +
  plot_layout(widths = c(8, 1))
```

## Overlap of unified gene sets
```{r overlap-unified-genes}
acute_gene_union <- readRDS(here(output_path, "union_acute_geneset.rds")) %>%
  mutate(class = "acute")
chronic_gene_union <- readRDS(here(
  output_path,
  "union_chronic_geneset.rds"
)) %>%
  mutate(class = "chronic")

a1 <- acute_gene_union %>% nrow()
a2 <- chronic_gene_union %>% nrow()
ca <- intersect(
  acute_gene_union %>% pull(gene),
  chronic_gene_union %>% pull(gene)
) %>%
  length()

grid.newpage()
v <- grid.grabExpr(draw.pairwise.venn(
  area1 = a1, area2 = a2, cross.area = ca,
  category = c("Acute", "Chronic"),
  # lty = "blank",
  cex = 1 / 12 * fz,
  fontfamily = rep("sans", 3),
  # fill = aachen_color(c("purple", "petrol")),
  # cat.col = aachen_color(c("purple", "petrol")),
  cat.cex = 1 / 12 * (fz + 1),
  cat.fontfamily = rep("sans", 2),
  cat.pos = c(340, 20),
  cat.just = list(c(0.5, 0), c(0.5, 0))
))
```

## Top exclusive chronic genes
### Heatmap
```{r exclusive-chronic-hmap}
keys <- key_mm %>%
  distinct(contrast, label2)

df <- readRDS(here(output_path, "ranked_exclusive_chronic_genes.rds"))
contrasts <- readRDS(here(output_path, "limma_result.rds"))

acute_contrasts <- c(
  "treat_vs_ctrl",
  "inLiver_lps_vs_ctrl",
  "ccl_8h_vs_0h", "ccl_24h_vs_0h", "ccl_48h_vs_0h",
  "apap_12h_vs_0h", "apap_24h_vs_0h", "apap_48h_vs_0h",
  "ph_0.5d", "ph_1d", "ph_2d",
  "bdl_vs_sham_1d"
)

mat_exclusive_chronic_genes <- df %>%
  filter(rank <= 100) %>%
  left_join(contrasts) %>%
  filter(contrast %in% acute_contrasts | treatment == "pure_ccl4") %>%
  inner_join(keys, by = "contrast") %>%
  mutate(gene = as_factor(gene)) %>%
  select(gene, label2, logFC) %>%
  untdy("gene", "label2", "logFC") %>%
  as.matrix()

exclusive_chronic_hmap <- Heatmap(
  t(mat_exclusive_chronic_genes),
  col = col_fun,
  cluster_rows = T, cluster_columns = T,
  show_row_dend = F,
  row_names_gp = gpar(fontsize = fz), column_names_gp = gpar(fontsize = fz - 2),
  name = "logFC",
  heatmap_legend_param = list(
    title_gp = gpar(
      fontface = "plain",
      fontsize = fz + 1
    ),
    labels_gp = gpar(fontsize = fz)
  ),
  row_gap = unit(2.5, "mm"),
  border = T,
  row_split = c(rep("Acute", 12), rep("Chronic", 3)),
  row_title_gp = gpar(fontsize = fz + 1)
)

exclusive_chronic_hmap
```

### Characteriatzion via TFs
```{r exclusive-chronic-tfs}
ora_res <- readRDS(here(
  output_path,
  "exclusive_genes_characterization.rds"
)) %>%
  filter(group %in% c("dorothea") & fdr <= 0.2) %>%
  filter(class == "chronic") %>%
  mutate(regulation = str_to_title(regulation))

chronic_tfs_up <- ora_res %>%
  filter(regulation == "Up") %>%
  filter(geneset %in% c("Hif1a", "Klf5")) %>%
  ggplot(aes(x = fct_reorder(geneset, p.value), y = -log10(p.value))) +
  geom_col() +
  geom_text(aes(label = geneset),
    color = "white", angle = 90, hjust = 1.5,
    size = fz / (14 / 5)
  ) +
  theme(axis.text.x = element_blank()) +
  labs(x = "TF", y = expression(-log["10"] * "(p-value)")) +
  my_theme(grid = "y", fsize = fz) +
  facet_rep_wrap(~regulation)

chronic_tfs_down <- ora_res %>%
  filter(regulation == "Down") %>%
  slice_min(order_by = p.value, n = 3) %>%
  ggplot(aes(x = fct_reorder(geneset, p.value), y = -log10(p.value))) +
  geom_col() +
  geom_text(aes(label = geneset),
    color = "white", angle = 90, hjust = 1.5,
    size = fz / (14 / 5)
  ) +
  theme(axis.text.x = element_blank()) +
  labs(x = "TF", y = expression(-log["10"] * "(p-value)")) +
  my_theme(grid = "y", fsize = fz) +
  facet_rep_wrap(~regulation)

chronic_tfs_up + chronic_tfs_down
```

### Characteriatzion via pathways
```{r exclusive-chronic-pws}
ora_res <- readRDS(here(
  output_path,
  "exclusive_genes_characterization.rds"
)) %>%
  filter(group %in% c("progeny") & fdr <= 0.2) %>%
  filter(class == "chronic") %>%
  mutate(regulation = str_to_title(regulation))

chronic_pw_up <- ora_res %>%
  filter(regulation == "Up") %>%
  filter(geneset %in% c("TGFb")) %>%
  ggplot(aes(x = fct_reorder(geneset, p.value), y = -log10(p.value))) +
  geom_col() +
  geom_text(aes(label = geneset),
    color = "white", angle = 90, hjust = 1.5,
    size = fz / (14 / 5)
  ) +
  theme(axis.text.x = element_blank()) +
  labs(x = "Pathway", y = expression(-log["10"] * "(p-value)")) +
  my_theme(grid = "y", fsize = fz) +
  facet_rep_wrap(~regulation)

chronic_pw_up
```

### Characterization via GO terms
```{r exclusive-chronic-go}
# get wordcounts for go terms
wordcounts <- readRDS(here(output_path, "go_wordcounts.rds")) %>%
  filter(class == "chronic")

# get ranking of go clusters
cluster_ranking <- readRDS(here(output_path, "go_cluster_ranking.rds")) %>%
  filter(class == "chronic") %>%
  group_by(cluster, class) %>%
  mutate(label = str_c(n(), " ", "GO terms")) %>%
  ungroup() %>%
  mutate(regulation = fct_rev(str_to_title(regulation))) %>%
  arrange(desc(cluster)) %>%
  mutate(description = fct_rev(as_factor(description)))

cluster_anno <- cluster_ranking %>%
  nest(data = c(rank, term)) %>%
  # find max peak for each cluster
  mutate(peak = data %>% map(function(data) {
    max(density(data$rank)$y)
  })) %>%
  unnest(c(peak)) %>%
  group_by(regulation) %>%
  mutate(max_peak = max(peak)) %>%
  ungroup() %>%
  distinct(
    cluster, regulation, description, peak, max_peak, max_rank,
    label
  ) %>%
  arrange(description) %>%
  group_by(regulation) %>%
  mutate(n_clusters = row_number()) %>%
  mutate(x_coord = case_when(
    n_clusters == 1 ~ 0,
    n_clusters == 2 ~ 1 * max_rank
  )) %>%
  ungroup()

# up-regulated genes
up_dens <- cluster_ranking %>%
  filter(regulation == "Up") %>%
  plot_go_rank_density() +
  my_theme(grid = "no", fsize = fz) +
  geom_text(
    data = filter(cluster_anno, regulation == "Up"),
    aes(x = x_coord, y = 1.1 * max_peak, label = label, color = description),
    inherit.aes = F, size = fz / (14 / 5), hjust = "inward",
    show.legend = F
  ) +
  scale_color_manual(values = aachen_color(c("petrol", "orange"))) +
  guides(color = guide_legend(nrow = 2)) +
  theme(legend.box.margin = margin(-10, 0, -20, 0))

up_cloud <- wordcounts %>%
  filter(regulation == "up") %>%
  slice_max(order_by = n, n = 15) %>%
  plot_wordcloud(fontsize = fz)

exclusive_chronic_up <- up_dens +
  inset_element(up_cloud,
    left = 0, bottom = 0, right = 1, top = 1,
    align_to = "panel"
  )

exclusive_chronic_up
```

## Top common genes
### Heatmap
```{r common-hmap}
keys <- key_mm %>%
  distinct(contrast, label2)

df <- readRDS(here(output_path, "ranked_common_genes.rds"))
contrasts <- readRDS(here(output_path, "limma_result.rds"))

acute_contrasts <- c(
  "treat_vs_ctrl",
  "inLiver_lps_vs_ctrl",
  "ccl_8h_vs_0h", "ccl_24h_vs_0h", "ccl_48h_vs_0h",
  "apap_12h_vs_0h", "apap_24h_vs_0h", "apap_48h_vs_0h",
  "ph_0.5d", "ph_1d", "ph_2d",
  "bdl_vs_sham_1d"
)

mat_common_genes <- df %>%
  filter(rank <= 100) %>%
  left_join(contrasts) %>%
  filter(contrast %in% acute_contrasts | treatment == "pure_ccl4") %>%
  inner_join(keys, by = "contrast") %>%
  mutate(gene = as_factor(gene)) %>%
  select(gene, label2, logFC) %>%
  untdy("gene", "label2", "logFC") %>%
  as.matrix()

common_hmap <- Heatmap(
  t(mat_common_genes),
  col = col_fun,
  cluster_rows = T, cluster_columns = T,
  show_row_dend = F,
  row_names_gp = gpar(fontsize = fz), column_names_gp = gpar(fontsize = fz - 2),
  name = "logFC",
  heatmap_legend_param = list(
    title_gp = gpar(
      fontface = "plain",
      fontsize = fz + 1
    ),
    labels_gp = gpar(fontsize = fz)
  ),
  row_gap = unit(2.5, "mm"),
  border = T,
  row_split = c(rep("Acute", 12), rep("Chronic", 3)),
  row_title_gp = gpar(fontsize = fz + 1)
)

common_hmap
```

### Characteriatzion via TFs
```{r common-tfs}
ora_res <- readRDS(here(
  output_path,
  "exclusive_genes_characterization.rds"
)) %>%
  filter(group %in% c("dorothea") & fdr <= 0.2) %>%
  filter(class == "common") %>%
  mutate(regulation = str_to_title(regulation))

common_tfs_up <- ora_res %>%
  filter(regulation == "Up") %>%
  filter(geneset %in% c("Klf5")) %>%
  ggplot(aes(x = fct_reorder(geneset, p.value), y = -log10(p.value))) +
  geom_col() +
  geom_text(aes(label = geneset),
    color = "white", angle = 90, hjust = 1.5,
    size = fz / (14 / 5)
  ) +
  theme(axis.text.x = element_blank()) +
  labs(x = "TF", y = expression(-log["10"] * "(p-value)")) +
  my_theme(grid = "y", fsize = fz) +
  facet_rep_wrap(~regulation)

common_tfs_down <- ora_res %>%
  filter(regulation == "Down") %>%
  slice_min(order_by = p.value, n = 3) %>%
  ggplot(aes(x = fct_reorder(geneset, p.value), y = -log10(p.value))) +
  geom_col() +
  geom_text(aes(label = geneset),
    color = "white", angle = 90, hjust = 1.5,
    size = fz / (14 / 5)
  ) +
  theme(axis.text.x = element_blank()) +
  labs(x = "TF", y = expression(-log["10"] * "(p-value)")) +
  my_theme(grid = "y", fsize = fz) +
  facet_rep_wrap(~regulation)

common_tfs_up + common_tfs_down
```

### Characteriatzion via pathways
```{r common-pws}
ora_res <- readRDS(here(
  output_path,
  "exclusive_genes_characterization.rds"
)) %>%
  filter(group %in% c("progeny") & fdr <= 0.2) %>%
  filter(class == "common") %>%
  mutate(regulation = str_to_title(regulation))

common_pw_up <- ora_res %>%
  filter(regulation == "Up") %>%
  filter(geneset %in% c("NFkB", "TNFa")) %>%
  ggplot(aes(x = fct_reorder(geneset, p.value), y = -log10(p.value))) +
  geom_col() +
  geom_text(aes(label = geneset),
    color = "white", angle = 90, hjust = 1.5,
    size = fz / (14 / 5)
  ) +
  theme(axis.text.x = element_blank()) +
  labs(x = "Pathway", y = expression(-log["10"] * "(p-value)")) +
  my_theme(grid = "y", fsize = fz) +
  facet_rep_wrap(~regulation)

# no common_pw_up

common_pw_up
```

### Characterization via GO terms
```{r common-go}
# get wordcounts for go terms
wordcounts <- readRDS(here(output_path, "go_wordcounts.rds")) %>%
  filter(class == "common")

# get ranking of go clusters
cluster_ranking <- readRDS(here(output_path, "go_cluster_ranking.rds")) %>%
  filter(class == "common") %>%
  group_by(cluster, class) %>%
  mutate(label = str_c(n(), " ", "GO terms")) %>%
  ungroup() %>%
  mutate(regulation = fct_rev(str_to_title(regulation))) %>%
  arrange(desc(cluster)) %>%
  mutate(description = fct_rev(as_factor(description)))

cluster_anno <- cluster_ranking %>%
  nest(data = c(rank, term)) %>%
  # find max peak for each cluster
  mutate(peak = data %>% map(function(data) {
    max(density(data$rank)$y)
  })) %>%
  unnest(c(peak)) %>%
  group_by(regulation) %>%
  mutate(max_peak = max(peak)) %>%
  ungroup() %>%
  distinct(
    cluster, regulation, description, peak, max_peak, max_rank,
    label
  ) %>%
  arrange(description) %>%
  group_by(regulation) %>%
  mutate(n_clusters = row_number()) %>%
  mutate(x_coord = case_when(
    n_clusters == 1 ~ 0,
    n_clusters == 2 ~ 1 * max_rank
  )) %>%
  ungroup()

# up-regulated genes
up_dens <- cluster_ranking %>%
  filter(regulation == "Up") %>%
  plot_go_rank_density() +
  my_theme(grid = "no", fsize = fz) +
  geom_text(
    data = filter(cluster_anno, regulation == "Up"),
    aes(x = x_coord, y = 1.1 * max_peak, label = label, color = description),
    inherit.aes = F, size = fz / (14 / 5), hjust = "inward",
    show.legend = F
  ) +
  scale_color_manual(values = aachen_color("turquoise")) +
  guides(color = guide_legend(nrow = 2)) +
  theme(legend.box.margin = margin(-10, 0, -20, 0))

up_cloud <- wordcounts %>%
  filter(regulation == "up") %>%
  slice_max(order_by = n, n = 15) %>%
  plot_wordcloud(fontsize = fz)

common_up <- up_dens +
  inset_element(up_cloud,
    left = 0, bottom = 0, right = 1, top = 1,
    align_to = "panel"
  )

# down-regulated genes
down_dens <- cluster_ranking %>%
  filter(regulation == "Down") %>%
  plot_go_rank_density() +
  my_theme(grid = "no", fsize = fz) +
  geom_text(
    data = filter(cluster_anno, regulation == "Down"),
    aes(x = x_coord, y = 1.1 * max_peak, label = label, color = description),
    inherit.aes = F, size = fz / (14 / 5), hjust = "inward",
    show.legend = F
  ) +
  scale_color_manual(values = aachen_color("purple")) +
  guides(color = guide_legend(nrow = 2)) +
  theme(legend.box.margin = margin(-10, 0, -20, 0))

down_cloud <- wordcounts %>%
  filter(regulation == "down") %>%
  slice_max(order_by = n, n = 15) %>%
  plot_wordcloud(fontsize = fz)

common_down <- down_dens +
  inset_element(down_cloud,
    left = 0, bottom = 0, right = 1, top = 1,
    align_to = "panel"
  )

common_up + common_down
```

## Top exclusive acute genes
### Heatmap
```{r exclusive-acute-hmap}
keys <- key_mm %>%
  distinct(contrast, label2)

df <- readRDS(here(output_path, "ranked_exclusive_acute_genes.rds"))
contrasts <- readRDS(here(output_path, "limma_result.rds"))

acute_contrasts <- c(
  "treat_vs_ctrl",
  "inLiver_lps_vs_ctrl",
  "ccl_8h_vs_0h", "ccl_24h_vs_0h", "ccl_48h_vs_0h",
  "apap_12h_vs_0h", "apap_24h_vs_0h", "apap_48h_vs_0h",
  "ph_0.5d", "ph_1d", "ph_2d",
  "bdl_vs_sham_1d"
)

mat_exclusive_acute_genes <- df %>%
  filter(rank <= 100) %>%
  left_join(contrasts) %>%
  filter(contrast %in% acute_contrasts | treatment == "pure_ccl4") %>%
  inner_join(keys, by = "contrast") %>%
  mutate(gene = as_factor(gene)) %>%
  select(gene, label2, logFC) %>%
  untdy("gene", "label2", "logFC") %>%
  as.matrix()

exclusive_acute_hmap <- Heatmap(
  t(mat_exclusive_acute_genes),
  col = col_fun,
  cluster_rows = T, cluster_columns = T,
  show_row_dend = F,
  row_names_gp = gpar(fontsize = fz), column_names_gp = gpar(fontsize = fz - 2),
  name = "logFC",
  heatmap_legend_param = list(
    title_gp = gpar(
      fontface = "plain",
      fontsize = fz + 1
    ),
    labels_gp = gpar(fontsize = fz)
  ),
  row_gap = unit(2.5, "mm"),
  border = T,
  row_split = c(rep("Acute", 12), rep("Chronic", 3)),
  row_title_gp = gpar(fontsize = fz + 1)
)

exclusive_acute_hmap
```

### Characteriatzion via TFs
```{r exclusive-acute-tfs}
ora_res <- readRDS(here(
  output_path,
  "exclusive_genes_characterization.rds"
)) %>%
  filter(group %in% c("dorothea") & fdr <= 0.2) %>%
  filter(class == "acute") %>%
  mutate(regulation = str_to_title(regulation))

acute_tfs_up <- ora_res %>%
  filter(regulation == "Up") %>%
  filter(geneset %in% c("Myc", "Trp53", "Stat3")) %>%
  ggplot(aes(x = fct_reorder(geneset, p.value), y = -log10(p.value))) +
  geom_col() +
  geom_text(aes(label = geneset),
    color = "white", angle = 90, hjust = 1.5,
    size = fz / (14 / 5)
  ) +
  theme(axis.text.x = element_blank()) +
  labs(x = "TF", y = expression(-log["10"] * "(p-value)")) +
  my_theme(grid = "y", fsize = fz) +
  facet_rep_wrap(~regulation)

acute_tfs_down <- ora_res %>%
  filter(regulation == "Down") %>%
  slice_min(order_by = p.value, n = 3) %>%
  ggplot(aes(x = fct_reorder(geneset, p.value), y = -log10(p.value))) +
  geom_col() +
  geom_text(aes(label = geneset),
    color = "white", angle = 90, hjust = 1.5,
    size = fz / (14 / 5)
  ) +
  theme(axis.text.x = element_blank()) +
  labs(x = "TF", y = expression(-log["10"] * "(p-value)")) +
  my_theme(grid = "y", fsize = fz) +
  facet_rep_wrap(~regulation)

acute_tfs_up + acute_tfs_down
```

### Characteriatzion via pathways
```{r exclusive-acute-pws}
ora_res <- readRDS(here(
  output_path,
  "exclusive_genes_characterization.rds"
)) %>%
  filter(group %in% c("progeny") & fdr <= 0.2) %>%
  filter(class == "acute") %>%
  mutate(regulation = str_to_title(regulation))

acute_pw_up <- ora_res %>%
  filter(regulation == "Up") %>%
  filter(geneset %in% c("MAPK", "EGFR", "TNFa")) %>%
  ggplot(aes(x = fct_reorder(geneset, p.value), y = -log10(p.value))) +
  geom_col() +
  geom_text(aes(label = geneset),
    color = "white", angle = 90, hjust = 1.5,
    size = fz / (14 / 5)
  ) +
  theme(axis.text.x = element_blank()) +
  labs(x = "Pathway", y = expression(-log["10"] * "(p-value)")) +
  my_theme(grid = "y", fsize = fz) +
  facet_rep_wrap(~regulation)

acute_pw_down <- ora_res %>%
  filter(regulation == "Down") %>%
  slice_min(order_by = p.value, n = 3) %>%
  ggplot(aes(x = fct_reorder(geneset, p.value), y = -log10(p.value))) +
  geom_col() +
  geom_text(aes(label = geneset),
    color = "white", angle = 90, hjust = 1.5,
    size = fz / (14 / 5)
  ) +
  theme(axis.text.x = element_blank()) +
  labs(x = "Pathway", y = expression(-log["10"] * "(p-value)")) +
  my_theme(grid = "y", fsize = fz) +
  facet_rep_wrap(~regulation)

acute_pw_up + acute_pw_down
```

### Characterization via GO terms
```{r exclusive-acute-go}
# get wordcounts for go terms
wordcounts <- readRDS(here(output_path, "go_wordcounts.rds")) %>%
  filter(class == "acute")

# get ranking of go clusters
cluster_ranking <- readRDS(here(output_path, "go_cluster_ranking.rds")) %>%
  filter(class == "acute") %>%
  group_by(cluster, class) %>%
  mutate(label = str_c(n(), " ", "GO terms")) %>%
  ungroup() %>%
  mutate(regulation = fct_rev(str_to_title(regulation))) %>%
  arrange(desc(cluster)) %>%
  mutate(description = fct_rev(as_factor(description)))

cluster_anno <- cluster_ranking %>%
  nest(data = c(rank, term)) %>%
  # find max peak for each cluster
  mutate(peak = data %>% map(function(data) {
    max(density(data$rank)$y)
  })) %>%
  unnest(c(peak)) %>%
  group_by(regulation) %>%
  mutate(max_peak = max(peak)) %>%
  ungroup() %>%
  distinct(
    cluster, regulation, description, peak, max_peak, max_rank,
    label
  ) %>%
  arrange(description) %>%
  group_by(regulation) %>%
  mutate(n_clusters = row_number()) %>%
  mutate(x_coord = case_when(
    n_clusters == 1 ~ 0,
    n_clusters == 2 ~ 1 * max_rank
  )) %>%
  ungroup()

# up-regulated genes
up_dens <- cluster_ranking %>%
  filter(regulation == "Up") %>%
  plot_go_rank_density() +
  my_theme(grid = "no", fsize = fz) +
  geom_text(
    data = filter(cluster_anno, regulation == "Up"),
    aes(x = x_coord, y = 1.1 * max_peak, label = label, color = description),
    inherit.aes = F, size = fz / (14 / 5), hjust = "inward",
    show.legend = F
  ) +
  scale_color_manual(values = aachen_color(c("blue", "maygreen"))) +
  guides(color = guide_legend(nrow = 2)) +
  theme(legend.box.margin = margin(-10, 0, -20, 0))

up_cloud <- wordcounts %>%
  filter(regulation == "up") %>%
  slice_max(order_by = n, n = 15) %>%
  plot_wordcloud(fontsize = fz)

exclusive_acute_up <- up_dens +
  inset_element(up_cloud,
    left = 0, bottom = 0, right = 1, top = 1,
    align_to = "panel"
  )

# down-regulated genes
down_dens <- cluster_ranking %>%
  filter(regulation == "Down") %>%
  plot_go_rank_density() +
  my_theme(grid = "no", fsize = fz) +
  geom_text(
    data = filter(cluster_anno, regulation == "Down"),
    aes(x = x_coord, y = 1.1 * max_peak, label = label, color = description),
    inherit.aes = F, size = fz / (14 / 5), hjust = "inward",
    show.legend = F
  ) +
  scale_color_manual(values = aachen_color("purple")) +
  guides(color = guide_legend(nrow = 2)) +
  theme(legend.box.margin = margin(-10, 0, -20, 0))

down_cloud <- wordcounts %>%
  filter(regulation == "down") %>%
  slice_max(order_by = n, n = 15) %>%
  plot_wordcloud(fontsize = fz)

exclusive_acute_down <- down_dens +
  inset_element(down_cloud,
    left = 0, bottom = 0, right = 1, top = 1,
    align_to = "panel"
  )

exclusive_acute_up + exclusive_acute_down
```

## Collage
### Figure 3
Main Figure.
```{r fig-3, fig.width=8.27, fig.height=11.69}
fig3 <- (((acute_gs_sim / v) + plot_layout(height = c(4, 1))) |
  acute_interstudy_enrichment) /
  grid.grabExpr(draw(exclusive_chronic_hmap)) /
  ((chronic_tfs_up | chronic_pw_up | exclusive_chronic_up | plot_spacer()) +
    plot_layout(width = c(1, 0.5, 2, 1))) /
  grid.grabExpr(draw(common_hmap)) /
  ((common_tfs_up | common_pw_up | common_up | plot_spacer()) +
    plot_layout(width = c(0.5, 1, 2, 1))) +
  plot_layout(height = c(1.5, 1.75, 0.75, 1.75, 0.75)) +
  plot_annotation(tag_levels = list(c(
    "A", "C", "B", "D", "E", "F", "G", "", "H", "I", "J", "K"
  ))) &
  theme(
    plot.tag = element_text(size = fz + 3, face = "bold"),
    legend.key.height = unit(11.5, "pt"),
    legend.key.width = unit(12.5, "pt")
  )

fig3

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

### Supplementary Figure 3.1
Plot to determine the time point of maximal liver damage.
```{r s-fig-3-1}
sfig3_1 <- max_liver_damage +
  theme(
    legend.key.height = unit(11.5, "pt"),
    legend.key.width = unit(12.5, "pt")
  )

sfig3_1

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

### Supplementary Figure 3.2
Direction of regulation for union of differential expressed genes for acute and chronic.
```{r s-fig-3-2}
sfig3_2 <- union_a + union_c +
  plot_layout(widths = c(8, 1)) +
  theme(
    legend.key.height = unit(11.5, "pt"),
    legend.key.width = unit(12.5, "pt")
  )

sfig3_2

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

### Supplementary Figure 3.3
Heatmap of exclusive acute genes and characterization of all up- and down regulated exclusive acute genes.
```{r s-fig-3-3, fig.width=7.87, fig.height=7.87}
sfig3_3 <- plot_spacer() /
  grid.grabExpr(draw(exclusive_acute_hmap)) /
  ((acute_tfs_up | acute_pw_up | exclusive_acute_up) +
    plot_layout(width = c(0.5, 0.5, 2))) /
  ((acute_tfs_down | acute_pw_down | exclusive_acute_down) +
    plot_layout(width = c(0.5, 0.5, 2))) +
  plot_layout(height = c(0, 2, 1, 1)) +
  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")
  )

sfig3_3

ggsave(here("figures/Supplementary Figure 3.3.pdf"), sfig3_3,
  width = 21, height = 20, units = c("cm")
)
ggsave(here("figures/Supplementary Figure 3.3.png"), sfig3_3,
  width = 21, height = 20, units = c("cm")
)
```

### Supplementary Figure 3.4
Characterization of down-regulated common genes.
```{r s-fig-3-4}
sfig3_4 <- (common_tfs_down | common_down) +
  plot_layout(width = c(0.5, 2)) +
  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")
  )

sfig3_4

ggsave(here("figures/Supplementary Figure 3.4.pdf"), sfig3_4,
  width = 21, height = 7, units = c("cm")
)
ggsave(here("figures/Supplementary Figure 3.4.png"), sfig3_4,
  width = 21, height = 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.