Bio Applied Scatac Chromatin

Single-Cell ATAC-seq: Chromatin Accessibility

References: Signac | SnapATAC2 | ArchR | chromVAR

scATAC-seq measures chromatin accessibility at single-cell resolution — a direct readout of where transcription factors bind and regulatory decisions are made, upstream of expression changes.

When to Use

• Processing a 10x Genomics scATAC-seq or multiome fragments.tsv.gz file into a cell-by-peak matrix.
• Normalizing sparse binary accessibility data (library-size normalization doesn't apply) and clustering cells.
• Diagnosing why PC1/LSI1 correlates with sequencing depth instead of biology.
• Linking distal peaks to genes (co-accessibility, peak-gene links) or scoring TF motif activity (chromVAR-style).
• QC'ing an ATAC library with TSS enrichment score or nucleosome fragment-length banding.

Version Compatibility

Python ≥3.10, SnapATAC2 ≥2.6, anndata ≥0.10, scikit-learn ≥1.3, umap-learn ≥0.5. R ≥4.3, Signac ≥1.13, Seurat ≥5.0. Cicero (R/monocle3) for co-accessibility. chromVAR ≥1.24 (Bioconductor) for motif deviations.

Prerequisites

pip install snapatac2 anndata scikit-learn umap-learn scipy pandas matplotlib or in R install.packages("Signac"); BiocManager::install("chromVAR"). Assumes familiarity with bio-single-cell-preprocessing (AnnData basics) and bio-genome-intervals-bed-file-basics (BED coordinates). A pseudo-bulk peak set (e.g. from MACS3) is a prerequisite input for the peak-by-cell matrix.

Technology

Tn5 tagmentation: Hyperactive Tn5 inserts into nucleosome-free (accessible) regions, fragmenting DNA and ligating sequencing adapters. For 10x: isolate nuclei first (removes mtDNA), then add Tn5 inside droplets.

Fragment file format (TSV.gz, Tabix-indexed):

chr1    10000    10200    ACGTCAGTACGT-1    1

Columns: chrom, start, end, cell barcode, read count. Compact and cell-indexed — faster than BAM for per-cell operations.

QC metrics:

1. TSS enrichment score — ratio of coverage at TSS vs. flanking background. Score > 4 acceptable; > 8 high quality.
2. Nucleosome banding — fragment length histogram shows mono- (147 bp), di- (294 bp), tri-nucleosome peaks + sub-nucleosomal (<147 bp) free DNA.
3. Fragments per cell — 10,000–50,000 unique fragments typical for 10x.

Peak Calling and TF-IDF Normalization

Why not call peaks per cell: each cell has only 0–1 reads per position — insufficient signal. Solution: pseudo-bulk peak calling (cluster cells → merge fragments per cluster → run MACS3 → build a master peak set).

TF-IDF handles sparse binary ATAC data where library-size normalization is meaningless:

• TF = fragments in peak / total fragments in cell (depth normalization)
• IDF = log(1 + n_cells / cells_with_peak_accessible) (downweights ubiquitous peaks)
• Final = TF × IDF — emphasizes cell-type-specific peaks

SnapATAC2 applies TF-IDF to raw fragment counts; ArchR binarizes (0/1) by default. Both work; binary is more robust to PCR duplicates.

Goal: build a cell-by-peak matrix and cluster cells by chromatin accessibility. Approach: load fragments, bin the genome into tiles, select variable features, run TF-IDF + spectral (LSI) embedding, then UMAP/Leiden.

import snapatac2 as snap

def cluster_scatac(fragment_file: str, bin_size: int = 500, n_comps: int = 30):
    """Build a tile matrix from a fragments.tsv.gz, normalize, and cluster.

    Returns an AnnData with obsm['X_spectral'] (LSI) and obs['leiden'].
    """
    data = snap.read(fragment_file, backed="r")
    snap.pp.make_tile_matrix(data, bin_size=bin_size)
    snap.pp.select_features(data, n_features=50000)
    snap.tl.spectral(data, n_comps=n_comps)  # TF-IDF + SVD (LSI) internally
    snap.pp.knn(data)
    snap.tl.leiden(data)
    snap.tl.umap(data)
    return data

Signac (R)

library(Signac); library(Seurat)

chrom_assay <- CreateChromatinAssay(
  counts = peak_counts, fragments = "fragments.tsv.gz",
  genome = "hg38", min.cells = 10
)
seurat_obj <- CreateSeuratObject(counts = chrom_assay, assay = "ATAC")
seurat_obj <- RunTFIDF(seurat_obj)
seurat_obj <- FindTopFeatures(seurat_obj, min.cutoff = "q75")
seurat_obj <- RunSVD(seurat_obj)
DepthCor(seurat_obj)  # LSI_1 usually shows high depth correlation -> drop it
seurat_obj <- RunUMAP(seurat_obj, reduction = "lsi", dims = 2:30)

LSI Dimensionality Reduction

Why PCA fails: scATAC matrices are binary/sparse; PC1 correlates with read depth (technical artifact), not biology.

LSI (SVD on the TF-IDF matrix): compute TF-IDF, apply truncated SVD, discard component 1 (almost always depth-correlated), then use components 2–30 for UMAP/clustering.

Goal: confirm LSI1 is a depth artifact and remove it before downstream analysis. Approach: correlate each LSI component against total fragment count per cell; drop any component with |r| > 0.5 (normally just component 1).

import numpy as np
from sklearn.decomposition import TruncatedSVD
from scipy.stats import pearsonr


def tfidf_lsi(counts: np.ndarray, n_components: int = 30):
    """TF-IDF normalize a binary cell-by-peak matrix and run LSI (SVD).

    counts: (n_cells, n_peaks) binary/count matrix.
    Returns (X_lsi, depth_corr) — depth_corr[i] is Pearson r of LSI comp i
    with per-cell fragment depth; drop components with |r| > 0.5.
    """
    n_cells = counts.shape[0]
    cell_totals = counts.sum(axis=1, keepdims=True) + 1e-6
    tf = counts / cell_totals
    cells_with_peak = (counts > 0).sum(axis=0) + 1
    idf = np.log1p(n_cells / cells_with_peak)
    tfidf = tf * idf[np.newaxis, :]

    svd = TruncatedSVD(n_components=n_components, random_state=42)
    x_lsi = svd.fit_transform(tfidf)

    depth = counts.sum(axis=1)
    depth_corr = np.array([pearsonr(x_lsi[:, i], depth)[0] for i in range(n_components)])
    return x_lsi, depth_corr


# Example: drop components correlated with depth, keep the rest
# x_lsi, corr = tfidf_lsi(counts_matrix)
# keep = np.where(np.abs(corr) <= 0.5)[0]
# x_lsi_clean = x_lsi[:, keep]

Co-accessibility and Peak-Gene Links

Cicero (Pliner et al. 2018): trains graphical LASSO on pseudo-cells to compute co-accessibility scores (0–1) within 500 kb. Score > 0.25 suggests a regulatory connection.

Peak-gene linking (Signac)

seurat_obj <- LinkPeaks(
  object = seurat_obj,
  peak.assay = "ATAC",
  expression.assay = "RNA",
  distance = 5e5
)

Validation: Hi-C/HiChIP (3D proximity), eQTL overlap, ChromHMM enhancer states, ENCODE cCREs.

Limitations: co-accessibility requires >500 cells; correlation ≠ causation; cell-type mixing creates spurious co-accessibility.

Pitfalls

• Coordinate systems: BED uses 0-based half-open; VCF/GFF use 1-based inclusive — mixing causes off-by-one errors.
• Depth confound: never use raw LSI1/PC1 for clustering without checking DepthCor/depth correlation first.
• Batch effects: always check for batch confounding (e.g. by sample/lane) before interpreting biological signal.
• Multiple testing: apply FDR correction (Benjamini-Hochberg) when testing thousands of peaks or peak-gene links simultaneously.

See Also

• bio-single-cell-scatac-analysis — broader single-cell ATAC preprocessing workflow.
• bio-atac-seq-atac-peak-calling — MACS3 pseudo-bulk peak calling.
• bio-atac-seq-motif-deviation — chromVAR TF motif deviation scores.
• bio-workflows-multiome-pipeline — joint RNA+ATAC multiome pipeline.