Bootcamp / Session 4 / 4.1
Session 4.1 · Manuel Garber, UMass Chan
Why differential expression uses models.
DESeq2, edgeR, and limma do not run a t-test on every gene. Four short notebooks show why: how expression measurements are distributed, how a t-test is already a linear model, what changes when the data are counts, and why a gene measured in three replicates needs to borrow information from all the others.
Recording
Manuel's lecture: a brief introduction to RNA-seq. Parts about unpublished data were cut from the recording and the slides.
The session in one sentence
The shape of one gene’s histogram is a weak reason to use a count model. The strong reasons are the mean–variance relationship, sharing information across genes, and the design matrix.
The four notebooks
Run each with uvx marimo edit --sandbox <file> from the notebook folder, Documents/bootcamp/session4/bootcamp-bulk-rnaseq-main. In a new terminal, go back to it first. Every code cell has a short explanation above it.
How are expression measurements distributed? 01_distributions.py
GAPDH and ACTB on microarrays look roughly normal. The same genes in RNA-seq come out as raw counts clumped by sequencing depth, until you normalize for depth. Stable genes at lower expression show where a normal distribution stops fitting the counts, below about 10 reads. The notebook ends with mean against variance within replicates.
The t-test is a linear model 02_ttest_is_lm.py
A two-sample t-test and the model y ~ group give the same t statistic and p-value. Writing it as a model pays off as soon as the experiment has structure: adding a batch term separates the treatment effect from batch, and in the simulation the batch model detects an effect the plain t-test almost always misses.
Count data and GLMs 03_counts_and_glm.py
The same simulated biology measured three ways: qPCR, microarray, and RNA-seq. Biological variation between samples plus Poisson sampling of reads gives the negative binomial distribution. A negative binomial GLM is then compared with ordinary regression on log counts, gene by gene and across many genes for false positives and power.
Estimating each gene’s variance 04_variance_estimation.py
With three replicates, a gene’s variance estimate is so noisy that its 95% range spans about 150-fold. Replicate-against-replicate scatters show a constant band for arrays and a funnel for RNA-seq, where Poisson noise dominates at low counts. Mean expression explains most of the gene-to-gene variation in RNA-seq variance, so a curve fitted across all genes, and shrinkage toward it (the idea behind limma’s eBayes and DESeq2’s dispersion prior), gives better estimates, checked against held-out donors.
The RNA-seq panels in notebook 1 read count tables that are not public, so Manuel shows them live. Everything else runs on your computer. Notebook 4 needs BULK_RNASEQ_DIR set, as in 4.0.
Next
In 4.2 an AI agent in VS Code sets up an RNA-seq run on Foundry, and you check it and launch it.
Further reading
- Law, Chen, Shi, Smyth (2014). voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biology 15:R29. Why limma with voom weights behaves much like a negative binomial GLM.
- Love, Huber, Anders (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15:550.
Materials
After the session
Two minutes of feedback, so the next session is pitched better. It also records your attendance, so fill it in whether you were in the room, on Zoom, or catching up on the recording.