Chapter 7: Preliminary Statistical Tests & Spreadsheets: t-Tests, FDR & Volcano Screening

1. Preliminary Statistical Testing in Gene Expression

Before executing generalized linear models (GLMs), biomedical researchers must master the statistical mechanics of hypothesis testing across thousands of genomic features simultaneously. For any individual gene \(i\), we test the null hypothesis (\(H_0\)) that mean expression in the Control group equals mean expression in the Treated group:

$$\mu_{\text{Control}, i} = \mu_{\text{Treated}, i}$$

Parametric vs Non-Parametric Tests

• Student's \(t\)-test: Assumes normality and equal variance (\(\sigma_1^2 = \sigma_2^2\)). Rarely satisfied in small-sample biological data.
• Welch's \(t\)-test: Does not assume equal variance. Much more robust for biological replicates with heteroskedasticity.
• Mann-Whitney \(U\) test (Wilcoxon rank-sum): Non-parametric rank-based test. Free from normality assumptions, but suffers from low statistical power when sample sizes are small (\(n < 8\)).

2. The Multiplicity Problem & Benjamini-Hochberg FDR

If you perform 20,000 independent statistical tests at a standard significance threshold of \(\alpha = 0.05\), purely random chance will generate:

$$20,000 \times 0.05 = 1,000 \text{ false positives!}$$

To control false discoveries without the draconian conservatism of Bonferroni correction (\(\alpha / N\)), the Benjamini-Hochberg (BH) False Discovery Rate (FDR) algorithm ranks all unadjusted \(p\)-values in ascending order \(p_{(1)} \le p_{(2)} \le \dots \le p_{(m)}\) and computes the adjusted \(q\)-value:

$$q_{(i)} = p_{(i)} \times \frac{m}{i}$$

3. Building a Complete Differential Expression Spreadsheet in R

# Vectorized Differential Screening Spreadsheet in R
set.seed(101)
n_genes <- 1000
sample_names <- paste0("Sample_", 1:12)
group <- factor(rep(c("Control", "Treated"), each = 6))

# Normalized log2 expression matrix
log2_expr <- matrix(rnorm(n_genes * 12, mean = 8, sd = 1.5), nrow = n_genes)
rownames(log2_expr) <- paste0("GENE_", 1:n_genes)
colnames(log2_expr) <- sample_names

# Introduce artificial differential expression in first 50 genes
log2_expr[1:50, group == "Treated"] <- log2_expr[1:50, group == "Treated"] + 2.0

# 1. Compute Mean, Log2 Fold Change, and Welch's t-test per gene
ctrl_idx <- which(group == "Control")
trt_idx  <- which(group == "Treated")

mean_ctrl <- rowMeans(log2_expr[, ctrl_idx])
mean_trt  <- rowMeans(log2_expr[, trt_idx])
log2_fc   <- mean_trt - mean_ctrl

# Apply vectorized t-test
p_values <- apply(log2_expr, 1, function(row) {
    t.test(row[trt_idx], row[ctrl_idx], var.equal = FALSE)$p.value
})

# 2. Benjamini-Hochberg Multiple Testing Correction
adj_p_values <- p.adjust(p_values, method = "BH")

# 3. Assemble Clinical Spreadsheet Dataframe
deg_table <- data.frame(
    Gene_Symbol = rownames(log2_expr),
    Mean_Control = round(mean_ctrl, 3),
    Mean_Treated = round(mean_trt, 3),
    Log2FoldChange = round(log2_fc, 3),
    P_Value = p_values,
    Adj_P_Value = adj_p_values,
    Significance = case_when(
        adj_p_values < 0.05 & log2_fc >= 1.0  ~ "Upregulated",
        adj_p_values < 0.05 & log2_fc <= -1.0 ~ "Downregulated",
        TRUE ~ "Not_Significant"
    )
)

# Sort by adjusted p-value
deg_table <- deg_table[order(deg_table$Adj_P_Value), ]

# Export to CSV for Excel review
write.csv(deg_table, "Preliminary_Differential_Expression_Spreadsheet.csv", row.names = FALSE)
cat("Summary of preliminary screening:\n")
print(table(deg_table$Significance))