Bio Atac Seq Nucleosome Positioning

Version Compatibility

Reference examples tested with: NucleoATAC 0.3.4+, ATACseqQC 1.26+, DANPOS 3.1+, samtools 1.19+, pysam 0.22+, pyBigWig 0.3+, BSgenome.Hsapiens.UCSC.hg38 1.4+, TxDb.Hsapiens.UCSC.hg38.knownGene 3.18+.

NucleoATAC is unmaintained since 2018 but remains the canonical ATAC-specific nucleosome caller; ATACseqQC, DANPOS3, and scprinter are actively developed alternatives. Verify versions before use:

• Python: pip show <package> then help(module.function) to check signatures
• R: packageVersion('<pkg>') then ?function_name to verify parameters
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws unexpected errors, introspect the installed package and adapt rather than retrying.

Nucleosome Positioning

"Where are the nucleosomes in my ATAC-seq data?" -> Use fragment-size classes (Tn5 cuts twice through naked DNA generating short fragments; once on each side of a single nucleosome generating ~147+linker fragments) to call nucleosome centers, occupancy scores, and the spacing pattern around regulatory elements.

• CLI: nucleoatac run --bed regions.bed --bam sample.bam --fasta genome.fa
• R: ATACseqQC::splitGAlignmentsByCut() -> fragment classes; factorFootprints() -> per-TF flanking nuc analysis
• CLI: python danpos.py dpos sample.bam (alternative; supports MNase, ATAC, DNase)
• Python: scprinter for multi-scale nucleosome inference

Nucleosome Physics for ATAC

A nucleosome wraps ~147 bp DNA in 1.65 turns. Adjacent nucleosomes are separated by 20-50 bp linker; mean nucleosome repeat length (NRL) is species-dependent:

NRL determines fragment-size peak positions. ATAC mono-nucleosome peak is at NRL (NOT 147 bp -- that's the protected length; ATAC fragments span the full nucleosome+linker). Di-nuc is at 2x NRL minus a small overlap.

Fragment-Size Classes (Buenrostro 2013, refined)

Mono-nucleosome window 180-247 bp is the Buenrostro 2013 convention; ATACseqQC uses 180-250. Adjust for the target organism's NRL.

V-Plot Interpretation

V-plots (fragment-size vs position) are diagnostic. X-axis is position relative to a feature (TSS, motif center); Y-axis is fragment size. Aggregate density forms characteristic patterns:

V-plots are the primary diagnostic for whether nucleosome-positioning analysis will succeed. Flat-band patterns mean no positioning information; classic V/W patterns mean positioning is recoverable.

Algorithmic Taxonomy

Methodology evolves; verify against current Schep 2015 (NucleoATAC), Chen 2013 (DANPOS), Hu 2025 (scPrinter) before locking pipelines.

+1 Nucleosome Calling

The +1 nucleosome (first nucleosome downstream of TSS, immediately bordering the NFR) is the most-studied positioning feature. Its position relative to TSS determines transcription initiation kinetics.

Canonical +1 position: +50 to +60 bp from TSS in metazoa; -100 to -120 bp from TATA in yeast; varies by gene type (Pol II vs Pol III, housekeeping vs developmental).

Calling strategy:

Goal: Identify each gene's +1 nucleosome, the first nucleosome downstream of the TSS that flanks the NFR.

Approach: Build gene-body intervals slopped around TSSs, run NucleoATAC over them to call per-base nucleosome positions, then pick the most-downstream-of-TSS nucleosome per gene.

# 1. Define gene-body intervals
bedtools slop -i genes.bed -g chrom.sizes -l 200 -r 1000 > gene_bodies.bed

# 2. Run NucleoATAC
nucleoatac run --bed gene_bodies.bed --bam sample.dedup.bam --fasta genome.fa \
    --out tss_nuc/ --cores 8

# 3. The first nucleosome downstream of each TSS in nucpos.bed is +1

A failure to detect a clear +1 peak in aggregate V-plot suggests TSS annotation is wrong or library is over-transposed.

Per-Tool Failure Modes

NucleoATAC -- Region size and depth dependence

Trigger: Short region BED (< 1 kb per region); shallow library (< 25M nuclear reads).

Mechanism: NucleoATAC fits an idealized V-plot template per region. Short regions provide too few fragments for stable correlation; shallow data provides noisy templates.

Symptom: No nucleosome calls in shallow regions; "occupancy" track is flat at zero.

Fix: Use regions >= 500 bp; merge adjacent peaks via bedtools to ensure region size; require >= 30M nuclear reads.

NucleoATAC -- Maintenance status

Trigger: Installing NucleoATAC in 2025+.

Mechanism: Last release 2018; pegs Python 3.6 in some installs; depends on outdated NumPy API.

Fix: Use a dedicated conda env (conda create -n nucleoatac python=3.7 numpy=1.18 scipy=1.5 pysam); accept it works but is no longer updated. Consider scprinter or DANPOS3 alternatives for new projects.

ATACseqQC factorFootprints -- Asymmetric nucleosome flanks

Trigger: Pioneer-factor binding sites where one face is on a nucleosome.

Mechanism: factorFootprints assumes symmetric flanking nucleosomes. Pioneer TFs (FOXA1, GATA) only have nucleosome on one side -> asymmetric output.

Symptom: Single shoulder in flanking signal; unbalanced V-plot.

Fix: Treat asymmetry as biological signal, not artefact. For pioneers, use stranded analysis.

DANPOS dpos with default parameters -- ATAC mismatch

Trigger: Running python danpos.py dpos with MNase defaults on ATAC.

Mechanism: DANPOS3's smoothing window and peak-calling defaults are tuned for MNase signal (smoother coverage). ATAC's sharper signal requires --smooth_width 80 --width 145 or similar; otherwise calls are over-smoothed.

Fix: Use ATAC-tuned parameters. See DANPOS docs for ATAC-specific recipe; or use NucleoATAC instead.

Mono-nucleosome filter window mis-set

Trigger: Using strict 147 bp filter for mono-nuc fraction; using 100-180 bp instead of 180-247.

Mechanism: Mono-nuc fragments are 180-247 bp because they span the nucleosome AND a linker. Filtering tighter excludes the legitimate signal.

Symptom: Mono-nuc count is much lower than expected (< 30% of NFR count).

Fix: Use Buenrostro 2013 windows: NFR < 100, mono 180-247, di 315-473.

Decision Tree by Goal

Estimating NRL from Fragment-Size Distribution

Goal: Estimate the nucleosome repeat length from ATAC fragment-size periodicity.

Approach: Collect proper-pair fragment lengths from the BAM, build a histogram, find density peaks via scipy find_peaks, and read off the mono-nucleosome peak position within the 150-250 bp window.

import numpy as np, pysam
from scipy.signal import find_peaks

bam = pysam.AlignmentFile('sample.bam', 'rb')
frag_lengths = [abs(r.template_length) for r in bam.fetch()
                if r.is_proper_pair and r.is_read1 and 0 < abs(r.template_length) < 1500]
hist, edges = np.histogram(frag_lengths, bins=300, range=(0, 1500))
peaks, _ = find_peaks(hist, distance=50, prominence=hist.max() * 0.05)
peak_positions = edges[peaks] + (edges[1] - edges[0]) / 2
# Mono peak should be ~NRL; di peak ~2*NRL
mono = peak_positions[(peak_positions > 150) & (peak_positions < 250)][0]
print(f'Estimated NRL: {mono:.0f} bp')

NRL inferred this way is approximate; for precision use autocorrelation on cumulative cleavage coverage instead.

V-Plot in Python

Goal: Build a fragment-size-by-position density plot to diagnose nucleosome positioning around a feature.

Approach: Iterate proper-pair fragments in a flank window around each feature center, accumulate counts into a (fragment_size x position) grid, and render the 2D density.

import numpy as np, pysam, matplotlib.pyplot as plt

def vplot(bam_path, regions_bed, max_size=600, flank=1000):
    bam = pysam.AlignmentFile(bam_path, 'rb')
    grid = np.zeros((max_size, 2 * flank))
    for line in open(regions_bed):
        chrom, start, *_ = line.strip().split('\t')
        center = int(start)
        for r in bam.fetch(chrom, max(0, center - flank), center + flank):
            if not r.is_proper_pair or not r.is_read1: continue
            size = abs(r.template_length)
            if size <= 0 or size >= max_size: continue
            frag_center = r.reference_start + size // 2
            x = frag_center - center + flank
            if 0 <= x < 2 * flank:
                grid[size, x] += 1
    return grid

g = vplot('sample.bam', 'tss.bed')
plt.imshow(g, aspect='auto', origin='lower', cmap='magma',
           extent=[-1000, 1000, 0, 600])
plt.xlabel('Distance from feature (bp)')
plt.ylabel('Fragment size (bp)')
plt.savefig('vplot.png', dpi=200, bbox_inches='tight')

V-plot quality is the most useful diagnostic before nucleosome calling. Classic V at TSS = positioning info recoverable; flat band = not.

Differential Nucleosome Positioning (DANPOS3 dpos)

# Compare control vs treatment nucleosome positions.
# The sample pair is the POSITIONAL argument (a:b means a minus b); -b is for background/input to
# subtract, and -c specifies a read-count to normalize to (an integer, NOT a control BAM path).
python danpos.py dpos condition2.bam:condition1.bam \
    -o danpos_diff/ \
    --paired 1 \
    --smooth_width 80

DANPOS reports four event types: shifted nucleosomes, gained, lost, fuzziness change. ENCODE has no official threshold; require >= 30 bp shift and FDR < 0.05 for nucleosome shift calls.

Full ATAC-tuned DANPOS3 recipe:

# --width 145: summit-scan window (DANPOS -jw/--width; default 40)
# --smooth_width 80: smoothing kernel width (DANPOS -z; default 20, widened for ATAC)
# -jd 145: min distance between adjacent nuc calls (single-dash short flag)
# --pheight 1e-5: occupancy P-value cutoff (DANPOS -p; dpos default 0). -q/--height is the separate density cutoff
# --frsz 200: fragment size used (mono-nuc)
python danpos.py dpos sample.bam \
    --paired 1 \
    --width 145 \
    --smooth_width 80 \
    -jd 145 \
    --pheight 1e-5 \
    --frsz 200 \
    --out danpos_out/

Verify exact flags with python danpos.py dpos --help; DANPOS3 (github.com/sklasfeld/DANPOS3) is invoked as python danpos.py, not a danpos3 executable, and installs from GitHub (the bioconda danpos package is DANPOS2). Its documentation has been spotty and flag names can drift across releases.

Adapted from DANPOS3 docs for ATAC; --smooth_width 80 widens the smoothing kernel to match ATAC's sharper signal vs MNase's broader cleavage. -jd 145 (single-dash short, alternative --distance 145) enforces nucleosome spacing >= 145 bp (one nucleosome footprint).

Histone Variant Detection from Fragment Size

Trigger: Suspected H2A.Z- or H3.3-containing nucleosomes; differential nucleosome composition between conditions.

Mechanism: H2A.Z replacement of H2A produces nucleosomes with weaker DNA-histone interaction (lower thermal stability); the H2A.Z population tends toward shorter fragment sizes than canonical H2A nucleosomes. H3.3 replacement is more subtle, but H3.3-H2A.Z double-variant nucleosomes are particularly destabilized at active promoters (Jin 2009 Nat Genet 41:941-945).

Detection: Aggregate fragment-size distribution at H2A.Z ChIP-seq peaks vs H3K4me3-only peaks; the H2A.Z population shows mean fragment size ~10 bp shorter. ATAC alone CANNOT definitively call H2A.Z; H2A.Z ChIP-seq is needed for ground truth. ATAC fragment-size analysis is a hypothesis generator.

# Per-region fragment-size mean as H2A.Z indicator
def region_frag_size(bam, region):
    sizes = [abs(r.template_length) for r in bam.fetch(*region)
             if r.is_proper_pair and r.is_read1 and 100 < abs(r.template_length) < 300]
    return np.mean(sizes) if sizes else np.nan

# Compare H2A.Z-positive vs H2A.Z-negative TSSs

Long-Read Single-Molecule Chromatin (Fiber-seq, NanoNOMe)

Alternative to short-read ATAC for nucleosome positioning:

Fiber-seq can detect nucleosome occupancy directly per single chromatin molecule (no aggregation needed). Resolves cell-cycle-dependent and stochastic positioning that bulk ATAC averages out. Preferred for fine-structure analysis of regulatory elements.

For most labs, short-read ATAC + NucleoATAC remains primary; Fiber-seq is special-purpose when single-molecule resolution is essential.

Nucleosome Fuzziness

Fuzziness measures how sharply positioned a nucleosome is across cells. Defined as the standard deviation of per-cell nucleosome center positions.

These ranges are field-convention bands (drawn from NucleoATAC / DANPOS practice); no single primary paper prescribes them -- verify against tool-specific documentation when reporting.

NucleoATAC reports per-nucleosome fuzziness in the fuzziness column (column 13) of .nucpos.bed -- a measure of how wide the signal peak is; well-positioned nucleosomes show low fuzziness (~20-50 bp), consistent with the table above.

Common Errors

References

• Schep AN et al 2015 Genome Res 25:1757 (NucleoATAC)
• Chen K et al 2013 Genome Res 23:341 (DANPOS)
• Buenrostro JD et al 2013 Nat Methods 10:1213 (ATAC fragment-size classes)
• Ou J et al 2018 BMC Genomics 19:169 (ATACseqQC)
• Hu Y et al 2025 Nature 638:779 (scPrinter/PRINT; multiscale footprints)
• Mavrich TN et al 2008 Nature 453:358 (+1 nucleosome positioning)
• Voong LN et al 2016 Cell 167:1555-1570 (high-resolution chemical nucleosome mapping)
• Jin C et al 2009 Nat Genet 41:941 (H3.3/H2A.Z double-variant nucleosome instability at active regions)
• Teif VB et al 2012 Nat Struct Mol Biol 19:1185 (NRL variation across cell types)

Related Skills

• atac-seq/atac-qc - Fragment-size periodicity QC
• atac-seq/atac-peak-calling - Nucleosome-aware MACS3 hmmratac
• atac-seq/footprinting - Per-TF flanking nucleosome analysis
• atac-seq/single-cell-atac - scprinter for sc nucleosome positioning
• chip-seq/peak-annotation - Annotate nucleosome positions to genes
• alignment-files/bam-statistics - Insert-size statistics upstream