parapeak 0.4.0

NGS peak caller with parallel computation and GC-corrected statistics

parapeak

An NGS peak caller supporting shallow (MiSeq) to deep coverage across ATAC-seq, genome-editing assays (CIRCLE-seq, GUIDE-seq), ChIP-seq, CUT&TAG, and similar experiments.

parapeak is inspired by MACS3 but differs in several key aspects:

Feature	MACS3	parapeak
Parallelism	Single-threaded	Multi-process, per chromosome
Genome size	User-supplied --gsize	Inferred from BAM headers
Statistical model	Poisson	Negative binomial and GC-corrected Z-score (both required)
GC correction	Separate tool	Built-in (no reference FASTA needed)
Blacklist	Not built-in	Applied at pileup construction time
Low-coverage support	Limited	Signal filters for MiSeq-depth data
Paired-end fragments	Fixed extension	Actual insert size from TLEN field
Unsorted BAM support	Requires sorted + indexed	Single-pass; no sorting or index required
Read QC filters	Basic	MAPQ threshold, QC-fail flag, duplicate/secondary/supplementary
Input	BAM/SAM/BED	BAM/SAM only (broken headers auto-repaired)

---

Installation

pip install -r requirements.txt
pip install -e .

Dependencies: pysam, numpy, scipy, numba

numba accelerates the negative-binomial survival function and the Benjamini–Hochberg isotonic step with JIT compilation. If numba is unavailable, parapeak falls back transparently to scipy.

---

Quick start

# Paired-end ATAC-seq, unsorted BAM
parapeak -t atac.bam --blacklist hg38-blacklist.v2.bed.gz \
         -o results/ -n ATAC -p 8

# CIRCLE-seq / genome-editing with NTC control (MiSeq depth)
parapeak -t treated.bam -c ntc.bam \
         --min-count 3 --min-fold 3.0 --pseudocount 2.0 \
         -o results/ -n edit -p 8

# Multiple replicates, stricter MAPQ filter
parapeak -t rep1.bam rep2.bam --min-mapq 30 -o results/ -p 4

BAM files do not need to be sorted or indexed. parapeak reads each file in a single sequential pass.

---

Options

Input:
  -t, --treated BAM [BAM ...]   Treated BAM/SAM files (sorting and indexing not required)
  -c, --control BAM [BAM ...]   Control/input BAM/SAM files
  --blacklist BED               Blacklist BED file (plain or gzip-compressed)

Output:
  -o, --output DIR              Output directory (default: parapeak_output)
  -n, --name NAME               Output file prefix (default: parapeak)
  --bedgraph-value TYPE         Value written to the bedGraph output.
                                fold_enrichment (default): (treated + pc) / (scaled control + pc)
                                pvalue: -log10 of the minimum NB / Z-score p-value per bin
                                pileup: raw treated read count per bin

Algorithm parameters:
  --local-window BP             Local background window size in bp (default: 1000)
  --fragment-size BP            Fallback extension for single-end or discordant
                                paired-end reads; estimated from paired-end insert
                                sizes if omitted, falls back to 200 bp
  --bin-size BP                 Pileup bin size in bp (default: 10)
  -q, --qvalue FLOAT            Q-value threshold for both methods (default: 0.05)
  --min-length BP               Minimum peak length in bp (default: 200)
  --max-gap BP                  Maximum gap to merge between significant bins (default: 30)

Read QC filters:
  --min-mapq INT                Minimum mapping quality (MAPQ). Reads below this
                                threshold are discarded (default: 20)
  --min-fragment BP             Minimum paired-end insert size accepted from TLEN.
                                Smaller inserts are treated as single-end (default: 10)
  --max-fragment BP             Maximum paired-end insert size accepted from TLEN.
                                Larger inserts are treated as single-end (default: 2000)

Signal filters:
  --min-count FLOAT             Minimum pileup at the peak summit bin (default: 5)
  --min-fold FLOAT              Minimum fold enrichment over background (default: 2.0)
  --pseudocount FLOAT           Background floor per local window (default: 1.0)

  -p, --threads INT             Number of parallel worker processes (default: 1)

---

Read QC filtering

The following filters are applied to each read in order before it contributes to the pileup. Counts for each filter are reported in the JSON run report.

Filter	SAM flag / field	Default behaviour
Unmapped	0x4	Discarded
Duplicate	0x400	Discarded
Secondary alignment	0x100	Discarded
Supplementary alignment	0x800	Discarded
QC fail	0x200	Discarded
Low mapping quality	MAPQ < --min-mapq	Discarded (default MAPQ < 20)
Read 2 of a pair	0x80	Skipped (not a quality filter; avoids double-counting)

Paired-end reads that pass all filters but whose TLEN falls outside [--min-fragment, --max-fragment] are not discarded; they are treated as single-end reads and extended by --fragment-size. The count of such reads is reported as insert_size_fallback in the JSON run report.

---

Single-end and paired-end handling

Single-end

Each mapped read is extended from its 5′ end by --fragment-size bp toward the sequenced fragment's expected end. --fragment-size is estimated automatically as the median insert size if the BAM contains paired reads; otherwise it defaults to 200 bp.

5'─────read──────3'──────────────extension──────────────►
|←────────────── fragment_size ────────────────────────→|

Paired-end

Only R1 is processed; R2 is skipped to avoid double-counting fragments.

For proper pairs whose TLEN (template length from the SAM field) falls within [--min-fragment, --max-fragment] (default 10–2000 bp), the actual fragment interval is used directly:

R1  5'──────────►
                    ◄──────────5'  R2
|←───────── TLEN (actual insert size) ─────────────────→|

This is important for experiments where fragment size varies:

• ATAC-seq: each pair of Tn5 cut sites defines an accessible region; using actual insert sizes correctly captures nucleosome-free regions (NFR, < 200 bp) and mono-nucleosomal fragments (∼200 bp).
• CIRCLE-seq / GUIDE-seq: read pairs span the circularised DNA around cut sites; TLEN reflects the circle diameter.
• ChIP-seq paired-end: sonication produces a range of fragment sizes; actual TLEN avoids over-smoothing short fragments.

For discordant pairs (TLEN outside [--min-fragment, --max-fragment]) and unpaired reads, the read is treated as single-end and extended by --fragment-size.

---

Signal filters for low-coverage data

At shallow sequencing depth (MiSeq, ~1–5 M reads), stochastic sampling means that ~60% of 10 bp bins in the control have zero reads. Without a signal floor, NB and Z-score tests can declare bins with only 1–2 reads significant wherever the control happens to have zero coverage, producing thousands of artefactual peaks.

Three filters address this, all applied after the BH q-value test:

--pseudocount (default: 1.0)

Adds a pseudocount (in units of reads per --local-window) to the scaled control pileup before NB background fitting and fold-enrichment calculation. Distributing the pseudocount across the local window gives each 10 bp bin a minimum background floor of pseudocount / (local_window / bin_size).

Fold enrichment is computed as:

fold = (mean_treatment + pseudocount) / (mean_scaled_control + pseudocount)

This eliminates infinite fold values when the control has zero reads and makes the metric interpretable at all coverage depths.

--min-count (default: 5.0)

Discards peaks whose summit bin has a pileup value below this threshold, regardless of the q-value. At typical fragment extension (200 bp), a summit pileup of 5 corresponds to approximately 2–3 overlapping read fragments; for sparser data where even 2–3 reads can be statistically significant against zero control, this filter prevents low-signal artefacts.

--min-fold (default: 2.0)

Discards peaks with fold enrichment (pseudocount-corrected) below this value. Combined with --pseudocount, a minimum fold of 2.0 ensures that the signal is at least twice the effective background, not just twice zero. Set to values < 1 to disable.

---

Recommended parameters by experiment

Experiment	SE/PE	Typical depth	--min-count	--min-fold	--pseudocount	--local-window	Notes
ATAC-seq (deep)	PE	> 20 M	10	2.0	0.5	1000	Nucleosome-free regions; actual insert sizes used
ATAC-seq (MiSeq)	PE	2–5 M	3	2.0	2.0	2000	Wider local window compensates for sparse control
CIRCLE-seq	PE	1–5 M	3	3.0	2.0	1000	Sparse library; higher fold reduces false positives
GUIDE-seq	PE	5–20 M	5	2.0	1.0	1000	On-target site may need blacklisting
ChIP-seq narrow (deep)	SE/PE	> 10 M	10	2.0	0.5	1000	TF binding sites
ChIP-seq narrow (MiSeq)	SE	1–5 M	3	2.0	2.0	5000	Wider local window for stable background
ChIP-seq broad	SE/PE	> 10 M	5	1.5	1.0	5000	Histone marks; lower fold threshold
CUT&TAG	PE	5–20 M	5	2.0	1.0	500	High signal-to-noise; narrow local window
FAIRE-seq	SE/PE	> 5 M	5	2.0	1.0	1000	Similar to ATAC-seq

Coverage guide

Sequencer	Typical read pairs	Recommended --min-count	Recommended --pseudocount
NovaSeq / HiSeq (deep)	> 30 M	10–20	0.5
NextSeq / HiSeq (standard)	10–30 M	5–10	1.0
MiSeq (V3 2×300)	3–10 M	3–5	1.0–2.0
MiSeq (V2 2×150)	1–3 M	3	2.0–5.0

---

Output files

File	Format	Description
<name>_peaks.narrowPeak	ENCODE narrowPeak (BED6+4)	Final peak calls
<name>_summits.bed	BED	Single-base summit positions
<name>_peaks.tsv	TSV	All score columns for downstream analysis
<name>.bedGraph	UCSC bedGraph	Genome-wide signal track (see --bedgraph-value)
<name>_run.json	JSON	Settings and read QC statistics

The narrowPeak columns follow the ENCODE specification:

chrom  start  end  name  score  strand  signalValue  -log10(p)  -log10(q)  summit_offset

The TSV file includes separate NB and Z-score p-values and q-values for inspection of which statistical test is the driver for each peak.

bedGraph signal track

The bedGraph file covers the full genome at the same bin resolution as the pileup (--bin-size). Bins with zero signal and blacklisted bins are omitted; consecutive bins with identical values are merged into a single record to reduce file size. The default value is fold enrichment:

fold = (treated_bin + pseudocount) / (scaled_control_bin + pseudocount)

Use --bedgraph-value pvalue to write −log₁₀(p) (minimum of the NB and Z-score p-values per bin), or --bedgraph-value pileup to write the raw treated read count per bin.

The track is compatible with IGV, the UCSC Genome Browser, and any tool that accepts bedGraph format.

JSON run report

The JSON file records all settings and per-file read statistics for reproducibility and QC review. Example:

{
  "run": {
    "timestamp": "2026-06-15T10:00:00Z",
    "tool": "parapeak"
  },
  "settings": {
    "treated": ["sample.bam"],
    "control": null,
    "blacklist": null,
    "bin_size": 10,
    "fragment_size": null,
    "min_mapq": 20,
    "min_fragment": 10,
    "max_fragment": 2000,
    "qvalue": 0.05,
    "min_count": 5.0,
    "min_fold": 2.0,
    "pseudocount": 1.0,
    "threads": 8,
    "bedgraph_value": "fold_enrichment"
  },
  "statistics": {
    "fragment_size_used": 185,
    "scale_factor": 1.0,
    "treated": {
      "total_reads": 42000000,
      "filtered_unmapped": 1200000,
      "filtered_duplicate": 3500000,
      "filtered_secondary": 80000,
      "filtered_supplementary": 12000,
      "filtered_qcfail": 45000,
      "filtered_low_mapq": 2100000,
      "skipped_read2": 17500000,
      "insert_size_fallback": 320000,
      "reads_used": 17243000
    },
    "control": null,
    "peaks_called": 48312
  }
}

---

Algorithm

1. Read pileup (single-pass)

Each BAM/SAM file is read in a single sequential pass using pysam, with no region argument. BAM files do not need to be sorted or indexed.

For each file, numpy arrays (one per chromosome, length = chromosome size / bin size) are pre-allocated and held in memory. As each read is consumed from the stream, QC filters are applied (see Read QC filtering), the fragment interval is determined (see Single-end and paired-end handling), and the corresponding bins are incremented by 1.

This approach has O(N_reads) I/O cost regardless of sort order. The previous per-chromosome fetch() strategy required N_chromosomes scans of an unsorted file; on a typical human genome (25 chromosomes) this is a ~25× reduction in I/O for unsorted input.

Blacklisted bins are zeroed after the pass completes.

If pysam rejects a file due to a malformed header (e.g. duplicate @HD lines produced by some aligners or merge tools), parapeak automatically repairs the header: for SAM files the raw text is rewritten with duplicate @HD lines removed; for BAM files the same repair is applied after dumping the file to SAM text via samtools view. The fixed content is written to a temporary file, read, and then deleted. If repair fails, the file is skipped with a warning and the run continues with the remaining inputs.

2. GC correction model

Without requiring a reference genome, parapeak estimates per-bin GC content from the sequences of mapped reads (a reliable proxy for reference GC%).

The genome is divided into non-overlapping windows (--local-window). For each window, the mean GC fraction of overlapping reads and the observed read depth are recorded. All windows across all chromosomes are pooled, binned into 20 equal GC% intervals (0–5%, 5–10%, …, 95–100%), and the mean and variance of depth are computed per interval.

A piecewise-linear interpolation curve f(GC%) → (expected_depth, variance) is fitted to these bin statistics. This curve is used in Step 4 to derive GC-corrected expected values and their variances for every genomic bin. When too few data points are available (sparse libraries), the model falls back to the global mean depth.

3. Negative-binomial p-values

The scaled control pileup is floored by adding pseudocount / (local_window / bin_size) per bin before fitting.

For each bin, the background is estimated at two scales:

• Global λ: genome-wide mean of the pseudocount-adjusted scaled control.
• Local λ: sliding-window mean over --local-window bp.

The larger lambda (more conservative, higher p-value) is used, so enrichment is only called when signal exceeds both local and global baselines.

A Negative Binomial NB(r, p) is fitted by the method of moments:

r = μ² / (σ² − μ)     p = μ / σ²

When σ² ≤ μ the distribution degenerates to Poisson (r → ∞). The survival function is evaluated with a numba-JIT log-PMF recurrence.

4. GC-corrected Z-score p-values

Each bin receives a GC-matched expected depth E and variance V from Step 2. Expected values are scaled to match the total treatment read count.

Z = (observed − E) / √V
p_Z = P(Z_standard > Z)   [one-sided upper tail]

5. Genome-wide Benjamini–Hochberg correction

P-values from all chromosomes are pooled and corrected independently for each method. Steps 2–5 are parallelised across chromosomes (-p workers).

q_NB = BH(p_NB)     q_Z = BH(p_Z)

A bin is significant only when both q_NB < threshold and q_Z < threshold.

6. Peak region assembly and signal filtering

Contiguous significant bins are merged across gaps of up to --max-gap bp. Regions shorter than --min-length bp are discarded.

Signal filters are then applied in order:

1. --min-count: summit pileup ≥ threshold
2. --min-fold: fold enrichment (pseudocount-corrected) ≥ threshold

The summit is the bin with the highest treatment pileup. The score in the output is −log₁₀ of the minimum q-value across both methods.

---

Comparison with MACS3

MACS3 uses a linked-list data structure internally, which limits its ability to parallelise across chromosomes. parapeak replaces this with per-chromosome NumPy arrays and Python multiprocessing, with inter-process communication only at the genome-wide BH correction step.

MACS3 requires sorted and indexed BAM files for random-access region queries. parapeak reads each file in a single forward pass, which avoids this requirement and is more I/O-efficient for unsorted input.

The true runtime bottleneck is BAM decompression (via the C library underlying pysam) and gzip decompression of the blacklist. These I/O costs dominate over Python-level arithmetic, so there is no need for a compiled extension language; numba JIT is sufficient for the compute-intensive statistical routines.

---

Author

Takaho A. Endo

License

MIT