Bio single cell scatac analysis skills scatac analysis

Version Compatibility

Reference examples tested with: Signac 1.13+, Seurat 5.0+, ArchR 1.0+

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name to verify parameters
• Python (SnapATAC2 alternative): pip 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.

scATAC-seq Analysis

"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.

• R: Signac::CreateChromatinAssay() -> RunTFIDF() -> FindTopFeatures() -> RunSVD() -> RunChromVAR()
• R (large data, on-disk): ArchR::createArrowFiles() -> addIterativeLSI() -> addReproduciblePeakSet()
• Python (scverse, >1M cells): snapatac2.pp.add_tile_matrix() -> tl.spectral() -> tl.macs3()

Governing Principle

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 Decision Table

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 Type: Tile vs Peak vs Gene Activity

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

TF-IDF + LSI: Diagnose the Depth Component

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.

Clustering on LSI

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

Consensus Peak Calling

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.

Differential Accessibility

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')

chromVAR Motif Deviations

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".

Gene Activity (Cluster-Level Only)

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))

Doublet Detection: Homotypic vs Heterotypic

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.

QC Thresholds

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.

Common Errors

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

Related Skills

• single-cell/multimodal-integration - joining the ATAC modality with RNA (Multiome WNN/MultiVI)
• single-cell/preprocessing - shared QC and filtering concepts from scRNA-seq
• single-cell/clustering - clustering and UMAP shared with scRNA-seq
• single-cell/doublet-detection - doublet concepts and rate expectations
• atac-seq/atac-peak-calling - bulk ATAC peak-calling background (MACS shift/extend)
• atac-seq/motif-deviation - chromVAR deviation scoring in depth
• chip-seq/motif-analysis - motif databases (JASPAR/cisBP) and enrichment testing

References

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).