The statistics underneath omics

A ground-up look at the methods every omics pipeline runs on your behalf — count distributions, dispersion, normalisation, multiple testing, and what a p-value means when you computed twenty thousand of them.

NoteThis walkthrough is being written

The outline below is the plan. If there is a piece of statistical machinery you have run a hundred times without being sure what it does, say which one and it will be prioritised: aaron.mitchell.dick@duke.edu.

Every omics pipeline makes statistical decisions for you. DESeq() is one function call, and inside it are a normalisation, a dispersion model, a shrinkage step, a hypothesis test, an outlier rule, and a filtering rule — six choices, any of which changes your gene list.

This series opens that box. It is not a statistics course and it does not assume one; it assumes you have run these pipelines and want to know what they did. The goal is that you can defend your analysis in a lab meeting and recognise when a default is wrong for your experiment.

Every section pairs the mathematics with runnable R, using the same six-sample mouse dataset from the bulk RNA-Seq series, so you can watch each choice change the answer.

Planned outline

Part 1 — Why count data breaks ordinary statistics

• Counts are not measurements: discreteness, non-negativity, and mean–variance coupling
• The Poisson model, and the exact experiment that shows it fails on real RNA-Seq
• Overdispersion: what biological variability does to the variance
• The negative binomial as a Poisson–gamma mixture, built up rather than asserted
• Why t.test() on log-CPM is not as wrong as people say, and when it is

Part 2 — Normalisation

• What normalisation is actually correcting for, and what it cannot correct for
• Why total-count scaling (CPM) fails when a few genes dominate the library
• Median-of-ratios (DESeq2) and TMM (edgeR), derived and compared side by side
• The composition-bias assumption both methods rest on, and how to check it
• Length correction: TPM, FPKM, and why tximport passes an offset to DESeq2 instead of using TPM
• Spike-ins, and the case where nothing else will save you

Part 3 — Dispersion estimation and shrinkage

• Why you cannot estimate variance from three replicates
• Information sharing across genes: the central trick that makes small-n omics work
• The dispersion-by-mean trend, fitted and interpreted
• Empirical Bayes shrinkage towards the trend, and dispersion outliers
• Effect size shrinkage: apeglm, ashr, normal, and what each assumes
• Reading a dispersion plot and knowing when it is telling you something is wrong

Part 4 — Testing twenty thousand hypotheses

• The Wald test and the likelihood ratio test, and when to choose which
• Family-wise error rate versus false discovery rate: two different questions
• Bonferroni, Benjamini–Hochberg, and Storey’s q-value
• Independent filtering: how DESeq2 buys power by discarding low-count genes, and why that is legitimate
• What an adjusted p-value of 0.05 licenses you to say, and what it does not
• The garden of forking paths in an omics pipeline, and pre-registration as a partial answer

Part 5 — Experimental design and multi-factor models

• Power analysis for RNA-Seq: replicates beat depth, almost always
• Biological versus technical replication, and the pseudoreplication trap
• Reading and writing model formulas: ~ batch + condition and what the terms mean
• Interaction terms — asking whether an effect differs between groups
• Paired and time-course designs
• Batch effects: modelling them, removing them, and the difference between the two
• Confounding you cannot fix after the fact, and the design meeting that prevents it

Who this is for

Anyone who has run DESeq2, limma, edgeR, or an equivalent, and has had the uncomfortable feeling of not being able to explain the output. Familiarity with means, variances, and p-values is assumed. Nothing beyond that is — the mathematics is built up from there, with every step written out.