Atac seq analysis

ATAC-seq Analysis: Fragment QC and TF Footprinting

When to Use

• Running ATAC-seq quality control: fragment-size distribution, nucleosome-free-region (NFR) fraction, FRiP
• Correcting Tn5 transposase insertion bias (+4/-5 bp offset) before any per-base signal analysis
• Computing footprint scores to detect TF occupancy from chromatin accessibility around motif sites
• Intersecting ATAC peaks with motif calls or blacklist regions using pybedtools
• Building accumulation (meta-profile) plots averaged across many genomic loci

Version Compatibility

pysam ≥0.22, pybedtools ≥0.10 (requires bedtools ≥2.30 on PATH), numpy ≥1.26, scipy ≥1.13, matplotlib ≥3.8, Python ≥3.10.

Prerequisites

• pip install pysam pybedtools numpy scipy matplotlib
• Coordinate-sorted, indexed BAM (see bio-alignment-files-bam-statistics); peak calls in BED (see bio-atac-seq-atac-peak-calling)
• Familiarity with BED interval semantics (bio-genome-intervals-bed-file-basics)

Quick Reference

Metric	Threshold	Notes
NFR fraction (< 150 bp)	≥ 40%	key ATAC-seq quality indicator
Mono-nucleosome fraction	~35% (150–300 bp)	characteristic ladder band
Di-nucleosome fraction	~15% (300–500 bp)
Tn5 forward offset	+4 bp	shift forward read starts +4
Tn5 reverse offset	-5 bp	shift reverse read ends -5
Footprint score > 1.5	strong footprint	flanking/central insertion ratio
FRiP (Fraction of Reads in Peaks)	≥ 0.20	enrichment quality check

Core Patterns

Goal: Quantify the nucleosomal ladder and check the NFR fraction QC gate. Approach: Bin fragment sizes from samtools stats/pysam template lengths into NFR/mono/di/tri classes and report proportions.

import numpy as np
import matplotlib.pyplot as plt

def classify_fragments(sizes):
    """Classify ATAC fragment sizes into nucleosomal ladder bins."""
    sizes = np.asarray(sizes)
    nfr   = sizes[(sizes > 50) & (sizes < 150)]
    mono  = sizes[(sizes >= 150) & (sizes < 300)]
    di    = sizes[(sizes >= 300) & (sizes < 500)]
    tri   = sizes[sizes >= 500]
    total = len(sizes)
    return {"NFR": len(nfr)/total, "Mono": len(mono)/total,
            "Di": len(di)/total, "Tri+": len(tri)/total}

fracs = classify_fragments(fragment_sizes)  # fragment_sizes: array of |TLEN| from BAM
print(f"NFR fraction: {fracs['NFR']:.1%}")  # target >= 40%

plt.hist(fragment_sizes, bins=200, range=(0, 800), density=True, color="steelblue")
plt.xlabel("Fragment size (bp)"); plt.ylabel("Density")
plt.title("ATAC-seq fragment size distribution (nucleosomal ladder)")

Goal: Extract Tn5-offset-corrected, per-base insertion counts from a BAM. Approach: Fetch reads in a region, shift forward reads +4 and reverse read ends -5, and tally insertion sites into a per-base array — the mandatory correction before any footprinting.

def tn5_corrected_insertions(bam_file, chrom, start, end):
    """Per-base Tn5 insertion counts with +4/-5 offset correction.

    Forward reads: insertion at reference_start + 4
    Reverse reads: insertion at reference_end - 5
    """
    import pysam
    insertions = np.zeros(end - start)
    with pysam.AlignmentFile(bam_file, "rb") as bam:
        for read in bam.fetch(chrom, start, end):
            if read.is_unmapped or read.is_secondary:
                continue
            pos = read.reference_start + 4 if not read.is_reverse else read.reference_end - 5
            if start <= pos < end:
                insertions[pos - start] += 1
    return insertions

Goal: Score how strongly a TF is bound at a set of motif sites from an averaged insertion profile. Approach: Compare depletion in a narrow central window (the footprint) against flanking accessibility — a high ratio indicates protection by a bound factor.

def footprint_score(profile, center_window=8, flank_window=(20, 60)):
    """Footprint score = mean flanking signal / mean central signal.

    profile: 1D array of insertion counts, centered on the motif midpoint.
    Higher score = deeper central depletion = stronger TF occupancy signal.
    """
    mid = len(profile) // 2
    central  = profile[mid - center_window : mid + center_window].mean()
    left     = profile[mid - flank_window[1] : mid - flank_window[0]].mean()
    right    = profile[mid + flank_window[0] : mid + flank_window[1]].mean()
    flanking = (left + right) / 2
    return flanking / (central + 1e-9)

Goal: Restrict motif sites to open chromatin and build a footprinting window, removing blacklisted regions. Approach: Use pybedtools to intersect, extend (slop), and subtract BED files without writing manual interval-overlap code.

import pybedtools

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

motifs_in_peaks = motifs.intersect(peaks, u=True)          # motifs overlapping ATAC peaks
motif_windows   = motifs_in_peaks.slop(b=200, genome="hg38")  # +/-200 bp footprint window

blacklist = pybedtools.BedTool("hg38_blacklist.bed")
clean = motif_windows.subtract(blacklist)                  # drop ENCODE blacklist regions

clean.sequence(fi="hg38.fa", fo="windows.fa")               # write windows as FASTA
print(f"Motifs in peaks: {len(motifs_in_peaks)}")

Goal: Run the full footprinting pipeline end-to-end and plot the aggregate meta-profile. Approach: Filter motifs to peak-overlapping sites, pull a Tn5-corrected profile per site, average, smooth, and score.

from scipy.ndimage import gaussian_filter1d

def footprinting_pipeline(bam_file, peaks_bed, motifs_bed, window=200):
    """End-to-end footprinting: filter motifs -> extract profiles -> score."""
    peaks  = pybedtools.BedTool(peaks_bed)
    motifs = pybedtools.BedTool(motifs_bed).intersect(peaks, u=True)

    profiles = []
    for site in motifs:
        chrom, center = site.chrom, (site.start + site.end) // 2
        profiles.append(tn5_corrected_insertions(
            bam_file, chrom, center - window, center + window))

    if not profiles:
        return None, None

    profiles = np.array(profiles)
    mean_profile = gaussian_filter1d(profiles.mean(axis=0), sigma=3)
    score = footprint_score(mean_profile)
    return mean_profile, score

mean_profile, fp_score = footprinting_pipeline("sample.bam", "peaks.bed", "ctcf.bed")

def accumulation_plot(profiles, positions, title="Meta-profile", sigma=5):
    """Plot mean +/- SEM insertion profile across many sites (profiles: n_sites x window)."""
    mean_p = profiles.mean(axis=0)
    sem_p  = profiles.std(axis=0) / np.sqrt(len(profiles))
    smoothed = gaussian_filter1d(mean_p, sigma=sigma)

    fig, ax = plt.subplots(figsize=(8, 4))
    ax.plot(positions, smoothed, color="steelblue", lw=2)
    ax.fill_between(positions, smoothed - sem_p, smoothed + sem_p, alpha=0.3, color="steelblue")
    ax.axvline(0, color="red", lw=1, ls="--", label="Motif center")
    ax.set_xlabel("Position relative to motif (bp)"); ax.set_ylabel("Mean Tn5 insertions")
    ax.set_title(title); ax.legend(frameon=False)
    return fig, ax

Pitfalls

• Missing Tn5 offset: uncorrected reads shift the apparent footprint by 4-5 bp; always apply +4/-5 before aggregating.
• Not filtering NFR vs. nucleosomal reads: footprinting should use only NFR (< 150 bp) fragments — nucleosomal reads dilute the signal.
• Genome blacklist regions: always subtract the ENCODE blacklist before analysis; these regions have artifactually high depth.
• Motif strand: CTCF and other factors bind asymmetrically; flip reverse-strand sites before averaging profiles.
• Footprint score vs. occupancy: a high score only indicates accessibility asymmetry; validate with ChIP-seq/CUT&RUN where possible.
• pybedtools slop genome argument: requires a genome file or dict, e.g. {'chr1': (0, 248956422), ...}; pybedtools.genome_registry.hg38 is a convenient built-in.

See Also

• bio-atac-seq-footprinting — dedicated deeper TF footprinting workflows (TOBIAS/HINT-ATAC style)
• bio-atac-seq-atac-peak-calling — MACS2/Genrich peak calling upstream of this QC
• bio-atac-seq-nucleosome-positioning — NucleoATAC-style nucleosome calling from fragment sizes
• bio-genome-intervals-bed-file-basics — BED/pybedtools interval fundamentals