plot-mouse-vs-human.Rmd

---
title: "Mouse vs. Human (Fig. 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 for the comparison of the chronic CCl4 model and patient cohorts.

## 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(gridExtra)
library(ggpubr)
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-mouse-vs-human"
output_path <- "output/meta-mouse-vs-human"

# 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"))
key_hs <- readRDS(here("data/meta-mouse-vs-human/contrast_annotation.rds"))
```

## Similarity of patient cohort gene sets
```{r sim-patient-gene-sets}
keys <- key_hs %>%
  distinct(contrast, label, source, phenotype)

j <- readRDS(here(output_path, "gene_set_similarity.rds")) %>%
  separate(set1, into = c("source", "phenotype", "contrast"), sep = "-") %>%
  inner_join(keys, by = c("contrast", "source", "phenotype")) %>%
  select(-source, -phenotype, -contrast) %>%
  rename(set1 = label) %>%
  separate(set2, into = c("source", "phenotype", "contrast"), sep = "-") %>%
  inner_join(keys, by = c("contrast", "source", "phenotype")) %>%
  rename(set2 = label) %>%
  select(set1, set2, similarity)

patient_gs_sim <- j %>%
  mutate(set1 = fct_rev(set1)) %>%
  ggplot(aes(
    x = set1, y = set2, fill = similarity,
    label = round(similarity, 3)
  )) +
  geom_tile(color = "black") +
  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()
  ) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  geom_text(size = (fz - 2) / (14 / 5)) +
  my_theme(fsize = fz, grid = "no")

patient_gs_sim
```

## Interstudy enrichment of patient cohorts
```{r interstudy-enrichment-patient}
keys <- key_hs %>%
  distinct(contrast, label, source, phenotype)

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

patient_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 = 2) +
  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"))

patient_interstudy_enrichment
```

## Mouse enrichment in patients
```{r mouse-enrichment-in-patients}
contrast_keys_mouse <- key_mm %>%
  distinct(signature = contrast, label2)
contrast_keys_human <- key_hs %>%
  unite(geneset, source, phenotype, contrast, sep = "-") %>%
  distinct(geneset, label)

gsea_res <- readRDS(here(output_path, "gsea_res.rds")) %>%
  inner_join(contrast_keys_mouse, by = "signature") %>%
  select(-signature) %>%
  rename(signature = label2) %>%
  inner_join(contrast_keys_human, by = "geneset") %>%
  select(-geneset) %>%
  rename(geneset = label) %>%
  mutate(
    label = gtools::stars.pval(padj),
    direction = fct_rev(str_to_title(direction))
  )

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

mm_enrichment_in_hs
```

## Overlap of leading-edge genes
```{r leading-edge-overlap}
keys <- key_mm %>%
  distinct(signature = contrast, class = time_label2)

le <- readRDS(here(output_path, "leading_edges.rds")) %>%
  inner_join(keys, by = "signature") %>%
  rename(regulation = direction) %>%
  arrange(class)

tables <- le %>%
  group_split(class)

# extract labellers
c1 <- tables[[1]] %>%
  distinct(class) %>%
  pull() %>%
  as.character()
c2 <- tables[[2]] %>%
  distinct(class) %>%
  pull() %>%
  as.character()
c3 <- tables[[3]] %>%
  distinct(class) %>%
  pull() %>%
  as.character()

t1 <- tables[[1]] %>% count(regulation)
t2 <- tables[[2]] %>% count(regulation)
t3 <- tables[[3]] %>% count(regulation)

le_overlap_plots <- c("up", "down") %>%
  map(function(r) {
    # set sizes of regulated genes
    a1 <- t1 %>%
      filter(regulation == r) %>%
      pull(n)
    a2 <- t2 %>%
      filter(regulation == r) %>%
      pull(n)
    a3 <- t3 %>%
      filter(regulation == r) %>%
      pull(n)

    a12 <- purrr::reduce(
      list(
        tables[[1]] %>% filter(regulation == r) %>% pull(gene),
        tables[[2]] %>% filter(regulation == r) %>% pull(gene)
      ),
      intersect
    ) %>%
      length()

    a23 <- purrr::reduce(
      list(
        tables[[2]] %>% filter(regulation == r) %>% pull(gene),
        tables[[3]] %>% filter(regulation == r) %>% pull(gene)
      ),
      intersect
    ) %>%
      length()

    a13 <- purrr::reduce(
      list(
        tables[[1]] %>% filter(regulation == r) %>% pull(gene),
        tables[[3]] %>% filter(regulation == r) %>% pull(gene)
      ),
      intersect
    ) %>%
      length()

    a123 <- purrr::reduce(
      list(
        tables[[1]] %>% filter(regulation == r) %>% pull(gene),
        tables[[2]] %>% filter(regulation == r) %>% pull(gene),
        tables[[3]] %>% filter(regulation == r) %>% pull(gene)
      ),
      intersect
    ) %>%
      length()



    grid.newpage()
    grid.grabExpr(draw.triple.venn(
      area1 = a1, area2 = a2, area3 = a3,
      n12 = a12, n23 = a23, n13 = a13,
      n123 = a123,
      category = c(c1, c2, c3),
      # lty = "blank",
      cex = 1 / 12 * fz,
      fontfamily = rep("sans", 7),
      # fill = aachen_color(c("purple", "petrol", "red")),
      # cat.col = aachen_color(c("purple", "petrol", "red")),
      cat.cex = 1 / 12 * (fz + 1),
      cat.fontfamily = rep("sans", 3),
      cat.pos = c(0, 0, 180),
      cat.prompts = T
    )) %>%
      as_ggplot() %>%
      grid.arrange(top = textGrob(str_to_title(r), gp = gpar(
        fontsize = fz,
        fontface = "bold"
      )))
  })

le_gene_overlap <- wrap_plots(le_overlap_plots)
```

## Top human-mouse consistent genes
### Heatmap
```{r top-human-mouse-consistent-genes}
contrast_keys_human <- key_hs %>%
  unite(contrast, source, phenotype, contrast, sep = "-") %>%
  distinct(contrast, label)

contrast_keys_mouse <- key_mm %>%
  filter(str_detect(contrast, "pure")) %>%
  distinct(contrast, label)

keys <- bind_rows(contrast_keys_mouse, contrast_keys_human)

# get degs from single cell data
sc_degs <- readRDS(here(data_path, "single_cell_degs.rds")) %>%
  enframe("celltype") %>%
  unnest(value) %>%
  filter(p_val_adj <= 0.05) %>%
  mutate(adjusted_logfc = case_when(
    avg_logFC >= 0 & cluster == "Uninjured" ~ -avg_logFC,
    avg_logFC < 0 & cluster == "Uninjured" ~ -avg_logFC,
    TRUE ~ avg_logFC
  ))

# load consistent genes
df <- readRDS(here(output_path, "consistent_genes.rds")) %>%
  left_join(sc_degs, by = "gene") %>%
  mutate(celltype = fct_explicit_na(celltype, na_level = "Unknown"))

# extract top 100 consistent genes
top_genes_df <- df %>%
  filter(rank <= 100) %>%
  inner_join(keys, by = "contrast") %>%
  mutate(label = fct_drop(label))

# build matrix for heatmap
m <- top_genes_df %>%
  distinct(gene, label, logFC) %>%
  spread(label, logFC) %>%
  data.frame(row.names = 1, check.names = F) %>%
  as.matrix()

# build cell type annotation
celltype_anno <- top_genes_df %>%
  distinct(gene, celltype, adjusted_logfc) %>%
  spread(celltype, adjusted_logfc, fill = 0) %>%
  data.frame(row.names = 1, check.names = F) %>%
  select(-Unknown)

col_fun_anno <- col_fun
ha <- HeatmapAnnotation(
  df = celltype_anno,
  show_legend = F,
  border = F,
  gap = unit(0.25, "mm"),
  col = list(
    MPs = col_fun_anno, pDCs = col_fun_anno, ILCs = col_fun_anno,
    Tcells = col_fun_anno, Bcells = col_fun_anno,
    `Plasma Bcells` = col_fun_anno, `Mast cells` = col_fun_anno,
    Endothelia = col_fun_anno, Mesenchyme = col_fun_anno,
    Mesothelia = col_fun_anno, Hepatocytes = col_fun_anno,
    Cholangiocytes = col_fun_anno
  ),
  annotation_name_gp = gpar(fontsize = fz - 1),
  simple_anno_size = unit(0.15, "cm")
)

# build heatmap
h_hmap <- Heatmap(
  t(as.matrix(m)),
  col = col_fun,
  cluster_rows = F,
  cluster_columns = T,
  show_row_dend = F,
  show_column_dend = T,
  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,
  top_annotation = ha,
  row_split = c(rep("Mouse", 3), rep("Human", 15)),
  row_title_gp = gpar(fontsize = fz + 1),
)

h_hmap
```

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

tfs_up <- ora_res %>%
  filter(regulation == "Up") %>%
  filter(geneset %in% c("SP1", "RELA", "NFKB1")) %>%
  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)

# no significant hits
# 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)

tfs_up
```

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

pw_up <- ora_res %>%
  filter(regulation == "Up") %>%
  filter(geneset %in% c("TGFb", "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)

pw_down <- ora_res %>%
  filter(regulation == "Down") %>%
  filter(geneset %in% c("Androgen")) %>%
  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)

pw_up + pw_down
```

### Characterization via GO terms
```{r consistent-go}
# get wordcounts for go terms
wordcounts <- readRDS(here(output_path, "go_wordcounts.rds"))

# get ranking of go clusters
cluster_ranking <- readRDS(here(output_path, "go_cluster_ranking.rds")) %>%
  group_by(cluster) %>%
  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("green", "purple"))) +
  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)

consistent_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)

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

consistent_up + consistent_down
```

## Expression of consistent genes in mouse
```{r expression-consistent-genes-in-mouse}
keys <- key_mm %>%
  distinct(contrast, label = time_label2)

le <- readRDS(here(output_path, "leading_edges_mgi.rds")) %>%
  rename(contrast = signature)

c <- readRDS(here("output/mouse-chronic-ccl4/limma_result.rds")) %>%
  filter(contrast_reference == "pure_ccl4") %>%
  left_join(le) %>%
  select(-regulation) %>%
  rename(regulation = direction) %>%
  replace_na(list(regulation = "ns")) %>%
  mutate(regulation = factor(regulation, levels = c("up", "down", "ns"))) %>%
  mutate(regulation = fct_recode(regulation,
    Up = "up", Down = "down",
    n.s. = "ns"
  )) %>%
  inner_join(keys, by = "contrast") %>%
  arrange(desc(regulation))

consistent_volcano_main <- c %>%
  filter(label == "Month 12") %>%
  ggplot(aes(
    x = logFC, y = -log10(pval), color = regulation,
    alpha = regulation
  )) +
  geom_point() +
  facet_rep_wrap(~label, scales = "free") +
  my_theme(grid = "y", fsize = fz) +
  scale_alpha_manual(values = c(0.7, 0.7, 0.2), guide = "none", drop = F) +
  scale_color_manual(
    values = c(aachen_color(c("red", "blue", "black50"))),
    breaks = c("Up", "Down"), labels = c("Up", "Down"),
    drop = F
  ) +
  labs(
    x = "logFC", y = expression(-log["10"] * "(p-value)"),
    color = str_wrap("Regulation", width = 15)
  )

consistent_volcano_supp <- c %>%
  filter(label != "Month 12") %>%
  ggplot(aes(
    x = logFC, y = -log10(pval), color = regulation,
    alpha = regulation
  )) +
  geom_point() +
  facet_rep_wrap(~label, scales = "free") +
  my_theme(grid = "y", fsize = fz) +
  scale_alpha_manual(values = c(0.7, 0.7, 0.2), guide = "none", drop = F) +
  scale_color_manual(
    values = c(aachen_color(c("red", "blue", "black50"))),
    breaks = c("Up", "Down"), labels = c("Up", "Down"),
    drop = F
  ) +
  labs(
    x = "logFC", y = expression(-log["10"] * "(p-value)"),
    color = str_wrap("Regulation", width = 15)
  )

consistent_volcano_main + consistent_volcano_supp
```

## Collage
### Figure 4
Main Figure.
```{r fig-4, fig.width=8.27, fig.height=11.69}

fig4 <- patient_interstudy_enrichment /
  ((mm_enrichment_in_hs | le_overlap_plots[[1]] | le_overlap_plots[[2]]) +
    plot_layout(widths = c(1, 2.25, 1.75))) /
  grid.grabExpr(draw(h_hmap)) /
  ((pw_up | tfs_up | consistent_up | consistent_volcano_main) +
    plot_layout(width = c(0.2, 0.3, 1.5, 1.5))) +
  plot_layout(height = c(0.75, 0.75, 1.75, 0.75)) +
  plot_annotation(tag_levels = list(c(
    "A", "B", "C", "", "D", "E", "F", "G",
    "", "H"
  ))) &
  theme(
    plot.tag = element_text(size = fz + 3, face = "bold"),
    legend.key.height = unit(11.5, "pt"),
    legend.key.width = unit(12.5, "pt")
  )

fig4

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

### Supplementary Figure 4.1
Gene similarity of patient cohorts.
```{r s-fig-4-1}
sfig4_1 <- patient_gs_sim +
  theme(
    legend.key.height = unit(11.5, "pt"),
    legend.key.width = unit(12.5, "pt")
  )

sfig4_1

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

### Supplementary Figure 4.2
Characterization of down-regulated consistent genes.
```{r s-fig-4-2}
sfig4_2 <- pw_down + consistent_down +
  plot_layout(widths = c(1, 8)) +
  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")
  )

sfig4_2

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