bio-single-cell-scatac-analysisAnalyze single-cell ATAC-seq with Signac/ArchR (R) and SnapATAC2 (Python alternative). Use when processing scATAC fragments, choosing a framework, calling consensus peaks, running TF-IDF/LSI while diagnosing the depth component, scoring chromVAR motif deviations against GC-matched backgrounds, detecting homotypic vs heterotypic doublets, or deciding whether…
Independent GitHub snapshot · Automated classification · Original documentation governs
Reference examples tested with: Signac 1.13+, Seurat 5.0+, ArchR 1.0+
Before using code patterns, verify installed versions match. If versions differ:
packageVersion('<pkg>') then ?function_name to verify parameterspip show snapatac2 then help(module.function)If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.
"Analyze my single-cell ATAC-seq data" -> Process fragments, QC on chromatin signal, reduce dimensions with TF-IDF/LSI, cluster, call consensus peaks per cell type, and score TF motif activity.
Signac::CreateChromatinAssay() -> RunTFIDF() -> FindTopFeatures() -> RunSVD() -> RunChromVAR()ArchR::createArrowFiles() -> addIterativeLSI() -> addReproduciblePeakSet()snapatac2.pp.add_tile_matrix() -> tl.spectral() -> tl.macs3()A zero in the cell-by-peak matrix is epistemically ambiguous: it can mean "closed in this cell" (biology) or "accessible but no Tn5 fragment captured here" (sampling). With ~2 DNA copies per diploid locus and shallow per-cell coverage, sampling dominates the zeros. The matrix is near-binary by sampling statistics, not by biology; underlying accessibility is continuous but observed as a Bernoulli-like draw.
Binarization is now disfavored. Among non-zero entries the count (1 vs 2 vs >2) is informative, and collapsing to 1 discards it (Martens 2024). Model fragment counts with a count likelihood (Paired-Insertion Counting, SnapATAC2; PoissonVI), never read counts (PCR noise). Caveat: the extra information lives in the count=2 tier, so the benefit scales with sequencing depth; binarized analyses of shallow data are leaving little on the table, deep data more.
Per-cell signal is near-binary by sampling, so single-cell single-gene quantitative claims are unreliable; aggregate to cluster/pseudobulk for graded signal.
Gene-activity scores are a weak cluster-level proxy, structurally not just empirically: (1) enhancer-to-promoter assignment is unknown, and any fixed-distance heuristic (Signac gene body + 2 kb, ArchR exponential decay to 100 kb) is wrong for genes whose enhancers sit outside the window or loop differently by cell type; (2) poised/bivalent promoters are accessible while the gene is silent, so accessibility-to-expression is not monotone. Use gene activity for cluster-level annotation and scRNA-integration anchoring only, never as a single-cell transcriptome surrogate.
The peak set depends on which cells called the peaks: peaks are called on cells already grouped, but the grouping used a feature matrix that depends on a peak/tile choice. This is a circularity. Rare populations unresolved in the first pass never get their peaks called, so their defining elements stay invisible, a self-reinforcing blind spot. This is why iterative per-cluster peak calling (ArchR iterative LSI) exists, and why testing differential accessibility on a peak set called from the same clustering is double-dipping.
Framework choice is an infrastructure decision (language, memory, multimodal needs), not a statistics decision. Scalability numbers describe each tool's most-optimized path, not every operation.
| Framework | Language / storage | Use when | Fails when |
|---|---|---|---|
| Signac | R, in-memory Seurat ChromatinAssay | Seurat-integrated multimodal (WNN), familiar Seurat API | ~10^5+ cells (RAM-bound; future workers copy the object) |
| ArchR | R, on-disk HDF5 Arrow files | Large R workflows (~1M cells), built-in iterative LSI/peak/GRN suite | Networked filesystems (HDF5 file-locking); not a portable matrix |
| SnapATAC2 | Python+Rust, AnnData backed | >1M cells (matrix-free spectral), scverse/scvi-tools stack | Less turnkey footprinting/GRN; faster-moving 2.x API |
| muon + scanpy | Python, MuData | Multimodal Python container (RNA+ATAC) | Not ATAC-optimized for the heaviest steps |
R<->Python interop (reticulate, zellkonverter, sceasy) loses information (ChromatinAssay slots, ArchR HDF5 do not round-trip); plan to stay in one ecosystem. Verify the current best-practice default against installed docs before committing.
| Matrix | What it is | Use when | Caveat |
|---|---|---|---|
| Tile/bin (500 bp) | Genome binned, no prior peaks | Initial LSI/clustering before peaks exist | Not biology-aware; 500 bp tiles vs 501 bp peaks (off-by-one feature bugs) |
| Peak (consensus) | Per-cluster MACS peaks merged to fixed width | Final accessibility quantification, DA testing | Requires peaks first; circular with clustering |
| Gene activity | Accessibility folded to per-gene scalar | Cluster annotation, scRNA-integration anchors | Weak proxy; repressed/bivalent genes fail; distal enhancers misassigned |
Goal: Reduce the sparse, near-binary, depth-confounded matrix without letting technical depth dominate.
Approach: Reweight peaks with TF-IDF, reduce with truncated SVD, then drop components that correlate with depth, diagnosed by DepthCor, not blindly dropping component 1.
obj <- RunTFIDF(obj) # method 1 (default) = log(TF x IDF), Stuart & Butler
obj <- FindTopFeatures(obj, min.cutoff = 'q0')
obj <- RunSVD(obj) # writes the 'lsi' reduction
DepthCor(obj, n = 10) # per-component Pearson correlation with nCount
# LSI_1 usually has |corr| > 0.95 with depth, but verify; occasionally it is component 2/3, or none
Component 1 captures depth ~90% of the time but the rule is symptom-based: compute each component's depth correlation (DepthCor, or ArchR corCutOff = 0.75) and drop whichever exceed the threshold. ArchR addIterativeLSI() recomputes LSI on variable features across clustering passes to reduce depth/batch artifacts. A reviewer flags blind dims = 2:30 with no depth-correlation diagnostic.
Goal: Cluster cells from the depth-cleaned LSI embedding.
Approach: Build the neighbor graph and UMAP on the retained LSI dimensions, then cluster.
dims_use <- 2:30 # set from DepthCor, not assumed
obj <- RunUMAP(obj, reduction = 'lsi', dims = dims_use)
obj <- FindNeighbors(obj, reduction = 'lsi', dims = dims_use)
obj <- FindClusters(obj, algorithm = 3, resolution = 0.5) # algorithm 3 = SLM
Goal: Call peaks per cell type and merge into a non-overlapping, reusable feature set, avoiding bias toward abundant cell types.
Approach: Pooled bulk calling misses rare-population elements; call per cluster on pseudobulk, then merge. ArchR's fixed-width iterative-overlap set is the most reproducible; Signac's CallPeaks is simpler but uses a variable-width union that drops significance metadata.
peaks <- CallPeaks(obj, group.by = 'seurat_clusters') # per-group MACS, then GRanges::reduce() union
peak_counts <- FeatureMatrix(fragments = Fragments(obj), features = peaks, cells = colnames(obj))
obj[['peaks']] <- CreateChromatinAssay(counts = peak_counts, fragments = Fragments(obj), annotation = Annotation(obj))
Fixed-width peaks (ArchR's 501 bp) remove per-peak length normalization and give a stable reusable feature space. ArchR ranks fixed-width candidates by significance, keeps the best, removes overlappers, and requires a peak in >=2 pseudobulk replicates (reproducibility, orthogonal to MACS q-value). Wrapper parameters differ (ArchR shift -75/extsize 150 with --nolambda; Signac/SnapATAC2 shift -100/extsize 200), which changes which weak peaks survive. Comparing peak sets across datasets requires re-quantifying against a unified set; peak boundaries are not portable.
Goal: Find peaks more accessible in one group, controlling for the depth confounder.
Approach: Use a logistic-regression test with total fragments as a latent variable; do not test on a peak set called from the same clustering being compared (double-dipping).
DefaultAssay(obj) <- 'peaks'
da <- FindMarkers(obj, ident.1 = 'cluster1', ident.2 = 'cluster2',
test.use = 'LR', latent.vars = 'nCount_peaks')
Goal: Find which TF motifs vary in accessibility across cells, corrected for GC content and depth.
Approach: Attach motif matches, then compute deviations against a GC- and accessibility-matched background; rank with z-scores, never raw deviations.
library(JASPAR2020); library(TFBSTools); library(motifmatchr)
library(BSgenome.Hsapiens.UCSC.hg38)
pfm <- getMatrixSet(JASPAR2020, opts = list(collection = 'CORE', tax_group = 'vertebrates', all_versions = FALSE))
obj <- AddMotifs(obj, genome = BSgenome.Hsapiens.UCSC.hg38, pfm = pfm)
obj <- RunChromVAR(obj, genome = BSgenome.Hsapiens.UCSC.hg38) # GC-matched background internally
DefaultAssay(obj) <- 'chromvar'
diff_motifs <- FindMarkers(obj, ident.1 = 'cluster1', ident.2 = 'cluster2',
mean.fxn = rowMeans, fc.name = 'avg_diff')
chromVAR's deviation is meaningful only against a GC- and accessibility-matched background; an unmatched background manufactures apparent enrichment for GC-rich motifs (most TF motifs are GC-rich). Use z-scores (background-normalized) for cross-motif ranking, raw deviations are not comparable across motifs. Motif != TF: paralogous TFs share near-identical motifs, so an enriched motif implicates a family, not a factor; motif presence != occupancy; and a footprint (TOBIAS, needs pseudobulk) is stronger occupancy evidence than motif-in-peak. Disambiguate with TF expression (Multiome) before claiming "TF X drives this program".
Goal: Approximate per-gene accessibility for marker-based annotation and scRNA anchoring.
Approach: Sum fragments over the gene body plus a promoter window; treat the output as a cluster-level aid, not measured RNA.
gene_act <- GeneActivity(obj) # gene body + 2 kb upstream, flat count, no distance weighting
obj[['ACT']] <- CreateAssayObject(counts = gene_act)
obj <- NormalizeData(obj, assay = 'ACT', scale.factor = median(obj$nCount_ACT))
Two strategies catch different doublet classes; run both and combine. Doublet callers are separate from QC metrics (TSS/nucleosome gate debris, not doublets).
| Tool | Principle | Catches | Key dependency |
|---|---|---|---|
| AMULET | >2 fragments overlapping a diploid locus -> Poisson + BH | Homotypic (same-type) | ~25k read pairs/cell for full recall |
ArchR addDoubletScores | Simulate doublets -> LSI/UMAP -> kNN; use DoubletEnrichment | Heterotypic (different-type) | LSI/UMAP quality; structurally blind to homotypic |
| scDblFinder ATAC | Simulate on nfeatures=25 aggregated meta-features | Heterotypic | Embedding quality |
AMULET silently under-calls below ~25k coverage; CNV/aneuploidy breaks its diploid null (amplified loci exceed 2 copies in true singlet cancer cells -> false positives); multinucleate/S-G2-M cells violate the <=2-copies assumption. ArchR prefers DoubletEnrichment over DoubletScore. scDblFinder uses nfeatures=25 (not 1000) and its authors recommend against clamulet.
| Metric | Signac column | Threshold | Basis |
|---|---|---|---|
| TSS enrichment | TSS.enrichment | >2-3 (Signac); >4 (ArchR human) | ENCODE signal/noise; threshold is annotation-dependent, not portable |
| Total fragments | nCount_peaks / nFrags | >1000 (often >3000) | removes empties/debris |
| Nucleosome signal | nucleosome_signal | <4 | banding quality; very low can mean over-transposition |
| FRiP | FRiP | >0.15-0.40 (study-dependent) | signal in peaks; depends on peak set and counting convention |
TSS scores are not comparable across pipelines/annotations; never port thresholds. TSSEnrichment(fast=TRUE) blocks later TSSPlot(). FRiP depends on the peak set (circular if the same cells) and counting convention (Signac counts fragments, CellRanger-ATAC counts insertions). Threshold from the joint distributions of the actual data, not copied defaults.
| Symptom | Cause | Fix |
|---|---|---|
| UMAP separates by depth, not biology | Did not drop the depth-correlated LSI component | Run DepthCor; drop components above threshold (often #1, verify) |
| Long flat run of zeros read as "closed" | Zeros are sampling-dominated, ambiguous | Interpret at cluster/pseudobulk level; check effective coverage before structural claims |
| Gene activity disagrees with RNA for a marker | Repressed/bivalent promoter is open but silent; distal enhancer outside window | Use gene activity for cluster annotation only; validate with multiome RNA |
| "Everything is GC-rich enriched" in chromVAR | Unmatched background | Use getBackgroundPeaks/RunChromVAR GC+accessibility-matched background; report z-scores |
| Rare cell type never appears | Peaks called from a coarse single-pass clustering missed its elements | Iterative per-cluster peak calling + re-clustering (ArchR iterative LSI) |
| DA peaks look inflated | Tested on a peak set called from the same clustering (double-dipping) | Call peaks independently of the comparison; treat as ranking |
| Doublets pass QC | TSS/nucleosome gate debris, not doublets | Run AMULET (homotypic) and ArchR/scDblFinder (heterotypic) and combine |
| AMULET finds few doublets in cancer | CNV breaks the diploid null; or coverage <25k pairs/cell | Use heterotypic callers in aneuploid samples; check per-cell coverage |
| "TF X drives this" from a motif | Motif implicates a family; presence != occupancy | Confirm with TF expression (multiome) and/or footprint (TOBIAS, pseudobulk) |
| QC, gene activity, and motifs all run but look wrong | Peaks, fragments, EnsDb annotation, and BSgenome are on different genome builds; coordinate mismatch is silently wrong (no crash) | Pin every reference to one build (e.g. all hg38); verify TSS enrichment and a known marker before trusting downstream |
Buenrostro JD, Giresi PG, Zaba LC, et al. Transposition of native chromatin for fast and sensitive epigenomic profiling (ATAC-seq). Nat Methods 10(12):1213-1218 (2013). Cusanovich DA, Daza R, Adey A, et al. Multiplex single-cell profiling of chromatin accessibility (TF-IDF/LSI). Science 348(6237):910-914 (2015). Stuart T, Srivastava A, Madad S, Lareau CA, Satija R. Single-cell chromatin state analysis with Signac. Nat Methods 18:1333-1341 (2021). Granja JM, Corces MR, Pierce SE, et al. ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis. Nat Genet 53:403-411 (2021). Zhang K, Zemke NR, Armand EJ, Ren B. A fast, scalable and versatile tool for analysis of single-cell omics data (SnapATAC2). Nat Methods 21(2):217-227 (2024). Schep AN, Wu B, Buenrostro JD, Greenleaf WJ. chromVAR: inferring transcription-factor-associated accessibility from single-cell epigenomic data. Nat Methods 14(10):975-978 (2017). Martens LD, Fischer DS, Theis FJ, Buettner F. Modeling fragment counts improves single-cell ATAC-seq analysis. Nat Methods 21(1):28-31 (2024). Miao Z, Kim J. Uniform quantification of single-nucleus ATAC-seq data with Paired-Insertion Counting (PIC) and a model-based insertion rate estimator. Nat Methods 21:32-36 (2024). Thibodeau A, Eroglu A, McGinnis CS, et al. AMULET: a novel read count-based method for effective multiplet detection from single-nucleus ATAC-seq data. Genome Biol 22:252 (2021). Germain P-L, Lun A, Garcia Meixide C, Macnair W, Robinson MD. Doublet identification in single-cell sequencing data using scDblFinder. F1000Research 10:979 (2022). Bentsen M, Goymann P, Schultheis H, et al. ATAC-seq footprinting unravels kinetics of transcription factor binding during zygotic genome activation (TOBIAS). Nat Commun 11:4267 (2020).