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.