# Single-cell RNA-seq differential expression between conditions: pseudobulk, not cells as replicates

> To compare conditions in single-cell RNA-seq, sum counts per sample within each cell type and test samples, not cells; a cell-level Wilcoxon test finds dozens of false positives. Plan enough animals. Python and R code.

Canonical: https://www.licklider.ai/guides/single-cell-pseudobulk/
Section: Guides › Nested measurements
Author: Licklider · Published 2026-11-12 · Last reviewed 2026-11-12 (UTC)

> **Quick answer.** To find genes that change between conditions (treated vs control, knockout vs wild type) in single-cell RNA-seq, the **animal or donor is the replicate, not the cell**. Within each cell type, **sum the counts of each sample's cells (pseudobulk)** and test the samples with a bulk RNA-seq method such as edgeR, DESeq2 or limma-voom, or, with dozens of samples, a mixed model with a random effect for the sample. Testing cells as replicates, as a default `FindMarkers` Wilcoxon test does, called on average 67 of 95 unchanged genes "significant" in our simulated studies (after Benjamini–Hochberg correction; Seurat's own default adjustment is the stricter Bonferroni, which reduces but does not fix the problem). Pseudobulk kept false discoveries to the share the correction allows, but with 4 mice per group and a gene-by-gene Welch's t-test it detected each doubled gene only 3.3% of the time on average, against 62.9% with 8 mice per group (edgeR or DESeq2 would do somewhat better): **plan the number of animals**, not the number of cells.

**Terms used below.** *Cell type (cluster)*: cells grouped by their expression profile. *Differential expression (DE)*: a gene's expression differs between conditions. *Pseudobulk*: the counts of all cells of one cell type from one sample, added together into one profile, like a bulk RNA-seq sample. *Biological replicate*: an independent sample, here a mouse. *UMI count*: the number of molecules detected per gene and cell (raw counts are non-negative whole numbers; decimals or negatives mean the values were normalised or logged). *edgeR, DESeq2, limma-voom*: R packages for bulk RNA-seq differential expression. *Marker gene*: a gene that distinguishes one cluster from the others. *Batch*: samples processed together; `~ batch + condition` is a model formula that adjusts for batch. *Quasi-likelihood F-test*: edgeR's recommended test for differential expression. *Counts per million (CPM)*: counts scaled by the sample's total counts. *Log2 fold change*: log2 of the ratio between conditions; 1 means doubled. *Wilcoxon rank-sum test*: compares two groups by ranks. *Welch's t-test*: compares two groups' means without assuming equal spreads. *Mixed model*: a regression with a random effect for each sample, so that cells from one sample are not treated as independent. *Benjamini–Hochberg (BH)*: a correction that controls the false discovery rate, the expected share of false positives among the genes called; *q* is the adjusted p-value. *Power*: the chance of detecting an effect that is really there. *False positive*: a "significant" result when there is no real effect.

## Which analysis for your question

| Question | Analysis | Report |
| --- | --- | --- |
| Which genes change between conditions in a cell type? | Pseudobulk per sample and cell type, then edgeR, DESeq2 or limma-voom; a mixed model with sample as a random effect is an alternative when there are dozens of samples, since its z-tests can be too liberal with a handful | Log2 fold change, q, n samples per condition |
| Which genes mark a cluster (cluster vs all other cells)? | Cell-level marker tests (for example the Wilcoxon test) are acceptable for describing clusters | Marker genes with effect sizes, labelled as descriptive |
| Does a cell type's proportion change? | Proportion per sample, Welch's t-test between conditions with samples as n ([comparing percentages](/guides/comparing-percentages-between-groups/)) | Difference in percentage points with 95% CI |
| Several batches or sequencing runs | Batch in the pseudobulk model (~ batch + condition), with both conditions in every batch; if each batch holds one condition, no analysis can separate them | As above, adjusted for batch |
| Only one sample per condition | No valid test of the condition; report as descriptive and add samples | — |

Squair and colleagues found that methods ignoring variation between biological replicates "are biased and prone to false discoveries. Indeed, the most widely used methods can discover hundreds of differentially expressed genes in the absence of biological differences", and that "pseudobulk methods abolished the false detection of DE genes" ([Squair et al. 2021](https://doi.org/10.1038/s41467-021-25960-2)). A benchmark of 18 methods concluded that "the pseudobulk methods performed generally best", with mixed models also "superior compared with the naïve single-cell methods" ([Junttila et al. 2022](https://doi.org/10.1093/bib/bbac286)). Zimmerman and colleagues propose mixed models "with a random effect for individual" and report little gain in power "after 100 cells per individual in most scenarios", recommending more individuals instead ([Zimmerman et al. 2021](https://doi.org/10.1038/s41467-021-21038-1)); their view that pseudobulk is underpowered was challenged by Murphy and Skene ([Murphy & Skene 2022](https://doi.org/10.1038/s41467-022-35519-4)), and the authors replied ([Zimmerman et al. 2022](https://doi.org/10.1038/s41467-022-35520-x)). Seurat's documentation for `FindMarkers` notes that "differential expression between groups (e.g. treatment conditions) should account for biological variation across samples, e.g. by using pseudobulk expression profiles" ([Seurat documentation](https://satijalab.org/seurat/reference/findmarkers)).

## Worked example: one cell type, 4 control and 4 treated mice

Constructed data: UMI counts for 100 genes in one cell type, 250 cells from each of 4 control and 4 treated mice. By design, **genes 001 to 005 double** in treated mice (log2 fold change 1) and **the other 95 stay the same**. Each mouse has its own level for every gene (a standard deviation of 0.3 on the natural-log scale between mice), as real animals do, and cells vary around their mouse's level. The random seed was the first one tried and was not chosen for its result. Stored code and data: [cells.csv](/guides/single-cell-pseudobulk/cells.csv), [model.py](/guides/single-cell-pseudobulk/model.py), [make_data.py](/guides/single-cell-pseudobulk/make_data.py), [analysis.py](/guides/single-cell-pseudobulk/analysis.py), [analysis.R](/guides/single-cell-pseudobulk/analysis.R) and [the analysis output](/guides/single-cell-pseudobulk/analysis_output.txt); analysis.py and analysis.R print the same output.

| Analysis | Genes with q < 0.05 |
| --- | --- |
| Cells as replicates (1,000 vs 1,000 cells), Wilcoxon test per gene | 63: all 5 changed genes **and 58 of the 95 unchanged** |
| **Pseudobulk (4 vs 4 mice), Welch's t-test on log2 CPM per gene** | **0** |
| Pseudobulk, five smallest p-values | genes 003, 002, 012, 001, 004 (p = 0.025 to 0.040) |
| Pseudobulk log2 fold change, genes 001 to 005 | 0.51, 1.09, 0.73, 0.84, 0.82 (largest among unchanged genes: 0.77) |

In this example the cell-level test called 58 unchanged genes (the simulation average is 67). It turned mouse-to-mouse differences into "treatment effects" for more than half of the unchanged genes: with 1,000 cells per group, any difference between the 4 treated and 4 control mice looks certain. The pseudobulk analysis asks the right question, whether the treated mice differ from the control mice, and with 4 mice per group the honest answer is "not clearly": four of the five changed genes are at the top of the ranking, but none passes the false discovery rate, and their estimated fold changes (0.51 to 1.09 around a true 1) scatter as much as the noise between 4 mice allows. In 1,000 simulated studies of this design ([simulation.py](/guides/single-cell-pseudobulk/simulation.py), [the simulation output](/guides/single-cell-pseudobulk/simulation_output.txt); fewer runs than on other pages because each study has at least 2,000 cells):

| Scenario | Cells as replicates: unchanged genes called (mean of 95) | Cells as replicates: at least one called | Pseudobulk: at least one unchanged gene called | Pseudobulk: share of the 5 doubled genes detected (mean) |
| --- | --- | --- | --- | --- |
| No true differences, 4 mice per group | **66.8** | 100.0% | 1.8% | — |
| 5 genes doubled, 4 mice per group | **67.4** | 100.0% | 3.0% | 3.3% (power) |
| 5 genes doubled, 8 mice per group | **67.6** | 100.0% | 22.3% | **62.9%** (power) |

More cells do not fix the cell-level test; more mice fix the pseudobulk test. These numbers come from a constructed 100-gene example; with thousands of genes, the false discovery correction is stricter and the number of mice needed is larger. With 8 mice per group, "at least one unchanged gene called" rises to 22.3% because the false discovery rate limits the share of false calls among all calls, not the chance of any false call, and allows a few alongside many true ones. Bulk RNA-seq tools such as edgeR and DESeq2 share information across genes and usually have more power than the gene-by-gene t-test used here for transparency.

The code expects one row per cell with columns cell, mouse, condition (control or treated) and one column of raw counts per gene (gene_…), for one cell type; run it once per cell type. From Seurat, combine `t(as.matrix(GetAssayData(obj, layer = "counts")))` with the sample and condition columns of `obj@meta.data`; from AnnData, `adata.X` (raw counts) with `adata.obs`. Rename the columns and labels to match, and run the code in the folder that holds the CSV file. Python needs NumPy, pandas and SciPy; the R code uses base R only, and the edgeR lines are an example to adapt (not run on this page).

```python
# Condensed: analysis.py also runs the cell-level Wilcoxon comparison
import numpy as np, pandas as pd
from scipy import stats

d = pd.read_csv("cells.csv")
genes = [c for c in d.columns if c.startswith("gene_")]
pb = d.groupby(["mouse", "condition"])[genes].sum().reset_index()   # pseudobulk: one profile per mouse
a = pb[genes].to_numpy(float); cpm = np.log2(a / a.sum(1, keepdims=True) * 1e6 + 1)
tr = (pb.condition == "treated").to_numpy()
q = stats.false_discovery_control(stats.ttest_ind(cpm[tr], cpm[~tr], equal_var=False).pvalue)
print((q < 0.05).sum(), pd.Series(q, index=genes).sort_values().head())   # 0 genes at q < 0.05; smallest q-values
```

```r
d <- read.csv("cells.csv"); genes <- grep("^gene_", names(d), value = TRUE)
pb <- aggregate(d[, genes], by = list(mouse = d$mouse, condition = d$condition), FUN = sum)   # one profile per mouse
a <- as.matrix(pb[, genes]); cpm <- log2(a / rowSums(a) * 1e6 + 1); tr <- pb$condition == "treated"
q <- p.adjust(apply(cpm, 2, function(v) t.test(v[tr], v[!tr])$p.value), "BH")
sum(q < 0.05); head(sort(q))                                                                # 0 genes at q < 0.05; smallest q-values
# Recommended in practice (edgeR, Bioconductor; not run here):
# library(edgeR); y <- DGEList(counts = t(a), group = pb$condition); y <- y[filterByExpr(y), , keep.lib.sizes = FALSE]
# y <- normLibSizes(y); design <- model.matrix(~ factor(condition, levels = c("control", "treated")), data = pb)   # control first, so coef 2 is treated vs control
# fit <- glmQLFit(y, design); topTags(glmQLFTest(fit, coef = 2))
```

## Before you analyse

1. **Plan the number of animals or donors per condition**, not the number of cells. Three per condition is a technical minimum, not a plan: the site's RNA-seq guidance is at least six ([how many biological replicates](/guides/how-many-biological-replicates/)), and our simulation needed about 8 to detect doubled genes most of the time. With 2 or 3 per condition, report fold changes as exploratory.
2. **Process conditions together.** Prepare, capture and sequence samples from both conditions in each batch, so that batch can be modelled rather than mixed up with condition.
3. **Aggregate within each cell type and sample.** Sum raw counts (not normalised values), drop sample–cell type combinations with very few cells by a rule set in advance (for example fewer than 10 cells), and test each cell type separately. Summed counts grow with the number of cells; normalising each pseudobulk sample (as edgeR and DESeq2 do) handles this, but report the number of cells per sample and check that it does not differ systematically between conditions.
4. **Use a bulk method on the pseudobulk counts** (edgeR, DESeq2 or limma-voom), with batch or sex in the design if needed ([nested data](/guides/nested-data-mixed-model/)).
5. **Keep cell-level tests for describing clusters**, such as marker genes, and say that those p-values describe cells, not conditions.

## Check what an AI assistant ran

1. **The assistant ran `FindMarkers` (or a Wilcoxon test) between conditions on all cells.** Ask for pseudobulk per sample and cell type, or a mixed model with sample as a random effect.
2. **The assistant reported hundreds of significant genes from 2 or 3 samples per condition.** Ask how many samples were the n.
3. **The assistant summed normalised or log values instead of raw counts.** Ask for raw counts summed per sample.
4. **The assistant did not include batch when samples were processed in batches.** Ask for batch in the design.
5. **The assistant summed all cells of a condition into one profile.** That leaves one sample per condition; ask for one profile per sample.

## Report it

> Differential expression between conditions was tested within [cell type] using pseudobulk profiles: raw UMI counts were summed per mouse ([4] control and [4] treated mice, [250] cells per mouse on average) and analysed with [edgeR quasi-likelihood F-tests]; p-values were adjusted across [100] genes with the Benjamini–Hochberg procedure (R [4.3.3]).

> [No gene] reached q < 0.05 with [4] mice per condition; the largest log2 fold changes were [1.09] ([gene 002]) and [0.84] ([gene 004]). Cell-level tests were used only to describe cluster markers.

Replace the bracketed values with your own.

## References

- Squair JW, Gautier M, Kathe C, et al. (2021). Confronting false discoveries in single-cell differential expression. *Nature Communications* 12:5692. [doi:10.1038/s41467-021-25960-2](https://doi.org/10.1038/s41467-021-25960-2)
- Junttila S, Smolander J, Elo LL (2022). Benchmarking methods for detecting differential states between conditions from multi-subject single-cell RNA-seq data. *Briefings in Bioinformatics* 23:bbac286. [doi:10.1093/bib/bbac286](https://doi.org/10.1093/bib/bbac286)
- Zimmerman KD, Espeland MA, Langefeld CD (2021). A practical solution to pseudoreplication bias in single-cell studies. *Nature Communications* 12:738. [doi:10.1038/s41467-021-21038-1](https://doi.org/10.1038/s41467-021-21038-1)
- Murphy AE, Skene NG (2022). A balanced measure shows superior performance of pseudobulk methods in single-cell RNA-sequencing analysis. *Nature Communications* 13:7851. [doi:10.1038/s41467-022-35519-4](https://doi.org/10.1038/s41467-022-35519-4)
- Zimmerman KD, Evans C, Langefeld CD (2022). Reply to: A balanced measure shows superior performance of pseudobulk methods in single-cell RNA-sequencing analysis. *Nature Communications* 13:7852. [doi:10.1038/s41467-022-35520-x](https://doi.org/10.1038/s41467-022-35520-x)
- Satija Lab. FindMarkers: gene expression markers of identity classes. *Seurat documentation* (checked 12 November 2026). [satijalab.org](https://satijalab.org/seurat/reference/findmarkers)

## Related

- [What counts as n?](https://www.licklider.ai/guides/what-counts-as-n/)
- [Reviewer says n is cells, not animals](https://www.licklider.ai/guides/reviewer-says-n-is-cells-not-animals/)
- [Nested data: average or mixed model?](https://www.licklider.ai/guides/nested-data-mixed-model/)
- [How many biological replicates?](https://www.licklider.ai/guides/how-many-biological-replicates/)

## How this page was made

- Written and published by Licklider. Recommendations cite peer-reviewed methods papers and official software documentation listed under References; cited passages were checked against the source text.
- The worked example uses constructed data. Its numbers were produced by stored code (Python 3.13 with SciPy 1.18.1; R 4.3.3, base R), and the two languages print identical output; the simulation was run in Python.
- Separate review passes checked the statistics, the citations and readability before publication.