13-meta-chronic-vs-acute.Rmd

---
title: "Meta analysis comparing chronic and acute mouse models"
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 integrate various acute liver damage mouse models with the chronic CCl4 mouse model to identify chronic exclusively and commonly regulated genes.

## 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(tidytext)
library(readxl)

library(fgsea)
library(dorothea)
library(progeny)
library(biobroom)

library(circlize)
library(AachenColorPalette)
library(lemon)
library(VennDiagram)
library(ComplexHeatmap)
library(ggwordcloud)

library(msigdf) # remotes::install_github("ToledoEM/msigdf@v7.1")

options("tidylog.display" = list(print))
source(here("code/utils-utils.R"))
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/meta-chronic-vs-acute"
output_path <- "output/meta-chronic-vs-acute"

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

## Merging data of all mouse models 
### Merging contrasts
Contrasts from all available mouse models are merged into a single object.
```{r "merging-contrast"}
# acute
tun <- readRDS(here("output/mouse-acute-tunicamycin/limma_result.rds")) %>%
  mutate(treatment = "tunicamycin", source = "wek", class = "acute")

lps <- readRDS(here("output/mouse-acute-lps/limma_result.rds")) %>%
  filter(contrast == "inLiver_lps_vs_ctrl") %>%
  mutate(treatment = "lps", source = "godoy", class = "acute")

acute_ccl4 <- readRDS(here("output/mouse-acute-ccl4/limma_result.rds")) %>%
  filter(contrast_reference == "ccl4") %>%
  select(-contrast_reference) %>%
  mutate(treatment = "ccl4", source = "godoy", class = "acute")

ph <- readRDS(here("output/mouse-acute-ph/limma_result.rds")) %>%
  filter(contrast_reference == "hepatec") %>%
  select(-contrast_reference) %>%
  mutate(treatment = "ph", source = "godoy", class = "acute")

apap <- readRDS(here("output/mouse-acute-apap/limma_result.rds")) %>%
  filter(contrast_reference == "apap") %>%
  select(-contrast_reference) %>%
  mutate(treatment = "apap", source = "ghallab", class = "acute")

bdl <- readRDS(here("output/mouse-acute-bdl/limma_result.rds")) %>%
  filter(contrast_reference == "bdl") %>%
  select(-contrast_reference) %>%
  mutate(treatment = "bdl", source = "ghallab", class = "acute")

# chronic
chronic_ccl4 <- readRDS(here("output/mouse-chronic-ccl4/limma_result.rds")) %>%
  filter(contrast_reference == "pure_ccl4") %>%
  select(-contrast_reference) %>%
  mutate(
    source = "ghallab", class = "chronic",
    treatment = "pure_ccl4"
  )

combined_contrasts <- bind_rows(
  tun, lps, acute_ccl4, apap,
  ph, bdl, chronic_ccl4
) %>%
  mutate(contrast = as_factor(contrast))

saveRDS(combined_contrasts, here(output_path, "limma_result.rds"))
```

### Merging meta data
Meta data from all available mouse models are merged into a single object.
```{r "merging-meta-data"}
# acute
tun <- readRDS(here("output/mouse-acute-tunicamycin/meta_data.rds")) %>%
  mutate(treatment = "tunicamycin", source = "wek", class = "acute")

lps <- readRDS(here("output/mouse-acute-lps/meta_data.rds")) %>%
  filter(origin == "liver") %>%
  transmute(sample, group = case_when(
    treatment == "control" ~ "control",
    TRUE ~ "treated"
  )) %>%
  mutate(
    treatment = "lps",
    source = "godoy",
    class = "acute"
  )

acute_ccl4 <- readRDS(here("output/mouse-acute-ccl4/meta_data.rds")) %>%
  mutate(
    group = case_when(
      time == 0 ~ "control",
      TRUE ~ "treated"
    ),
    treatment = "ccl4",
    source = "godoy"
  ) %>%
  select(sample, group, time, treatment, source) %>%
  mutate(
    time = as.character(time),
    class = "acute"
  )

ph <- readRDS(here("output/mouse-acute-ph/meta_data.rds")) %>%
  mutate(
    group = case_when(
      time == 0 ~ "control",
      TRUE ~ "treated"
    ),
    treatment = "ph",
    source = "godoy"
  ) %>%
  select(-c(surgeon, year, mouse)) %>%
  mutate(
    time = as.character(time),
    class = "acute"
  )

apap <- readRDS(here("output/mouse-acute-apap/meta_data.rds")) %>%
  mutate(
    group = case_when(
      time == 0 ~ "control",
      TRUE ~ "treated"
    ),
    treatment = "apap",
    source = "ghallab"
  ) %>%
  select(sample, group, time, treatment, source) %>%
  mutate(
    time = as.character(time),
    class = "acute"
  )

bdl <- readRDS(here("output/mouse-acute-bdl/meta_data.rds")) %>%
  mutate(
    group = case_when(
      treatment == "sham" ~ "control",
      TRUE ~ "treated"
    ),
    treatment = "bdl",
    source = "ghallab"
  ) %>%
  select(-mouse) %>%
  mutate(
    time = as.character(time),
    class = "acute"
  )

# chronic
chronic_ccl4 <- readRDS(here("data/mouse-chronic-ccl4/meta_data.rds")) %>%
  mutate(
    group = case_when(
      treatment == "ccl4" ~ "treated",
      treatment %in% c("oil", "ctrl") ~ "control"
    ),
    treatment = "pure_ccl4",
    source = "ghallab"
  ) %>%
  select(sample, group, time, treatment, source) %>%
  mutate(
    time = as.character(time),
    class = "chronic"
  )

combined_meta_data <- bind_rows(
  tun, lps, acute_ccl4, ph, apap, bdl,
  chronic_ccl4
)

saveRDS(combined_meta_data, here(output_path, "meta_data.rds"))
```

## Number of mice
Barplot showing number of mice per study.
```{r mouse-number}
df <- readRDS(here(output_path, "meta_data.rds")) %>%
  count(treatment, source, group, class)

df %>%
  ggplot(aes(
    x = n, fct_reorder(interaction(treatment, class, sep = "_"), n),
    group = group, fill = group
  )) +
  geom_col(position = "dodge") +
  labs(x = "Number of mice", y = "Study") +
  my_theme(grid = "x", fsize = fz)
```

## Gene coverage
Barplot showing the gene coverage of the acute and chronic mouse models.
```{r gene-coverage}
contrasts <- readRDS(here(output_path, "limma_result.rds")) %>%
  distinct(gene, treatment, source, class) %>%
  count(treatment, source, class)

contrasts %>%
  ggplot(aes(
    x = n, fct_reorder(interaction(class, treatment, sep = "_"), n),
    group = class
  )) +
  geom_col() +
  labs(x = "Gene coverage", y = NULL) +
  my_theme(grid = "x", fsize = fz)
```

## PCA of z-scores
To check for batch effects we check the clustering of perturbed samples in PCA space (based on z-scores)
```{r z-score-pca}
# load z-scores and combine into a single data frame
z_scores_mat <- list.files(here("output"),
  pattern = "z_scores.rds",
  recursive = TRUE, full.names = TRUE
) %>%
  # discard human experiments
  keep(str_detect(., pattern = "mouse-")) %>%
  map(function(path) {
    readRDS(path) %>%
      rownames_to_column("gene")
  }) %>%
  reduce(.f = full_join, by = "gene") %>%
  column_to_rownames("gene")

# remove genes with infinity or NA
z_scores_mat <- z_scores_mat[is.finite(rowSums(z_scores_mat)), ]

# load meta data
meta <- readRDS(here(output_path, "meta_data.rds")) %>%
  filter(group == "treated")

pca_result <- do_pca(z_scores_mat, meta, top_n_var_genes = NULL)

saveRDS(pca_result, here(output_path, "z_score_pca.rds"))

plot_pca(pca_result, feature = "treatment") +
  my_theme()
```


## Time point of maximal liver damage
Here we identify for the acute CCl4, APAP and PH mouse model the time point with the strongest deregulated gene expression profile, based on PCA space. This time point is considered the time point of maximal liver damage.
```{r time-point-max-liver-damage}
pca_result_ccl4 <- readRDS(here("output/mouse-acute-ccl4/pca_result.rds")) %>%
  pluck("coords") %>%
  select(sample, PC1, time) %>%
  mutate(
    treatment_abbr = "CCl4",
    time = as.character(time)
  )

pca_result_apap <- readRDS(here("output/mouse-acute-apap/pca_result.rds")) %>%
  pluck("coords") %>%
  select(sample, PC1, time) %>%
  mutate(
    treatment_abbr = "APAP",
    time = as.character(time)
  )

pca_result_ph <- readRDS(here("output/mouse-acute-ph/pca_result.rds")) %>%
  pluck("coords") %>%
  select(sample, PC1, time) %>%
  mutate(
    treatment_abbr = "PH",
    time = as.character(round(24 * as.numeric(as.character(time)), 0))
  )

pca_result <- bind_rows(pca_result_ccl4, pca_result_apap, pca_result_ph) %>%
  mutate(time = ordered(as.numeric(time))) %>%
  arrange(treatment_abbr, time)

pca_dist <- pca_result %>%
  group_by(time, treatment_abbr) %>%
  summarise(mean_pc1_coord = mean(PC1)) %>%
  ungroup() %>%
  arrange(treatment_abbr) %>%
  mutate(control_coord = case_when(time == 0 ~ mean_pc1_coord)) %>%
  fill(control_coord) %>%
  mutate(dist = abs(mean_pc1_coord - control_coord)) %>%
  group_by(treatment_abbr) %>%
  mutate(max = case_when(
    dist == max(dist) ~ "yes",
    TRUE ~ "no"
  )) %>%
  ungroup() %>%
  filter(time != 0)

saveRDS(pca_dist, here(output_path, "pca_dist.rds"))

pca_dist %>%
  ggplot(aes(x = time, y = dist, fill = max)) +
  geom_col() +
  my_theme(grid = "y", fsize = fz) +
  scale_fill_manual(values = aachen_color(c("black50", "green"))) +
  labs(
    x = "Time in hours",
    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")
```

## Interstudy analysis of acute mouse models
### Mutual similarity of differential expressed genes
This analysis computes the similarity of differential expressed genes for specific contrast of the acute mouse models. Similarity is measured with the Jaccard Index.
```{r deg-similarity}
# select specific contrasts from the acute mouse models
contrast_of_interest <- c(
  "treat_vs_ctrl", "inLiver_lps_vs_ctrl",
  "ccl_24h_vs_0h", "apap_24h_vs_0h", "ph_0.5d",
  "bdl_vs_sham_1d"
)

contrasts <- readRDS(here(output_path, "limma_result.rds")) %>%
  filter(contrast %in% contrast_of_interest)

top_n_genes = seq(100, 3000, 100)

j = map_dfr(top_n_genes, function(n) {
  # populate gene sets with a fixed size selected by effect size (t-value)
  mat_top <- contrasts %>%
    group_by(contrast, treatment, source, class) %>%
    top_n(n, abs(statistic)) %>%
    mutate(key = row_number()) %>%
    ungroup() %>%
    unite(geneset, class, source, treatment, contrast, sep = "-") %>%
    mutate(geneset = as_factor(geneset)) %>%
    select(geneset, gene, key) %>%
    untdy(key, geneset, gene)
  
  # usage of jaccard index for balanced set sizes
  set_similarity(mat_top, measure = "jaccard", tidy = T) %>%
    mutate(top_n_genes = n)
  
})

saveRDS(j, here(output_path, "gene_set_similarity.rds"))

# distribution of mean jaccard index with varying gene set size
j %>%
  group_by(top_n_genes) %>%
  summarise(mean_similarity = mean(similarity),
            sd_similarity = sd(similarity)) %>%
  ggplot(aes(x=top_n_genes, y=mean_similarity)) +
  geom_pointline() +
  geom_errorbar(aes(ymin = mean_similarity - sd_similarity,
                    ymax = mean_similarity + sd_similarity)) +
  my_theme() +
  labs(x = "Gene set size", y="Mean Jaccard Index")

# figure with fixed gene set size
j %>%
  filter(top_n_genes == 500) %>%
  ggplot(aes(x = set1, y = set2, fill = similarity)) +
  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)) +
  my_theme(fsize = fz, grid = "no")
```

### Mutual enrichment of differential expressed genes
This analysis explores whether the top differential expressed genes of specific contrasts of the acute mouse models are consistently regulated across the acute mouse models.
```{r mututal-enrichment-degs, fig.height=6}
# select specific contrasts from the acute mouse models
contrast_of_interest <- c(
  "treat_vs_ctrl", "inLiver_lps_vs_ctrl",
  "ccl_24h_vs_0h", "apap_24h_vs_0h", "ph_0.5d",
  "bdl_vs_sham_1d"
)

contrasts <- readRDS(here(output_path, "limma_result.rds")) %>%
  filter(contrast %in% contrast_of_interest)

# construct signature matrix/data frame
signature_df <- contrasts %>%
  unite(signature, class, source, treatment, contrast, sep = "-") %>%
  mutate(signature = as_factor(signature)) %>%
  untdy("gene", "signature", "statistic")

top_n_genes = c(seq(100, 1000, 100), 1500, 2000, 3000, 5000)

gsea_res_top <-  map_dfr(top_n_genes, function(n) {

  # populate gene sets with a fixed size selected by effect size (t-value)
  genesets_top <- contrasts %>%
    mutate(direction = case_when(
      sign(statistic) >= 0 ~ "up",
      sign(statistic) < 0 ~ "down"
    )) %>%
    group_by(class, source, treatment, contrast, direction) %>%
    top_n(n, abs(statistic)) %>%
    ungroup() %>%
    unite(geneset, class, source, treatment, contrast, sep = "-") %>%
    unite(geneset, geneset, direction, sep = "|") %>%
    mutate(geneset = as_factor(geneset)) %>%
    select(geneset, gene)

  # run gsea
  set.seed(123)
  run_gsea(signature_df, genesets_top, tidy = T) %>%
    separate(geneset, into = c("geneset", "direction"), sep = "[|]") %>%
    mutate(
      signature = as_factor(signature),
      geneset = as_factor(geneset),
      top_n_genes = n
    )
})

saveRDS(gsea_res_top, here(output_path, "interstudy_enrichment.rds"))


# distribution of mean enrichment score with varying gene set size
gsea_res_top %>%
  filter(abs(ES) != 1) %>%
  group_by(top_n_genes, direction) %>%
  summarise(mean_es = mean(NES),
            sd_es = sd(NES)) %>%
  ggplot(aes(x=top_n_genes, y=mean_es)) +
  geom_pointline() +
  geom_errorbar(aes(ymin = mean_es - sd_es,
                    ymax = mean_es + sd_es)) +
  geom_hline(yintercept = 0, lty = "dashed") +
  my_theme() +
  labs(x = "Gene set size", y="Mean ES") +
  facet_rep_wrap(~direction)


# # individual gsea plot
# make_gsea_plot(
#   signature = signature_df %>%
#     rownames_to_column("gene") %>%
#     select(gene, t = `acute-wek-tunicamycin-treat_vs_ctrl`) %>%
#     drop_na(),
#   geneset = genesets_top %>%
#     filter(geneset == "acute-godoy-ccl4-ccl_24h_vs_0h|up"),
#   gene_level_stat = t
# )

# summary gsea plot with a fixed gene set size
gsea_res_top %>%
  filter(top_n_genes == 500) %>%
  mutate(label = gtools::stars.pval(padj)) %>%
  ggplot(aes(x = signature, y = geneset, fill = ES)) +
  geom_tile() +
  geom_text(aes(label = label)) +
  facet_wrap(~direction) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_fill_gradient2() +
  my_theme(fsize = fz, grid = "no") +
  labs(x = "Signature", y = "Gene set")
```

## Construction of unified acute and chronic genes
In this analysis we want to identify genes specific for  chronic or acute liver damage and genes that are deregulated in both scenarios.

### Acute
#### Pool of differential expressed genes in acute mouse models
Filter for differential expressed genes of the acute mouse models. Visual inspection suggest that the differential expressed genes are consistently regulated across the acute mouse models.
```{r acute-pool}
contrasts <- readRDS(here(output_path, "limma_result.rds")) %>%
  assign_deg(fdr_cutoff = 1e-4)

acute_gene_pool <- contrasts %>%
  filter(class == "acute") %>%
  filter(regulation != "ns") %>%
  # remove late bdl time point as this could already be a chronic damage
  filter(contrast != "bdl_vs_sham_21d")

acute_gene_pool %>%
  mutate(statistic = case_when(
    statistic >= 25 ~ 25,
    TRUE ~ statistic
  )) %>%
  ggplot(aes(
    x = contrast, y = fct_reorder(gene, statistic, mean),
    fill = statistic
  )) +
  geom_tile() +
  scale_fill_gradient2() +
  theme(
    axis.text.x = element_text(angle = 90),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) +
  labs(y = "Genes", x = NULL) +
  my_theme(grid = "no", fsize = fz)

saveRDS(acute_gene_pool, here(output_path, "acute_gene_pool.rds"))
```

#### Unify acute genes
Acute genes are unified and a median t-statistic is computed for each gene.
```{r acute-union}
acute_gene_pool <- readRDS(here(output_path, "acute_gene_pool.rds"))

acute_gene_union <- acute_gene_pool %>%
  group_by(gene) %>%
  summarise(
    m = mean(sign(statistic)), n = n(),
    median_statistic = median(statistic),
    median_logFC = median(logFC)
  ) %>%
  distinct(gene, median_statistic, median_logFC)

saveRDS(acute_gene_union, here(output_path, "union_acute_geneset.rds"))
```

### Chronic
#### Pool of differential expressed genes in chronic mouse model
Filter for differential expressed genes of the chronic mouse model. Visual inspection suggest that the differential expressed genes are consistently regulated across the chronic contrasts.
```{r chronic-pool}
contrasts <- readRDS(here(output_path, "limma_result.rds")) %>%
  assign_deg(fdr_cutoff = 1e-4)

chronic_gene_pool <- contrasts %>%
  filter(class == "chronic") %>%
  filter(regulation != "ns")

chronic_gene_pool %>%
  mutate(statistic = case_when(
    statistic >= 25 ~ 25,
    TRUE ~ statistic
  )) %>%
  ggplot(aes(
    x = contrast, y = fct_reorder(gene, statistic, mean),
    fill = statistic
  )) +
  geom_tile() +
  scale_fill_gradient2() +
  theme(
    axis.text.x = element_text(angle = 90),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) +
  labs(y = "Genes", x = NULL) +
  my_theme(grid = "no", fsize = fz)

saveRDS(chronic_gene_pool, here(output_path, "chronic_gene_pool.rds"))
```

#### Unify chronic genes
Chronic genes are unified and a median t-statistic is computed for each gene.
```{r chronic-union}
chronic_gene_pool <- readRDS(here(output_path, "chronic_gene_pool.rds"))

chronic_gene_union <- chronic_gene_pool %>%
  group_by(gene) %>%
  summarise(
    m = mean(sign(statistic)), n = n(),
    median_statistic = median(statistic),
    median_logFC = median(logFC)
  ) %>%
  distinct(gene, median_statistic, median_logFC)

saveRDS(chronic_gene_union, here(output_path, "union_chronic_geneset.rds"))
```

### Overlap of unified gene sets
Venn diagram showing the overlap of unified chronic and acute genes.
```{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()

v <- 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)
)
```

### Assign membership to each gene
Genes are assigned a membership: i) exclusive chronic, ii) exclusive acute, iii) common. 
```{r assign-membership}
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")

# assign membership to the genes
m <- bind_rows(acute_gene_union, chronic_gene_union) %>%
  add_count(gene) %>%
  mutate(membership = case_when(
    n == 2 ~ "common",
    n == 1 & class == "acute" ~ "acute",
    n == 1 & class == "chronic" ~ "chronic"
  )) %>%
  dplyr::select(-n)

saveRDS(m, here(output_path, "gene_membership.rds"))
```

## Extraction of top exclusive and common genes
For each membership class the the genes are ranked and the top genes are visualized.

### Exclusive chronic
Exclusive chronic genes are ranked based on a metric that prioritizes genes that have a high consensus chronic gene-level statistic and at the same time are consistently not deregulated in selected acute contrasts.
```{r rank-exclusive-chronic-genes}
exclusive_chronic_genes <- readRDS(here(output_path, "gene_membership.rds")) %>%
  filter(membership == "chronic" & class == "chronic") %>%
  distinct(gene, chronic_statistic = median_statistic)

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

df <- contrasts %>%
  filter(contrast %in% acute_contrasts) %>%
  inner_join(exclusive_chronic_genes, by = "gene") %>%
  group_by(gene, chronic_statistic) %>%
  summarise(acute_statistic = median(statistic), var = var(statistic), n = n()) %>%
  ungroup()

ranked_exclusive_chronic_genes <- df %>%
  # consider only genes that are available in at least in 5 acute contrasts
  filter(n >= 5) %>%
  # compute empirical metric that maximizes if the chronic statistic is high,
  # and the acute statistic and variance is low
  mutate(importance = chronic_statistic * (1 / acute_statistic) * sqrt(1 / var)) %>%
  arrange(-abs(importance), -chronic_statistic) %>%
  mutate(rank = row_number())

saveRDS(
  ranked_exclusive_chronic_genes,
  here(output_path, "ranked_exclusive_chronic_genes.rds")
)
```

Extraction of the top 100 exclusive chronic genes. Their expression in acute and chronic mouse models is visualized in a heatmap.
```{r vis-exclusive-chronic-genes}
ranked_exclusive_chronic_genes <- readRDS(
  here(output_path, "ranked_exclusive_chronic_genes.rds")
)

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

Heatmap(t(mat_exclusive_chronic_genes),
  col = col_fun,
  cluster_rows = F, cluster_columns = T,
  row_names_gp = gpar(fontsize = fz), column_names_gp = gpar(fontsize = fz - 4),
  name = "logFC",
  row_gap = unit(2.5, "mm"),
  border = T,
  row_split = c(rep("Acute", 12), rep("Chronic", 3))
)
```

### Common
Common genes are ranked based on a metric that prioritize genes that have a high consensus chronic gene-level statistic and at the same are time consistently regulated in the same direction as in the chronic scenario in selected acute contrasts. 
```{r rank-common-genes}
common_genes <- readRDS(here(output_path, "gene_membership.rds")) %>%
  filter(membership == "common" & class == "chronic") %>%
  distinct(gene, chronic_statistic = median_statistic)

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

df <- contrasts %>%
  filter(contrast %in% acute_contrasts) %>%
  inner_join(common_genes, by = "gene") %>%
  group_by(gene, chronic_statistic) %>%
  summarise(acute_statistic = median(statistic), var = var(statistic), n = n()) %>%
  ungroup()

ranked_common_genes <- df %>%
  # consider only genes that are available in at least in 5 acute contrasts
  filter(n >= 5) %>%
  # compute empirical metric that maximizes if the chronic and acute statistic
  # is high and acute variance is low
  mutate(importance = chronic_statistic * acute_statistic * sqrt(1 / var)) %>%
  arrange(-abs(importance), -chronic_statistic) %>%
  mutate(rank = row_number(-importance))

saveRDS(ranked_common_genes, here(output_path, "ranked_common_genes.rds"))
```

Extraction of the top 100 common genes. Their expression in acute and chronic mouse models is visualized in a heatmap.
```{r vis-common-genes}
ranked_common_genes <- readRDS(
  here(output_path, "ranked_common_genes.rds")
)

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

Heatmap(t(mat_common_genes),
  col = col_fun,
  cluster_rows = F, cluster_columns = T,
  row_names_gp = gpar(fontsize = fz), column_names_gp = gpar(fontsize = fz - 4),
  name = "logFC",
  row_gap = unit(2.5, "mm"),
  border = T,
  row_split = c(rep("Acute", 12), rep("Chronic", 3))
)
```

### Exclusive acute
Exclusive acute genes are ranked based on a metric that prioritize genes that have a high consensus acute gene-level statistic and at the same time are consistently not deregulated in the chronic contrasts.

```{r rank-exclusive-acute-genes}
exclusive_acute_genes <- readRDS(here(output_path, "gene_membership.rds")) %>%
  filter(membership == "acute" & class == "acute") %>%
  distinct(gene, acute_statistic = median_statistic)

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

df <- contrasts %>%
  filter(treatment == "pure_ccl4") %>%
  inner_join(exclusive_acute_genes, by = "gene") %>%
  group_by(gene, acute_statistic) %>%
  summarise(chronic_statistic = median(statistic), var = var(statistic), n = n()) %>%
  ungroup()

ranked_exclusive_acute_genes <- df %>%
  # compute empirical metric that maximizes if the acute statistic is high,
  # and the chronic statistic and variance is low
  mutate(importance = (1 / chronic_statistic) * acute_statistic * sqrt(1 / var)) %>%
  arrange(-abs(importance), -acute_statistic) %>%
  mutate(rank = row_number())

saveRDS(
  ranked_exclusive_acute_genes,
  here(output_path, "ranked_exclusive_acute_genes.rds")
)
```

Extraction of the top 100 exclusive acute genes. Their expression in acute and chronic mouse models is visualized in a heatmap.
```{r vis-exclusive-acute-genes}
ranked_exclusive_acute_genes <- readRDS(
  here(output_path, "ranked_exclusive_acute_genes.rds")
)

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

Heatmap(t(mat_exclusive_acute_genes),
  col = col_fun,
  cluster_rows = F, cluster_columns = T,
  row_names_gp = gpar(fontsize = fz), column_names_gp = gpar(fontsize = fz - 4),
  name = "logFC",
  row_gap = unit(2.5, "mm"),
  border = T,
  row_split = c(rep("Acute", 12), rep("Chronic", 3))
)
```

## Characterization of exclusive and common genes
All exclusive chronic, exclusive acute and common genes are characterized GO terms, [PROGENy's](http://saezlab.github.io/progeny/) pathways and [DoRothEA's](http://saezlab.github.io/dorothea/) TFs. As statistic over-representation analysis is used.
```{r exclusive-common-gene-characterizations}
exclusive_chronic_genes <- readRDS(
  here(output_path, "ranked_exclusive_chronic_genes.rds")
) %>%
  select(gene, statistic = chronic_statistic, rank) %>%
  mutate(class = "chronic")

exclusive_acute_genes <- readRDS(
  here(output_path, "ranked_exclusive_acute_genes.rds")
) %>%
  select(gene, statistic = acute_statistic, rank) %>%
  mutate(class = "acute")

common_genes <- readRDS(here(output_path, "ranked_common_genes.rds")) %>%
  filter(sign(chronic_statistic) == sign(acute_statistic)) %>%
  select(gene, statistic = chronic_statistic, rank) %>%
  mutate(class = "common")

signatures <- bind_rows(
  exclusive_chronic_genes, exclusive_acute_genes,
  common_genes
) %>%
  mutate(regulation = ifelse(statistic >= 0, "up", "down"))

# load gene sets of GO terms, pathways, tfs)
genesets <- load_genesets() %>%
  filter(confidence %in% c(NA, "A", "B", "C"))

# run over-representation analysis
ora_res <- signatures %>%
  nest(sig = c(-class, -regulation)) %>%
  dplyr::mutate(ora = sig %>% map(run_ora,
    sets = genesets, min_size = 10,
    options = list(alternative = "greater"),
    background_n = 20000
  )) %>%
  select(-sig) %>%
  unnest(ora)

saveRDS(ora_res, here(output_path, "exclusive_genes_characterization.rds"))
```

## GO downstream analysis
In this section significant GO terms are summarized in two different ways. i) Which words appear the most among the GO terms and ii) Position of members of manually created GO cluster in a ranked list of significant GO terms (based on p-value)

### Text analysis
GO terms are splitted into words and their frequency is counted.
```{r go-text-analysis}
# list of words that will be ignored
stop_go_words <- readRDS(here("data/annotation/stop_go_words.rds"))

go_terms <- readRDS(here(output_path, "exclusive_genes_characterization.rds")) %>%
  filter(group == "go" & fdr <= 0.05) %>%
  mutate(
    term = str_remove(geneset, "GO_"),
    term = str_replace_all(term, "_", " "),
    term = str_to_lower(term)
  ) %>%
  select(regulation, term, class)

go_wordcounts <- go_terms %>%
  # glue words that should be treated as one
  mutate(
    term = str_replace(term, "smooth muscle", "smoothmuscle"),
    term = str_replace(term, "amino acid", "aminoacid"),
    term = str_replace(term, "cell cycle", "cellcycle"),
    term = str_replace(term, "endoplasmic reticulum stress", "endoplasmicreticulumstress"),
    term = str_replace(term, "endoplasmic reticulum", "endoplasmicreticulum")
  ) %>%
  unnest_tokens(word, term) %>%
  anti_join(stop_words, by = "word") %>%
  # correct for abbreviations
  mutate(word = case_when(
    word == "er" ~ "endoplasmicreticulum",
    TRUE ~ word
  )) %>%
  # remove words which are pure numbers
  filter(!str_detect(word, "^[0-9]+")) %>%
  # count word frequency
  count(class, regulation, word, sort = T) %>%
  # dissect prior glued words
  mutate(word = case_when(
    word == "smoothmuscle" ~ "smooth-muscle",
    word == "aminoacid" ~ "amino-acid",
    word == "cellcycle" ~ "cell-cycle",
    word == "endoplasmicreticulumstress" ~ "endoplasmic-reticulum-stress",
    word == "endoplasmicreticulum" ~ "endoplasmic-reticulum",
    TRUE ~ word
  )) %>%
  # remove meaning less words
  anti_join(stop_go_words, by = "word")

saveRDS(go_wordcounts, here(output_path, "go_wordcounts.rds"))


go_wordcounts %>%
  group_by(class, regulation) %>%
  slice_max(order_by = n, n = 15) %>%
  ungroup() %>%
  ggplot(aes(label = word, size = n)) +
  geom_text_wordcloud() +
  scale_size_area(max_size = fz / (14 / 5)) +
  facet_grid(class ~ regulation) +
  theme(axis.line = element_blank()) +
  my_theme(grid = "no", fsize = fz)
```

### Cluster ranking
Distribution of manually created GO-cluster.
```{r go-cluster-ranking}
go_cluster_mapping <- tribble(
  ~cluster, ~description,
  4, "Proliferation",
  5, "Metabolism",
  6, "Migration",
  7, "Development and Morphogenesis",
  9, "ER stress",
  10, "Proliferation",
  11, "Metabolism"
)

go_cluster_chronic <- read_excel(here(data_path, "manual_go_cluster_anno.xlsx"),
  sheet = "chronic"
) %>%
  mutate(class = "chronic")
go_cluster_acute <- read_excel(here(data_path, "manual_go_cluster_anno.xlsx"),
  sheet = "acute"
) %>%
  mutate(class = "acute")
go_cluster_common <- read_excel(here(data_path, "manual_go_cluster_anno.xlsx"),
  sheet = "common"
) %>%
  mutate(class = "common")

go_cluster <- bind_rows(go_cluster_chronic, go_cluster_acute, go_cluster_common) %>%
  distinct(class, regulation, term, cluster) %>%
  drop_na() %>%
  inner_join(go_cluster_mapping)

go_terms_ranking <- readRDS(
  here(output_path, "exclusive_genes_characterization.rds")
) %>%
  mutate(
    term = str_remove(geneset, "GO_"),
    term = str_replace_all(term, "_", " "),
    term = str_to_lower(term)
  ) %>%
  filter(fdr <= 0.05) %>%
  group_by(regulation, class) %>%
  mutate(rank = row_number(p.value)) %>%
  mutate(max_rank = max(rank)) %>%
  ungroup() %>%
  select(class, regulation, term, rank, max_rank)

go_cluster_ranking <- go_terms_ranking %>%
  inner_join(go_cluster)

saveRDS(go_cluster_ranking, here(output_path, "go_cluster_ranking.rds"))

go_cluster_ranking %>%
  mutate(regulation = as_factor(regulation)) %>%
  ggplot(aes(x = rank, colour = description)) +
  stat_density(geom = "line", position = "identity") +
  geom_rug() +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    legend.position = "top"
  ) +
  labs(
    x = str_wrap("GO term ranking based on p-value"), y = "Density",
    color = NULL
  ) +
  facet_rep_wrap(~ class + regulation, scales = "free", drop = T) +
  scale_x_continuous(breaks = function(x) {
    as.numeric(
      gsub("^0", 1, unique(floor(pretty(seq(0, (max(x) + 1) * 1.1)))))
    )
  }) +
  my_theme(grid = "no", fsize = fz)
```

```{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.