Chapter 15 — Bulk and Single-Cell ATAC-seq
Where this chapter sits. Companion to Lecture 15 and Tutorial 15. This is the chapter where the workshop pivots away from RNA. Prerequisites: Chapter 1 (general QC concepts; the ATAC QC metrics are different but the approach is the same).
15.1 What ATAC-seq measures, and why it’s different from RNA-seq
ATAC-seq (Assay for Transposase-Accessible Chromatin with sequencing) measures chromatin accessibility — which regions of DNA are physically open and accessible to a hyperactive Tn5 transposase1. The raw measurement is binary at each genomic position: a fragment was found there, or it wasn’t. Accessibility is the regulatory layer that gates transcription, which is why it complements RNA so well2; the single-cell version was pioneered by combinatorial-indexing atlases of accessibility across whole organisms3.
This is fundamentally different from RNA-seq, which measures transcription output (which mRNAs are present). The two layers are correlated but not redundant:
- RNA tells you what the cell is currently producing — its working transcriptome
- ATAC tells you what the cell is currently able to produce — its regulatory potential
ATAC also reveals transcription factor binding sites indirectly: at sub-nucleosomal resolution, an open chromatin region’s exact width tells you whether a TF is bound there (footprinting). This is a layer of biology entirely invisible to RNA-seq.
The Tn5 transposase used in ATAC-seq binds DNA as a homodimer and inserts sequencing adapters at both cut sites simultaneously. Because Tn5 is large (~55 kDa per monomer), it physically cannot insert into regions wrapped around nucleosomes — only nucleosome-free regions (NFRs) are accessible. The fragment length distribution therefore carries structural information: fragments shorter than ~150 bp come from NFRs (the “sub-nucleosomal” or “nucleosome-free” peak), fragments of ~200 bp come from mononucleosomal DNA, and longer fragments come from di- and tri-nucleosomal regions. A healthy ATAC library shows a characteristic ladder pattern in the fragment size distribution with peaks at multiples of ~200 bp. This ladder is used as a quality control metric (the nucleosome signal in §15.3): a flat distribution or absent ladder indicates poor Tn5 activity or a degraded sample.
Unlike scRNA-seq, scATAC-seq does not directly measure gene expression. It measures which genomic regions are accessible, and those regions include promoters, enhancers, and other regulatory elements. Signac’s GeneActivity() converts this chromatin accessibility signal into a gene-level proxy: it sums the fragment counts in a window around each gene’s transcription start site (by default the 2 kb upstream of the TSS plus the gene body). The result looks like a gene × cell expression matrix but is derived entirely from accessibility data. This “gene activity score” is an approximation — accessibility at a promoter is correlated with, but not identical to, transcription. The score is sufficient for cluster annotation and label transfer from RNA, but it should not be interpreted as quantitative expression. A cell with high gene activity for CD14 is likely a monocyte; quantitative comparisons of CD14 activity levels between clusters are less reliable than the same comparison made with RNA.
15.2 The 10x scATAC-seq chemistry
Single-cell ATAC works similarly to scRNA but with one crucial difference: the unit of capture is a fragment (a pair of Tn5 cut sites at distance ~50–500 bp), not a polyA-tailed mRNA. The 10x scATAC bead carries primers that hybridize to Tn5-introduced ME (mosaic-end) sequences rather than polyT.
A run produces:
- A fragments file (
fragments.tsv.gz) — the key output. Each line ischrom start end barcode duplicate-count - A peaks-by-cell counts matrix — peaks (genomic intervals where many fragments concentrate) × cells, integer fragment counts
- A summary
web_summary.htmlwith QC
The fragments file is what Signac and ArchR consume. The peaks-by-cell matrix is what most analyses use.
15.3 Per-cell QC metrics in scATAC-seq
Six metrics characterize scATAC cell quality; none has a universal threshold, but all have clear directional interpretations. Table 1 summarizes all six and their expected ranges for healthy PBMC data.
| Metric | What it measures | Healthy range (PBMC) | High value means | Low value means |
|---|---|---|---|---|
nCount_peaks |
Total fragments per cell | ≥ 5,000 | Very active chromatin or doublet | Empty or lysed nucleus |
nFeature_peaks |
Distinct peaks with ≥ 1 fragment | ≥ 1,000 | Doublet (two cells’ accessibility) | Low-complexity or dead cell |
| TSS enrichment score | Enrichment of fragments at ±2 kb of TSS | > 4 | Strong signal, high-quality cell | Poor cell quality; ambient chromatin dominates |
| Nucleosome signal | Ratio of mono-nucleosomal to sub-nucleosomal fragments | < 4 (lower is better) | Over-digestion or poor library quality | Normal; near-zero is expected for NFR-dominant libraries |
| Blacklist ratio | Fragments overlapping ENCODE blacklist regions | < 0.05 | Artifactual signal from repetitive elements | Expected; very low is good |
| FRiP (Reads in Peaks) | Fraction of fragments inside called peaks | > 0.4 | High signal-to-background | Background noise dominates; poor enrichment |
The four standard metrics are:
nCount_peaks — total fragments per cell. Like nCount_RNA, just at peak resolution. 5,000 fragments per cell is a healthy minimum for PBMC data; tissues with smaller or more fragile nuclei may require adjusting this threshold. Very high fragment counts (> 100,000) can indicate doublets or multiplets, analogous to the nCount_RNA upper filter in scRNA-seq QC.
nFeature_peaks — number of distinct peaks with at least one fragment. Like nFeature_RNA. A cell with many fragments but few distinct peaks has low complexity — it is repeatedly sequencing the same small set of regions, which suggests a degraded or contaminated nucleus.
TSS enrichment score — fraction of fragments at transcription start sites (within ±2 kb). Active TSSs are reliably open in essentially all cell types, so a healthy cell should show a sharp accessibility peak at TSSs. Healthy cells have TSS enrichment > 4 (high signal, low background). Low TSS = high background, suggesting poor cell quality or extensive ambient chromatin. This metric is conceptually similar to nCount_RNA being correlated with library quality in scRNA-seq — it is a signal-to-noise ratio, not a biological measurement.
Nucleosome signal / fragment size distribution — the characteristic ATAC profile shows peaks at <100 bp (nucleosome-free), ~200 bp (mononucleosome), ~400 bp (dinucleosome). A flat distribution suggests poor digestion or contamination. Signac’s NucleosomeSignal() computes the ratio of mononucleosomal to sub-nucleosomal fragments, which should be low in a healthy library (most fragments come from open, nucleosome-free regions). A high nucleosome signal in a cell means most of its detected fragments are from closed (nucleosomal) DNA — the Tn5 is cutting into nucleosome-wrapped regions, which typically means the nucleus was not properly permeabilized or was of poor quality.
Two more advanced metrics:
- Blacklist ratio — fragments overlapping ENCODE blacklist regions (repetitive sequences that produce mapping artifacts). High blacklist ratio = artifactual signal, suggesting the cell is dominated by repeats, not biology.
- Reads in peaks (FRiP) — fraction of fragments inside called peaks. >0.4 is healthy; <0.2 suggests background dominates over signal.
Tutorial 15 walks all of these in Signac.
15.4 The Signac workflow
Signac4 is the Seurat-flavored R toolkit for scATAC; ArchR5 is the main scalable alternative, and a recent benchmark compares the available computational pipelines head-to-head6. The workflow:
- Load the fragments file + peak matrix →
ChromatinAssay - Compute QC metrics (TSS enrichment, nucleosome signal, blacklist ratio)
- Filter cells by these metrics
- Normalize with TF-IDF (term-frequency–inverse-document-frequency) — analogous to log-normalize but for binary peak presence
- Dimensionality reduction with LSI (latent semantic indexing) — analogous to PCA but on TF-IDF–normalized binary data
- Skip the first LSI component (it usually correlates with depth, not biology)
- UMAP on LSI components 2–30
- Cluster with Seurat’s standard graph-based methods
The conceptual shapes match RNA: get a QC’d, normalized, low-dim representation, cluster, visualize. The math is different (binary peaks not counts, TF-IDF not log-normalize, LSI not PCA) but the workflow scaffolding is the same.
Why TF-IDF, not log-normalization? ATAC data is binary or near-binary: each peak in each cell is either accessible (1 or a small count) or not (0). The data is extremely sparse — in a 100,000-peak × 10,000-cell matrix, more than 95% of entries are zero. Log-normalization, which divides by total counts and takes a log, is designed for count data with a roughly Poisson-to-negative-binomial distribution. For binary accessibility data, TF-IDF is more appropriate: it upweights rare peaks (high IDF — these are informative markers of specific cell types) and downweights ubiquitously-accessible peaks (low IDF — these are promoters active in all cells and carry little cell-type information). After TF-IDF, SVD (singular value decomposition) produces LSI components, which play the same role as PCs in RNA-seq.
FindTopFeatures(pbmc, min.cutoff = "q0") selects peaks above a minimum accessibility cutoff before running TF-IDF and SVD. Default with min.cutoff = "q0": keep all peaks (the 0th percentile, meaning no filtering). Setting min.cutoff = "q75" retains only the top 25% most accessible peaks, reducing noise from very rarely accessible peaks. For most analyses "q0" (all peaks) or "q5" (discard bottom 5%) is appropriate; raising the cutoff aggressively can discard peaks specific to rare cell types. The value "q0" is explicitly used in the tutorial to match the Signac vignette and is appropriate for the PBMC 10k dataset; for noisier data or datasets with many cell types consider "q5" or "q10".
After RunSVD(pbmc), downstream functions use dims = 2:30 rather than dims = 1:30. Default if unspecified: dims = 1:30, which includes the depth-confounded first LSI component. LSI component 1 almost universally tracks sequencing depth (total fragment count) rather than biology: DepthCor(pbmc) shows the Pearson correlation between each LSI component and nCount_ATAC, and component 1 typically has |r| > 0.9. Including it drives UMAP layout and cluster assignments by library size rather than chromatin state. Always run DepthCor and exclude any early components with |depth correlation| > 0.75.
CreateChromatinAssay(counts = counts, min.cells = 10, min.features = 200) filters peaks and cells at object creation. min.cells = 10 default: a peak must be accessible in at least 10 cells to be retained. Lowering this (e.g. to 5) retains more peaks and may recover rare-cell-type–specific elements, but increases the matrix size and noise. Raising it (e.g. to 20–50) removes rare peaks aggressively — appropriate for atlas-scale data where rare peaks dominate the sparse background. The min.features = 200 filter removes cells with fewer than 200 accessible peaks — analogous to nFeature_RNA > 200 in scRNA-seq QC.
15.5 Gene-activity scores — translating ATAC to RNA-like
Signac’s GeneActivity() translates ATAC into pseudo-RNA: for each gene, sum the fragment counts in a window around the TSS (e.g. ±2 kb upstream + the gene body). The result is a genes × cells matrix that looks like RNA expression but is derived from ATAC.
Why bother?
- It lets you use scRNA-derived markers to annotate scATAC clusters. A monocyte cluster in scATAC will have high gene activity for
CD14, LYZeven though those genes aren’t directly measured. - It enables label transfer from scRNA-seq (with cell types) to scATAC (without). The
FindTransferAnchorsworkflow projects an scATAC dataset against a labeled scRNA reference and copies the labels over.
Gene-activity scores are not as good as actual RNA expression — they’re an approximation — but they’re sufficient for cluster annotation and many downstream questions.
The window used by GeneActivity() matters: the default ±2 kb from the TSS captures the proximal promoter and a bit of the gene body. For genes regulated primarily by distal enhancers (common in immune genes), accessibility at the TSS may be low even when the gene is actively transcribed. Widening the window to ±5 kb or including a separate enhancer set can improve gene activity scores for such genes, though at the cost of noisier scores from overlap between adjacent gene windows. Always check that known markers for your expected cell types show appropriately elevated gene activity scores after computing them.
The dimensionality-reduction pipeline for scATAC-seq uses TF-IDF normalization followed by latent semantic indexing (LSI, implemented as SVD) rather than the log-normalization + PCA pipeline used for RNA. The reason is the nature of the data: ATAC peak counts are near-binary (each peak is either accessible or not in a given cell), severely sparse (>95% zeros), and lack the continuous count distribution that makes Poisson-based methods appropriate. TF-IDF (term frequency × inverse document frequency) is a normalization from information retrieval that weighs features (peaks) by their rarity across cells: a peak accessible in only 2% of cells (high IDF) is upweighted relative to a ubiquitously open promoter (low IDF). LSI / SVD then finds the major axes of variance in this TF-IDF–normalized binary matrix, analogous to how PCA operates on a scaled RNA matrix. The resulting “LSI components” play the same conceptual role as PCs and are used for UMAP, clustering, and label transfer. The practical warning is that LSI component 1 almost always captures sequencing depth rather than biology — it must be discarded before downstream analysis.
15.6 Peak calling vs pre-called peaks
The scATAC workflow assumes you have peaks. There are two ways:
Pre-called peaks (10x default): Cell Ranger ATAC calls peaks across the bulk-equivalent fragment pile and gives you that peak set. Fast; uses the same peaks for all cells.
De novo peak calling per cluster: cluster the cells first (using the pre-called peaks), then re-call peaks within each cluster, then re-cluster. This finds cell-type-specific peaks that the pan-cell-type call missed. Slower; produces a per-cluster peak set that you can use for downstream cell-type-specific motif analysis.
For a quick analysis the pre-called peaks are fine. For publication-quality cell-type-specific regulatory analysis, the de novo per-cluster recall is worth it.
15.7 Motif enrichment and TF activity
ChromVAR7 computes per-cell motif accessibility scores: for each known TF binding motif, how enriched is its presence in this cell’s peaks vs the genome-wide expectation? The score is per-cell, so you can plot it on UMAP, do DE on it, etc.
An alternative for motif enrichment is HOMER (findMotifsGenome.pl) — a command-line tool that takes a BED file of peaks (or a subset of peaks such as a cluster’s differential peaks) and reports known motif enrichment against a background. HOMER is not per-cell (it works on a peak set, not individual cells), making it better suited to the question “which TF motifs are enriched in this cluster’s peaks relative to all other peaks?” rather than the per-cell question chromVAR answers. Input to HOMER is a tab-separated BED-format file with chr / start / end / strand columns; the -size flag controls the scan window — use -size given to score each peak’s actual width, or a fixed window such as -size 200 centered on each peak summit (use one or the other, not both). The genome assembly specified to HOMER must match the alignment build of your peaks.
The ATAC pipeline often ends with:
- Cluster cells by chromatin accessibility (LSI + Leiden)
- Annotate clusters via gene activity / RNA label transfer
- Compute per-cluster differential motif accessibility — which TFs are most “active” in each cluster?
- Identify positive TF regulators by joining motif accessibility with gene expression (TFs whose motif is up and whose gene expression is up)
For tutorial-level worked examples of all of this, see scNotebooks Module 11 — the ArchR walkthrough on Kumegawa et al. 2022 breast-cancer data is more complete than our Signac one.
15.7b Clustering algorithm choice in scATAC
The tutorial uses FindClusters(algorithm = 3) — the SLM (Smart Local Moving) algorithm — rather than the default Louvain (algorithm = 1). SLM allows within-cluster re-optimization after each pass, producing more stable partitions on large, noisy graphs such as scATAC data. Leiden (algorithm = 4) is the modern default in many pipelines and generally preferred for its convergence guarantees. Either SLM or Leiden is a reasonable choice for scATAC; the key decision is to avoid the basic Louvain on noisy data where small parameter changes produce unstable clusterings.
15.7c Coordinate conversion: liftOver
If you need to bring a third-party peak set, motif reference, or annotation onto a different genome build than your data (for example, overlaying an hg38 ENCODE peak set onto an hg19-aligned experiment), use coordinate liftOver with rtracklayer and a UCSC chain file. Do not attempt to remap by gene name — coordinate offsets between hg19 and hg38 are nontrivial (insertions, deletions, chromosome rearrangements). The chain files (hg19ToHg38.over.chain, hg38ToHg19.over.chain, mouse equivalents mm9ToMm10, mm10ToMm39) are available from UCSC. liftOver() returns a GRangesList; always check lengths(result) == 1 to confirm each interval mapped to exactly one target interval before downstream use.
15.7d QC threshold design consideration
The thresholds in the Signac PBMC vignette (and echoed in the tutorial) are a reasonable starting point for healthy PBMC data but should not be applied mechanically to other tissues or datasets. An adaptive alternative is MAD-based outlier detection via scater::isOutlier(): flag cells whose metrics are more than 3 MADs below (or above, for nucleosome signal) the sample median. This approach scales automatically as overall library quality shifts and is less sensitive to the exact distribution than fixed thresholds.
15.8 hg19 vs hg38 — a non-trivial pitfall
The workshop’s PBMC 10k scATAC-seq dataset (used in Tutorial 15) is the v1 / hg19 release. The annotation package paired with it is EnsDb.Hsapiens.v75 (hg19).
The “Next GEM v1.1” re-release of the same sample is on hg38 with filenames prefixed atac_v1_pbmc_10k_nextgem_*. Mixing them produces bizarre coordinate errors that take hours to debug.
The rule: match the EnsDb annotation version to the genome the data is aligned to. v75 = hg19; v86 / v98 / v110 = various hg38 builds. Cell Ranger ATAC’s web_summary.html reports the genome.
15.9 Common errors
Table 2 collects the errors most commonly encountered in a Signac scATAC workflow, their causes, and the fix.
| You see | What’s wrong | What to do |
|---|---|---|
EnsDb.Hsapiens.v75 doesn’t match coordinates |
hg19 vs hg38 mismatch | Check the dataset’s reference genome; use the matching EnsDb |
fragments.tsv.gz is missing .tbi |
No tabix index | tabix fragments.tsv.gz to regenerate |
CreateChromatinAssay warns “duplicate fragments” |
Same fragment recorded under different barcodes | Usually safe (≤ 5%); high rates suggest sample-mixing |
| LSI component 1 dominates the UMAP | LSI1 correlates with sequencing depth, not biology | Skip dimension 1: use dims = 2:30 in RunUMAP |
| Gene activity scores are uniformly low | Promoter window too small | Default is ± 2 kb; widen to ±5 kb if your tissue has distal regulatory elements |
15.10 Where this material is also discussed
- Lecture: Lecture 15
- Tutorial: Tutorial 15
- Talapas raw-data processing: Tutorial Talapas RawData ATAC
- Alternative toolkit (ArchR): scNotebooks Module 11
15.11 Going further
The assay and its biology come from the ATAC-seq method paper1 and the chromatin-accessibility review2; single-cell ATAC was established by combinatorial-indexing atlases3. The analysis toolkits are Signac4 and ArchR5, with ChromVAR7 for TF-motif activity. To choose a protocol and pipeline, see the two benchmarks — of scATAC protocols8 and of computational methods6. For joint RNA+ATAC, see Cell Ranger ARC and Seurat’s WNN workflow9. The full curated list is on Key Papers, Reviews & Benchmarks.