# Bio Applied Tf Footprinting

> Detect TF footprints in ATAC-seq via Tn5 offset correction, insertion-profile aggregation, footprint scoring, and pybedtools intersect/slop/closest. Use when doing TF footprinting, ATAC-seq Tn5 bias correction, footprint scoring, or motif-site meta-profile plots.

- Skill: `pavel-kravchenko/bio-applied-tf-footprinting` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-tf-footprinting`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-tf-footprinting/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/bio-applied-tf-footprinting

---


# TF Footprinting & Chromatin Accessibility

## When to Use

- Determining whether a specific TF (e.g. CTCF) is actively bound at accessible chromatin sites, not just "open" but occupied.
- Computing Tn5 insertion profiles around motif centers from ATAC-seq BAMs and scoring the depth of the central footprint.
- Intersecting ATAC peaks, TF motif calls, blacklist regions, and TSS annotations with pybedtools before footprinting.
- Building accumulation (meta-profile) plots comparing high- vs low-affinity binding sites, or bound vs unbound sites.
- Setting up inputs for production footprinting tools (TOBIAS, HINT-ATAC) or interpreting their output.

## Version Compatibility

- Python ≥3.10, NumPy ≥1.24, SciPy ≥1.11, pandas ≥2.0
- pybedtools ≥0.9 (requires BEDTools ≥2.30 in `PATH`)
- pysam ≥0.22 for BAM parsing
- TOBIAS ≥0.16 (optional, for production bias-corrected footprinting)

## Prerequisites

- `pip install numpy scipy pandas pysam pybedtools matplotlib`
- BEDTools installed system-wide (`conda install -c bioconda bedtools`)
- Familiarity with BED coordinate conventions (0-based half-open) and basic ATAC-seq QC (fragment size distribution)
- Related concepts: PWM/motif scanning (see `bio-core-motif-discovery`), ATAC-seq peak calling (see `bio-applied-scatac-chromatin`)

## Background

**ATAC-seq**: Tn5 transposase cuts and tags open chromatin. Nucleosome-free regions (NFRs) are cut frequently; nucleosome-occupied DNA is protected.

**Footprinting**: A bound TF blocks Tn5 at the binding site → local depletion of cuts at the motif center flanked by elevated cuts on either side.

| Feature | TF-bound sites | Unbound accessible sites |
|---|---|---|
| Cuts at motif center | Depleted (footprint) | Elevated |
| Cuts flanking motif | Elevated | Uniform |
| Footprint score | High (>1.5) | ~1.0 |

**Fragment size ladder** (QC on the BAM before footprinting):

| Peak | Size | Interpretation |
|---|---|---|
| NFR | <150 bp | Nucleosome-free — most informative for TF binding |
| Mono | ~200 bp | Mononucleosome |
| Di | ~400 bp | Dinucleosome |
| Tri+ | ~600 bp | Trinucleosome and higher |

NFR fraction >40% of total fragments = good library quality.

**Analysis workflow**: map reads → BAM; extract Tn5 insertion positions (5' end, +4/-5 bp offset); define TF motif sites (PWM scan or ChIP-seq peaks); compute insertion profiles in a ±200 bp window around each site; average across sites → aggregate profile; compute footprint score = flanking signal / central signal.

## Tn5 Insertion Profile from a BAM

**Goal:** turn raw ATAC-seq reads into a per-position insertion count profile around a set of motif centers, with the mandatory Tn5 offset correction.

**Approach:** for each read overlapping a window around a site, take its 5' end, shift +4 bp (forward strand) or -5 bp (reverse strand), then bin the corrected position relative to the site center.

```python
import numpy as np
import pysam


def accumulation_plot(bam, sites, window=200):
    """
    Compute the average Tn5 insertion profile across all sites.

    bam: pysam.AlignmentFile (opened in 'rb' mode)
    sites: list of (chrom, center) tuples (center = motif midpoint, 0-based)
    window: bp on each side of center to accumulate
    Returns: np.ndarray of shape (2*window+1,), mean insertions per site.
    """
    profile = np.zeros(2 * window + 1)
    for chrom, center in sites:
        for read in bam.fetch(chrom, max(0, center - window), center + window):
            if read.is_unmapped:
                continue
            # Tn5 dimer offset correction: +4 bp fwd strand, -5 bp rev strand
            if not read.is_reverse:
                pos = read.reference_start + 4
            else:
                pos = read.reference_end - 5 - 1
            offset = pos - center
            if -window <= offset <= window:
                profile[offset + window] += 1
    return profile / max(len(sites), 1)


with pysam.AlignmentFile("atac.bam", "rb") as bam:
    sites = [("chr1", 1_000_500), ("chr1", 2_004_210)]  # motif centers
    raw_profile = accumulation_plot(bam, sites, window=200)
```

## Footprint Score

$$\text{Footprint Score} = \frac{\text{mean(flanking signal)}}{\text{mean(central signal) + } \varepsilon}$$

Score >1 means the center is depleted relative to flanks. Well-footprinted TFs: 1.5–3.0.

```python
from scipy.ndimage import gaussian_filter1d


def footprint_score(profile, center_window=8, flank_window=(20, 60)):
    """
    Compute the footprint score from an insertion profile centered on a motif.

    profile: 1D array, position 0 (index len//2) is the motif center.
    center_window: +/- bp around the center counted as the footprint.
    flank_window: (inner, outer) bp range on each side used as flanking signal.
    """
    mid = len(profile) // 2
    central = profile[mid - center_window: mid + center_window].mean()
    left_flank = profile[mid - flank_window[1]: mid - flank_window[0]].mean()
    right_flank = profile[mid + flank_window[0]: mid + flank_window[1]].mean()
    flanking = (left_flank + right_flank) / 2
    return flanking / (central + 1e-6)


# Smooth before scoring — raw Tn5 counts are noisy at single-bp resolution
profile_smooth = gaussian_filter1d(raw_profile, sigma=2)
score = footprint_score(profile_smooth)
print(f"Footprint score: {score:.2f}")
```

Statistically compare bound vs. background site scores with a one-tailed Mann-Whitney U test (`scipy.stats.mannwhitneyu(scores_bound, scores_unbound, alternative="greater")`).

## pybedtools Interval Operations

**Goal:** go from raw ATAC peaks + motif calls to a clean, footprinting-ready site list.

**Approach:** intersect motifs with peaks to keep only accessible sites, extend to a footprinting window, drop blacklist regions, and annotate with nearest TSS.

```python
import pybedtools

peaks = pybedtools.BedTool("atac_peaks.bed")
motifs = pybedtools.BedTool("ctcf_motifs.bed")

# Motifs overlapping ATAC peaks (accessible, candidate bound sites)
motifs_in_peaks = motifs.intersect(peaks, u=True)

# Negative control: motif sites NOT in accessible chromatin
background = motifs.intersect(peaks, v=True)

# Extend to +/-200 bp window for footprinting (requires a genome file or dict)
motifs_extended = motifs_in_peaks.slop(b=200, genome="hg38")

# Remove ENCODE blacklist regions
blacklist = pybedtools.BedTool("hg38_blacklist.bed")
clean = motifs_extended.subtract(blacklist)

# Annotate with closest gene TSS
tss = pybedtools.BedTool("hg38_tss.bed")
nearest = motifs_in_peaks.closest(tss, d=True)
```

Install: `pip install pybedtools` (requires BEDTools on `PATH`; `conda install -c bioconda bedtools` is the easiest route).

## Production Footprinting (TOBIAS)

For real (non-simulated) footprinting, use a dedicated bias-corrected tool rather than raw Tn5 counts — Tn5 has strong sequence preference that a naive footprint score does not correct for.

```bash
# 1. Correct Tn5 sequence bias
TOBIAS ATACorrect --bam atac.bam --genome hg38.fa --peaks peaks.bed --outdir atac_corrected/

# 2. Score footprints (per-position footprint score, bigWig output)
TOBIAS FootprintScores --signal atac_corrected/atac_corrected.bw --regions peaks.bed --output footprints.bw

# 3. Differential footprinting across conditions
TOBIAS BINDetect --motifs motifs.jaspar --signals cond1.bw cond2.bw --genome hg38.fa --peaks peaks.bed
```

## Pitfalls

- **Tn5 offset correction is mandatory**: the Tn5 dimer inserts with a +4 bp (forward strand) / -5 bp (reverse strand) offset from the read 5' end. Skip this and footprint centers will be shifted and scores diluted.
- **Fragment size selection**: use only NFR fragments (<150 bp) for footprinting. Mono/dinucleosomal fragments dilute the footprint signal.
- **Low-occupancy TFs don't footprint**: only highly occupied (>50% of sites bound) TFs produce detectable footprints. CTCF, cohesin, and pioneer factors work well; transient binders (MYC) rarely do.
- **Sequence bias**: Tn5 has a strong insertion sequence preference (~10 bp periodic). Always compare to an expected (background/Tn5-bias-corrected) insertion profile, not just raw observed counts — this is what TOBIAS ATACorrect does for you.
- **Coordinate systems**: BED is 0-based half-open; VCF/GFF/GTF are 1-based inclusive — mixing them causes off-by-one errors when merging peak/motif/TSS files.
- **`slop`/`genome` argument**: pybedtools `slop` needs a genome file or `{chrom: length}` dict matching your reference build; a mismatched build silently truncates or drops intervals near chromosome ends.
- **Batch effects**: always check for batch confounding (sequencing run, extraction date) before interpreting differential footprinting signal between conditions.

## See Also

- `bio-core-motif-discovery` — PWM scanning to define motif site positions before footprinting
- `bio-applied-scatac-chromatin` — ATAC-seq/scATAC peak calling and chromatin accessibility upstream of footprinting
- `bio-applied-chipseq-pipeline` — ChIP-seq peak calling as an alternative source of TF binding sites
- `atac-seq-analysis` — broader ATAC-seq processing pipeline (alignment through peak calling)