Pseudobulk differential expression in single-cell RNA-seq: why cells are not replicates
Why cell-level tests produce false positives in scRNA-seq, how pseudobulk aggregation works, when to use mixed models, and a practical DESeq2/edgeR workflow.
Key takeaways
- In multi-sample single-cell studies, the donor or animal is the statistical unit, not the cell.
- Pseudobulk: sum counts per sample and cell type, then use bulk RNA-seq methods (DESeq2, edgeR, limma-voom).
- Pseudobulk controls false discoveries better than cell-level tests in benchmarks.
- You need at least 3, and preferably 4 or more, biological replicates per condition.
"Differential expression between conditions" is the question behind most single-cell experiments: what changes in macrophages from treated mice, or in T cells from responders? It is also where the most consequential statistical mistake in single-cell analysis happens.
The problem: thousands of cells, a handful of replicates
Suppose you profile three knockout and three wild-type mice and get 6,000 microglia per group. A cell-level test (Wilcoxon, t-test, or MAST without random effects) compares 6,000 cells to 6,000 cells. It assumes each cell is an independent observation.
It is not. Cells from one mouse share its genotype, age, cage, dissection time and handling. The real sample size for the question "does the knockout change microglial expression?" is three versus three. Treating cells as replicates confuses technical precision within an animal with biological reproducibility across animals.
The result is predictable: thousands of "significant" genes, many driven by one unusual animal. Squair and colleagues compared many methods against ground truth and found that methods ignoring between-replicate variation systematically produce false discoveries. Pseudobulk methods performed best. Zimmerman and colleagues reached the same conclusion with mixed models.
The fix: aggregate to the replicate level
Pseudobulk analysis restores the right unit of inference:
- Annotate cell types across all samples (after QC, ambient correction and doublet removal).
- Sum raw counts for every gene, across all cells of a cell type, within each sample. Use raw counts, not normalized or integrated values.
- Filter sample-cell type combinations with too few cells (we typically require at least 10 to 20 cells) and lowly expressed genes.
- Test with bulk RNA-seq methods that model biological variance: DESeq2, edgeR (quasi-likelihood) or limma-voom.
- Shrink fold changes (apeglm or ashr) and control the false discovery rate (Benjamini-Hochberg) within each cell type.
Your design matrix can then include everything bulk RNA-seq supports: paired samples (before and after treatment in the same patient), batch as a covariate, sex, age and interaction terms.
What our experts avoid: running pseudobulk on integrated or batch-corrected values. Integration outputs are for embeddings and clustering. Differential testing needs raw counts, with batch handled in the design matrix.
When to use mixed models instead
Mixed models with a random effect for donor (for example NEBULA, or MAST with a random intercept) model cell-level data while accounting for donor structure. They can help when:
- the number of cells per sample varies widely,
- you want to model cell-level covariates, or
- cell types are rare and pseudobulk profiles would be very noisy.
They are slower and more sensitive to model choices. For most designs, pseudobulk is the robust default, and we often run both as a sensitivity check.
How many replicates do you need?
Pseudobulk makes the replicate count visible, which is uncomfortable but honest. Our general guidance:
| Replicates per group | What you can claim |
|---|---|
| 1 | Descriptive differences only, no statistical test |
| 2 | Very low power; only large, consistent effects |
| 3 | Minimum for a valid test; moderate power |
| 4-6 | Good power for typical effects |
| 8+ | Subtle effects, covariates, heterogeneous human cohorts |
Human cohorts vary more between individuals than inbred mice, so they need more replicates. If you are planning an experiment, see our single-cell experimental design guide.
Interpreting the results
Pseudobulk typically gives far fewer significant genes than cell-level tests, often tens or hundreds rather than thousands. That is not a loss of sensitivity. The removed genes were mostly false positives. In our experience, the remaining genes form more coherent pathways, replicate better in follow-up experiments and survive review.
Two practical points:
- Rank genes by shrunken effect size, not by p-value alone. Tiny fold changes with tiny p-values are rarely the biology you care about.
- Check robustness: does the result hold if you drop the most extreme replicate? We run leave-one-out sensitivity analyses for key findings.
Differential abundance is a separate question
Pseudobulk DE asks whether expression changes within a cell type. Whether the amount of a cell type changes is a different question, and it needs compositional methods (scCODA, propeller) or neighborhood-based tests (Milo), again with replicates as the unit.
How OmicsQuest handles this
Every condition comparison in our single-cell RNA-seq service uses pseudobulk or mixed models by default, with sensitivity analyses for key claims. If a reviewer has asked you to "account for biological replicates", see how we resolved exactly that in seven working days.
Frequently asked questions
What is pseudobulk analysis?
Pseudobulk analysis aggregates (usually sums) the raw counts of all cells of one cell type within each sample, producing one expression profile per sample and cell type. These profiles are then tested between conditions with bulk RNA-seq methods such as DESeq2, edgeR or limma-voom, which model variation between biological replicates.
Can I do pseudobulk with only one sample per condition?
No. With one sample per condition there is no estimate of biological variability, so no valid test of a condition effect exists. You can describe differences, but statistical claims need replicates.
Is FindMarkers in Seurat wrong?
FindMarkers with default settings is fine for finding cluster markers within a dataset, as an exploratory step. It is not appropriate for comparing conditions across donors, because it treats cells as independent. Use pseudobulk or mixed models for condition comparisons.
Work with OmicsQuest. Related services:
References
- Squair JW, Gautier M, Kathe C, et al. Confronting false discoveries in single-cell differential expression. Nat Commun. 2021;12:5692.
- Zimmerman KD, Espeland MA, Langefeld CD. A practical solution to pseudoreplication bias in single-cell studies. Nat Commun. 2021;12:738.
- Crowell HL, Soneson C, Germain PL, et al. muscat detects subpopulation-specific state transitions from multi-sample multi-condition single-cell transcriptomics data. Nat Commun. 2020;11:6077.
- Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
- Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139-140.
- Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47.
- Zhu A, Ibrahim JG, Love MI. Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences. Bioinformatics. 2019;35:2084-2092.
- Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc B. 1995;57:289-300.