Sanity check · Single-Nucleus RNA-seq
How to Sanity-Check Marker Genes and Cell Type Labels in Single-Nucleus RNA-seq
FindMarkers will happily hand you a stress signature, an ambient RNA gradient, or a pseudoreplication artifact and let you name a cluster after it.
By Ming "Tommy" Tang, Director of Bioinformatics in Big Pharma · Reviewed September 2026 · 5 min read
You clustered your nuclei, ran FindAllMarkers() or sc.tl.rank_genes_groups(), and now you're staring at a table sorted by adjusted p-value, about to write "Cluster 4 = oligodendrocytes" into your methods section. The top genes look plausible enough that you almost don't check further. That's the moment this page is for.
The problem is that FindMarkers doesn't know what a cell type is. It returns whatever genes separate one group of cells from the rest, and in snRNA-seq that's just as likely to be a dissociation-stress signature, a sequencing-depth artifact, or ambient RNA leaking into low-count droplets as it is a real identity gene. If you inherited your QC instincts from whole-cell scRNA-seq, half of them actively mislead you here: mitochondrial percentage means something different, ribosomal genes behave differently, and canonical markers built from secreted or cytoplasmic proteins can go missing for reasons that have nothing to do with biology.
Get this wrong and the error doesn't stay contained. Every downstream differential expression test, every trajectory inference, every cross-study comparison inherits the mislabeled cluster. This page gives you an ordered set of checks, cheapest first, to run against a marker list before you trust it, plus the fixes for when it's broken and the cases where the "weird" marker is actually correct.
What it looks like when it's happening
- Top marker for a cluster is FOS, JUN, EGR1, or another immediate-early/stress gene instead of anything resembling cell identity
- A ribosomal (RPS/RPL) gene sits near the top of a marker list even though ribosomes should be largely excluded from isolated nuclei
- The same gene shows up in the top 10 markers of three or four different clusters
- A canonical marker you expect from scRNA-seq literature (often a secreted protein) is barely detected anywhere in the dataset
- FindMarkers with only.pos = TRUE returns hundreds to thousands of genes at p_val_adj near zero for a single cluster
- A marker gene has high pct.1 but pct.2 is nearly as high, so the fold change is driven by a small numeric gap rather than real specificity
- One cluster has conspicuously low MALAT1 expression relative to every other cluster
- The marker list for the same annotated cell type looks completely different between your snRNA-seq run and a matched scRNA-seq dataset
Why it happens
Nuclei are not cells. When you lyse frozen tissue and pull out nuclei, you get mostly unspliced pre-mRNA rather than mature cytoplasmic transcripts, and that pre-mRNA maps disproportionately to introns: 34 to 58% of UMIs in a typical nuclei sample are intronic. That's not contamination, it's the signal, and if your pipeline doesn't count intronic reads you lose roughly half your sensitivity before you ever call a marker. Nuclei also carry almost no mitochondrial or ribosomal transcript, since both organelles are excluded during nuclear extraction, and they carry a lot of MALAT1, a nuclear lncRNA that's abundant in essentially every healthy nucleus. None of these baseline facts match what you'd expect from whole-cell data, and every one of them can masquerade as a "marker" if you don't account for it first.
The statistical side compounds this. The default Wilcoxon or t-test in FindMarkers and rank_genes_groups treats every cell (or nucleus) as an independent observation. It isn't: cells from the same animal or donor are correlated, so per-cell tests systematically underestimate variance and overstate significance. A widely cited reanalysis of an Alzheimer's disease snRNA-seq dataset found the original per-cell approach called roughly 14,000 "significant" DEGs; proper pseudobulk analysis at the sample level dropped that to 26. A marker list built the naive way isn't just noisy, it's often mostly noise dressed up with tiny p-values.
There's also a subtler biological cause specific to nuclei: apparent differential expression between clusters can reflect differences in how fast a gene's spliced transcript gets exported from the nucleus, not differences in transcription. A gene that looks "upregulated" in one cluster's nuclear RNA might just be retained in the nucleus longer there, while a gene that looks silent might be transcribed at a normal rate but exported and degraded before you ever counted it. High nuclear abundance in snRNA-seq is not a stand-in for high biological activity the way it often is in scRNA-seq.
Finally, ambient RNA contributes 3 to 35% of total counts per droplet, and this fraction isn't uniform: it concentrates in low-input or damaged nuclei, where genuinely abundant transcripts from lysed cells bleed into every neighboring droplet's barcode. A cluster with low RNA content will pick up the ambient profile of whatever cell type is most abundant in your prep, and that shows up as a "marker" that's really just background.
The checks
Run them in order. Each one tells you what healthy looks like and what the problem looks like.
0/9 checked · saved in this browser
Before plotting anything, sort the FindMarkers/rank_genes_groups output by avg_log2FC and look at pct.1 (fraction of cluster expressing the gene) against pct.2 (fraction of everything else expressing it) for every candidate marker.
rmarkers <- FindMarkers(seurat_obj, ident.1 = "4", test.use = "wilcox") markers[order(-markers$avg_log2FC), c("avg_log2FC", "pct.1", "pct.2", "p_val_adj")]- Healthy
- pct.1 well above pct.2 (e.g. 0.7 vs 0.1) with a clear gap for the top candidates.
- Red flag
- pct.1 is low (under 0.3) even for the top-ranked gene, or pct.1 and pct.2 are both high or both low. The gene either isn't expressed broadly in the cluster or isn't specific to it, regardless of what the p-value says.
Plot MALAT1 expression across all clusters. MALAT1 is a nuclear lncRNA that should be near-ubiquitous and roughly uniform across healthy nuclei.
pythonsc.pl.violin(adata, ['MALAT1'], groupby='leiden', rotation=90)- Healthy
- Broadly similar MALAT1 expression across every cluster.
- Red flag
- One cluster sits far below the dataset median (roughly under 1%) for MALAT1. That cluster is damaged nuclei or debris, and its 'markers' are QC artifacts, not identity genes.
Check the alignment/counting command (Cell Ranger --include-introns, or the STARsolo/alevin equivalent) and the reported intronic fraction in the run summary.
- Healthy
- An intronic fraction in the 34-58% range for a nuclei prep, confirming pre-mRNA is being captured.
- Red flag
- Near-zero intronic fraction. The data was processed in exonic-only 'cells' mode, genes-per-nucleus and total counts are deflated for every nucleus, and your clustering and marker list are built on roughly half the real signal.
Compute the mitochondrial UMI fraction per nucleus and compare across clusters.
rseurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, pattern = "^MT-") VlnPlot(seurat_obj, features = "percent.mt", group.by = "seurat_clusters")- Healthy
- Near-zero mitochondrial percentage (typically under 1-5%) uniformly across all clusters, since mitochondria are excluded during nuclear extraction.
- Red flag
- A cluster or subset of nuclei with elevated mito% (above ~10%). That's whole-cell or cytoplasmic contamination in those droplets, not a distinct cell type, and its markers should not be trusted.
Plot the top markers for every cluster together in a single dot plot spanning all clusters, so you see specificity or the lack of it in one view.
pythonsc.pl.dotplot(adata, var_names=top_markers, groupby='leiden')- Healthy
- Each marker lights up (large, dark dots) in essentially one cluster and is near-absent elsewhere.
- Red flag
- A marker shows meaningful expression (dot size above ~20-30%) in three or more clusters. It's not specific, and is likely driven by ambient RNA, cell-cycle phase, or sequencing depth rather than identity.
Build a short list of known stress/immediate-early genes (FOS, FOSB, JUN, JUNB, EGR1, IER2, HSPA1A, HSPA1B and similar) and intersect it with each cluster's top 10-20 markers.
- Healthy
- Little to no overlap between the stress gene list and top markers.
- Red flag
- Three or more of a cluster's top 10 markers fall in the stress gene list. The cluster is separating on processing stress, not cell identity, even in nuclei prepared from frozen tissue.
Run CellBender, SoupX, or DecontX on the raw (unfiltered) count matrix and recompute log fold change for your candidate markers on the corrected counts.
- Healthy
- Marker fold change stays roughly similar in magnitude and direction after correction.
- Red flag
- A marker's log fold change collapses toward zero once ambient RNA is removed. The apparent difference was background contamination, most severe in low-count or damaged-nuclei clusters (which can carry 3-35% ambient signal).
Aggregate raw counts to the sample level within each cluster and rerun differential expression with DESeq2 or edgeR, treating each donor/sample as one observation.
rpb <- AggregateExpression(seurat_obj, group.by = c("sample", "cluster"), assays = "RNA", slot = "counts") # feed pb$RNA into DESeq2 with sample as the experimental unit- Healthy
- The marker remains significant at a reasonable FDR (e.g. under 0.05) and points the same direction across most samples, using only a handful of donors.
- Red flag
- Per-cell Wilcoxon gives p_val_adj near zero for hundreds to thousands of genes, but pseudobulk collapses that list to a handful of genes or none. That gap is the signature of pseudoreplication, and it has been documented shrinking a ~14,000-gene per-cell DEG list down to 26 genes in a published reanalysis.
Pull transcript or gene length (or total intron length) for each of your top markers across all clusters and look for a systematic skew toward the longest genes in the genome.
- Healthy
- Marker genes span a normal range of lengths, similar to the background gene set.
- Red flag
- Top markers are consistently among the longest, most intron-rich genes genome-wide. Since intronic UMI counts scale with gene length, this is a length artifact riding on top of (or instead of) real identity signal.
What to do about it
Remove ambient RNA before calling markers
When: Fold change collapses after correction in the ambient RNA check, especially in low-count or debris-adjacent clusters.
Run CellBender (most precise background estimate per the available benchmarks), SoupX, or DecontX on the raw unfiltered matrix, then recompute FindMarkers/rank_genes_groups on the corrected counts instead of raw counts.
Caveat: These tools need the raw, unfiltered barcode matrix including empty droplets, so you can't run them after you've already discarded that data. An overestimated ambient fraction can also strip out real low-expression markers, so sanity-check a few known positive controls before and after correction.
Exclude or regress out stress genes as marker candidates
When: Three or more top markers in a cluster overlap your curated stress/immediate-early gene list.
Either drop the stress gene list from the marker candidate pool before ranking, or compute a stress module score per nucleus and regress it out prior to differential expression.
Caveat: Don't do this globally across the whole count matrix if you're comparing disease vs. control conditions, since stress-response genes can be genuine biology there. Only exclude them when you've confirmed the signal tracks with processing batch or dissociation time, not with your experimental variable.
Switch to pseudobulk for the confirmatory statistic
When: Per-cell Wilcoxon returns an implausibly long list of 'significant' genes, or you're about to report a p-value in a paper.
Aggregate counts to sample level within each cluster and run DESeq2 or edgeR. Keep the per-cell Wilcoxon test for screening and ranking candidates, but require pseudobulk confirmation before you call anything a validated marker.
Caveat: Needs at least three biological replicates (donors/samples) per group to have any power. With one or two samples per condition, pseudobulk DE is underpowered and you're stuck relying on effect size and literature concordance instead.
Drop and recluster damaged/debris nuclei rather than annotating them
When: A cluster shows low MALAT1, low intronic fraction, and low total counts together.
Filter out that population before re-running clustering, using MALAT1 and intronic fraction jointly as QC covariates, not just total UMI count.
Caveat: Don't threshold on UMI/gene count alone. Some legitimate cell types are naturally low-RNA-content, and you'll silently delete them if MALAT1 and intronic fraction both look normal but counts are just low.
Replace missing scRNA-seq markers with nuclear-detectable surrogates
When: A canonical marker from scRNA-seq literature is essentially undetected in your nuclei data.
Check whether the gene product is secreted or predominantly cytoplasmic once spliced; if so, look for a nuclear-detectable marker from a nuclei-specific reference (e.g. a published single-nucleus atlas for your tissue) instead of forcing the scRNA-seq marker to work.
Caveat: Nuclei-specific reference atlases are sparser than whole-cell ones and may not exist for your tissue or species, which can mean building and citing your own literature-supported panel rather than relying on an off-the-shelf list.
When not to "fix" it
Don't treat the biology-driven quirks of nuclei as bugs to correct. A 34-58% intronic read fraction is expected and desirable in a nuclei prep; stripping or down-weighting intronic reads throws away the sensitivity gain (roughly +59% genes per cell, +48% total UMIs) that counting introns was supposed to buy you. Near-zero mitochondrial and ribosomal signal across every cluster uniformly is correct, not a QC failure, since both organelles are excluded from nuclear extraction. High NEAT1 expression is real paraspeckle biology co-localized with MALAT1, not a QC flag, and shouldn't be filtered the way a low-MALAT1 flag is. Similarly, if a canonical scRNA-seq marker is genuinely absent because its transcript is exported and degraded in the cytoplasm before nuclear capture, that's real nuclear-vs-cytoplasmic localization, not a failed experiment, and forcing that marker back in by loosening detection thresholds will just reintroduce noise. And if a cluster's markers skew toward long genes because that population is legitimately biased toward high expression of long transcripts in a context your intronic fraction and MALAT1 checks both look normal for, don't apply a blanket length correction that could erase a real identity signal along with the artifact.
Five things experienced analysts do here
- Always plot a candidate marker on a dot plot spanning every cluster before you accept it, never just a violin plot restricted to the one cluster you're trying to name.
- Keep a standing stress/dissociation gene list (FOS, JUN, EGR1, HSPA1A/B and similar) in your annotation script and check overlap with every new marker table automatically, before you look at the biology.
- Use per-cell Wilcoxon or the Presto-accelerated version only for screening and ranking candidates; require a pseudobulk DESeq2/edgeR result before you write a marker into a figure or paper.
- Carry MALAT1 and intronic fraction as per-nucleus QC covariates from the raw matrix all the way through to your final cluster annotations, not just at the initial filtering step.
- When you borrow a marker panel from a whole-cell scRNA-seq atlas, check each gene's protein product for secretion or fast cytoplasmic export before trusting it to show up in nuclei; swap in a nuclear-detectable surrogate rather than lowering your detection threshold to force it.
Questions people ask
- Why do ribosomal genes show up as top markers in my snRNA-seq clusters if ribosomes are supposed to be excluded from nuclei?
Ribosomal protein transcripts should be largely depleted in isolated nuclei since ribosomes are cytoplasmic. If RPS/RPL genes dominate a marker list anyway, it usually points to ambient RNA contamination or incomplete removal of cytoplasmic debris in that cluster. Check the ambient RNA fraction with CellBender, SoupX, or DecontX and see whether the ribosomal signal drops after correction before trusting those markers.
- Should I use FindAllMarkers with only.pos = TRUE for snRNA-seq data?
It's fine for a first-pass screen since it keeps the output readable, but treat the resulting p-values with caution. The default per-cell Wilcoxon test inflates significance through pseudoreplication, so a huge only.pos list is a starting point for candidates, not a validated marker set. Confirm anything you plan to publish with pseudobulk DESeq2 or edgeR at the sample level.
- Do FOS and JUN show up as dissociation-stress artifacts in snRNA-seq the way they do in scRNA-seq?
Yes. Even though nuclei are extracted from frozen tissue rather than freshly dissociated live cells, the mechanical and enzymatic steps of nuclear isolation can still induce a stress-response signature. Screen every marker list against a curated immediate-early/stress gene set before naming a cluster, the same way you would for whole-cell data.
- Why doesn't FindMarkers detect a canonical marker gene I expected to see?
Two common causes: the gene's protein product is secreted or predominantly cytoplasmic once spliced, so its mature transcript is exported and degraded before nuclear capture, or intronic reads weren't counted during alignment, which deflates sensitivity for every gene including that one. Check your Cell Ranger or STARsolo settings for intron inclusion first, then check the gene's known subcellular transcript localization.
Related pages
- Guide · How to Find and Remove Doublets in Single-Nucleus RNA-seq
- Guide · How to Sanity-Check Marker Genes and Cell Type Labels in Single-Cell RNA-seq
- Guide · How to Sanity-Check Marker Genes and Cell Type Labels in Spatial Transcriptomics
- Guide · How to Tell If You Overclustered in Single-Nucleus RNA-seq
- Guide · How to Avoid Pseudoreplication in Single-Nucleus RNA-seq
- Glossary · Marker gene
- Glossary · Count matrix
- Glossary · Effect size
Related reading on the blog
Sources
- Chapter 19 Single-nuclei RNA-seq processing | Orchestrating Single-Cell Analysis with Bioconductor (OSCA) — Nuclear export rate differences and apparent DE, nuclear vs cytoplasmic localization
- Chapter 11 Single-nuclei RNA-seq processing | Advanced Single-Cell Analysis with Bioconductor — scRNA-seq markers for secreted proteins not transferring to snRNA-seq
- Avoiding false discoveries in single-cell RNA-seq by revisiting the first Alzheimer's disease dataset — Pseudoreplication inflating DEGs from ~14,000 to 26 after pseudobulk correction; QC thresholds
- Gene expression markers of identity classes, FindMarkers • Seurat — FindMarkers default Wilcoxon test and output columns
- Scanpy rank_genes_groups documentation — rank_genes_groups methods and inflated p-value limitation
- Perspectives on single-nucleus RNA sequencing in different cell types and tissues — Absence of mitochondrial and ribosomal transcripts in nuclei
Part of the Marker gene sanity series.