Batch in the Design

The first remedy costs one term. Adding the batch to the design formula, ~ type + condition, asks the negative binomial model to estimate a coefficient for the library type and to report the condition effect with the batch held fixed. The counts are never touched. The model that tests the biology is the same model that carries the nuisance, which is why this is the default answer whenever the batch is recorded and the design is not fully confounded.

The idea is the linear model one from introductory statistics, the same one the Practical Bioinformatics book teaches for two groups:

expression = b0 + b1 * condition + b2 * batch + error

The coefficient b1 is the treatment effect controlling for batch, which is the number the experiment exists to estimate. DESeq2 fits that model, with a negative binomial error, one gene at a time.

Two designs, one matrix

Both designs run on the same filtered matrix. The first ignores the batch. The second carries it.

library(DESeq2)
library(pasilla)

# Counts and annotation, loaded the way the DESeq2 vignette does.
pas_anno <- system.file("extdata", "pasilla_sample_annotation.csv",
                        package = "pasilla", mustWork = TRUE)
coldata <- read.csv(pas_anno, row.names = 1)
coldata <- coldata[, c("condition", "type")]
rownames(coldata) <- sub("fb", "", rownames(coldata))
coldata$condition <- relevel(factor(coldata$condition), ref = "untreated")
coldata$type <- factor(coldata$type)

pas_cts <- system.file("extdata", "pasilla_gene_counts.tsv",
                       package = "pasilla", mustWork = TRUE)
cts <- as.matrix(read.csv(pas_cts, sep = "\t", row.names = "gene_id"))
cts <- cts[, rownames(coldata)]

# Keep genes with at least 10 reads total.
keep <- rowSums(cts) >= 10
cts <- cts[keep, ]

# Design without the batch term.
dds_simple <- DESeqDataSetFromMatrix(cts, coldata, design = ~ condition)
dds_simple <- DESeq(dds_simple, quiet = TRUE)

res_simple <- results(dds_simple,
                      contrast = c("condition", "treated", "untreated"),
                      alpha = 0.05)
summary(res_simple)

out of 9921 with nonzero total read count
adjusted p-value < 0.05
LFC > 0 (up)       : 407, 4.1%
LFC < 0 (down)     : 431, 4.3%
outliers [1]       : 1, 0.01%
low counts [2]     : 1347, 14%
(mean count < 5)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

# Design with the batch term, batch first.
dds_batch <- DESeqDataSetFromMatrix(cts, coldata,
                                    design = ~ type + condition)
dds_batch <- DESeq(dds_batch, quiet = TRUE)

res_batch <- results(dds_batch,
                     contrast = c("condition", "treated", "untreated"),
                     alpha = 0.05)
summary(res_batch)

out of 9921 with nonzero total read count
adjusted p-value < 0.05
LFC > 0 (up)       : 493, 5%
LFC < 0 (down)     : 593, 6%
outliers [1]       : 0, 0%
low counts [2]     : 1539, 16%
(mean count < 6)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

The contrast argument names the comparison explicitly, and with two variables in the design that is not a formality. results() without it returns the last variable in the formula, which here would be the condition effect but with the batch first is the library type. Naming the contrast keeps the table about the biology.

What the extra term did

The two summaries differ in three places: the significant gene counts, the outlier count, and nothing else that matters. The next blocks count both lists and compare them.

sig_simple <- rownames(res_simple)[!is.na(res_simple$padj) &
                                   res_simple$padj < 0.05]
sig_batch <- rownames(res_batch)[!is.na(res_batch$padj) &
                                 res_batch$padj < 0.05]

length(sig_simple)
length(sig_batch)
length(intersect(sig_simple, sig_batch))

[1] 838
[1] 1086
[1] 826

# Log fold changes of the shared genes, one design against the other.
shared <- intersect(rownames(res_simple), rownames(res_batch))
lfc_simple <- res_simple[shared, "log2FoldChange"]
lfc_batch <- res_batch[shared, "log2FoldChange"]

median(abs(lfc_simple - lfc_batch), na.rm = TRUE)
cor(lfc_simple, lfc_batch, use = "complete.obs")

[1] 0.03129861
[1] 0.9759554

The direction of the first result surprises people the first time. Modeling the batch returned more significant genes, not fewer. The batch variance left the residual and the standard errors shrank, so 1086 genes cleared the threshold where 838 did before, with 826 of the original 838 still on the list. The log fold changes themselves barely moved, a median shift of three hundredths of a log2 unit across the shared genes, because the coefficient the design estimates for the condition is barely changed by carrying one more nuisance term. What changed is the uncertainty around it.

The outlier count moved too. The simple design flagged one gene as a Cooks distance outlier, meaning one sample sat far from the model for that gene. With the library type in the design, that distance is explained and the flag disappears. A batch term regularly cleans up exactly the kind of sample a first analysis suspects of being bad data.

The same idea in a pipeline

The nf-core/differentialabundance pipeline takes this route through the blocking column of its contrasts file, where a contrast declares the variables to control for. The site runs that pipeline for real in the differential abundance guide, and its contrasts file carries the same logic this page wrote into a design formula.

One limit closes the page. This route worked because every cell of the design table was occupied, so the model could tell the two effects apart. When the batch is fully confounded with the condition, the design term cannot be estimated, the fit fails or silently absorbs the biology, and the only honest answer is that the experiment cannot support the comparison. The next page covers the route that rewrites the data instead, which is sometimes the only option left, and measures what it costs.