What is pseudoreplication in single-nucleus RNA-seq?

It's treating individual nuclei as independent statistical replicates when they actually come from a small number of donors. Nuclei from the same donor share genotype, dissociation batch, and ambient RNA background, so they're correlated, not independent. The correct unit of replication for a condition-level comparison is the donor, not the nucleus.

Why does FindMarkers give p-values like 1e-100 or smaller?

The default Wilcoxon test computes its standard error across cells, so with thousands of nuclei it implicitly assumes thousands of independent samples. If the true replication unit is 2 or 3 donors, that standard error is far too small, and any nonzero mean difference, even one driven by a single donor, comes out as statistically certain even though it wouldn't replicate.

How many donors do I need for snRNA-seq differential expression?

There's no settled minimum for every tissue type, but the pattern is consistent: brain-region snRNA-seq studies with only 2 to 3 donors per condition see per-nucleus tests report thousands of 'significant' genes collapse to a handful once tested at the donor level. More donors buys you more real power than more nuclei per donor does; power scales with donor count, not with total cell count alone.

Should I use pseudobulk or a mixed model for snRNA-seq DE?

Pseudobulk aggregation followed by DESeq2 is the simpler default and controls the false positive rate well when nuclei counts per donor are roughly balanced. A generalized linear mixed model, or MAST with a random effect for donor, outperforms pseudobulk when cell sampling is unbalanced across donors, but it's more involved to set up and fit correctly.

Can I just run Harmony or ComBat to remove donor-to-donor variation before testing?

No. Treating donor identity as a batch to correct away worsens pseudoreplication and increases false positives, because it can erase the between-donor variance your test needs to tell real, consistent effects apart from one donor's idiosyncrasy. Use integration for aligning cell types across donors for clustering or annotation, and keep donor as an explicit covariate or random effect in the actual DE model.

Sanity check · Single-Nucleus RNA-seq

How to Avoid Pseudoreplication in Single-Nucleus RNA-seq

5,000 nuclei from two donors is not 5,000 replicates: pick pseudobulk or a mixed model before you trust a single p-value out of FindMarkers.

By Ming "Tommy" Tang, Director of Bioinformatics in Big Pharma · Reviewed September 2026 · 4 min read

You ran FindMarkers() on 5,000 nuclei pooled across two or three donors per condition, and the output has hundreds of genes with p-values around 1e-100 to 1e-200. It looks like a clean hit list. It is actually the most common statistical error in single-nucleus work: treating every nucleus as an independent replicate when the real unit of replication is the donor the nucleus came from.

The stakes are not abstract. A gene list built on pseudoreplicated p-values will not reproduce in a third donor, a validation cohort, or a reviewer's reanalysis. You'll spend weeks chasing qPCR or IHC validation for genes that were never really significant, and the two or three genes that would have replicated get buried under thousands of false positives.

This page gives you a fast way to tell whether your current DE table has this problem, and a concrete path to redo it correctly with pseudobulk aggregation or a mixed model, using data you already have, in under an hour.

What it looks like when it's happening

FindMarkers or a per-nucleus Wilcoxon test returns hundreds to thousands of DE genes with p-values below 1e-100
Volcano plot shows a cloud of points with tiny p-values paired with log2FC near zero (under 0.1)
UMAP or PCA shows nuclei clustering by donor within a cell type, and the 'DE' gene set tracks that donor split rather than the condition split
Dropping one donor and rerunning the exact same test produces a substantially different top-gene list
Subsampling to equal nuclei per donor changes which genes come out significant
A pseudobulk or donor-level test on the same contrast returns a handful of genes instead of thousands

Why it happens

Nuclei from the same donor are correlated with each other in a way nuclei from different donors are not: same genotype, same dissociation batch, same ambient RNA background, same time-of-death or freeze-thaw history. A cell-level test (Wilcoxon, t-test, the Seurat default) computes its standard error across cells, so it silently assumes those thousands of nuclei are independent draws from the population. They are not. The real sample size for a donor-level biological question is the number of donors, often 2 to 3 per condition in brain-region snRNA-seq studies. Feed a test n=5,000 when the true n is 2, and the standard error collapses toward zero, so any nonzero mean difference, however small or however driven by a single donor, reads as astronomically significant.

snRNA-seq amplifies this. Roughly half the signal per nucleus is unspliced pre-mRNA, so counting pipelines have to include intronic reads or library sizes drop and cell-type calls get inconsistent, and ambient RNA from lysed cytoplasm makes up a larger share of each droplet than in whole-cell scRNA-seq. Both of these introduce donor-specific technical structure on top of donor-specific biology, so "donor" ends up explaining even more of the per-nucleus variance in snRNA-seq than it does in scRNA-seq, which makes the pseudoreplication problem worse, not milder, in nuclei data.

There's a second, separate inflation source worth naming so you don't conflate it with this one: clustering cells and then testing for DE genes between the same clusters is double dipping, and it also produces implausibly small p-values even with correct replication units. Pseudoreplication and double dipping stack on top of each other in a typical scRNA-seq/snRNA-seq workflow, which is part of why per-cell DE tables so often look too good to be true.

Droplet count data also follow a negative binomial distribution, not a normal one, so tests built on normal-approximation assumptions already misestimate variance at the cell level before the replication-unit problem is even factored in. Fixing the replication unit (pseudobulk, mixed model) and fixing the distributional assumption (negative binomial GLM, as in DESeq2) usually need to happen together.

The checks

Run them in order. Each one tells you what healthy looks like and what the problem looks like.

0/6 checked · saved in this browser

01Count donors per condition
Tabulate donor_id against condition in your metadata before running any DE test. This costs one line and tells you immediately whether you have any real replication at all.
r
```
table(seurat_obj$donor_id, seurat_obj$condition)
```
Healthy
At least 2, ideally 3+ donors per condition, with roughly similar nuclei counts contributed by each donor.
Red flag
One donor per condition (zero replication, no test can fix this), or wildly unbalanced counts like 4,000 nuclei from one donor against 200 from another.
02Identify what unit your current DE test actually used
Check the exact call: was it FindMarkers() on the per-cell object with the default Wilcoxon test, or was it run on a pseudobulk/aggregated object with test.use = 'DESeq2'?
Healthy
DE testing for any condition-level claim runs on donor-aggregated pseudobulk counts, or on a mixed model with donor as a random effect.
Red flag
Default FindMarkers() Wilcoxon test run directly on thousands of individual nuclei to answer a condition-level question.
03Plot effect size against p-value
Make a volcano plot: log2FC on x, -log10(p) on y, for your current DE table.
r
```
library(ggplot2)
ggplot(de_table, aes(avg_log2FC, -log10(p_val_adj))) + geom_point(alpha = 0.3)
```
Healthy
A spread of p-values, and the most significant genes also carry meaningful log2FC (commonly above 0.25 to 1, though the right cutoff depends on your biology).
Red flag
A wall of points with p-values below 1e-50 but log2FC clustered near zero. Statistically certain, biologically nothing.
04Leave-one-donor-out sensitivity check
Drop one donor, rerun the same DE test on the remaining donors, and compare the top 50 to 100 genes against the full-donor result. Repeat for each donor.
Healthy
The core signature is substantially preserved (most top genes still appear) when any single donor is removed.
Red flag
The top-gene list changes almost entirely depending on which donor you drop, meaning one donor was driving the whole result.
05Aggregate to pseudobulk and re-run at the donor level
Sum raw counts per donor per cell type in Seurat, then hand the aggregated matrix to DESeq2 with median-of-ratios normalization and a negative binomial GLM.
r
```
pb <- AggregateExpression(
  seurat_obj,
  group.by = c("condition", "donor_id", "seurat_annotations"),
  return.seurat = TRUE
)

dds <- DESeqDataSetFromMatrix(
  countData = pseudobulk_counts,
  colData = sample_metadata,
  design = ~ condition + batch
)
dds <- DESeq(dds)
```
Healthy
A much shorter DE gene list than the per-cell test. One published reanalysis of an snRNA-seq Alzheimer's dataset found a 549-fold reduction in significant genes at FDR 0.05 after correcting for pseudoreplication.
Red flag
The pseudobulk gene list is roughly the same size as the per-cell list (check your group.by columns, you may have aggregated across donors instead of within them), or zero genes survive because you only have 2 donors per condition and DESeq2 has no residual degrees of freedom left to test with.
06Check the design formula only contains sample-level variables
Print the DESeq2 design and confirm every term is invariant within a donor: condition, batch, sex, age. Nothing that varies cell-to-cell.
r
```
design(dds)
```
Healthy
~ condition, optionally plus batch/sex/age as covariates that take one value per donor.
Red flag
A per-cell covariate like nCount_RNA or cluster identity leaking into a formula that's supposed to be run on aggregated, donor-level data.

What to do about it

Pseudobulk aggregation with DESeq2

When: You have 3 or more donors per condition and nuclei counts per donor are roughly balanced.

Aggregate raw counts by donor and cell type (Seurat's AggregateExpression(), or scuttle::aggregateAcrossCells(), or presto::collapse_counts()), then run DESeq2 with median-of-ratios normalization and a design formula restricted to sample-level variables.

Caveat: Pseudobulk is conservative and underpowered when cell counts per donor are highly unbalanced, and with only 2 donors per condition you have almost no residual degrees of freedom, so the resulting p-values are close to meaningless regardless of how significant they look.

Mixed model with donor as a random effect

When: Cell counts are unbalanced across donors, or pseudobulk leaves you underpowered and you still need to model within-donor structure directly.

Fit a GLMM, or MAST with a random effect for individual, with condition as the fixed effect and donor as the random effect, modeling correlation within donor instead of aggregating it away.

Caveat: Harder to set up and slower to fit than pseudobulk, and it cannot invent statistical power that the experiment doesn't have: with only 2 to 3 donors per condition, power scales with donor count, not with how many extra nuclei you profiled per donor.

Keep donor identity as a modeled covariate, not something you integrate away

When: Someone proposes running ComBat or Harmony on the count matrix to 'remove donor effects' before differential expression.

Use integration methods only to align shared cell types for annotation or clustering across donors. For the DE contrast itself, keep donor as an explicit term in the pseudobulk or mixed model instead of pre-correcting the matrix.

Caveat: Treating donor as a batch to be corrected away worsens the pseudoreplication problem and inflates false positives further, because it can erase the exact between-donor variance the test needs to distinguish real signal from a single donor's idiosyncrasy.

Add donors, not nuclei, in the next experiment

When: Your pseudobulk or mixed-model result comes back with zero genes at a reasonable FDR and you have reason to believe there's real biology to find.

Budget statistical power around donor count. 5,000 nuclei from 2 donors gives you no more power to detect a replicable effect than a smaller, well-designed experiment with more donors. If tissue or nuclei are the bottleneck, weigh whether bulk RNA-seq on a sorted population across more donors answers the question more cheaply than adding single-nucleus depth to the same two donors.

Caveat: Costs more tissue, more sequencing budget, and more time than reanalyzing what you already have, so it's a next-experiment fix, not a today fix.

When not to "fix" it

If your question is descriptive rather than inferential, don't force pseudobulk. Counting how many nuclei fall into each cluster, annotating cell types on a shared reference atlas, or finding marker genes that define a cluster within one dataset are all cell-level tasks where per-cell tests are the right tool, because you're not making a donor-level claim about a condition effect. Likewise, if your "cells" really are technical replicates from the same donor and library (assessing library-prep noise, not biology), the per-cell variance is the variance you want to model, and pseudobulk would throw away the signal you're trying to measure. And if you truly have only one donor per condition, no amount of pseudobulk or mixed modeling manufactures a second replicate: say so plainly to collaborators rather than reporting a p-value that implies a comparison you can't make.

Five things experienced analysts do here

State your n as a sentence before writing any DE code: 'N donors per condition,' not 'N nuclei.' If you can't say it cleanly, stop and fix the experimental framing first.
Keep two objects: the full per-nucleus object for clustering and annotation, and a separate donor-aggregated pseudobulk object for any condition-level statistical claim. Don't reuse the clustering object's test for a replication question.
Run a leave-one-donor-out check on every DE table before trusting it. It's a five-minute rerun and it catches the case where one donor is quietly driving the whole signature.
Audit your design formula for covariates that don't vary within a donor. The moment a per-cell metric like nCount_RNA or a cluster label sneaks into a formula meant to run on aggregated samples, you're back to pseudoreplication in disguise.
Treat 'thousands of significant genes' from a 2 to 3 donor comparison as a warning sign, not a result. Real, replicable biology at that donor count usually survives correction as tens of genes, not thousands, so price out how many donors you actually need before you run the experiment, not after you see the p-values.

Questions people ask

What is pseudoreplication in single-nucleus RNA-seq?: It's treating individual nuclei as independent statistical replicates when they actually come from a small number of donors. Nuclei from the same donor share genotype, dissociation batch, and ambient RNA background, so they're correlated, not independent. The correct unit of replication for a condition-level comparison is the donor, not the nucleus.
Why does FindMarkers give p-values like 1e-100 or smaller?: The default Wilcoxon test computes its standard error across cells, so with thousands of nuclei it implicitly assumes thousands of independent samples. If the true replication unit is 2 or 3 donors, that standard error is far too small, and any nonzero mean difference, even one driven by a single donor, comes out as statistically certain even though it wouldn't replicate.
How many donors do I need for snRNA-seq differential expression?: There's no settled minimum for every tissue type, but the pattern is consistent: brain-region snRNA-seq studies with only 2 to 3 donors per condition see per-nucleus tests report thousands of 'significant' genes collapse to a handful once tested at the donor level. More donors buys you more real power than more nuclei per donor does; power scales with donor count, not with total cell count alone.
Should I use pseudobulk or a mixed model for snRNA-seq DE?: Pseudobulk aggregation followed by DESeq2 is the simpler default and controls the false positive rate well when nuclei counts per donor are roughly balanced. A generalized linear mixed model, or MAST with a random effect for donor, outperforms pseudobulk when cell sampling is unbalanced across donors, but it's more involved to set up and fit correctly.
Can I just run Harmony or ComBat to remove donor-to-donor variation before testing?: No. Treating donor identity as a batch to correct away worsens pseudoreplication and increases false positives, because it can erase the between-donor variance your test needs to tell real, consistent effects apart from one donor's idiosyncrasy. Use integration for aligning cell types across donors for clustering or annotation, and keep donor as an explicit covariate or random effect in the actual DE model.

Guide · How to Avoid Pseudoreplication in Single-Cell RNA-seq
Guide · How to Avoid Pseudoreplication in Spatial Transcriptomics
Guide · How to Avoid Pseudoreplication in Single-Cell ATAC-seq
Compare · Seurat vs DESeq2: Which One Should You Use?
Glossary · Effect size
Glossary · P-value
Glossary · Log fold change (log2FC)

Sources

Single-cell RNA-seq: Pseudobulk differential expression analysis — Supports the definition of pseudoreplication, the donor-as-replicate framing, and the pseudobulk/DESeq2 workflow and design formula guidance.
A practical solution to pseudoreplication bias in single-cell studies — Supports the GLMM vs pseudobulk comparison, the MAST random-effects approach, the ComBat warning, and the sqrt(cell count) power argument.
Quantification and statistical modeling of droplet-based single-nucleus RNA-sequencing data — Supports the intronic-read counting requirement and the negative binomial distribution of droplet-based snRNA-seq counts.
Differential expression testing — Source for the AggregateExpression() and FindMarkers(test.use = 'DESeq2') commands and the donor-level grouping guidance.
Avoiding false discoveries in single-cell RNA-seq by revisiting the first Alzheimer's disease dataset — Source for the 549-fold reduction in DE genes and the brain-region 2-3 donor pattern.
scRNA-seq as pseudobulk for DEseq2: does AggregateExpression() normalize counts? — Supports the note that AggregateExpression() sums raw counts without normalizing, requiring downstream DESeq2/edgeR normalization.

Part of the Pseudoreplication series.