# Reading a differential expression table: raw p, BH-adjusted p, a fold-change filter applied after testing, and a test
# against a fold-change threshold (TREAT, McCarthy & Smyth 2009). Student's t-test per gene on log2 expression, 4 vs 4.
# R 4.3.3, base R.
d <- read.csv("log2_expression.csv"); truth <- read.csv("truth.csv")$true_log2fc
ctl <- as.matrix(d[, grep("^control_", names(d))]); trt <- as.matrix(d[, grep("^treated_", names(d))])
n1 <- ncol(ctl); n2 <- ncol(trt); df <- n1 + n2 - 2
lfc <- rowMeans(trt) - rowMeans(ctl)
s2 <- ((n1 - 1) * apply(ctl, 1, var) + (n2 - 1) * apply(trt, 1, var)) / df; se <- sqrt(s2 * (1 / n1 + 1 / n2))
p <- 2 * pt(abs(lfc) / se, df, lower.tail = FALSE); q <- p.adjust(p, "BH")
tau <- 1                                              # threshold: a 2-fold change
p_treat <- pt((abs(lfc) - tau) / se, df, lower.tail = FALSE) + pt((abs(lfc) + tau) / se, df, lower.tail = FALSE)
q_treat <- p.adjust(p_treat, "BH")
show <- function(lab, called) cat(sprintf("%s: %d genes; %d truly unchanged; %d with a true change below 2-fold (including unchanged)\n",
                                          lab, sum(called), sum(called & truth == 0), sum(called & abs(truth) < tau)))
show("Raw p < 0.05", p < 0.05)
show("BH-adjusted p (padj) < 0.05", q < 0.05)
show("padj < 0.05 and |log2 fold change| > 1, filtered after testing", q < 0.05 & abs(lfc) > tau)
show("Test against |log2 fold change| > 1 (TREAT), adjusted p < 0.05", q_treat < 0.05)
cat(sprintf("Truly changed genes: %d, of which %d by more than 2-fold\n", sum(truth != 0), sum(abs(truth) > tau)))