Bio Atac Seq Motif Deviation

Version Compatibility

Reference examples tested with: chromVAR 1.24+, motifmatchr 1.24+, JASPAR2024 0.99+, TFBSTools 1.40+, BSgenome.Hsapiens.UCSC.hg38 1.4+, SummarizedExperiment 1.32+, limma 3.58+, ggplot2 3.5+, Matrix 1.6+, ArchR 1.0.2+, Signac 1.13+.

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

• R: packageVersion('<pkg>') then ?function_name to verify parameters

If code throws unexpected errors, introspect the installed package and adapt rather than retrying.

Motif Deviation (chromVAR)

"Which TF motifs explain accessibility variation across my samples or cells?" -> Compute per-sample (or per-cell) deviation z-scores: how many standard deviations above expectation each TF motif's accessibility falls, controlling for GC content and overall accessibility via matched background peak sets.

• R: chromVAR::computeDeviations(counts, motifs) -> per-sample z-scores
• R: chromVAR::computeVariability(dev) -> per-motif variance ranking
• Single-cell alternative: Signac::RunChromVAR() (wrapper with matched defaults) or ArchR::addDeviationsMatrix()

chromVAR answers a different question than footprinting: footprinting asks "is this specific motif site bound?", chromVAR asks "do peaks containing this motif have systematically more or less accessibility than expected?" The two are complementary.

What chromVAR Computes

For each (motif, sample) pair:

• Raw deviation = Sum of accessibility at peaks containing the motif - expected from a matched-GC, matched-accessibility background.
• Bias-corrected deviation = Raw deviation / SD of background deviations.
• Z-score = (corrected deviation - mean across cells) / SD across cells. Reported as the principal output.

Z-scores are signed: positive = motif more accessible in this sample than population average; negative = less. Magnitudes 2-5 are typical for biologically interesting motifs; >5 indicates strong covariation with sample state.

Algorithmic Taxonomy

Methodology evolves; verify against current chromVAR (Schep 2017), ArchR (Granja 2021), and Signac (Stuart 2021) benchmarks before locking pipelines.

chromVAR vs Footprinting -- Different Questions

chromVAR is fundamentally a summary statistic over many motif sites. Footprinting is per-site classification. Use chromVAR when motif site count >> 100; use footprinting when specific sites matter.

Per-Tool Failure Modes

chromVAR -- Too few peaks or too few reads

Trigger: ATAC peakset < 5000 peaks; per-sample read depth < 1500 in peaks.

Mechanism: chromVAR's background sampling requires enough peaks to find matched GC + accessibility partners. Sparse sampling at low peak count creates correlated null distributions, inflating both positive and negative z-scores.

Symptom: Variability scores all > 5 (suspiciously high); top variable motifs are dominated by AT-rich or GC-rich sequences regardless of biology.

Fix: Verify peakset is at full ATAC scale (typically 50k-200k peaks). For sc ATAC, aggregate cells to clusters of >= 500 cells before running.

chromVAR background peaks -- Default is good, custom requires care

Trigger: Calling getBackgroundPeaks() with non-default niterations or bias.

Mechanism: Default niterations=50 yields 50 matched background peaks per foreground peak. Reducing niterations increases noise; increasing slows linearly without much accuracy gain.

Symptom: Custom backgrounds inflate variability when niterations < 30.

Fix: Stick to defaults unless benchmarking. If running on huge cell counts, test on subsample first.

chromVAR on broadly accessible cell types -- Z-scores compressed

Trigger: Multiple cell types in dataset have very different overall accessibility magnitudes.

Mechanism: chromVAR's correction normalizes for total accessibility; cell types with high background accessibility have compressed z-scores even if their motif-specific signal is strong.

Symptom: PCA on z-scores does not separate cell types as cleanly as raw counts.

Fix: Run chromVAR per-cell-type-cluster (separate runs) when global accessibility differs by > 5x. Alternatively use ArchR's per-cluster background.

chromVAR on bulk samples without enough variation -- All z-scores near zero

Trigger: All bulk samples are technical replicates or very similar.

Mechanism: Z-scores normalize across the sample population; if there is no across-sample variability, all z-scores collapse to zero.

Symptom: Variability ranking is unstable across runs; top motifs change.

Fix: chromVAR is designed for variability; if the dataset has only one biological condition replicated, use footprinting or differential accessibility instead. chromVAR needs 6+ samples with biological variation to be informative.

Signac::RunChromVAR -- Motif matching mismatch

Trigger: Motif assay added before peak set finalized; peak coordinates change.

Mechanism: RunChromVAR matches motifs to peaks at the time it's called; if peaks change downstream (e.g., after merge), the motif annotations become stale.

Symptom: Some peaks have NA motif annotations; deviation matrix has missing entries.

Fix: Run AddMotifs() -> RunChromVAR() AFTER finalizing peakset. Re-run if peaks change.

ArchR::addDeviationsMatrix -- TileMatrix vs PeakMatrix

Trigger: Calling on tile matrix when peak matrix is more appropriate.

Mechanism: ArchR can compute deviations on either tiles (regular bins) or peaks. Peaks are biologically meaningful; tiles add noise from intergenic background.

Fix: Use matrixName='PeakMatrix' after addReproduciblePeakSet. Tile-based deviations are mainly for embedding, not biology.

Decision Tree by Setting

chromVAR Workflow (Bulk)

Goal: Compute per-sample TF-motif accessibility z-scores corrected for GC bias and total signal.

Approach: Build a SummarizedExperiment from peak counts, add GC bias, filter sparse samples and peaks, match JASPAR motifs to peaks, sample matched background peaks, then compute deviations and per-motif variability.

library(chromVAR); library(motifmatchr); library(BSgenome.Hsapiens.UCSC.hg38)
library(JASPAR2024); library(TFBSTools); library(SummarizedExperiment)

peaks <- rtracklayer::import('consensus_peaks.bed')           # GRanges
counts <- as.matrix(read.delim('peak_counts.tsv', row.names=1))    # rows = peaks, cols = samples
se <- SummarizedExperiment(assays=list(counts=counts), rowRanges=peaks)
se <- addGCBias(se, genome=BSgenome.Hsapiens.UCSC.hg38)

# Filter: depth >= 1500 reads/sample, FRiP >= 0.15; drop peaks with < 10 total fragments
se <- filterSamples(se, min_depth=1500, min_in_peaks=0.15, shiny=FALSE)
se <- filterPeaks(se, non_overlapping=TRUE, min_fragments_per_peak=10)

# Motifs: JASPAR vertebrate CORE (default for human/mouse)
# JASPAR2024 + TFBSTools incompatibility (TFBSTools issue #39): getMatrixSet does not dispatch on the
# JASPAR2024 object directly. Open the SQLite handle and pass that to getMatrixSet instead.
library(RSQLite)
jaspar2024 <- JASPAR2024::JASPAR2024()
sq <- dbConnect(SQLite(), db(jaspar2024))
pfm <- getMatrixSet(sq, opts=list(collection='CORE', tax_group='vertebrates'))
# JASPAR2020 (older) accepts the package object directly: getMatrixSet(JASPAR2020, opts=...)
motif_ix <- matchMotifs(pfm, se, genome=BSgenome.Hsapiens.UCSC.hg38, p.cutoff=5e-05)

# Background peaks: matched GC + accessibility (default 50 iterations is fine)
bg <- getBackgroundPeaks(object=se, niterations=50)

# Compute deviations
dev <- computeDeviations(object=se, annotations=motif_ix, background_peaks=bg)
zscores <- deviationScores(dev)         # motif x sample matrix of z-scores (deviations() returns raw bias-corrected deviations)
variability <- computeVariability(dev) # per-motif variability ranking

Differential Motif Activity (limma on z-scores)

Goal: Identify TF motifs whose chromVAR z-scores differ significantly between conditions.

Approach: Build a contrast design matrix, fit limma's linear model on the motif-x-sample z-score matrix with empirical Bayes moderation, and pull motifs at adjusted p < 0.05.

library(limma)
groups <- factor(colData(se)$condition, levels=c('control', 'treated'))
design <- model.matrix(~groups)
fit <- lmFit(zscores, design); fit <- eBayes(fit)
diff_motifs <- topTable(fit, coef=2, number=Inf, p.value=0.05)   # adj.P.Val column

Use adj.P.Val (limma's BH FDR), not FDR (which limma does not return). logFC is the difference in z-scores between groups; magnitudes 0.5-2 typical.

chromVAR for Single-Cell ATAC (Signac)

Goal: Compute per-cell TF-motif z-scores in a Seurat scATAC workflow and call cluster-marker motifs.

Approach: Open the JASPAR2024 SQLite handle, attach motifs to the Seurat object via AddMotifs, run RunChromVAR to build the chromvar assay, then call FindAllMarkers with mean.fxn=rowMeans for z-score-appropriate differential.

library(Signac); library(Seurat); library(JASPAR2024); library(TFBSTools)
library(BSgenome.Hsapiens.UCSC.hg38); library(RSQLite)

# Assume `seurat_obj` has an ATAC assay with consensus peaks
# JASPAR2024 + TFBSTools workaround (see TFBSTools issue #39):
jaspar2024 <- JASPAR2024::JASPAR2024()
sq <- dbConnect(SQLite(), db(jaspar2024))
pfm <- getMatrixSet(sq, opts=list(collection='CORE', tax_group='vertebrates'))
seurat_obj <- AddMotifs(seurat_obj, genome=BSgenome.Hsapiens.UCSC.hg38, pfm=pfm)
seurat_obj <- RunChromVAR(seurat_obj,
                          genome=BSgenome.Hsapiens.UCSC.hg38,
                          new.assay.name='chromvar')
DefaultAssay(seurat_obj) <- 'chromvar'

# Per-cluster differential motifs.
# `mean.fxn` is the standard FindAllMarkers/FindMarkers control for the per-feature summary.
# `fc.name` controls the output column name and is accepted by Seurat 4.x/5.x; if it errors,
# fall back to renaming the output column post-hoc.
markers <- FindAllMarkers(seurat_obj, only.pos=TRUE, mean.fxn=rowMeans, fc.name='avg_diff')

mean.fxn=rowMeans is required for z-score-style data; the default fold-change function (designed for log-counts) makes no sense on chromVAR z-scores.