# Module 8: designs with more than one factor.  Video: T10.
# Working directory: the "workshop" folder.
# Note: the 2020 version of this exercise built the model matrix but then still fitted ~conds.
# Here the design really is ~ day + MI.
suppressPackageStartupMessages(library(DESeq2))
dir.create("results", showWarnings = FALSE)

meta <- read.delim("data/metadata_multi.txt", row.names = 1)
meta$MI <- factor(meta$MI)      # 0 = sham, 1 = infarct  (as factors, not numbers!)
meta$day <- factor(meta$day)    # 1 or 3
meta
count <- round(as.matrix(read.delim("data/count_gene.txt", row.names = 1, check.names = FALSE)[, rownames(meta)]))

# ---- A. look at the model matrix ------------------------------------------------------------------------
model.matrix(~ day + MI, meta)
# EXERCISE A1: what does each column mean? Which samples have a 1 in "MI1"?

# ---- B. fit ~ day + MI ------------------------------------------------------------------------------------
dds <- DESeqDataSetFromMatrix(count, meta, design = ~ day + MI)
dds <- DESeq(dds)
resultsNames(dds)
res_mi <- results(dds, name = "MI_1_vs_0")          # the MI effect, adjusted for day
summary(res_mi, alpha = 0.05)
write.table(as.data.frame(res_mi), "results/deseq_multifactor.txt", sep = "\t", quote = FALSE, col.names = NA)
res_day <- results(dds, name = "day_3_vs_1")        # the day effect, adjusted for MI
summary(res_day, alpha = 0.05)

# ---- C. compare with the single-factor results from module 5 ------------------------------------------------
r1 <- read.delim("results/deseq_1D.txt", row.names = 1)
r3 <- read.delim("results/deseq_3D.txt", row.names = 1)
rm_ <- as.data.frame(res_mi)
sig <- function(r) rownames(r)[!is.na(r$padj) & r$padj < 0.05]
c(day1 = length(sig(r1)), day3 = length(sig(r3)), both_days = length(intersect(sig(r1), sig(r3))),
  multifactor = length(sig(rm_)))
plot(r1$log2FoldChange, rm_[rownames(r1), "log2FoldChange"], pch = 16, cex = 0.3,
     xlab = "LFC, day 1 only", ylab = "LFC, ~ day + MI")
abline(0, 1, col = "red")

# ---- D. an interaction: does MI act differently at day 3? ------------------------------------------------------
ddi <- DESeqDataSetFromMatrix(count, meta, design = ~ day + MI + day:MI)
ddi <- DESeq(ddi, quiet = TRUE)
resultsNames(ddi)
res_int <- results(ddi, name = "day3.MI1")
summary(res_int, alpha = 0.05)
# EXERCISE D1: in words, what does a significant day3.MI1 coefficient mean for a gene?
# EXERCISE D2: with only 2 replicates per group, would you trust these interaction results? Why (not)?

# ---- E. functional analysis on the multifactor result ------------------------------------------------------------
# EXERCISE E1: rerun module 7, section C (fgsea) on res_mi. Which pathways are shared with the day 1 result?