# Atac Seq Analysis

> Analyze ATAC-seq BAM/BED data with pysam and pybedtools — fragment-size QC, NFR fraction, Tn5 +4/-5 offset correction, and TF footprint scoring around motif sites. Use when doing ATAC-seq QC, computing nucleosome-free-region fraction, correcting Tn5 insertion bias, or scoring transcription-factor footprints from chromatin accessibility data.

- Skill: `pavel-kravchenko/atac-seq-analysis` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/atac-seq-analysis`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/atac-seq-analysis/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Web & Frontend
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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.

```python
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.

```python
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.

```python
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.

```python
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.

```python
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