# Module 6: multiple testing and the false discovery rate.  Video: T8.
# Working directory: the "workshop" folder. Needs results/deseq_1D.txt from module 5.
suppressPackageStartupMessages(library(DESeq2))
set.seed(7)

# ---- A. what happens when nothing changes? A null experiment -------------------------------------
# Take the 4 SHAM samples and split them into two fake groups, A and B, so that each group has one
# day-1 and one day-3 sample. There is no real difference between A and B.
meta <- read.delim("data/metadata.txt", row.names = 1)
count <- round(as.matrix(read.delim("data/count_gene.txt", row.names = 1, check.names = FALSE)[, rownames(meta)]))
sham <- rownames(meta)[grepl("^SHAM", meta$condition)]
data.frame(sample = sham, condition = meta[sham, "condition"])
fake <- factor(c("A", "B", "B", "A"))           # A = 02 (day 1) + 07 (day 3), B = 03 (day 1) + 06 (day 3)
ddn <- DESeqDataSetFromMatrix(count[, sham], data.frame(row.names = sham, grp = fake), ~grp)
ddn <- DESeq(ddn, quiet = TRUE)
rn <- results(ddn)
hist(rn$pvalue, breaks = 20, col = "grey70", main = "p-values: fake groups", xlab = "p-value")
c(raw_p_below_0.05 = sum(rn$pvalue < 0.05, na.rm = TRUE), padj_below_0.05 = sum(rn$padj < 0.05, na.rm = TRUE))
# EXERCISE A1: now set fake <- factor(c("A", "A", "B", "B")), so that A = day 1 and B = day 3.
#              Rerun. What happens to the histogram and to padj, and why? (Hint: is there still
#              "no difference" between the groups?)

# ---- B. the real comparison --------------------------------------------------------------------------
r1 <- read.delim("results/deseq_1D.txt", row.names = 1)
hist(r1$pvalue, breaks = 20, col = "#3f7cac", main = "p-values: MI vs SHAM, day 1", xlab = "p-value")
m <- sum(!is.na(r1$pvalue))
c(tests = m, expected_false_if_nothing_changed = 0.05 * m, raw_p_below_0.05 = sum(r1$pvalue < 0.05, na.rm = TRUE))

# ---- C. Benjamini-Hochberg by hand -----------------------------------------------------------------------
p <- c(0.001, 0.008, 0.012, 0.035, 0.041, 0.06, 0.2, 0.45, 0.62, 0.9)    # the toy example from the video
bh_by_hand <- function(p) {
  n <- length(p); o <- order(p)
  adj <- p[o] * n / seq_len(n)                      # step 2: p x tests / rank
  adj <- rev(cummin(rev(adj)))                      # step 3: never bigger than the value below
  out <- numeric(n); out[o] <- pmin(adj, 1); out
}
cbind(p, by_hand = round(bh_by_hand(p), 3), p.adjust = round(p.adjust(p, "BH"), 3), bonferroni = p.adjust(p, "bonferroni"))

# ---- D. why p.adjust() on all p-values does not match DESeq2's padj ---------------------------------------
pa_all <- p.adjust(r1$pvalue, "BH")
c(DESeq2_padj_below_0.05 = sum(r1$padj < 0.05, na.rm = TRUE), p.adjust_all_below_0.05 = sum(pa_all < 0.05, na.rm = TRUE))
# DESeq2 first removes genes with a low mean count (independent filtering), then adjusts:
tested <- !is.na(r1$padj)
pa_f <- p.adjust(r1$pvalue[tested], "BH")
all.equal(pa_f, r1$padj[tested])
# EXERCISE D1: count the genes with NA in pvalue, and those with a p-value but NA in padj. What is
#              different about them? (Hint: look at baseMean.)
c(na_pvalue = sum(is.na(r1$pvalue)), na_padj_only = sum(!is.na(r1$pvalue) & is.na(r1$padj)))