mouse-acute-apap.Rmd

---
title: "Acute APAP 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 APAP induced acute liver damage. The transcriptomic profiles were measured at different time points ranging from 1 hour 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(janitor)

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

## Global variables for this analysis
Definition of global variables that are used throughout the entire analysis.
```{r analysis-specific-params, cache=FALSE}
# i/o
data_path <- "data/mouse-acute-apap"
output_path <- "output/mouse-acute-apap"

# 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)
colnames(expr) <- str_c("sample",
  str_remove(colnames(expr), "_\\(Mouse430_2\\).CEL"),
  sep = "_"
)

# 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", "rep"), remove = F) %>%
  dplyr::select(-tmp) %>%
  mutate(
    rep = parse_number(rep),
    time = parse_number(key)
  ) %>%
  mutate(time = case_when(
    str_detect(key, "d") ~ time * 24,
    str_detect(key, "h") ~ time,
    str_detect(key, "control") ~ 0
  )) %>%
  mutate(time = ordered(time)) %>%
  mutate(group = case_when(
    str_detect(key, "d") ~ str_c("d", parse_number(key)),
    str_detect(key, "h") ~ str_c("h", parse_number(key)),
    str_detect(key, "control") ~ key
  )) %>%
  mutate(group = factor(group, levels = c(
    "control", "h1", "h6", "h12", "d1", "d2",
    "d4", "d6", "d8", "d16"
  )))

# 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 APAP 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 apap treatment
  apap_1h_vs_0h = h1 - control,
  apap_6h_vs_0h = h6 - control,
  apap_12h_vs_0h = h12 - control,
  apap_24h_vs_0h = d1 - control,
  apap_48h_vs_0h = d2 - control,
  apap_96h_vs_0h = d4 - control,
  apap_144h_vs_0h = d6 - control,
  apap_192h_vs_0h = d8 - control,
  apap_384h_vs_0h = d16 - control,

  # consecutive time point comparison
  consec_1h_vs_0h = h1 - control,
  consec_6h_vs_1h = h6 - h1,
  consec_12h_vs_6h = h12 - h6,
  consec_24h_vs_12h = d1 - h12,
  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 = factor(contrast, levels = c(
    "apap_1h_vs_0h", "apap_6h_vs_0h", "apap_12h_vs_0h", "apap_24h_vs_0h",
    "apap_48h_vs_0h", "apap_96h_vs_0h", "apap_144h_vs_0h", "apap_192h_vs_0h",
    "apap_384h_vs_0h",
    "consec_1h_vs_0h", "consec_6h_vs_1h", "consec_12h_vs_6h", 
    "consec_24h_vs_12h", "consec_48h_vs_24h", "consec_96h_vs_48h", 
    "consec_144h_vs_96h", "consec_192h_vs_144h", "consec_384h_vs_192h"
  ))) %>%
  mutate(contrast_reference = case_when(
    str_detect(contrast, "apap") ~ "apap",
    str_detect(contrast, "consec") ~ "consec"
  ))

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

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

df %>%
  filter(contrast_reference == "apap") %>%
  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 == "apap") %>%
  select(-contrast_reference) %>%
  pivot_wider(names_from = class, values_from = logFC) %>% 
  write_delim(here(output_path, "stem/input/apap.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 == "apap") %>%
  distinct() %>%
  plot_stem_profiles(model_profile = F, ncol = 2) +
  labs(x = "Time in Hours", y = "logFC") +
  my_theme(grid = "y", fsize = fz)
```