pku-yuangroup 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>thenhelp(module.function)to check signatures - CLI:
<tool> --versionthen<tool> --helpto 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(formerlyFootprintScores) ->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
| Tool | Bias model | Scoring | Min depth | Strength | Fails when |
|---|---|---|---|---|---|
| TOBIAS (BINDetect) | +/-12 bp k-mer window (--k_flank 12), dinucleotide weight matrix (DWM) |
Two-step: continuous footprint score then motif-anchored bound/unbound classification | >= 50M nuclear reads | Mature, peer-reviewed (Bentsen 2020), differential support, modular pipeline | Below 50M reads; sequencing errors near motif inflate background |
| HINT-ATAC | Hidden-Markov + dinucleotide bias correction | HMM emits open/footprint/closed states; calls ranked footprints | >= 50M | Single-step; integrates motif matching; handles DNase too | Less control over individual stages; HMM occasionally over-segments |
| Wellington (pyDNase) | DNase-developed; ATAC adaptation by post-shift | Cleavage-rate Poisson Z-score | >= 50M (DNase >= 80M) | Original footprinting framework; well-validated for DNase | Designed for DNase II; ATAC-specific bias not corrected as carefully |
| PIQ | Bayesian latent variable on cut sites | Genome-wide PWM scan + cleavage profile | >= 30M (lower because of model) | Per-TF posterior probabilities; works on lower depth | Outdated; not actively maintained; harder to install |
| scprinter | Multi-scale CNN-based footprint and TF activity | Resolves footprints at multiple TF size scales (CTCF vs nuclear receptors) | >= 1M cells (sc) or 50M (bulk) | Modern ML approach; single-cell; multi-scale resolves problematic TF families | Newer tool; benchmarks evolving; GPU recommended |
| TOBIAS + scprinter combination | TOBIAS bias correction + scprinter scoring | Two-step bridging | >= 50M | Combines the best bias model with multi-scale scoring | Manual pipeline, no single CLI |
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:
| Strand | Read 5' end correction |
|---|---|
| + strand | shift +4 bp downstream |
| - strand | shift -5 bp upstream |
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
| Method | Approach | When to use |
|---|---|---|
| TOBIAS ATACorrect | +/-12 bp k-mer window, dinucleotide weight matrix (DWM) | Default for most ATAC; fast |
| chromBPNet bias model (Pampari 2024) | CNN trained on naked-DNA control or k-mer baseline | Best when sequence context complex; handles low-complexity flanks |
| seqOutBias (Martins 2018) | Genome-wide k-mer frequency scaling (observed vs expected cut counts) | Independent of footprinting tool; works upstream |
| HINT-ATAC dinucleotide | HMM-integrated dinucleotide bias | Built into HINT pipeline; less control |
| Naked-DNA empirical | Sequence Tn5 on protein-free DNA | Gold standard for non-model organisms; expensive wet-lab |
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.
Decision Tree by Goal
| Goal | Recommended pipeline |
|---|---|
| Identify all TFs differentially bound between two conditions | TOBIAS ATACorrect (per condition) -> ScoreBigwig -> BINDetect with --cond_names |
| Find the strongest single-TF binding (e.g., CTCF) | TOBIAS PlotAggregate over JASPAR CTCF motif sites; verify V-shape |
| Per-cell footprinting (scATAC) | scprinter (single-cell mode); avoid TOBIAS unless pseudobulking by cluster |
| Multi-scale TF activity (handle short and long simultaneously) | scprinter or TOBIAS + custom multi-scale |
| Differential nuclear-receptor binding | TOBIAS pooled-replicate footprints + ChIP cross-validation; raw ATAC alone often misses transient binding |
| Plant / non-model organism | TOBIAS or HINT-ATAC with custom motifs; bias model retrained from genomic background |
TOBIAS Three-Step Pipeline
Goal: Call bound/unbound TF motif sites per condition and detect differential occupancy across two conditions.
Approach: Run ATACorrect to subtract Tn5 bias from cleavage counts, ScoreBigwig to compute a continuous per-base footprint score, then BINDetect to anchor footprints to motif positions and produce per-TF differential bound calls with p-values.
# Step 1: Bias correction
TOBIAS ATACorrect \
--bam cond1.bam --genome hg38.fa \
--peaks consensus.bed --blacklist hg38-blacklist.v2.bed \
--outdir cond1_corrected/ --cores 16
# Step 2: Per-base footprint scoring (continuous)
TOBIAS ScoreBigwig \
--signal cond1_corrected/cond1_corrected.bw \
--regions consensus.bed \
--output cond1_footprints.bw \
--cores 16
# Step 3: Motif-anchored bound/unbound calls + differential
TOBIAS BINDetect \
--motifs JASPAR2024_CORE_vertebrates.pfm \
--signals cond1_footprints.bw cond2_footprints.bw \
--genome hg38.fa --peaks consensus.bed \
--outdir bindetect/ \
--cond_names cond1 cond2 \
--cores 16
BINDetect output columns: output_prefix, motif info, condition counts (cond1_bound, cond2_bound), cond1_mean_score, cond2_mean_score, cond1_cond2_change (differential), cond1_cond2_pvalue (one-sided per direction).
Differential Reading
| BINDetect output | Interpretation |
|---|---|
cond1_cond2_change > 0, low pvalue |
TF more bound in cond1 |
cond1_cond2_change < 0, low pvalue |
TF more bound in cond2 |
Both cond1_bound and cond2_bound near 0 |
Motif present but no footprint either condition; TF likely not active |
cond1_bound >> cond2_bound but change small |
High dynamic range; differential per-site rather than aggregate |
The differential score is the difference in mean footprint score across motif sites, not a fold-change. Magnitudes around 0.1-0.5 are typical for biologically relevant changes (TOBIAS BINDetect tutorials / Bentsen 2020 examples; no formally published cutoff -- calibrate against positive controls in the current dataset).
Reconciling TOBIAS vs HINT-ATAC
| Pattern | Likely cause | Action |
|---|---|---|
| TOBIAS calls binding, HINT does not | TOBIAS more sensitive; HINT's HMM filters edge calls | Trust if motif is canonical; suspect for novel/weak motifs |
| HINT calls binding, TOBIAS does not | HINT's HMM occasionally over-segments and reports spurious | Verify by aggregate footprint at the called sites |
| Both call same TF as differential but opposite directions | Different bias correction model; different bound/unbound thresholds | Re-check ATACorrect output; one bias model may be miscalibrated |
| Both flat | Library too shallow; chromatin too closed at motif sites | Pool replicates; consider scprinter multi-scale |
Operational rule: For high-confidence reporting, require two-tool concordance (TOBIAS + HINT-ATAC OR TOBIAS + ChIP-seq overlap > 50%). Single-tool calls should be reported as exploratory.
NFR-Filtering Before Footprinting
Goal: Restrict footprinting input to nucleosome-free (sub-100 bp) fragments where TF binding signal lives.
Approach: Stream the BAM through awk, keep header lines and fragments whose insert size is between -100 and 100 bp, then re-index.
# Filter to fragments < 100 bp (NFR) -- TF binding lives here, not on nucleosomes
samtools view -h sample.bam | \
awk 'substr($0,1,1)=="@" || ($9 > 0 && $9 < 100) || ($9 < 0 && $9 > -100)' | \
samtools view -b > sample.nfr.bam
samtools index sample.nfr.bam
# Use this NFR BAM as input to TOBIAS ATACorrect
Filtering NFR strengthens footprint signal but discards di-nucleosome-borne information. Keep the unfiltered BAM for nucleosome-positioning analysis.
Motif Database Choice
| Database | Coverage | Format | Notes |
|---|---|---|---|
| JASPAR 2024 CORE vertebrates | ~880 vertebrate non-redundant motifs (curated, experimentally derived) | JASPAR PFM, MEME, etc. | Default for vertebrate ATAC |
| HOCOMOCO v12 | ~1443 curated motifs (v12 CORE) | JASPAR PFM | Best for resolving paralogues; provides secondary motif subtypes per TF |
| CIS-BP 2.0 | ~80,000 motifs across 1000+ species | PWM, .meme | Broadest coverage including non-model species |
| MEME-CHIP / homer | Custom from peaks | .meme, .motif | When de novo motif needed |
JASPAR motifs are conservatively curated; HOCOMOCO is comprehensive for human/mouse with quality scores per motif (A/B/C/D); CIS-BP excels for non-model organisms.
Common Errors
| Error / symptom | Cause | Solution |
|---|---|---|
| Aggregate footprint inverted (peak instead of dip) | No bias correction; or wrong genome FASTA | Run ATACorrect; verify FASTA matches BAM build |
| BINDetect reports zero bound sites | Default cutoff too stringent; or peakset too narrow | Raise --bound-pvalue (default 0.001) and inspect; verify peaks include where binding expected |
| TOBIAS ATACorrect out of memory | Genome FASTA huge or many cores | Reduce --cores; use samtools faidx to confirm FASTA index exists |
| Differential score noisy / random | Per-condition bias correction inconsistent | Re-run ATACorrect with identical peakset and blacklist for each |
| Empty motif file warning | JASPAR PFM format mismatch | Use MEME suite to convert; TOBIAS expects JASPAR format |
| HINT-ATAC reports many tiny footprints | Default HMM over-segments | Use --organism flag explicitly; check --region-file is consensus, not raw peaks |
| Wellington crashes on paired-end ATAC | Wellington was DNase-targeted, single-end model | Use TOBIAS instead, or convert paired-end to cuts-only BED |
| Per-site footprint is V-shape but aggregate is flat | Mixing strands; some motifs on - strand | Aggregate function should handle strand; verify input motif strand column |
References
- Buenrostro JD et al 2013 Nat Methods 10:1213 (ATAC-seq protocol)
- Karabacak Calviello A et al 2019 Genome Biol 20:42 (protocol-specific Tn5/DNase bias modeling for footprinting)
- Bentsen M et al 2020 Nat Commun 11:4267 (TOBIAS framework, benchmark)
- Li Z et al 2019 Genome Biol 20:45 (HINT-ATAC)
- Piper J et al 2013 NAR 41:e201 (Wellington / pyDNase)
- Sherwood RI et al 2014 Nat Biotechnol 32:171 (PIQ)
- Hu Y et al 2025 Nature 638:779 (scPrinter/PRINT; multi-scale single-cell footprinting)
- Martins AL et al 2018 NAR 46:e9 (seqOutBias bias correction alternative)
- Castro-Mondragon JA et al 2022 NAR 50:D165 (JASPAR 2022, recently 2024)
- Vorontsov IE et al 2024 NAR 52:D154 (HOCOMOCO v12)
Related Skills
- atac-seq/atac-peak-calling - Generate input peakset (NFR-only optional)
- atac-seq/atac-qc - Confirm depth >= 50M before footprinting
- atac-seq/single-cell-atac - scprinter for single-cell footprinting
- atac-seq/motif-deviation - Complementary chromVAR for accessibility variability
- atac-seq/deep-learning-atac - chromBPNet bias correction alternative; in silico variant effect at footprints
- atac-seq/allele-specific-accessibility - Observed allelic imbalance at TF-bound sites
- chip-seq/peak-annotation - Cross-validate footprints with ChIP peaks
- sequence-manipulation/motif-search - Underlying motif scanning patterns
- gene-regulatory-networks/scenic-regulons - Downstream regulatory network inference