Why do I get zero significant genes after FDR correction in Seurat?

Seurat's FindMarkers() defaults to Bonferroni, not Benjamini-Hochberg, and Bonferroni corrects across every gene in the object, which is extremely conservative for a genome-wide screen. Recompute p.adjust(p_val, method="BH") on the same raw p-values and see if the list reappears; if it does, you were Bonferroni-conservative, not truly null. If BH also returns nothing, check your replicate count with a pseudobulk test before assuming the comparison has no signal.

Should I use Bonferroni or Benjamini-Hochberg (FDR) for scRNA-seq differential expression?

Benjamini-Hochberg is the standard choice for exploratory marker discovery across thousands of genes, because it controls the expected proportion of false discoveries rather than the chance of any single false positive. Reserve Bonferroni for a small number of pre-specified, confirmatory hypotheses, using it as the default for whole-transcriptome marker screening, as Seurat's FindMarkers() does, routinely erases real signal.

What is independent filtering in DESeq2, and does it apply to scRNA-seq marker tests?

Independent filtering removes low-count genes based on mean normalized expression before FDR adjustment, choosing the threshold that maximizes rejections while keeping FDR controlled. DESeq2 applies it automatically when you run results(dds), but Seurat's FindMarkers() and Scanpy's rank_genes_groups() do not do this by default, so low-count, high-dispersion genes stay in the correction and either dilute power for real hits or surface as spurious noise unless you filter genes yourself beforehand.

How many biological replicates do I actually need for a valid FDR-corrected scRNA-seq comparison?

You need independent donors or samples, not cells, cells from the same donor are correlated and don't count as independent replicates. Aim for at least 3 samples per condition and run the comparison at the pseudobulk (sample) level with DESeq2, edgeR, or limma-voom; fewer than that and an honest test will often return few or no significant genes, which reflects real statistical power, not a broken pipeline.

Why does FindAllMarkers give a different adjusted p-value than running FindMarkers on the same two clusters directly?

FindAllMarkers() runs a separate FindMarkers() call per cluster and reports each call's own Bonferroni-corrected p_val_adj without additionally correcting for the number of clusters you tested. Apply `p_val_adj <- pmin(p_val_adj * length(x = idents.all), 1)` after the fact, or just rerun FindMarkers() directly on the pair you care about if you need a defensible number for a figure.

Sanity check · Single-Cell RNA-seq

How to Handle Multiple Testing and FDR in Single-Cell RNA-seq

Bonferroni zeroes out your gene list, Benjamini-Hochberg inflates it, and neither number means anything until you fix pseudoreplication first.

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

You ran FindMarkers() or rank_genes_groups() between two clusters or two conditions, got a table back, sorted by p_val_adj, and now you're staring at one of two bad outcomes: either every gene in the panel is "significant," or almost nothing survives correction and the list you need for the figure is empty. Both are the same underlying problem wearing different faces, you tested tens of thousands of genes, sometimes across dozens of cluster pairs, and the correction method you used (often without choosing it deliberately) is either too conservative or not accounting for what actually varies in your data.

This matters because the reviewer's first question on any scRNA-seq DE table is "how did you correct for multiple testing, and how many biological replicates is that based on." If the honest answer is "one donor per condition, tested cell by cell," the p-values in your table are not defensible regardless of which correction method produced them. Getting this wrong either buries real signal under an overly strict correction, or hands a collaborator a marker list that is mostly pseudoreplication noise dressed up as FDR < 0.05.

In the next hour you can work out which failure mode you're in: check which correction your tool actually ran, recompute it the other way, count your true biological replicates, and run the pseudobulk comparison that tells you whether your per-cell hits are real. That's enough to know whether your gene list needs a different test, a different filter, or is simply reporting an honest null.

What it looks like when it's happening

FindMarkers() between two clusters returns hundreds of genes with p_val_adj < 0.05, and nearly every gene in the panel clears the threshold
After Seurat's default Bonferroni correction across all ~20,000 genes, p_val_adj is 1 for almost every row and the significant gene list is empty
A gene shows up as significant in FindAllMarkers()'s per-cluster table but loses significance (or gets a very different p_val_adj) when you rerun FindMarkers() directly on that same cluster pair
scanpy.tl.rank_genes_groups() gives a very different count of significant genes depending on whether corr_method is 'benjamini-hochberg' or 'bonferroni'
The DE gene count explodes when you test at the per-cell level but collapses to near zero once you pseudobulk by sample and rerun DESeq2/edgeR
The same one or two highly expressed genes (often ribosomal or mitochondrial) top the marker list for almost every cluster in the dataset
Adjusted p-values barely change no matter how aggressively you pre-filter low-count genes, or they change dramatically depending on the filter threshold
A reviewer or collaborator asks how many donors/samples the comparison is based on and the honest answer is one or two per condition

Why it happens

Every DE call tests thousands of genes at once. At an unadjusted p < 0.05, testing 20,000 genes against pure noise yields roughly 1,000 "significant" hits by construction, that's the definition of a 5% false positive rate applied 20,000 times. Multiple testing correction exists to control that inflation, but Bonferroni and Benjamini-Hochberg control different things: Bonferroni controls the family-wise error rate (the chance of even one false positive across all tests) by dividing your significance threshold by the number of tests, which is brutally conservative when you're screening the whole transcriptome. Benjamini-Hochberg controls the false discovery rate (the expected proportion of false positives among your hits), which is far more forgiving and is the standard choice for exploratory marker discovery. Seurat's FindMarkers() defaults to Bonferroni across every feature in the object; Scanpy's rank_genes_groups() defaults to Benjamini-Hochberg. If you don't know which one ran, you don't know what your p_val_adj column actually promises.

FindAllMarkers() compounds this: it runs a separate FindMarkers() call per cluster and reports each call's own Bonferroni-corrected p_val_adj, but it does not additionally correct for the fact that you ran that test once per cluster. A gene can look like a solid marker for cluster 3 purely because you gave yourself many independent chances to find one, and the reported p_val_adj undersells how many tests actually happened.

The bigger distortion, and the one most analysts never check, is pseudoreplication. A DE test's power comes from its sample size, and in single-cell data the tempting sample size is "number of cells," which can be in the tens of thousands. But cells from the same donor are correlated: they share genotype, environment, ambient RNA background, and sequencing depth. The true unit of biological replication is the donor or sample, not the cell. Testing condition differences cell by cell treats those correlated cells as independent observations, which shrinks the standard error far below what the actual number of biological replicates supports and manufactures tiny p-values for effects that are really just donor-to-donor or batch-to-batch variation.

Low-count genes make this worse in the other direction. Genes with few UMIs anywhere in the matrix have high dispersion and are prone to zero-inflation and outlier-driven spurious calls that have nothing to do with true biology. DESeq2 handles this with independent filtering, dropping low-mean-count genes before adjustment based on an optimization curve, which improves power for the genes that remain. Most scRNA-seq DE workflows do not do this automatically, so those noisy low-count genes stay in the denominator of your correction, diluting power for the genes that matter, or surface directly as spurious hits. Layer on ambient RNA (3-35% of counts per cell are background, not the encapsulated cell) and depth differences that are confounded with condition or batch, and you get systematic technical signal that clears FDR thresholds looking exactly like biology.

The checks

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

0/8 checked · saved in this browser

01Inspect the raw p-value histogram before trusting p_val_adj
Find out which correction actually ran, then plot a histogram of the raw p_val column, not the adjusted one. Seurat's FindMarkers() defaults to Bonferroni across all genes in the object; Scanpy's rank_genes_groups() defaults to corr_method='benjamini-hochberg'.
r
```
markers <- FindMarkers(obj, ident.1 = "clusterA", ident.2 = "clusterB")
hist(markers$p_val, breaks = 50)
```
Healthy
A roughly flat distribution from 0 to 1 with a spike near 0, that spike is your real signal, the flat part is the null distribution behaving correctly.
Red flag
A histogram skewed toward 0 across almost every gene (classic sign of pseudoreplication inflating your effective sample size) or a histogram with no spike at all (no real signal present, so any correction returning few hits is correct, not broken).
02Recompute with the other correction method and compare hit counts
Bonferroni (FWER) and Benjamini-Hochberg (FDR) diverge most when thousands of genes are tested. Take the same raw p-values and adjust them both ways, then compare how many genes clear 0.05 under each.
r
```
markers$p_val_adj_BH <- p.adjust(markers$p_val, method = "BH")
sum(markers$p_val_adj < 0.05)      # Bonferroni
sum(markers$p_val_adj_BH < 0.05)   # BH
```
Healthy
The BH-significant list is a superset of the Bonferroni list, usually several-fold larger. If the two are nearly identical, either very few genes were tested or the signal is unusually strong.
Red flag
Bonferroni gives zero hits while BH gives hundreds. That is not "nothing significant", it means Bonferroni is controlling family-wise error across ~20,000 simultaneous tests, which is the wrong error rate for exploratory marker screening.
03For FindAllMarkers, check whether the per-group correction was applied
FindAllMarkers() reports each cluster's own Bonferroni p_val_adj without correcting for the number of clusters you tested against. Apply the extra correction and compare the result to a direct pairwise FindMarkers() call on the same two clusters.
r
```
all_markers$p_val_adj <- pmin(all_markers$p_val_adj * length(unique(Idents(obj))), 1)
```
Healthy
Genes flagged as significant in the corrected FindAllMarkers table are still significant when you run FindMarkers() directly on that specific cluster pair.
Red flag
A gene is a top marker in the FindAllMarkers table but is not significant (or has a very different p_val_adj) when tested directly pairwise, the reported significance didn't account for how many cluster comparisons you actually ran.
04Count true biological replicates, not cells
Cross-tabulate donors/samples against condition. The cell count is not your sample size for a between-condition DE test.
r
```
table(obj$sample, obj$condition)
```
Healthy
At least 3 independent donors or samples per condition, each contributing many cells.
Red flag
One or two donors per condition, or "replicates" that are really technical splits of the same biological sample, any per-cell p-value here is testing within-donor noise, not the condition effect.
05Cross-tab condition against batch and per-cell depth
Before trusting any FDR-significant gene list, check whether condition is confounded with the technical variables that drive DE on their own: sequencing batch and median UMI/genes per cell.
r
```
table(obj$condition, obj$batch)
aggregate(nCount_RNA ~ condition, data = obj@meta.data, FUN = median)
```
Healthy
Conditions are spread across multiple batches with comparable median depth per cell.
Red flag
One condition maps to exactly one batch, or one condition has systematically higher median UMI/cell, naive DE will over-call genes that are simply better sampled, not biologically different.
06Plot significant hits against mean expression
Make an MA-style plot (log fold-change vs. mean expression, or -log10(p_val_adj) vs. mean expression) and check where your significant genes sit on the expression axis.
Healthy
Significant genes span a range of expression levels, with most hits at moderate-to-high mean expression.
Red flag
Nearly all significant genes cluster at the lowest expression bin, the signature of dispersion-driven false positives that DESeq2's independent filtering step is specifically designed to remove before FDR adjustment is applied.
07Run the pseudobulk comparison and see if the hit list survives
Aggregate raw counts per sample x cluster (sum UMIs across cells within each sample-cluster combination), then run DESeq2, edgeR, or limma-voom on that pseudobulk matrix with donor as the unit of replication. Compare the significant gene count and list against your per-cell test.
r
```
pb <- AggregateExpression(obj, group.by = c("sample", "cluster"), assays = "RNA", slot = "counts")
# feed pb into DESeq2::DESeqDataSetFromMatrix() with sample-level design
```
Healthy
Pseudobulk returns fewer significant genes than the per-cell test, but the genes that do overlap agree in direction of fold-change.
Red flag
Per-cell testing returns hundreds of "significant" genes while pseudobulk on the exact same comparison returns none, that gap is pseudoreplication, not biology, and the per-cell number should not go in the paper.
08Check whether the same gene tops every cluster's marker list
Tabulate the top marker gene(s) across all clusters from a FindAllMarkers()/rank_genes_groups() run and count how often the same gene recurs.
Healthy
Marker lists are largely cluster-specific with limited overlap of top genes across unrelated clusters.
Red flag
The same highly expressed gene (often ribosomal or mitochondrial) tops nearly every cluster's list, a signature of ambient RNA contamination inflating apparent DE rather than true cluster-defining biology ([[lesson-43]]).

What to do about it

Switch from per-cell testing to pseudobulk DE for condition comparisons

When: You're comparing a condition or treatment effect across samples, not just annotating clusters, any test where donors, not cells, are the unit that should be replicated.

Aggregate UMI counts per sample x cluster and run DESeq2, edgeR, or limma-voom with sample as the replicate unit and donor/batch as covariates if needed. Treat the FDR from these tools as the number that goes in the paper.

Caveat: Loses power relative to per-cell testing, especially with fewer than 3-4 samples per condition, you may legitimately end up with few or zero significant genes, and that's the correct answer, not a failure of the method.

Apply the extra group correction after FindAllMarkers

When: You used FindAllMarkers() and are reporting per-cluster marker significance in a table or figure.

Multiply each cluster's p_val_adj by the number of clusters tested and cap at 1: p_val_adj <- pmin(p_val_adj * length(x = idents.all), 1), or rerun the specific pairwise comparisons you care about directly with FindMarkers().

Caveat: Makes an already conservative Bonferroni correction stricter still, which can zero out real markers in small or rare clusters, check whether that cost is acceptable before reporting it as your final list.

Switch from Bonferroni to Benjamini-Hochberg for exploratory marker discovery

When: You're screening for candidate markers or cluster-defining genes, not testing a single pre-specified hypothesis, and Bonferroni is returning an empty list.

In Seurat, recompute p_val_adj with p.adjust(markers$p_val, method = "BH"); in Scanpy, set corr_method='benjamini-hochberg' in rank_genes_groups (already the default there).

Caveat: BH is intentionally more permissive, genes passing BH at 5% still need downstream validation (fold-change magnitude, expression in an independent dataset) before you treat them as confirmed markers.

Pre-filter low-count genes before running the test

When: Your MA-style plot shows significant hits concentrated at the lowest expression bin, or the tested gene set includes thousands of genes detected in a handful of cells.

Filter to genes expressed in a minimum fraction of cells before calling FindMarkers()/rank_genes_groups() (Seurat's min.pct argument, default 0.1, is one lever) so the multiple testing budget isn't spent on genes with no power to detect anything.

Caveat: The filter threshold is itself a judgment call, set it too aggressively and you erase markers for rare cell populations that are only expressed in a small subset of cells.

Decontaminate ambient RNA before differential expression

When: The same one or two highly expressed genes dominate the marker list for nearly every cluster.

Run CellBender, DecontX, or SoupX on the raw (unfiltered) count matrix before downstream normalization and DE testing; CellBender gives the most precise background estimates and the largest improvement in marker specificity among the three ([[lesson-43]]).

Caveat: Adds a compute step and requires the raw unfiltered matrix, not just the filtered cell-by-gene matrix most pipelines start from; correction is an estimate, not a ground truth removal.

Redesign the comparison when condition is confounded with batch or depth

When: Your cross-tab shows condition mapping cleanly onto a single sequencing batch, or one condition has systematically different per-cell depth.

Include batch as a covariate in the pseudobulk model where samples exist in more than one batch per condition; where they don't, do not report an FDR value for that comparison at all, flag the confound and note that additional samples spanning batches are needed ([[lesson-365]]).

Caveat: If condition and batch are perfectly confounded (every case sample sequenced on one run, every control on another), no statistical correction fixes it, only a redesigned experiment or more samples does.

When not to "fix" it

When you only have two or three irreplaceable biological samples per arm, rare patient biopsies, a limited animal cohort, and the pseudobulk DESeq2 comparison returns zero significant genes after FDR correction, that zero is usually the honest, correct answer: the study is underpowered at the level of true biological replication. Reverting to per-cell testing to manufacture hits, or reporting raw p-values instead of adjusted ones, doesn't solve the power problem, it just launders pseudoreplication into a number you can't defend under review. Similarly, if the task is marker discovery for cluster annotation rather than a confirmatory between-condition test, chasing FDR < 0.05 for every top marker misapplies a confirmatory framework to an exploratory one, rank candidates by fold-change and detection rate (pct.1 vs pct.2), and treat the p-value as a filter, not the headline.

Five things experienced analysts do here

Report both the number of biological replicates and the number of cells in every methods section, reviewers now specifically ask for the replicate count because pseudoreplication is a known trap in scRNA-seq DE, and burying it in cell counts alone is a red flag.
Before trusting any per-cell DE p-value for a condition comparison, run the pseudobulk version on the same data, if the two disagree wildly on hit count, trust pseudobulk, not the per-cell number.
Check the correction method on every single tool call, not once per project, Seurat's FindMarkers() defaults to Bonferroni, FindAllMarkers() needs a manual extra correction for the number of groups tested, and Scanpy defaults to Benjamini-Hochberg; mixing tools without checking silently swaps FWER control for FDR control.
Treat "zero significant genes after correction" as informative rather than a bug to engineer around, it usually means the real effect size is small relative to the number of independent biological replicates you actually have, not that the test is broken.
Decouple statistical significance from biological importance in every marker table, sort and report by log fold-change and percent-expressing alongside p_val_adj, because with thousands of cells even a trivial 5% expression difference becomes formally "significant."

Questions people ask

Why do I get zero significant genes after FDR correction in Seurat?: Seurat's FindMarkers() defaults to Bonferroni, not Benjamini-Hochberg, and Bonferroni corrects across every gene in the object, which is extremely conservative for a genome-wide screen. Recompute p.adjust(p_val, method="BH") on the same raw p-values and see if the list reappears; if it does, you were Bonferroni-conservative, not truly null. If BH also returns nothing, check your replicate count with a pseudobulk test before assuming the comparison has no signal.
Should I use Bonferroni or Benjamini-Hochberg (FDR) for scRNA-seq differential expression?: Benjamini-Hochberg is the standard choice for exploratory marker discovery across thousands of genes, because it controls the expected proportion of false discoveries rather than the chance of any single false positive. Reserve Bonferroni for a small number of pre-specified, confirmatory hypotheses, using it as the default for whole-transcriptome marker screening, as Seurat's FindMarkers() does, routinely erases real signal.
What is independent filtering in DESeq2, and does it apply to scRNA-seq marker tests?: Independent filtering removes low-count genes based on mean normalized expression before FDR adjustment, choosing the threshold that maximizes rejections while keeping FDR controlled. DESeq2 applies it automatically when you run results(dds), but Seurat's FindMarkers() and Scanpy's rank_genes_groups() do not do this by default, so low-count, high-dispersion genes stay in the correction and either dilute power for real hits or surface as spurious noise unless you filter genes yourself beforehand.
How many biological replicates do I actually need for a valid FDR-corrected scRNA-seq comparison?: You need independent donors or samples, not cells, cells from the same donor are correlated and don't count as independent replicates. Aim for at least 3 samples per condition and run the comparison at the pseudobulk (sample) level with DESeq2, edgeR, or limma-voom; fewer than that and an honest test will often return few or no significant genes, which reflects real statistical power, not a broken pipeline.
Why does FindAllMarkers give a different adjusted p-value than running FindMarkers on the same two clusters directly?: FindAllMarkers() runs a separate FindMarkers() call per cluster and reports each call's own Bonferroni-corrected p_val_adj without additionally correcting for the number of clusters you tested. Apply p_val_adj <- pmin(p_val_adj * length(x = idents.all), 1) after the fact, or just rerun FindMarkers() directly on the pair you care about if you need a defensible number for a figure.

Guide · How to Choose Cell QC Thresholds in Single-Cell RNA-seq
Guide · How to Find and Remove Doublets in Single-Cell RNA-seq
Guide · How to Log-Transform Counts Without Fooling Yourself in Single-Cell RNA-seq
Guide · How to Sanity-Check Marker Genes and Cell Type Labels in Single-Cell RNA-seq
Guide · How to Tell If You Overclustered in Single-Cell RNA-seq
Glossary · Effect size
Glossary · P-value
Glossary · False discovery rate (FDR)

Sources

Understanding p value, multiple comparisons, FDR and q value — Source for the 20,000-test false positive math and the Bonferroni vs BH distinction used in why_it_happens
Seurat Differential Expression Vignette — Confirms FindMarkers() defaults to Bonferroni correction across all features
adjusted p-values in FindAllMarkers are incorrect · Issue #3384 — Source for the extra per-group correction formula applied to FindAllMarkers output
scanpy.tl.rank_genes_groups documentation — Confirms rank_genes_groups defaults to Benjamini-Hochberg and lists corr_method options
DESeq2: Clarification on independent filtering threshold values · Bioconductor Support — Source for the independent filtering mechanism referenced in checks and fixes
Multi-subject scRNA-seq Analysis (aggregateBioVar vignette) — Source for the pseudoreplication mechanism and pseudobulk aggregation fix
A comparison of methods accounting for batch effects in differential expression analysis of UMI count based single cell RNA sequencing — Source for the batch-confounded FDR inflation check and fix

Part of the Multiple testing and FDR series.