Is FindMarkers wrong to use on spatial transcriptomics data?

Not wrong, just answering a different question than you think. FindMarkers on pooled spots or cells is fine for finding cluster markers within a dataset, where the goal is description. It becomes pseudoreplication when you use its p-values to claim a condition effect across patients or sections, because it treats every spot or cell as an independent replicate when the real replicate is the section or donor.

How many biological replicates do I need for pseudobulk DESeq2 in spatial data?

There's no fixed minimum you should trust blindly here. The honest check is the same one used in bulk RNA-seq: run a leave-one-sample-out analysis and see whether your top gene list and fold-change direction survive dropping any single section or donor. If the result depends heavily on one sample, you don't have enough replicates yet, regardless of what number you were aiming for.

Does pseudobulk aggregation throw away the spatial information I paid for?

Pseudobulk and spatial neighborhood statistics answer different questions, so you need both, not one instead of the other. Use pseudobulk with DESeq2 or muscat to compare conditions across sections or donors. Use Moran's I or a spatial mixed model to ask whether expression varies with location within a section. Collapsing to pseudobulk for the first question doesn't stop you from running spatial statistics separately for the second.

Is muscat better than running DESeq2 on pseudobulk directly?

muscat is mainly a convenience wrapper: it aggregates single-cell or spatial counts to cluster-by-sample pseudobulk and then calls a backend like edgeR, DESeq2, or limma per cluster for you. Statistically you get the same pseudobulk-vs-GLMM equivalence either way; muscat just saves you from writing the aggregation and looping code for every cluster by hand.

Why do Visium spots need pseudobulk even though each spot has whole-transcriptome coverage?

Whole-transcriptome coverage per spot doesn't make spots independent. A 55-micron Visium spot already mixes 1 to 10 cells, and neighboring spots on the same section share diffused transcripts and local tissue environment, so they're spatially autocorrelated on top of being pseudoreplicates of the same section. Aggregating spots to the section level is what restores an independent unit of measurement.

Sanity check · Spatial Transcriptomics

How to Avoid Pseudoreplication in Spatial Transcriptomics

A wall of p-values near zero on pooled spots or cells almost always means you counted the wrong thing as your sample size.

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

You ran Seurat's FindMarkers (or a pooled Wilcoxon test in Scanpy) across two conditions in your Visium or Xenium dataset, feeding it every spot or cell from every section at once. The output has hundreds of genes with p-values so small R prints them as 0, and for a moment it feels like the cleanest result you've ever pulled out of a spatial experiment.

The catch: those p-values describe how confident you are that the spots or cells inside this dataset differ. They say nothing about whether a third patient or a fresh section would show the same thing. If that result lands in a grant, a paper, or a target-nomination memo, and a validation cohort flattens it, the failure is on the sample-size line, not the biology.

This page gives you checks to run in the next hour, on the analysis you already have, to tell whether you're looking at a real condition effect or a spot/cell count wearing a p-value's clothes, plus the pseudobulk and spatial-model fixes that put the test back on the correct unit.

What it looks like when it's happening

FindMarkers or a pooled Wilcoxon test between two conditions returns hundreds of genes with p-values reported as 0 or below 1e-100.
The volcano plot shows a flat wall of points pinned to the top of the y-axis because -log10(p) is off the chart for most 'hits'.
The adjusted p-value column is full of exact zeros after BH/Bonferroni correction, for genes with only a 1.2 to 1.5-fold change.
Re-running the same test after dropping one patient or one section changes the top gene list substantially.
Plotting a 'significant' gene's expression per sample_id shows two patients in the same condition group pointing in opposite directions.
Coloring the UMAP or spatial plot by sample or section id, instead of by condition, reproduces the same pattern you attributed to biology.

Why it happens

The mechanics are the same double-dipping problem that inflates cluster-vs-cluster marker p-values: a rank-based or t-like test estimates variance from how much cells or spots disagree with each other, then divides a mean difference by that variance scaled down by sample size n. When n is the cell or spot count, thousands of units, that denominator shrinks toward zero for any nonzero difference, technical or biological, and the p-value follows it to zero. You are not testing whether the condition effect is real, you are testing whether you sampled enough units of noise to make the standard error tiny, and in single-cell and spatial data you always have.

Spatial data adds a layer scRNA-seq doesn't carry. A Visium spot is not a cell, it's a 55-micron well holding roughly 1 to 10 cells, so a spot is already a pseudoreplicate of whatever cells landed in it. Neighboring spots on the same section then share diffusing transcripts, microenvironment, and processing history, so their values are spatially autocorrelated: the reading at one spot partially predicts the reading next to it. Imaging panels like Xenium give you true single cells, but those cells still sit inside the same autocorrelated tissue and the same section-level batch effects, so the independence problem doesn't disappear, it just moves down one level.

The variance that actually matters for a condition claim lives between sections and between donors, not inside one section. Two sections from the same animal processed on the same day differ from a third section from a different animal for reasons that have nothing to do with treatment: dissection angle, permeabilization time, RNA quality, exact tissue depth. Pool spots or cells across sections without telling the model which section each one came from, and it can't separate that between-section variance from your treatment effect, so it defaults to acting as if there is none. That is why comparing two control samples with no engineered difference can still throw off hundreds of "significant" genes under a naive per-cell test.

Collapsing to one count per biological sample before testing restores the correct denominator, and it is not a heuristic fudge: a pseudobulk model equipped with a proper size-factor offset has the same point estimates and standard errors as a full negative-binomial mixed model treating donor as a random effect, while running 10 to 600x faster and converging far more reliably on lowly-expressed genes. That is why DESeq2-on-pseudobulk and muscat, not a GLMM, are the standard first move.

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 independent biological units per condition
Before running any test, tabulate how many unique samples (patients, blocks, or sections) sit inside each condition group. This is the number you should be prepared to write as 'n' in a methods section.
r
```
table(seurat_obj$condition, seurat_obj$sample_id)
sapply(split(seurat_obj$sample_id, seurat_obj$condition), function(x) length(unique(x)))
```
Healthy
Each condition has multiple distinct sample_ids, and that count matches the n you plan to report.
Red flag
A handful of sample_ids each contribute thousands of cells or spots, but the 'n' you were about to report in the DE test is the cell or spot total.
02Inspect the colData/design matrix actually fed to the model
Print the object handed to DESeqDataSetFromMatrix, or the metadata slot FindMarkers is testing on, and check the row count and what varies within it.
r
```
nrow(cluster_metadata)   # should equal (# samples) x (# cell types), not (# cells)
colData(dds)
```
Healthy
One row per sample x cell-type combination, with the condition variable constant within each sample_id.
Red flag
The row count is in the thousands and matches the raw cell or spot count.
03Compare single-cell/spot p-values against pseudobulk p-values for the same genes
Run the pooled DE test as usual, then aggregate to pseudobulk per sample x cell type with AggregateExpression and rerun DESeq2 on the identical contrast; line up the p-values for the top 20 genes from each run.
r
```
pb <- AggregateExpression(seurat_obj, group.by = c("sample_id","cell_type"), return.seurat = FALSE)
dds <- DESeqDataSetFromMatrix(pb$RNA, colData = sample_metadata, design = ~ condition)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition","case","control"))
```
Healthy
Pseudobulk p-values land orders of magnitude higher than the pooled single-cell run for genuine effects, and some single-cell 'hits' drop out entirely.
Red flag
Single-cell p-values sit at 1e-100 or below while the same genes barely clear p<0.05, or fail to, at the sample level; that gap is the pseudoreplication signature.
04Negative-control permutation within each sample
Shuffle the condition label among cells or spots within each sample (keep sample_id fixed) and rerun the exact same pooled DE test. You destroyed the true grouping but kept everything else, so nothing should come out significant.
Healthy
A near-uniform p-value histogram, with roughly the number of 'significant' genes you'd expect by chance at your alpha (about 5% at p<0.05).
Red flag
Hundreds of genes still show tiny p-values after permutation, proving the test's null distribution is wrong regardless of biology.
05Check Moran's I / spatial autocorrelation of your top genes within a section
For spatial data specifically, run spatial-variability testing on the top DE genes within a single section to see if their expression already tracks physical location independent of condition.
r
```
object <- FindSpatiallyVariableFeatures(object, assay = "Spatial", features = top_de_genes, selection.method = "moransi")
```
Healthy
Some spatial autocorrelation is normal since tissue has structure; you flag it and either add a spatial covariate or confirm your pseudobulk aggregation (summing whole sections) already absorbs it.
Red flag
A 'condition-associated' gene shows strong Moran's I inside every section, meaning the effect may be distance-from-a-landmark (tissue edge, necrotic core) that happens to correlate with which section it is.
06Leave-one-sample-out stability check
Refit the pseudobulk DESeq2 model once per sample, each time excluding one biological sample, and compare the top DE gene list and log2FC direction across the refits.
Healthy
The core gene set and fold-change direction stay consistent no matter which sample is dropped.
Red flag
Dropping a single sample flips the sign of log2FC or wipes out most of the 'significant' gene list; one section or patient is driving the whole result.

What to do about it

Rebuild the test at sample level with pseudobulk + DESeq2

When: You're comparing two or more conditions across multiple sections or patients and want a single, defensible p-value per gene.

Sum raw counts per sample x cell type with Seurat's AggregateExpression (or aggregate.Matrix on a SingleCellExperiment), build a colData table with exactly one row per sample x cell-type combination, then run DESeqDataSetFromMatrix -> DESeq -> results -> lfcShrink with design = ~ condition.

Caveat: You trade cell/spot resolution for validity, and with only two or three samples per group the model is still fragile, so inspect per-sample values, not just the final p-value, before trusting the call.

Use an offset-pseudobulk model when you need GLMM-grade variance estimates fast

When: You have many subjects or sections (tens to hundreds) and need mixed-model statistical properties without the runtime, since a full NB-GLMM can take hours and fails to converge on lowly-expressed genes.

Aggregate counts per sample as above, then add an explicit offset equal to log(sum of cell-level size factors within that sample) to the model instead of relying on default library-size normalization alone.

Caveat: The equivalence to a GLMM only holds if the offset is computed correctly at the cell level before summing; get the offset wrong and you reintroduce the same bias pseudobulk was meant to fix.

Let muscat handle multi-cluster, multi-sample aggregation for you

When: You have several cell types or clusters and want differential state testing across samples and conditions without hand-rolling the aggregation and per-cluster DE calls.

Feed a SingleCellExperiment with sample_id, cluster_id and group_id colData columns into muscat's pseudobulk workflow; it aggregates to cluster-sample pseudobulk internally and calls your chosen backend (edgeR, DESeq2, limma) per cluster.

Caveat: muscat automates the bookkeeping, not the study design; it still needs real biological replicates per condition and won't rescue an experiment with one section per group.

Add a spatial covariance term when the question is about location, not condition

When: The comparison you actually care about is within a section, for example tumor core versus margin, and spatial autocorrelation between neighboring spots is what's distorting your variance, not a missing biological replicate.

Fit a spatial mixed model with an exponential covariance structure over spot distance instead of, or on top of, pseudobulk, so nearby spots are modeled as correlated rather than independent.

Caveat: This fixes within-section non-independence, it does not create biological replicates across subjects; if your real question is case versus control across patients, you still need pseudobulk across sections first.

Use per-cell/per-spot DE only as a hypothesis generator, confirm in pseudobulk

When: You're doing exploratory marker discovery or cell-type annotation and want maximum sensitivity to flag candidates, not a final claim.

Run FindMarkers as usual to rank candidates, then intersect the list with genes that are also nominally significant in the sample-level pseudobulk test before the list goes into a report, figure, or grant.

Caveat: This discards genuinely weak-but-real effects that pseudobulk lacks power to detect with few samples, so treat it as triage, not a substitute for adequate replicate numbers.

When not to "fix" it

Pseudobulk correction restores the biological replicate axis for between-condition inference; it's the wrong tool when that axis isn't your actual question. If you're annotating cell types or clusters within one section, or describing a spatial gradient (distance from a necrotic core, distance from a tumor-immune boundary) within a single sample, per-spot or per-cell resolution is exactly what you need, and collapsing to pseudobulk would erase the signal you're trying to see. Likewise, if you deliberately profiled multiple sections from the same block to estimate technical variance rather than biological variance, don't relabel those sections as biological replicates just to inflate n in a pseudobulk model; report the technical variance as what it is.

Five things experienced analysts do here

Write down your real n, unique donors or sections, before you run any DE test, not after a reviewer asks for it.
Treat any FindMarkers p-value below roughly 1e-50 on pooled cells or spots as a pseudoreplication smell, not a result, until pseudobulk confirms it.
Build the pseudobulk colData table by hand and check nrow matches the sample x cell-type count, since a silent aggregation bug quietly reintroduces cell-level n.
Run the within-sample label-permutation negative control on any new spatial DE pipeline once, so you know its baseline false-positive rate before you trust its real output.
Color every diagnostic plot by sample or section id before you color it by condition, so batch and section effects are visible before you mistake them for biology.

Questions people ask

Is FindMarkers wrong to use on spatial transcriptomics data?: Not wrong, just answering a different question than you think. FindMarkers on pooled spots or cells is fine for finding cluster markers within a dataset, where the goal is description. It becomes pseudoreplication when you use its p-values to claim a condition effect across patients or sections, because it treats every spot or cell as an independent replicate when the real replicate is the section or donor.
How many biological replicates do I need for pseudobulk DESeq2 in spatial data?: There's no fixed minimum you should trust blindly here. The honest check is the same one used in bulk RNA-seq: run a leave-one-sample-out analysis and see whether your top gene list and fold-change direction survive dropping any single section or donor. If the result depends heavily on one sample, you don't have enough replicates yet, regardless of what number you were aiming for.
Does pseudobulk aggregation throw away the spatial information I paid for?: Pseudobulk and spatial neighborhood statistics answer different questions, so you need both, not one instead of the other. Use pseudobulk with DESeq2 or muscat to compare conditions across sections or donors. Use Moran's I or a spatial mixed model to ask whether expression varies with location within a section. Collapsing to pseudobulk for the first question doesn't stop you from running spatial statistics separately for the second.
Is muscat better than running DESeq2 on pseudobulk directly?: muscat is mainly a convenience wrapper: it aggregates single-cell or spatial counts to cluster-by-sample pseudobulk and then calls a backend like edgeR, DESeq2, or limma per cluster for you. Statistically you get the same pseudobulk-vs-GLMM equivalence either way; muscat just saves you from writing the aggregation and looping code for every cluster by hand.
Why do Visium spots need pseudobulk even though each spot has whole-transcriptome coverage?: Whole-transcriptome coverage per spot doesn't make spots independent. A 55-micron Visium spot already mixes 1 to 10 cells, and neighboring spots on the same section share diffused transcripts and local tissue environment, so they're spatially autocorrelated on top of being pseudoreplicates of the same section. Aggregating spots to the section level is what restores an independent unit of measurement.

Guide · How to Choose Cell QC Thresholds in Spatial Transcriptomics
Guide · How to Detect Batch Effects in Spatial Transcriptomics
Guide · How to Detect Integration Over-Correction in Spatial Transcriptomics
Guide · How to Sanity-Check Marker Genes and Cell Type Labels in Spatial Transcriptomics
Guide · How to Tell If You Overclustered in Spatial Transcriptomics
Glossary · False discovery rate (FDR)
Glossary · Log fold change shrinkage
Glossary · Log fold change (log2FC)

Sources

Pseudobulk with proper offsets has the same statistical properties as generalized linear mixed models in single-cell case-control studies — Mathematical equivalence and speed/convergence advantages of offset-pseudobulk over GLMMs
Differential gene expression analysis of spatial transcriptomic experiments using spatial mixed models — Spatial mixed models with exponential covariance and how they control Type I error inflation from spatial autocorrelation
Confronting false discoveries in single-cell differential expression — Single-cell DE methods generating hundreds of false-positive genes even between two control samples
AggregateExpression function documentation — Seurat function for summing counts across cells into sample-level pseudobulk
PseudobulkExpression function documentation — Seurat's pseudobulk normalization function and its default parameters
Differential state analysis with muscat — muscat's cluster-sample pseudobulk aggregation workflow for differential state analysis
Single-cell and spatial transcriptomics data analysis: Pseudobulk differential expression analysis — Hands-on aggregation and DESeq2 code for pseudobulk differential expression
Analysis, visualization, and integration of spatial datasets with Seurat — Moran's I and FindSpatiallyVariableFeatures for detecting spatial autocorrelation

Part of the Pseudoreplication series.