# Module 5: differential expression with DESeq2.  Video: T7.
# Working directory: the "workshop" folder.
suppressPackageStartupMessages(library(DESeq2))
dir.create("results", showWarnings = FALSE)

meta <- read.delim("data/metadata.txt", row.names = 1)
# Why counts and not TPM? DESeq2 models counts with a negative binomial and does its own normalisation.
count <- read.delim("data/count_gene.txt", row.names = 1, check.names = FALSE)[, rownames(meta)]
count <- round(as.matrix(count))     # kallisto gives estimated (non-integer) counts

coldata <- data.frame(row.names = rownames(meta), conds = factor(meta$condition))
dds <- DESeqDataSetFromMatrix(countData = count, colData = coldata, design = ~conds)
dds <- DESeq(dds)                     # size factors, dispersions, GLM fit, Wald test

# ---- A. MI vs SHAM at day 1, and at day 3 ------------------------------------------------------
res1 <- results(dds, contrast = c("conds", "MI_1D", "SHAM_1D"))   # first = numerator, second = reference
res3 <- results(dds, contrast = c("conds", "MI_3D", "SHAM_3D"))
summary(res1, alpha = 0.05)
write.table(as.data.frame(res1), "results/deseq_1D.txt", sep = "\t", quote = FALSE, col.names = NA)
write.table(as.data.frame(res3), "results/deseq_3D.txt", sep = "\t", quote = FALSE, col.names = NA)

r1 <- as.data.frame(res1)
up <- sum(r1$padj < 0.05 & r1$log2FoldChange > 0, na.rm = TRUE)
down <- sum(r1$padj < 0.05 & r1$log2FoldChange < 0, na.rm = TRUE)
c(up = up, down = down)
head(r1[order(r1$padj), ], 10)

# ---- B. one gene, by hand: the Wald test -----------------------------------------------------------
g <- "Serpine1"
plotCounts(dds, gene = g, intgroup = "conds")
r1[g, ]
z <- r1[g, "log2FoldChange"] / r1[g, "lfcSE"]
p <- 2 * pnorm(-abs(z))
c(z = z, stat_from_DESeq2 = r1[g, "stat"], p = p, p_from_DESeq2 = r1[g, "pvalue"])
# the fold change is (almost exactly) the ratio of the group means of the normalised counts:
nc <- counts(dds, normalized = TRUE)
log2(mean(nc[g, meta$condition == "MI_1D"]) / mean(nc[g, meta$condition == "SHAM_1D"]))
# EXERCISE B1: do the same for "Nppa" and "Postn". Which is significant, and why?

# ---- C. plots --------------------------------------------------------------------------------------
plotMA(res1, ylim = c(-8, 8), alpha = 0.05)
with(r1, plot(log2FoldChange, -log10(padj), pch = 16, cex = 0.4,
              col = ifelse(!is.na(padj) & padj < 0.05, "#c2185b", "grey60"),
              main = "volcano: MI vs SHAM, day 1"))
abline(h = -log10(0.05), lty = 2)

# ---- D. shrunken fold changes, for ranking and plotting ---------------------------------------------
# Low-count genes can have huge, noisy fold changes. lfcShrink pulls them towards 0.
# ("apeglm" or "ashr" are better choices if installed; "normal" needs nothing extra.)
res1s <- lfcShrink(dds, contrast = c("conds", "MI_1D", "SHAM_1D"), type = "normal")
plotMA(res1s, ylim = c(-8, 8), alpha = 0.05, main = "shrunken")

# ---- E. questions ------------------------------------------------------------------------------------
# Q1  At padj < 0.05, how many genes go up and down at day 1? At day 3?
# Q2  What is the most significant gene at each day? Look it up: does it make sense in a heart after MI?
# Q3  Why do some genes have NA in padj? (Hint: summary(res1) and metadata(res1)$filterThreshold)
summary(res3, alpha = 0.05)