mouse-acute-ccl4.Rmd

---
title: "Acute CCl4 mouse model"
author: "Christian H. Holland"
date: "2020-12-19"
output: 
  workflowr::wflow_html:
    code_folding: hide
editor_options:
  chunk_output_type: console
---

```{r chunk-setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  autodep = TRUE,
  cache = TRUE,
  message = FALSE,
  warning = FALSE
)
```

## Introduction
Here we analysis a mouse model of CCl4 induced acute liver damage. The transcriptomic profiles were measured at 9 different time points ranging from 2 hours to 16 days.

## Libraries and sources
These libraries and sources are used for this analysis.
```{r libs-and-src, message=FALSE, warning=FALSE, cache=FALSE}
library(mouse4302.db)

library(tidyverse)
library(tidylog)
library(here)

library(oligo)
library(annotate)
library(limma)
library(biobroom)
library(progeny)
library(dorothea)

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

library(AachenColorPalette)
library(cowplot)
library(lemon)
library(patchwork)

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

# graphical parameters
# fontsize
fz <- 9
```

## Data processing
### Load .CEL files and quality control
The array quality is controlled based on the relative log expression values (RLE) and the normalized unscaled standard errors (NUSE).
```{r load-cel-files}
# load cel files and check quality
platforms <- readRDS(here("data/annotation/platforms.rds"))
raw_eset <- list.celfiles(here(data_path), listGzipped = T, full.names = T) %>%
  read.celfiles() %>%
  ma_qc() # Discarding in total 1 arrays: 6h_m3_(Mouse430_2).CEL
```

### Normalization and probe annotation
Probe intensities are normalized with the `rma()` function. Probes are annotated with MGI symbols.
```{r normalization-and-annotation}
eset <- rma(raw_eset)

# annotate microarray probes with mgi symbols
expr <- ma_annotate(eset, platforms)
pattern <- str_c(c(
  "Patricio_",
  "_\\(Mouse430_2\\).CEL",
  "_\\(Mouse430_2\\)_2.CEL"
), collapse = "|")

colnames(expr) <- str_c(
  "sample", str_remove_all(colnames(expr), pattern),
  sep = "_"
) %>%
  str_to_lower()

# save normalized expression
saveRDS(expr, here(output_path, "normalized_expression.rds"))
```

### Build meta data
Meta information are parsed from the sample names.
```{r build-meta-data}
# build meta data
meta <- colnames(expr) %>%
  enframe(name = NULL, value = "sample") %>%
  separate(sample, into = c("tmp", "key", "run"), remove = F) %>%
  dplyr::select(-tmp) %>%
  separate(key, into = c("group", "rep"), sep = "m", remove = T) %>%
  mutate(time = case_when(
    group == "k" ~ 0,
    str_detect(group, "h") ~ parse_number(group),
    str_detect(group, "d") ~ parse_number(group) * 24
  )) %>%
  mutate(time = ordered(time)) %>%
  mutate(group = case_when(
    group == "k" ~ "control",
    str_detect(group, "h") ~ str_c("h", parse_number(group)),
    TRUE ~ group
  )) %>%
  mutate(group = factor(group, levels = c(
    "control", "h2", "h8", "d1", "d2",
    "d4", "d6", "d8", "d16"
  ))) %>%
  arrange(sample)

# save meta data
saveRDS(meta, here(output_path, "meta_data.rds"))
```

## Exploratory analysis
### PCA of normalized data
PCA plot of normalized expression data contextualized based on the time point. Only the top 1000 most variable genes are used as features.
```{r pca-norm-data}
expr <- readRDS(here(output_path, "normalized_expression.rds"))
meta <- readRDS(here(output_path, "meta_data.rds"))

pca_result <- do_pca(expr, meta, top_n_var_genes = 1000)

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

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

## Differential gene expression analysis
### Running limma
Differential gene expression analysis via limma with the aim to identify the effect of CCl4 intoxication for the different time points.
```{r running-limma}
# load expression and meta data
expr <- readRDS(here(output_path, "normalized_expression.rds"))
meta <- readRDS(here(output_path, "meta_data.rds"))

stopifnot(colnames(expr) == meta$sample)

# build design matrix
design <- model.matrix(~ 0 + group, data = meta)
rownames(design) <- meta$sample
colnames(design) <- levels(meta$group)

# define contrasts
contrasts <- makeContrasts(
  # effect of ccl4 treatment
  ccl_2h_vs_0h = h2 - control,
  ccl_8h_vs_0h = h8 - control,
  ccl_24h_vs_0h = d1 - control,
  ccl_48h_vs_0h = d2 - control,
  ccl_96h_vs_0h = d4 - control,
  ccl_144h_vs_0h = d6 - control,
  ccl_192h_vs_0h = d8 - control,
  ccl_384h_vs_0h = d16 - control,

  # consecutive time point comparison
  consec_2h_vs_0h = h2 - control,
  consec_8h_vs_2h = h8 - h2,
  consec_24h_vs_8h = d1 - h8,
  consec_48h_vs_24h = d2 - d1,
  consec_96h_vs_48h = d4 - d2,
  consec_144h_vs_96h = d6 - d4,
  consec_192h_vs_144h = d8 - d6,
  consec_384h_vs_192h = d16 - d8,
  levels = design
)

limma_result <- run_limma(expr, design, contrasts) %>%
  assign_deg()

deg_df <- limma_result %>%
  mutate(contrast = fct_inorder(contrast)) %>%
  mutate(contrast_reference = case_when(
    str_detect(contrast, "ccl") ~ "ccl4", 
    str_detect(contrast, "consec") ~ "consec"
  ))

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

### Volcano plots
Volcano plots visualizing the effect of CCl4 on gene expression.
```{r volcano-plots}
df <- readRDS(here(output_path, "limma_result.rds"))

df %>%
  filter(contrast_reference == "ccl4") %>%
  plot_volcano() +
  my_theme(grid = "y", fsize = fz)
```

## Time series clustering
Gene expression trajectories are clustered using the [STEM](http://www.cs.cmu.edu/~jernst/stem/) software. The cluster algorithm is described [here](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-7-191).

### Prepare input
```{r prepare-stem-input}
# prepare input for stem analysis
df <- readRDS(here(output_path, "limma_result.rds"))

stem_inputs <- df %>%
  mutate(class = str_c("Hour ", parse_number(as.character(contrast)))) %>%
  mutate(class = factor(class, levels = unique(.$class))) %>%
  select(gene, class, logFC, contrast_reference)

stem_inputs %>%
  filter(contrast_reference == "ccl4") %>%
  select(-contrast_reference) %>%
  pivot_wider(names_from = class, values_from = logFC) %>%
  write_delim(here(output_path, "stem/input/ccl4.txt"), delim = "\t")
```

### Run STEM
STEM is implemented in Java. The .jar file is called from R. Only significant time series clusters are visualized.
```{r run-stem}
# execute stem
stem_res <- run_stem(here(output_path, "stem"), clear_output = T)

saveRDS(stem_res, here(output_path, "stem_result.rds"))

stem_res %>%
  filter(p <= 0.05) %>%
  filter(key == "ccl4") %>%
  distinct() %>%
  plot_stem_profiles(model_profile = F, ncol = 2) +
  labs(x = "Time in Hours", y = "logFC") +
  my_theme(grid = "y", fsize = fz)
```

### Cluster characterization
STEM clusters are characterized by 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 cluster-characterization}
stem_res = readRDS(here(output_path, "stem_result.rds"))

signatures = stem_res %>%
  filter(p <= 0.05) %>%
  distinct(profile, gene, p_profile = p)

genesets = load_genesets() %>%
  filter(confidence %in% c(NA,"A", "B", "C")) 

ora_res = signatures %>%
  nest(sig = c(-profile)) %>%
  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, "stem_characterization.rds"))
```