Sanity check · Single-Cell ATAC-seq
How to Sanity-Check Marker Genes and Cell Type Labels in Single-Cell ATAC-seq
A marker with a huge fold-change can still be sequencing depth or dissociation stress in disguise, here's how to tell before you write the label on the UMAP.
By Ming "Tommy" Tang, Director of Bioinformatics in Big Pharma · Reviewed September 2026 · 5 min read
You've clustered your scATAC-seq data and now you're running FindMarkers() or getMarkerFeatures() to put a name on each cluster. The output looks reassuringly familiar: gene names, log2 fold-changes, adjusted p-values, sorted top to bottom. It's tempting to take row one and call it done, the way you would with scRNA-seq marker tables you've read a hundred times.
The catch is that scATAC-seq marker genes are one inferential step removed from expression. You didn't measure transcripts; you measured Tn5 cut sites in a peak-by-cell matrix and then derived a "gene activity" score from accessibility near a promoter or gene body. That extra inference step is exactly where sequencing depth, dissociation stress, and noisy proxy scores sneak in and get ranked as if they were biology. Label a cluster from an artifact and every downstream figure, every claim about that cell type, inherits the mistake.
This page gives you an ordered set of checks to run against your own marker table in the next hour, starting with just reading the gene names and ending with motif-level and cross-modality confirmation. By the end you'll know which of your top markers survive scrutiny and which ones you should drop before anyone else builds on your annotation.
What it looks like when it's happening
- Top of your FindMarkers()/getMarkerFeatures() table is dominated by MT-, RPL/RPS, or FOS/JUN genes instead of expected lineage genes.
- The candidate marker gene lights up as a smooth gradient across the whole UMAP rather than being restricted to one cluster.
- A cluster's top markers change substantially when you switch from default Wilcoxon testing to a depth-corrected test.
- pct.1 for your top marker is low (well under 30%) even though the reported log2FC is large.
- Gene activity score for the marker doesn't correlate with paired RNA expression in Multiome data.
- A small, rare cluster returns almost no significant markers at all.
- The same 'marker' gene keeps turning up as a top hit in several unrelated clusters.
Why it happens
Gene activity in scATAC-seq is a derived score, not a direct measurement: Signac's GeneActivity() and ArchR's GeneScoreMatrix sum accessibility over a promoter window and gene body and treat that sum as a stand-in for expression. Because the underlying peak-by-cell matrix is sparse and many loci (especially housekeeping and ribosomal genes) are open in nearly every cell type, small cell-to-cell differences in library complexity get amplified into apparent "differential accessibility" at exactly those universally-open loci.
Marker tests rank by statistical separation between groups, not by biological plausibility. When the default test (Wilcoxon rank-sum, unweighted) is used without bias correction, and clusters differ even slightly in nFrags or TSS enrichment, which happens routinely because library complexity often tracks cell size or type, depth becomes a systemic confounder. Thousands of loci shift together with depth, so a depth-correlated gene can post a huge nominal significance despite carrying no cell-identity information.
Dissociation stress is a separate and biologically real mechanism, not a normalization bug. FOS, JUN, and other immediate-early genes activate transcriptionally within minutes of tissue dissociation, and the chromatin at these loci genuinely opens as part of that stress response. Unlike shot noise, this shows up as real differential accessibility, it reflects the cell's handling history, not its identity, and no amount of post-hoc normalization removes it because the chromatin state itself changed.
Finally, the choice of gene activity model and peak set shapes what you can even detect. Promoter-only, 2kb-upstream models can pick up a neighboring gene's regulatory elements and misattribute signal, while peaks called on the pooled dataset represent common accessible regions and often omit the specific regulatory elements that define rare or small cell populations. That starves small clusters of real markers while making broadly-open loci look artificially "specific" for large ones.
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
Before plotting anything, print the top 10-20 markers per cluster from your
FindMarkers()/getMarkerFeatures()output, ranked by fold-change. Scan the gene symbols against a standing denylist: mitochondrial genes (MT-), ribosomal proteins (RPL/RPS), histones (HIST1H,H2A,H2B), and immediate-early/dissociation-stress genes (FOS,JUN,JUNB,EGR1,HSPA1A,HSPA1B).- Healthy
- Top markers are recognizable lineage genes for the tissue you're profiling (e.g. CD3D, MS4A1, LYZ in PBMCs).
- Red flag
- Several of the top 10 hits are MT-, RPL/RPS, histone, or FOS/JUN-family genes instead of lineage markers.
Pull per-cell nFrags (or log10(nFrags)) and TSSEnrichment (ArchR's getCellColData(), or the equivalent metadata columns in your Signac object) and plot them as violins or boxplots grouped by cluster.
- Healthy
- QC metric distributions overlap substantially across clusters, with only the normal variation you'd expect from different cell sizes.
- Red flag
- The cluster you're trying to label also stands out with systematically higher or lower nFrags/TSSEnrichment than the rest, that's a depth axis, not a cell-identity axis.
Re-run the marker test with depth and quality modeled as covariates, then compare the new top-marker list to the original one gene by gene.
r# ArchR: bias-matched background controls for depth/quality markersGS <- getMarkerFeatures( ArchRProj = proj, groupBy = "Clusters", testMethod = "wilcoxon", bias = c("TSSEnrichment", "log10(nFrags)"), useMatrix = "GeneScoreMatrix" ) # Signac: logistic regression with fragment count as a latent variable FindMarkers(atac_obj, test.use = "LR", latent.vars = "nFrags")- Healthy
- The marker list is largely unchanged after bias correction, the same genes stay on top.
- Red flag
- The top markers shuffle substantially, or previously significant genes drop out entirely once depth is modeled as a covariate.
Look at the pct.1/pct.2 (or equivalent mean-difference) columns in the marker output. Require the marker to be accessible/active in most cells of the cluster (pct.1 high) and clearly higher than the rest (pct.1 >> pct.2), not just a large log2FC.
- Healthy
- pct.1 covers a majority of cells in the cluster and is well separated from pct.2.
- Red flag
- Huge log2FC but low pct.1 (under roughly 20-30%), the fold-change is being driven by a handful of outlier or doublet cells, not the cluster as a whole.
Violin plot or feature plot the gene activity score (or peak accessibility) for the candidate marker across every cluster, not only the one you're naming.
- Healthy
- A clean bimodal pattern: high and mostly on in the target cluster, low and mostly off everywhere else.
- Red flag
- The signal looks like a smooth gradient correlated with sequencing depth across the whole UMAP, or 'high' cells are scattered across multiple unrelated clusters.
Build a small panel of 3-5 literature-established markers for each candidate cell type, use meta-analytic, redundant marker sets rather than a single gene, and check that they co-vary in the same cluster.
- Healthy
- Several independent canonical markers for the same cell type all point to the same cluster.
- Red flag
- The cell-type call rests on one gene, and that gene isn't independently documented as a marker for that cell type in literature you'd actually cite.
Map the cluster's differentially accessible peaks to nearby genes, then test for motif enrichment against a GC/width-matched background rather than raw genomic background.
r# Signac: map DA peaks to nearest genes, then test motif enrichment ClosestFeature(atac_obj, regions = da_peaks) matched_bkgd <- MatchRegionStats( regions = da_peaks, bkgd = all_peaks, stats = c("GC", "width") )- Healthy
- DA peaks near the marker gene are enriched for motifs of transcription factors known to be active in that cell type.
- Red flag
- No motif enrichment survives GC/width-matched background testing, or the enriched motifs are ubiquitous factors unrelated to the expected cell identity.
For Multiome data, correlate the ATAC gene activity score with the actual RNA count for the marker gene in the same cells. For ATAC-only data, compare against a matched public scRNA-seq reference for the same tissue.
- Healthy
- Gene activity and RNA expression are positively correlated for the marker across cells or clusters.
- Red flag
- High gene activity/accessibility with little or no corresponding RNA expression, the locus is accessible ('poised') but not actually transcribed, so it isn't safe to call it a functional marker.
What to do about it
Re-run marker detection with bias correction
When: Check 2/3 shows the target cluster differs systematically in nFrags or TSSEnrichment from the rest, or the marker list changes substantially under a depth-corrected test.
Use ArchR's getMarkerFeatures(..., bias = c('TSSEnrichment','log10(nFrags)')) or Signac's FindMarkers(test.use='LR', latent.vars='nFrags'), and only keep markers that stay significant after correction.
Caveat: If a cell type genuinely differs in chromatin complexity or nucleus size for real biological reasons, forcing depth correction can wash out true signal along with the artifact, check the corrected list against known biology before trusting it blindly.
Strip known artifact gene classes before ranking
When: Mitochondrial, ribosomal, histone, or immediate-early stress genes are cluttering the top of the marker table across multiple clusters.
Exclude MT-, RPL/RPS, HIST1H/H2A/H2B, and FOS/JUN/JUNB/EGR1/HSPA1A/HSPA1B (or your organism's equivalents) from the candidate list before interpreting results, or compute and regress out a stress/cell-cycle score.
Caveat: These same genes can be legitimate biology in specific contexts, such as MKI67 in a genuinely proliferating population or a stress response you're deliberately studying, don't exclude them reflexively.
Require multi-marker, literature-matched confirmation
When: A cell-type call currently rests on a single gene.
Curate a small panel of 3-5 independent, well-established markers per expected cell type for your tissue, and only assign the label if several of them agree in the same cluster.
Caveat: Takes longer per cluster and needs a decent marker reference for the tissue; for novel or poorly characterized cell states you may not have 3-5 established markers to draw on.
Validate against matched RNA when it exists
When: You have Multiome (paired RNA+ATAC) data, or a matched scRNA-seq dataset from the same tissue and condition.
Correlate the ATAC gene activity score with actual RNA counts for candidate markers, and keep only markers that agree across both modalities in your final cell-type calls.
Caveat: Not available for ATAC-only experiments without a suitable reference; matched scRNA-seq from a different study introduces its own batch effects.
Re-call peaks per cluster before trusting small-cluster markers
When: A small or rare cluster returns few or no significant markers, or its markers look unusually sparse compared to large clusters.
Peaks called on the pooled dataset are dominated by common accessible regions and under-represent regulatory elements unique to rare populations. Re-call peaks per cluster (or per broad cell type), regenerate the peak-by-cell matrix, then re-test markers for the small cluster.
Caveat: Computationally more expensive, and per-cluster peak calling on very small clusters (a few dozen cells) can itself produce unreliable peaks.
Add motif-level evidence before finalizing a label
When: The candidate marker's cell-type identity is ambiguous, or the gene activity signal looks diffuse rather than sharply cluster-restricted.
Run motif enrichment on the cluster's differentially accessible peaks against a GC/width-matched background (MatchRegionStats() + hypergeometric test) and check whether the enriched motifs correspond to transcription factors known to define that cell type.
Caveat: GC-content and generically 'sticky' motifs (e.g. CTCF) can dominate enrichment results if the background isn't properly matched, producing a false sense of confirmation.
When not to "fix" it
Don't strip out or bias-correct a signal just because it looks like the usual suspects when the biology of your experiment predicts it directly. A cluster of proliferating cells (tumor, stem-cell niche, activated lymphocytes) should legitimately show cell-cycle gene accessibility as a top marker, that's the phenotype, not an artifact. If your experimental design compares dissociation protocols, timepoints, or tissue handling, elevated FOS/JUN/EGR1 accessibility in one condition may be the actual result you're reporting, not noise to regress away. And real cell types genuinely differ in chromatin complexity (large neurons versus small lymphocytes, for instance), if depth-bias correction erases a distinction that also holds up under matched RNA and motif evidence, trust the biology over the correction.
Five things experienced analysts do here
- Keep a standing denylist of noise genes (MT-, RPL/RPS, FOS/JUN/EGR1, histones) for your organism and check every new marker table against it before you look at anything else.
- Never rank markers by fold-change alone, always pull pct.1/pct.2 or the equivalent specificity columns, since a handful of outlier cells can produce a huge fold-change with almost no biological meaning.
- Treat gene activity scores as a hypothesis, not a measurement, they're a proxy built from accessibility windows, so confirm with matched RNA whenever it's available before publishing a label.
- When a marker list looks suspicious, re-run the test with depth as a covariate before you touch anything else, it's the cheapest way to separate technical confound from real signal.
- For rare or small clusters, check whether the peak set was called on the pooled dataset before concluding the cluster has 'no markers', its distinguishing regions may never have been called as peaks in the first place.
Questions people ask
- Why do ribosomal or mitochondrial genes show up as top markers in my scATAC-seq gene activity matrix?
They usually reflect background accessibility and sequencing-depth correlation rather than true chromatin state differences. These loci are highly accessible in essentially all cells, so any depth difference between clusters surfaces there first, and residual mitochondrial signal can dominate marker lists even after standard mtDNA% filtering (5-10%).
- Should I trust FindMarkers() with default settings in Signac for scATAC-seq?
Default Wilcoxon rank-sum testing on peaks or gene activity is a reasonable starting point, but it doesn't correct for sequencing depth. Re-run with test.use='LR' and latent.vars='nFrags' (or ArchR's bias-matched getMarkerFeatures) and see whether your marker list survives the correction.
- Is a high gene activity score enough to call a cell type?
No. Gene activity is a proxy computed from promoter and gene-body accessibility, and can be high at a locus that's open but not transcribed. Confirm with matched RNA when it's available, or with multiple independent markers plus motif-level evidence before trusting a single gene's score.
- My rare cluster has almost no significant markers, is my clustering wrong?
Not necessarily. Peaks called on the pooled dataset often miss accessible regions specific to rare cell types, which starves small clusters of usable features. Re-call peaks per cluster or per broad cell type before concluding the cluster genuinely lacks distinguishing markers.
- Can I reuse scRNA-seq marker gene lists to annotate scATAC-seq clusters?
You can use them as a starting reference for gene activity scoring, but validate each one against your ATAC data first. Gene activity noise, dissociation stress, and depth artifacts mean an scRNA-seq marker panel doesn't transfer one-to-one without checking specificity and motif support.
Related pages
- Guide · How to Avoid Pseudoreplication in Single-Cell ATAC-seq
- Guide · How to Call Peaks You Can Trust in Single-Cell ATAC-seq
- Guide · How to Choose Cell QC Thresholds in Single-Cell ATAC-seq
- Guide · How to Detect Batch Effects in Single-Cell ATAC-seq
- Guide · How to Find and Remove Doublets in Single-Cell ATAC-seq
- Glossary · Log fold change (log2FC)
- Glossary · Peak calling
- Glossary · scATAC-seq
Related reading on the blog
Sources
- ArchR: Robust and scalable analysis of single-cell chromatin accessibility data — getMarkerFeatures() bias-matched testing and FDR/Log2FC cutoff thresholds
- Analyzing PBMC scATAC-seq with Signac — FindMarkers(), GeneActivity(), and ClosestFeature() usage and defaults
- Motif analysis with Signac — MatchRegionStats() and GC/width-matched hypergeometric motif testing
- ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis — gene score model design and TSS enrichment as the primary QC metric
- Learning single-cell chromatin accessibility profiles using meta-analytic marker genes — redundant, co-expressed marker panels improving annotation accuracy
- Mitochondrial genes at top gene markers - GitHub Issue #7741 — mitochondrial genes surviving as top markers even after standard mtDNA% QC filtering
- Benchmarking automated cell type annotation tools for single-cell ATAC-seq data — annotation tool performance and degradation with reference mislabeling
Part of the Marker gene sanity series.