bio-applied-tf-footprinting

name	bio-applied-tf-footprinting
description	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.
tool_type	python
primary_tool	pybedtools

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.

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.

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.

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.

# 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)