Chapter 6 — DESeq2: Bulk Fundamentals + Pseudobulk DE
Where this chapter sits. Companion to Lecture 06 and Tutorial 06 — DESeq2: Bulk Primer + Pseudobulk DE. Prerequisites: Appendix B — Statistical Foundations, especially §§ 1–3 (NB / Poisson / over-dispersion), § 4 (mean-variance relationship), §§ 5–7 (testing), § 8 (why per-cell DE fails).
Part A — DESeq2 fundamentals on bulk RNA-seq
6.1 Why we learn DESeq2 on bulk first
Pseudobulk differential expression for scRNA-seq (Part B of this chapter) is exactly bulk DESeq2 applied to per-cell-type aggregated counts. So we learn the bulk model in isolation here, on the canonical bulk dataset (airway from Bioconductor — 8 samples of dexamethasone-treated airway smooth muscle cells), and then in Part B we apply it almost mechanically to the single-cell aggregation.
This is the most heavily used statistical model in functional genomics1. Its closest relatives — edgeR2 and limma-voom3,4 — share the same negative-binomial-or-linear-model framing and routinely agree; benchmarks of bulk DE methods find that with adequate replication these established tools all perform well, and that adding biological replicates buys far more than adding sequencing depth5. Earning a working understanding of this model is worth several lectures.
Differential expression (DE) asks: for cells of the same type, which genes are expressed at systematically different levels between two conditions? The question is answered at the level of estimated means — is the expected count of gene \(i\) higher in the treated group than the control group, by more than chance alone? The answer depends on having enough independent observations (biological replicates) to separate true between-condition signal from within-condition biological noise. Adding more cells or more sequencing depth does not substitute for adding more independent samples, because the limiting noise is biological variability across subjects, not counting noise within a single library.
Table 1 compares the three standard tools so you can navigate the literature and justify your choice.
| Method | Statistical model | Dispersion handling | LFC shrinkage | Pseudobulk support | Workshop use |
|---|---|---|---|---|---|
| DESeq21 | Negative-binomial GLM | Per-gene shrinkage toward a parametric trend | apeglm / ashr / normal priors |
Yes — standard approach | Primary tool (this chapter) |
| edgeR2 | Negative-binomial GLM (quasi-likelihood variant available) | Empirical Bayes across genes (common, trended, or tagwise) | Not built in; use limma::treat for threshold tests |
Yes | Alternative; see §6.21 |
| limma-voom3,4 | Linear model on variance-stabilized counts | Mean-variance weights via voom; eBayes shrinkage |
Moderated t-statistics; treat for threshold |
Yes | Alternative for large cohorts |
6.2 The negative-binomial GLM
DESeq2 fits, for each gene, a generalized linear model:
\[\log(\mu_{ij}) = X_j^\top \beta_i + \log(s_j)\]
where:
- \(\mu_{ij}\) is the expected count of gene \(i\) in sample \(j\)
- \(X_j\) is the design matrix row for sample \(j\) (encodes condition, batch, etc.)
- \(\beta_i\) is the per-gene coefficient vector (what we want)
- \(s_j\) is the size factor for sample \(j\) — depth normalization, fitted separately
The likelihood is negative binomial:
\[Y_{ij} \sim \mathrm{NB}(\mu_{ij}, \alpha_i)\]
with per-gene dispersion \(\alpha_i\) that controls how much variance exceeds the Poisson baseline. The variance is \(\mu_{ij}(1 + \alpha_i \mu_{ij})\).
The choice of NB rather than Poisson is non-negotiable for RNA-seq counts. Poisson assumes variance = mean; real RNA-seq has variance > mean (over-dispersion), and Poisson testing is anti-conservative — it produces falsely-significant p-values. The NB adds a per-gene parameter to absorb the extra variance; results are honest. See Appendix B § 1–2 for the underlying distributions.
When you test 20,000 genes simultaneously, the raw per-gene p-value is a misleading quantity: at \(\alpha = 0.05\), roughly 1,000 genes would appear significant by chance alone. The Benjamini–Hochberg false discovery rate (FDR) procedure6 ranks the raw p-values from smallest to largest, compares each to a scaled threshold \((\text{rank}/\text{number of tests}) \times q\), and adjusts the p-values so that the expected proportion of false positives among all declared significant genes equals \(q\) (typically 5%). The column padj in DESeq2 output is this BH-adjusted p-value; filtering on padj < 0.05 means that you accept a 5% false-discovery rate across your entire significant list, not a 5% per-test error rate. Never filter on the raw pvalue column in a multi-gene test, and never report raw p-values in a publication when adjusted values are available.
6.3 Size factors — depth normalization
Different samples have different sequencing depth. Without correction, a sample with twice as many total reads has twice the count for every gene, regardless of biology. The DESeq2 size factor \(s_j\) corrects this:
\[s_j = \mathrm{median}_i\left( \frac{y_{ij}}{(\prod_k y_{ik})^{1/n}} \right)\]
This is the median of ratios estimator. It computes a per-gene “geometric mean across samples” baseline, then for each sample takes the median of (sample \(j\)’s count / baseline) across genes. The result is robust — a few highly-expressed genes can’t dominate — which is why it performs far better than simple library-size normalization when a small number of genes are strongly upregulated in one condition. If you simply divided by total library size, those few up-regulated genes would make everything else appear artificially down-regulated in the treated samples.
DESeqDataSetFromMatrix(...) and the subsequent DESeq(dds) call do this automatically; you don’t typically interact with size factors directly. But knowing the formula clarifies why DESeq2 is robust to a few genes shifting wildly between samples. To inspect the estimated size factors, use sizeFactors(dds) after fitting; values far from 1.0 (e.g. > 3 or < 0.3) are a sign that samples have very different depths or that the library composition is extreme. In well-controlled RNA-seq experiments, size factors typically range from 0.5 to 2.0.
6.4 Dispersion estimation and shrinkage
Every gene has its own dispersion \(\alpha_i\), but with only a few samples per group, per-gene estimates are extremely noisy. DESeq2 borrows strength across genes by shrinking dispersion estimates toward a fitted trend (mean expression vs dispersion).
The procedure:
- Maximum-likelihood per-gene dispersion estimates
- Fit a parametric trend: log-dispersion as a function of log-mean
- Shrink each gene’s estimate toward the trend, weighted by its uncertainty
The result: high-confidence genes (high counts, many samples) keep their MLE estimates; low-confidence genes (low counts) are pulled to the trend. The dispersions(dds) slot holds the final shrunk values; plotDispEsts(dds) visualizes the trend.
This is a Bayesian-flavored procedure — the trend is the prior, the per-gene MLE is the likelihood, the shrunk value is the posterior. See Appendix B § 7 for the connection.
Reading the dispersion plot. plotDispEsts(dds) shows three things: the raw per-gene MLE estimates (black dots), the fitted trend (red line), and the final shrunken values (blue dots with circles). A healthy plot has: (a) a smooth trend that decreases monotonically as mean expression rises — high-count genes are less variable in relative terms; (b) black dots scattered around the red line without systematic structure (systematic departures suggest batch effects or outlier samples); (c) blue dots that fall on the trend for low-expressed genes but diverge toward the MLE for high-expressed genes where the data are informative. If the trend is flat or U-shaped, something is wrong with the dataset — usually an outlier sample, a confound in the design, or far too few samples for dispersion estimation to work. That is the time to go back and inspect the VST PCA plot (vst(dds) → plotPCA()), which often reveals the outlier or batch structure causing the problem.
DESeqDataSetFromMatrix(design = ~ condition) specifies which sample-level covariates enter the model. Default: no default — this argument is required. The simplest design ~ condition models only the treatment effect. Adding a blocking variable — ~ donor + condition for paired donors, or ~ batch + condition for multi-batch experiments — accounts for that covariate in the model so the condition coefficient reflects the treatment effect within each level of the covariate, substantially increasing power in paired or blocked designs. The design must match the columns of colData exactly; a mismatch produces an error at DESeq() time. Always call relevel() on all factor columns before DESeq() so the sign of every log-fold-change is intentional.
6.5 The design matrix and contrasts
The design = ~ ... argument specifies which sample-level covariates affect expression. The simplest case:
DESeqDataSetFromMatrix(countData = ..., colData = ..., design = ~ condition)Here condition is a factor with levels untreated and treated. The model has two coefficients per gene: an intercept (log expected expression in untreated, the reference level) and a conditiontreated coefficient (log fold-change of treated vs untreated).
For paired designs:
design = ~ donor + conditionThis adds donor as a covariate, soaking up donor-level variation. The condition coefficient now estimates the treatment effect within each donor — much more powerful when donors are paired.
For interactions:
design = ~ donor + group + condition + group:conditionThis tests whether the effect of condition differs between groups.
Always run relevel() to set your reference level explicitly:
dds$condition <- relevel(dds$condition, ref = "untreated")If you skip this, R’s alphabetical default may pick the wrong reference, flipping the sign of every fold-change in your output. Hard to debug after the fact.
For paired designs, adding donor as a fixed effect means that each subject’s measurements are modeled jointly and the condition effect is estimated within each donor, analogous to a paired t-test. The design formula ~ donor + condition does not literally mean “donor causes expression” — it means “account for donor-level mean differences when estimating the condition coefficient.” For real multi-donor single-cell data (such as the ifnb dataset’s eight lupus donors), this is almost always the correct choice over ~ condition alone, because inter-donor variability in baseline expression typically dwarfs the stimulus effect and its omission inflates within-group variance. In a balanced study where every donor appears in every condition, the paired design can reduce the number of significant genes to a fraction compared to an unpaired design — not because the biology changed, but because the unpaired analysis was confusing donor-specific expression differences with condition effects.
For interaction designs — asking whether the condition effect differs across groups — the formula becomes ~ group + condition + group:condition. The group:condition term captures the difference-in-differences: does the treatment response depend on which group a sample belongs to? Interactions are powerful but require substantially more replication because the interaction has fewer degrees of freedom than the main effects.
DESeqDataSetFromMatrix(design = ~ donor + condition) adds a blocking variable to the model. Default: required argument with no default. The simplest valid design ~ condition tests the condition effect ignoring all other structure. Adding ~ donor + condition for paired donors (each donor measured under each condition) or ~ batch + condition for multi-batch experiments is the standard upgrade: the blocking variable soaks up the covariate’s variance contribution so the condition coefficient estimates the treatment effect within each level of the covariate. This is directly analogous to a paired t-test. For the ifnb eight-donor dataset, upgrading from ~ condition to ~ donor + condition substantially increases power and produces far more reproducible gene lists. The blocking variable must be a factor column in colData and must not be collinear with condition (e.g. if every donor appears in only one condition, ~ donor + condition is unidentifiable).
AggregateExpression(seu, group.by = c("donor","stim","celltype")) followed by a manual count filter sets the minimum number of cells required for a pseudobulk column to be included. De-facto default: 10 cells. Below that threshold the pseudobulk profile is dominated by stochastic dropout rather than cell-type biology: one or two cells happen to express or not express any given gene by chance, and the resulting count column looks nothing like the true cell-type profile. Increasing the threshold to 20–50 is appropriate when libraries are shallow or the cell type is rare in some donors. There is no way to compensate for a too-small pseudobulk column with better statistics — the signal simply is not there. The consequence of including under-threshold columns is inflated dispersion estimates and reduced power for every other cell type in the same DESeq2 run.
6.6 The Wald test and p-values
DESeq2’s default per-coefficient test is the Wald test:
\[z = \frac{\hat\beta - 0}{\hat{SE}(\hat\beta)}\]
Under the null (\(\beta = 0\)), \(z\) is approximately standard normal. The two-sided p-value is \(2(1 - \Phi(|z|))\).
For testing whether any level of a multi-level factor matters, use the likelihood ratio test (LRT):
dds <- DESeq(dds, test = "LRT", reduced = ~ donor)This compares the full model (~ donor + condition) to the reduced model (~ donor); the LRT statistic is \(-2 \log(L_{reduced}/L_{full})\) and is \(\chi^2\)-distributed. Useful when condition has 3+ levels.
For both tests, the output table from results(dds) has columns described in Table 2.
results(dds) output table, their meaning, and how to interpret each. padj = NA indicates the gene was excluded by independent filtering or had a Cook’s-distance outlier — this is normal behavior, not a bug (§6.9).
| Column | Meaning | How to use it |
|---|---|---|
baseMean |
Mean of normalized counts across all samples | Low values (< 5) flag unreliable LFC estimates — apply LFC shrinkage |
log2FoldChange |
Effect size on base-2 log scale (unshrunk unless lfcShrink is called) |
Ranking and visualization: always shrink first with lfcShrink |
lfcSE |
Standard error of the log2 fold-change estimate | Use to assess confidence; wide SE → noisy estimate |
stat |
Wald z-statistic (or LRT \(\chi^2\)) | Directly from the GLM; positive = up in numerator condition |
pvalue |
Raw per-gene test p-value | Never filter on this alone — apply BH correction first |
padj |
Benjamini–Hochberg FDR6 | The primary filter: padj < 0.05 is the standard threshold |
6.7 LFC shrinkage — never skip
Genes with very low counts (baseMean < 5) have huge LFC standard errors. Their raw log2FoldChange values are noisy and dominated by sampling luck. Sort the unshrunk results by LFC and your “top hits” will be a list of garbage low-count genes.
The fix: shrinkage:
res_shrunk <- lfcShrink(dds, coef = "condition_treated_vs_untreated",
type = "apeglm")apeglm7 is the recommended estimator. It pulls noisy LFCs toward zero, weighted by their uncertainty. High-confidence genes barely change; low-count noisy genes shrink dramatically.
Coefficient naming. apeglm requires the exact coefficient name, not a contrast vector. The name is always of the form <variable>_<level>_vs_<reference> — e.g. dex_treated_vs_untreated for the airway dataset after renaming, or condition_STIM_vs_CTRL for the ifnb pseudobulk loop. Get the exact string with resultsNames(dds) before calling lfcShrink; a mismatch gives a confusing “coefficient not found” error. If you renamed factor levels (as Tutorial 06 does for airway’s trt→treated/untrt→untreated), the coefficient name reflects the new level names.
Always shrink before:
- Volcano plots
- Ranking by effect size
- Any “top genes” list
The unshrunk LFC is fine for hypothesis testing (the p-value cares about the SE-scaled LFC, not the LFC itself), but for ranking and visualization, shrunk values are far more honest.
The shrinkage is a Bayesian procedure with apeglm’s prior — see Appendix B § 7 for the connection to the dispersion shrinkage in §6.4. It’s the same idea: borrow strength from all genes to stabilize each gene’s estimate.
An important nuance: since approximately 2016, DESeq2 computes pvalue and padj from the unshrunken MLE log-fold-changes — shrinkage does not enter the Wald test. This keeps the FDR well-calibrated. The shrunken LFCs and their SEs are exclusively for ranking, volcano/MA plots, and effect-size reporting after testing. Older tutorials occasionally show shrinkage before the results() call — this was the workflow in early DESeq2 versions and is now incorrect.
lfcShrink(dds, coef = ..., type = "apeglm") controls the shrinkage estimator. lfcShrink is a separate function you call after DESeq() — it is not part of model fitting. Its current default is type = "apeglm" (the legacy default was "normal", before DESeq2 1.18 in 2017), and apeglm is the recommended estimator7. "apeglm" uses an adaptive Cauchy (heavy-tailed) prior that shrinks toward zero very aggressively for noisy low-count genes but barely touches high-confidence genes. "ashr" (adaptive shrinkage) is an alternative that produces similar results and also gives posterior probabilities; use it if you need compatibility with a contrast vector rather than a coefficient name. Avoid "normal" for standalone shrinkage calls — it was the legacy default and is strictly dominated by the other two for ranking purposes. The coef argument must match exactly one string in resultsNames(dds).
6.7b Reading the MA plot
Before moving to the workflow, it is worth knowing how to read the MA plot that plotMA(res_shrunk) produces, because it is the primary post-fit sanity check. The M axis (y) is the log-fold-change; the A axis (x) is the mean normalized count (baseMean). Every gene is a dot; red dots are statistically significant at the chosen alpha.
What a healthy MA plot shows:
- A horizontal cloud centered near zero. Most genes do not change between conditions; the cloud of grey dots should be narrow and symmetric around the x-axis. A cloud that drifts above or below zero for all genes is a sign of failed normalization.
- Narrowing toward the left (low counts). Without shrinkage, the cloud fans out dramatically at the left — the high variance of low-count estimates. With
lfcShrink(type = "apeglm"), these estimates are pulled toward zero and the cloud is tight even at low baseMean. - Significant genes (red) scattered left-to-right. The most strongly up- and down-regulated genes should span the range of baseMean values. If all red dots sit in the far-right (high-count) corner, the signal may be driven by a handful of highly expressed housekeeping genes shifting due to normalization failure.
- No banana-shaped distortion. If the cloud curves upward at one end of the baseMean axis, there is a systematic bias — often from improper sample pairing or a covariate not accounted for in the design.
6.8 The full bulk DESeq2 workflow, in one block
library(DESeq2); library(airway)
data(airway)
dds <- DESeqDataSetFromMatrix(
countData = assay(airway),
colData = colData(airway),
design = ~ cell + dex
)
# Pre-filter low-count genes (speed only — does not affect FDR)
dds <- dds[rowSums(counts(dds)) >= 10, ]
# Reference level (so log2FC sign means what you think)
dds$dex <- relevel(dds$dex, ref = "untrt")
# Fit
dds <- DESeq(dds)
# Get results for the contrast of interest
res <- results(dds, contrast = c("dex", "trt", "untrt"))
# Shrink LFCs for ranking and visualization
res <- lfcShrink(dds, coef = "dex_trt_vs_untrt", res = res, type = "apeglm")
# Filter to significant genes
sig <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
nrow(sig)Tutorial 06 Part A walks through this exact recipe with the airway dataset.
6.9 Independent filtering — what padj = NA means
You’ll see some genes have padj = NA even when their pvalue is small. This is independent filtering: DESeq2 automatically excludes very-low-count genes from BH adjustment, on the grounds that they have insufficient power to ever be significant and including them dilutes the FDR for the rest.
The filter is selected adaptively to maximize the number of detected genes at the chosen alpha. It is not a bug; it is making FDR more powerful. If you need to re-include filtered genes for a custom analysis, set independentFiltering = FALSE in results().
Cook’s-distance filtering also produces NAs: extreme outlier samples in a small experiment can break the dispersion estimate for one gene, and DESeq2 NAs it. Inspect with assays(dds)[["cooks"]].
Why padj = NA is actually helping you. When you test 15,000 genes, even just 2,000 of them having essentially zero power to be detected dilutes the BH budget for the 13,000 that do. Independent filtering removes those zero-power genes from the denominator of the BH correction, so the remaining genes get a less stringent adjusted threshold. The practical effect: you detect more true-positive genes at the same FDR level than you would if you included the un-testable ones. The genes excluded by independent filtering are low-count genes with baseMean below an adaptively chosen cutoff — not biologically important genes. You can see the cutoff with attr(res, "filterThreshold") and see what was removed with metadata(res)$filterThreshold. The summary message from summary(res) tells you how many genes fall into each category.
independentFiltering)
results(dds, independentFiltering = TRUE) is the default. It adaptively removes very-low-count genes from multiple testing correction to maximize discoveries at the chosen alpha. Default: TRUE. Setting independentFiltering = FALSE includes all genes in the BH denominator, which reduces power slightly but is sometimes needed when you must have a padj value for every tested gene (e.g. for downstream tools that require a complete result table). The filter threshold is chosen automatically and stored in attr(res, "filterThreshold"); inspect it to see how many genes were excluded and at what mean-expression level.
6.10 Common errors (bulk)
Table 3 collects the errors most frequently encountered when running the bulk DESeq2 workflow.
| You see | What’s wrong | What to do |
|---|---|---|
Wrong sign on log2FoldChange |
Reference level was alphabetical default | relevel() your factor before DESeq() |
apeglm errors with “name not found” |
You passed a contrast vector instead of coefficient name | Use coef = "<name>" from resultsNames(dds) |
Many padj = NA |
Independent filtering | This is normal; if you really want all genes, independentFiltering = FALSE |
| Dispersion plot is U-shaped | Mean-variance modeling failed; usually too few samples or batch confound | Inspect; consider model design |
| Significant gene count is 0 | Real biological effect is small, or you’re underpowered | Lower LFC threshold; check sample size; verify dispersion looks right |
Part B — Pseudobulk DE on single-cell data
6.10b Pseudobulk DE at a glance
Before diving into the implementation details, Table 4 summarizes the complete pseudobulk workflow in order. Each step is explained in the sections that follow.
| Step | Action | Key function | Pitfall to watch |
|---|---|---|---|
| 1. Aggregate | Sum UMI counts per (donor × condition × cell type) | AggregateExpression() |
Run after JoinLayers(); build colData from metadata, not parsed column names |
| 2. Filter | Drop pseudobulk columns with fewer than 10 cells | Manual count filter | A cell type absent from any donor is not a bug — handle NA levels downstream |
| 3. Build DDS | One DESeqDataSetFromMatrix per cell type |
DESeqDataSetFromMatrix() |
Include donor as blocking variable: design = ~ donor + condition |
| 4. Fit | Run the NB GLM, shrink dispersions, compute Wald test | DESeq() |
Check plotDispEsts(dds) per cell type |
| 5. Extract | Pull results for the condition contrast; shrink LFC | results() → lfcShrink() |
Use resultsNames(dds) to get the exact coefficient name for apeglm |
| 6. Combine | Bind all cell-type results into one table | map_dfr(cell_types, run_de) |
Drop padj = NA rows before downstream filtering |
In a multi-donor single-cell experiment, thousands of cells come from each donor, but the donor — not the cell — is the true unit of biological replication. When cells from the same donor are treated as independent observations in a statistical test, the effective sample size becomes artificially inflated by 100–10,000-fold, variances collapse to near zero, and the resulting p-values are anti-conservative: you detect thousands of “significant” genes that are really donor-specific idiosyncrasies, not condition effects8. Pseudobulk aggregation fixes this by summing (or averaging) raw UMI counts across all cells of a given type within each donor to produce one observation per (donor × condition × cell type). The resulting matrix has columns equal in number to the number of donor–condition pairs (typically tens, not thousands) and is analyzed with a standard bulk DE method that treats each column as one biological replicate. The unit of test is now the donor, not the cell, and inferences are statistically honest.
6.11 The cardinal sin
The single most-violated rule in scRNA-seq:
Cells from the same donor are not statistically independent. If your experimental unit is the donor, your statistical unit must be the donor too — not the cell.
Per-cell differential expression (Seurat::FindMarkers between conditions) treats every cell as an independent observation. With 5,000 cells per condition, the effective N becomes huge, and the variance estimates collapse to almost zero. The result is anti-conservative p-values — thousands of “significant” genes that wouldn’t replicate in a new cohort.
This is documented empirically: Squair et al. showed that per-cell DE methods report many more “significant” genes than honest pseudobulk methods on the same data, with no corresponding increase in real signal — most of the extra hits are artifacts of treating cells as replicates8. The statistical root cause — pseudoreplication — and the donor-level aggregation that fixes it have been laid out independently as well9, and the original single-cell DE benchmark already flagged that bulk-style methods on aggregated data behave better than naive per-cell tests10.
Pseudobulk DE fixes this. The fix is conceptually simple: aggregate per-cell counts to one observation per (donor, condition, cell type), then run a standard bulk DESeq2 model (Part A) on the aggregated matrix. The unit of test is now the donor, not the cell, and inferences are honest.
When is FindMarkers appropriate? Seurat::FindMarkers is not wrong in general — it is the right tool for identifying genes that distinguish clusters within a single sample or condition, where the question is “what genes mark this cell type?” rather than “what genes change between conditions?” It is also useful as a rapid exploratory tool when you have a single sample per condition and no pseudobulk is possible. The mistake is specifically using per-cell tests for condition-level contrasts in multi-donor experiments, where the cell-as-replicate fallacy inflates significance. Within a single subject, cells are not independent but the inflation is bounded; across subjects claiming cells as independent is the dangerous case.
Mixed-effects alternatives. A second family of cell-level methods addresses pseudoreplication by modeling each cell while explicitly accounting for donor-level non-independence through a random effect. Methods in this family — nebula, glmmTMB, and MAST with a random effect on sample — add a per-donor random intercept to a cell-level count model, so the effective unit of replication is correctly the donor even though the model fits at the cell level. The literature is genuinely divided: some benchmarks favor pseudobulk for simplicity and robustness8, while others find well-specified mixed models competitive9. The practical stance for this workshop: default to pseudobulk, which is simple, hard to fool, and interpretable. Reach for a mixed-effects cell-level model when you have many cells per donor and need the extra power the cell-level information can provide, and only after verifying the random-effect specification matches your data’s structure.
6.12 The aggregation step
For each (donor × cell type) combination, sum the UMI counts across all cells. The result:
- Rows: genes (same as the per-cell counts)
- Columns: (donor × cell type) pseudobulk samples — typically dozens, not thousands
- Values: integer total UMI counts
This donor-level aggregation is the heart of the muscat framework for multi-sample, multi-condition single-cell DE11, which formalizes pseudobulk and provides convenience wrappers around it. In Seurat:
pb <- AggregateExpression(seu,
assays = "RNA",
slot = "counts",
group.by = c("donor", "stim", "celltype"),
return.seurat = FALSE)$RNAThe output has columns named like donor1_CTRL_CD14Mono. Build a parallel colData data frame describing each pseudobulk column’s metadata:
meta_pb <- data.frame(
group_id = colnames(pb),
donor = sub("_.*", "", colnames(pb)),
stim = sub("_[^_]+$", "", sub("^[^_]+_", "", colnames(pb))),
celltype = sub(".*_", "", colnames(pb))
)(Or build it from the original seu@meta.data with distinct() — more robust than parsing column names if your sample/cell-type names contain underscores.)
6.12b Reading the pseudobulk aggregation output
The AggregateExpression call returns a matrix where columns correspond to (donor × condition × cell type) combinations. Before running any statistics, inspect the resulting dimensions and the table(meta_pb$celltype, meta_pb$condition) printout to understand your data structure. For ifnb with eight donors, two conditions (CTRL and STIM), and eight major cell types, a fully balanced dataset would produce 128 columns (8 × 2 × 8). In practice you will typically have fewer because some cell types are absent or too sparse in specific donors — this is biologically normal. The fraction of missing cells-per-group combinations is itself an important annotation: a cell type appearing in only one of eight donors under one condition is genuinely rare or condition-specific, not a data quality failure.
The column naming convention used by AggregateExpression follows the order of group.by and joins values with underscores. When any level of group.by contains underscores (common in cell-type annotations like "CD14+ Monocytes") or starts with a digit (the ifnb donor IDs), the resulting column names can be ambiguous or malformed. The tutorial pre-empts this by renaming donor IDs to include a letter prefix (d101, d107, …) and using spaces or hyphens only within a single grouping level. A robust defensive strategy is to build the metadata data frame directly from the original Seurat @meta.data using distinct() rather than parsing column names.
6.13 Filter low-cell groups
A pseudobulk column built from very few cells is noise. The de-facto field default is 10 cells minimum per (donor × condition × cell type) combination. Below that, drop the column.
cells_per_group <- seu@meta.data |>
count(donor, stim, celltype) |>
unite(group_id, donor, stim, celltype, sep = "_")
meta_pb <- meta_pb |>
left_join(cells_per_group, by = "group_id") |>
rename(n_cells = n)
keep <- meta_pb$n_cells >= 10
pb <- pb[, keep]
meta_pb <- meta_pb[keep, ]The threshold balances signal against noise. Pseudobulk built from 1–2 cells is dominated by stochastic dropout — even a real cell type’s pseudobulk profile is unrecognizable at that count. 10 is the floor; for very sparse cell types or shallow libraries, 20–50 is safer.
This is also where you discover missing data: a cell type may be present in some donors but not others. Pseudobulk gracefully handles this — the GLM is robust to missing combinations as long as the design is identifiable. But if a cell type appears in only one condition, you can’t test condition effects in it.
6.14 The per-cell-type DESeq2 loop
Run DESeq2 separately for each cell type:
run_de <- function(ct) {
keep <- meta_pb$celltype == ct
if (sum(keep) < 4 || length(unique(meta_pb$stim[keep])) < 2) {
return(NULL)
}
dds <- DESeqDataSetFromMatrix(
countData = pb[, keep],
colData = as.data.frame(meta_pb[keep, ]),
design = ~ donor + stim
)
dds <- DESeq(dds, quiet = TRUE)
res <- results(dds, contrast = c("stim", "STIM", "CTRL"))
res <- lfcShrink(dds, coef = "stim_STIM_vs_CTRL", res = res, type = "apeglm")
as.data.frame(res) |>
rownames_to_column("gene") |>
mutate(celltype = ct)
}
de_all <- map_dfr(unique(meta_pb$celltype), run_de) |>
filter(!is.na(padj))Note: Tutorial 06 uses the simpler design = ~ condition for portability (it works whether donors are real or synthetic), which is why the Going further section there shows upgrading to ~ donor + condition. When you have real paired donors (as the muscData::Kang18_8vs8() loader provides via the ind column — eight lupus donors, each measured in both CTRL and STIM), adding donor as a blocking variable is the right choice — it blocks on donor-level variance, the effect is estimated within each donor, and power increases substantially compared to ~ condition alone. The ifnb pseudobulk workflow starts from the ~23,000 post-QC cells retained after Tutorial 01’s filtering step, aggregated per (donor × condition × cell type).
6.15 Why per-cell-type, not a single joint test?
You could imagine running a single pseudobulk DESeq2 with design = ~ donor + celltype + stim + celltype:stim and reading off the celltype:stim interaction term. Why don’t we?
Three reasons argue against the joint model. First, different cell types have different dispersion structures: mean-variance trends differ between, say, CD4 T cells and CD14 monocytes, and a single pooled dispersion fit will over-shrink some cell types and under-shrink others. Second, different cell types have different tested gene sets — a gene expressed only in monocytes should not be tested at all in T cells, and the per-cell-type loop applies the right gene filter to each model automatically. Third, the questions are different: per-cell-type DE answers “in CD14 monocytes, which genes change?” while the joint interaction term answers “which genes change more in CD14 monocytes than in the average cell type?” — a valid but less common question.
The cost of per-cell-type DE is multiple tests per gene (one per cell type). FDR correction within each cell type is appropriate; cross-cell-type meta-correction is rarely worth it.
What to do when a cell type has no significant genes. Before concluding that a cell type shows no DE, check three things: (1) the dispersion plot for that cell type using plotDispEsts(dds) — if dispersion is poorly estimated (all dots far from the trend, or an unusual shape), the model fit may be unreliable; (2) the number of cells per pseudobulk column, especially whether any donors have very few cells for this cell type; (3) the raw LFC distribution — are the gene-level LFCs biologically implausible (near zero for everything) or do they suggest a real but underpowered effect? A cell type with a genuine null result shows an LFC distribution centered at zero with narrow spread; an underpowered result shows the expected biology in the LFC direction but no significance. These two patterns call for very different interpretations and follow-up actions.
6.16 What the output looks like
Each row of de_all is a (gene × cell type) per-cell-type DE result. Filter:
sig <- de_all |> filter(padj < 0.05, abs(log2FoldChange) > 1)
sig |> count(celltype)This gives you the per-cell-type “significant gene count”. For a strong contrast like ifnb STIM vs CTRL, expect 50–500 significant genes per major cell type (mostly ISGs and downstream). For weaker contrasts, expect a handful.
The resulting per-cell-type DE table is the input to:
- Chapter 7 — Functional Analysis: turn the gene list into biological themes
- Chapter 8 — Differential Abundance: the parallel question — do cell-type proportions change?
6.17 Replication requirements
Pseudobulk’s honesty depends on having enough donors. Power scales with the number of independent biological replicates, not cells. Table 5 shows approximate power at a moderate effect size as a function of replicate count; the estimates assume a standard NB model with typical overdispersion values from bulk RNA-seq.
| Design (donors per group) | Approx. power at LFC = 1 | What to expect |
|---|---|---|
| 1 vs 1 | 0% (no degrees of freedom) | No valid p-values — report fold-changes only |
| 2 vs 2 | ~20% | Detect only the largest effects; expect few significant genes |
| 4 vs 4 | ~55% | Reasonable for pilot experiments; report sample size prominently |
| 8 vs 8 | ~85% | Adequate for well-powered studies |
| 12 vs 12 | ~95% | Suitable for clinical cohort analyses |
If you have only 2 vs 2, you can still run pseudobulk — it’ll be honest — but you should expect few significant hits and report sample size openly. If you have only 1 vs 1, no DE method can give you valid p-values; report descriptive fold-changes only.
6.18 Donor IDs in the workshop ifnb dataset
The workshop loads ifnb via muscData::Kang18_8vs8() from Bioconductor (see Tutorial 01), which preserves the real per-cell donor IDs from the original Kang et al. (2017) study in the ind column — eight lupus donors, each measured in both CTRL and STIM. Tutorial 06 detects this column automatically and uses the eight real donors as biological replicates, producing biologically meaningful pseudobulk p-values.
A synthetic pseudo-donor fallback is included for objects that arrive without donor metadata (e.g. older SeuratData-sourced ifnb objects, which lost the ind column). In that fallback, four pseudo-donors per condition are created by random cell assignment — a teaching shortcut that makes the DESeq2 mechanics run correctly but produces p-values with no biological validity. If you trigger the fallback, the tutorial prints a clear warning; don’t report those results in a paper.
For your own data with real donor IDs, the workflow gives valid inferences without modification.
6.19 Common errors (pseudobulk)
Table 6 collects the errors most frequently encountered when running the pseudobulk workflow.
| You see | What’s wrong | What to do |
|---|---|---|
AggregateExpression returns a list, not a matrix |
API differs across Seurat versions | In v5, returned object is Assay5; pull $RNA$counts |
| Underscores in donor or cell-type names break parsing | Hard-to-debug join failures | Build meta_pb from the original Seurat metadata, don’t parse column names |
| Many cell types have no significant genes | Genuinely small effect, or low N | Inspect dispersion plots (plotDispEsts); check sample size |
apeglm errors with “coefficient not found” |
Coef name doesn’t match | Run resultsNames(dds), copy the exact name |
| Some donor × cell type combos missing | Real biology + filter | Tabulate n_cells before fitting; document explicitly |
JoinLayers required error in v5 |
Seurat v5 assay still has split layers | Run JoinLayers(seu) before AggregateExpression |
6.20 Where this material is also discussed
- Lecture: Lecture 06 — DESeq2: Bulk Fundamentals + Pseudobulk DE
- Tutorial: Tutorial 06 — DESeq2: Bulk Primer + Pseudobulk DE
- Statistical foundations: Appendix B §§ 1–8
- Functional follow-up: Chapter 7 — Functional Analysis
- Compositional follow-up: Chapter 8 — Differential Abundance
6.21 Going further
The model paper is DESeq21, with apeglm for LFC shrinkage7; its alternatives are edgeR2 and limma-voom3,4, and the bulk benchmark showing replicates beat depth is essential design reading5. For the single-cell side, the must-reads are the pseudoreplication papers8,9, the DE benchmark10, and the muscat pseudobulk framework11. The full curated list is on Key Papers, Reviews & Benchmarks.