Bio Atac Seq Footprinting

Version Compatibility

Reference examples tested with: TOBIAS 0.16+, RGT HINT-ATAC 1.0.2+, Wellington (pyDNase) 0.3+, scprinter 0.1+, samtools 1.19+, bedtools 2.31+, pyBigWig 0.3+, MEME suite 5.5+.

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

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

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

TF Footprinting

"Identify TF binding footprints in my ATAC-seq data" -> Detect short DNA stretches (typically 6-20 bp) of reduced Tn5 cleavage within accessible regions, where a bound TF physically protects DNA. Requires (1) Tn5 sequence-bias correction, (2) per-base footprint scoring, (3) motif-anchored detection.

• CLI: TOBIAS ATACorrect -> TOBIAS ScoreBigwig (formerly FootprintScores) -> TOBIAS BINDetect
• CLI: rgt-hint footprinting --atac-seq (HINT-ATAC, single-step)
• CLI: wellington_footprints.py (legacy DNase, adapted for ATAC)
• Python: scprinter (multi-scale, single-cell aware; Hu 2025 Nature)

Tn5 has a strong sequence preference (Karabacak Calviello 2019), reading approximately +/- 4 bp around the insertion site. Without bias correction, "footprints" reflect Tn5 sequence preference rather than TF binding. This is the single most important step.

Algorithmic Taxonomy

Methodology evolves; verify against the current Bentsen 2020, Karabacak Calviello 2019, and scPrinter (Hu 2025) benchmarks. ATAC footprinting power saturates above 100M nuclear reads; below 50M, weak-binding TFs (transient occupancy) cannot be reliably called.

Tn5 Bias and Why Correction Matters

Trigger: Tn5 inserts preferentially at certain k-mers (Karabacak Calviello 2019 measured the protocol-specific 6-mer insertion-bias model used for bias correction). The preference is reproducible and biologically uninteresting.

Mechanism: Without correction, every "TF footprint" near a high-bias k-mer reads as occupancy. Conversely, regions with low-bias flanks but real binding may show no footprint dip relative to corrected expectation.

Symptom: Aggregate footprint at random GC-rich motifs shows V-shape; aggregate at AT-rich motifs shows inverse-V (peak instead of dip).

Fix: Apply ATACorrect (TOBIAS), seqOutBias (Martins 2018), or HINT's dinucleotide model. Bias correction subtracts the Tn5-expected per-base profile from observed cleavage. After correction, V-shape is preserved only at TF-bound sites.

Goal: Subtract Tn5 sequence-bias from the per-base cleavage signal so residual footprints reflect TF binding rather than enzyme preference.

Approach: Run TOBIAS ATACorrect over the deduplicated BAM with the reference genome, consensus peaks, and ENCODE blacklist; it emits per-condition uncorrected, bias, expected, and corrected bigWigs for downstream scoring.

# TOBIAS ATACorrect: produces uncorrected, bias, expected, and corrected bigWigs
TOBIAS ATACorrect \
    --bam sample.dedup.bam \
    --genome hg38.fa \
    --peaks consensus_peaks.bed \
    --blacklist hg38-blacklist.v2.bed \
    --outdir corrected/ \
    --cores 8
# Output: sample_uncorrected.bw, sample_bias.bw, sample_expected.bw, sample_corrected.bw

Tn5 Cut Geometry: +4 / -5 Dual-Cut

Tn5 dimers cut both strands of DNA but with a 9 bp staggered offset. The cleavage event creates two free 5' ends: one shifted +4 bp from the binding center on the forward strand and -5 bp on the reverse strand. Footprinting tools must apply this shift before per-base counting:

Trigger: Computing per-base Tn5 cut signal manually.

Symptom: Footprint aggregates show ~9 bp asymmetry (apex shifted from motif center).

Fix: Apply +4/-5 shift before counting; TOBIAS, HINT-ATAC, and scprinter handle this internally. Custom analyses must apply explicitly. deepTools provides this via alignmentSieve --ATACshift (which applies the canonical +4 / -5 shift in one step).

Bias Correction Alternatives

For non-model organisms (no published Tn5 bias model), naked-DNA control is required. chromBPNet's bias model is the modern standard for human/mouse and outperforms TOBIAS at low-complexity sequence contexts (Pampari 2024). See atac-seq/deep-learning-atac.

In Silico Variant Effect at Footprinted TF Motifs

Trigger: A GWAS-fine-mapped or rare variant falls inside a TOBIAS-bound motif site.

Mechanism: Sequence-based DL models (chromBPNet, Enformer) predict per-base accessibility at ref vs alt allele; combined with footprint evidence (TOBIAS bound site overlap), this produces a mechanistic hypothesis: "variant disrupts binding of TF X at enhancer Y."

Workflow: Run TOBIAS BINDetect to identify bound motif sites; for variants in bound sites, score with chromBPNet (atac-seq/deep-learning-atac) for ref/alt log2FC; |log2FC| > 1 supports functional disruption. Cross-reference with allele-specific accessibility (atac-seq/allele-specific-accessibility) for observed evidence.

Per-TF Footprinting Failure Modes

Different TF families produce different footprint signatures. The same tool can report a clean V-shape for one family and noise for another.

CTCF -- The gold standard

Trigger: Footprinting CTCF.

Mechanism: CTCF has high ChIP-seq concordance, deep V-shaped footprint (~19-20 bp protected), strong sequence specificity. ChIP-seq overlap is typically >70%.

Symptom: Aggregate corrected footprint shows clean ~20 bp dip with bilateral cleavage shoulders.

Verification: Always validate footprinting output by checking CTCF first; if CTCF footprint is shallow, the bias correction or depth is the problem, not the biology.

Nuclear receptors (ER, AR, GR) -- Transient binding

Trigger: Glucocorticoid response, hormone-stimulated systems.

Mechanism: Steroid receptors bind transiently (residence time minutes vs hours for CTCF); average ATAC sample captures binding probability < 30% per allele.

Symptom: Aggregate footprint is shallow or absent despite ChIP-seq peaks at the same sites.

Fix: Use scprinter's multi-scale model OR limit to ChIP-validated sites OR pool replicates for higher effective depth. Do not interpret absence of footprint as absence of binding.

Pioneer TFs (FOXA1, GATA, OCT4) -- Half-site footprint

Trigger: Pioneer-factor binding to nucleosomal DNA.

Mechanism: Pioneer factors bind one DNA face; the back face is on the histone octamer. Footprint is asymmetric (one side protected, other side accessible).

Symptom: Aggregate plot shows asymmetric V; one shoulder is taller than the other.

Fix: Use single-stranded scoring; HINT-ATAC has stranded mode. Treat asymmetric footprints as biologically meaningful, not artefactual.

AP-1 family (FOS, JUN) -- Heterodimer composite footprint

Trigger: AP-1 enrichment; the JASPAR motif is composite of multiple heterodimer combinations.

Mechanism: Different AP-1 dimers (FOS+JUN, FOS+JUNB, JUNB+JUNB) bind slightly different motifs. JASPAR entries are degenerate; scoring averages over all.

Fix: Use specific HOCOMOCO motifs per heterodimer when distinguishing matters. Otherwise accept the composite call.

ZBTB family / BTB-zinc finger -- Dynamic / unfootprintable

Trigger: Footprinting ZBTB16, BCL6, others.

Mechanism: Dynamic binding kinetics + cofactor-mediated stabilization mean steady-state occupancy is highly variable. Some ZBTBs simply do not produce reliable ATAC footprints despite genuine ChIP-seq binding.

Fix: Document the failure; use ChIP-seq for these TFs. Footprinting cannot rescue everything.

Forkhead / homeobox (FOX, HOX) -- Short footprint < 8 bp

Trigger: Short-motif TFs.

Mechanism: Footprint extent matches motif length; <8 bp footprints are at the resolution limit of Tn5 (which has ~4 bp positional uncertainty).

Fix: Multi-scale scoring (scprinter); aggregate over thousands of sites; do not rely on per-site calls for short motifs.