Where this lecture fits

• Previous: Lec 08 — Differential Abundance
• You are here: Lec 15 — measuring chromatin accessibility, in bulk and per cell
• Next: Lec 13 — WGCNA
• Companion tutorial: Tutorial 15 — scATAC-seq with Signac
• Companion book chapter: Chapter 15 — Bulk and Single-Cell ATAC-seq
• Alternative toolkit: scNotebooks Module 11 walks the same workflow with ArchR instead of Signac, on Kumegawa et al. 2022 breast-cancer data — useful for breadth

Goals of this lecture

• Understand what ATAC-seq measures and how it complements RNA-seq
• See how the bulk ATAC-seq workflow turns reads into peaks
• Translate the bulk workflow to the single-cell setting (scATAC-seq, Signac)
• Read scATAC QC metrics — fragment counts, TSS enrichment, nucleosome signal
• Connect chromatin clusters to RNA-defined cell types via gene activity

The biology in one slide

• DNA in the nucleus is wrapped around histones into nucleosomes
• Regions where transcription factors and RNA polymerase need access (promoters, enhancers, insulators) are nucleosome-depleted → “open chromatin”
• ATAC-seq uses a hyperactive Tn5 transposase that simultaneously cuts open DNA and tags it with sequencing adapters — only accessible regions are cut
• Map reads back to the genome → peaks of read coverage = open regions

ATAC-seq — overview

ATAC-seq overview — chromatin accessibility from Tn5 to peaks

• The whole assay in one picture: open chromatin → Tn5 cuts → tagged fragments → sequenced reads → peaks
• Same data structure as ChIP-seq downstream, very different sample prep upstream

ATAC vs. other accessibility assays

ATAC tools comparison

ATAC tools comparison — bulk vs single-cell stacks

• The processing pipelines for bulk and scATAC share the same input/output shape but use different tooling at each step
• Where Cell Ranger ATAC, Signac, ArchR, and SnapATAC2 sit relative to each other

From reads to peaks

1. Sequence paired-end reads (50–75 bp), trim adapters
2. Align with bwa-mem or bowtie2 (--very-sensitive)
3. Filter: drop duplicates, mitochondrial reads (chrM is hugely over-represented), low-MAPQ reads
4. Shift Tn5 offsets — +4/-5 bp on the two strands (Tn5 binds as a dimer)
5. Call peaks with MACS2 callpeak --nomodel --shift -75 --extsize 150 -B --SPMR (or -f BAMPE for paired-end)
6. Annotate peaks with ChIPseeker or HOMER against a TSS database

Common QC metrics in bulk ATAC

scATAC vs. scRNA — same shape, different math

The 10x Multiome / scATAC fragment file

The primary input for Signac is a fragments.tsv.gz — every line is one Tn5 cut event:

chr1   10125   10570   AAACGAAAGAGGAGCC-1   2

chrom start end cell_barcode duplicate_count

A tabix index (.tbi) is built so Signac can pull fragments for any genomic region in O(log n).

Tutorial dataset — PBMC 10k v1, hg19

scATAC QC metrics — what to look at

# Standard Signac QC subset
pbmc <- subset(pbmc,
               subset = peak_region_fragments > 3000 &
                        peak_region_fragments < 20000 &
                        pct_reads_in_peaks    > 15 &
                        blacklist_ratio       < 0.05 &
                        nucleosome_signal     < 4 &
                        TSS.enrichment        > 2)

The Signac workflow

library(Signac); library(Seurat); library(EnsDb.Hsapiens.v75)

counts <- Read10X_h5("data/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5")
chrom_assay <- CreateChromatinAssay(
  counts    = counts,
  sep       = c(":", "-"),
  fragments = "data/atac_v1_pbmc_10k_fragments.tsv.gz",
  min.cells = 10, min.features = 200)

pbmc <- CreateSeuratObject(counts = chrom_assay, assay = "ATAC",
                           meta.data = read.csv("data/atac_v1_pbmc_10k_singlecell.csv",
                                                row.names = 1))

# ATAC-specific normalization + reduction
pbmc <- RunTFIDF(pbmc)
pbmc <- FindTopFeatures(pbmc, min.cutoff = "q5")
pbmc <- RunSVD(pbmc)

# Skip LSI component 1 — usually correlates with depth, not biology
pbmc <- RunUMAP(pbmc, reduction = "lsi", dims = 2:30)
pbmc <- FindNeighbors(pbmc, reduction = "lsi", dims = 2:30) |>
        FindClusters(algorithm = 3)

Gene-activity scores

To label scATAC clusters using familiar gene markers, build a gene activity matrix — sum fragments overlapping each gene body + 2 kb upstream:

gene.activities <- GeneActivity(pbmc)
pbmc[["RNA"]] <- CreateAssayObject(counts = gene.activities)
pbmc <- NormalizeData(pbmc, assay = "RNA",
                      normalization.method = "LogNormalize",
                      scale.factor = median(pbmc$nCount_RNA))

DefaultAssay(pbmc) <- "RNA"
FeaturePlot(pbmc, c("MS4A1","CD3D","LEF1","NKG7","TREM1","LYZ"),
            max.cutoff = "q95", ncol = 3)

Label transfer from scRNA to scATAC

If you also have a matched scRNA dataset (10x Multiome, or a published PBMC reference such as Signac’s pbmc_10k_v3.rds), use anchor-based label transfer:

pbmc_rna <- readRDS("data/pbmc_10k_v3.rds")
anchors <- FindTransferAnchors(reference = pbmc_rna, query = pbmc,
                               reduction = "cca")
pred    <- TransferData(anchorset = anchors, refdata = pbmc_rna$celltype,
                        weight.reduction = pbmc[["lsi"]], dims = 2:30)
pbmc <- AddMetaData(pbmc, pred)

Downstream you can reach for

• Differential accessibility — FindMarkers(test.use = "LR", latent.vars = "nCount_peaks")
• Motif enrichment — Signac::FindMotifs against JASPAR2020, or HOMER (findMotifsGenome.pl) from the command line
• Per-cell TF activity — chromVAR deviations
• Coverage plots — CoveragePlot() to visualize fragment pileups by cell type

Design considerations

• QC thresholds are dataset-specific. The thresholds in the QC code above (peak_region_fragments > 3000, TSS.enrichment > 2, etc.) are starting points — always inspect the violin plots and density scatters first before cutting.
• min.cutoff = "q0" vs "q5" in FindTopFeatures — "q0" keeps all peaks; "q5" drops the lowest-frequency 5%. Use "q0" if you have few cells or rare cell types.
• SLM (algorithm 3) vs Louvain/Leiden — SLM (algorithm 3) often gives more stable partitions on noisy scATAC graphs; Leiden (algorithm 4) is the modern default for RNA but either is fine.
• Gene activity is a proxy — open chromatin near a gene loosely correlates with expression; a gene can be open without being transcribed (poised) or transcribed via a distal enhancer. Use gene activity for labelling, not expression quantification.

What to remember from Lecture 15

• ATAC-seq measures chromatin accessibility via Tn5 tagmentation
• Bulk ATAC-seq goes reads → peaks → annotation; scATAC-seq adds a barcode layer per Tn5 event
• scATAC math is LSI, not PCA — and the first component is usually depth, drop it
• Per-cell QC has its own metrics: TSS enrichment, nucleosome signal, FRiP
• Gene activity scores + label transfer let you reuse RNA cell-type labels in chromatin space

Coming up next

• Lec 13 — WGCNA: co-expression networks for bulk and pseudobulk data
• Hands-on: Tutorial 15 — scATAC-seq with Signac

Further reading

• Signac vignette — PBMC scATAC
• ArchR book — alternative scATAC stack, fast on large datasets
• ENCODE ATAC-seq pipeline standards — bulk QC thresholds