GPTomics/bio-single-cell-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).