RapidMACS

Standalone C++ library implementing MACS3's paired-end narrow peak-calling path for FRAG, BAMPE, and BEDPE input, with optional treatment/control local lambda. The implementation has been extracted from Chromap-suite/src/peak_caller/, where it has been validated end-to-end inside libchromap's concurrent mode and against chromap_callpeaks on the 100K and full PBMC 3k multiome benchmarks.

RapidMACS is open source under the MIT license. It is a specialized paired-end narrow-peak implementation, not a complete replacement for MACS3.

Two consumers are intended:

1. libchromap — links RapidMACS to call peaks in-process during the ATAC mapping pass (concurrent / memory-source mode preserved via a thin FragmentIterator adapter over the in-memory fragment store).
2. rapidmacs standalone CLI — reads a Chromap/ARC 5-col fragments TSV, paired-end BAM, or BEDPE and emits narrowPeak plus summits. Replaces / supersedes Chromap-suite's chromap_callpeaks.

See:

• TECHNICAL-METHODOLOGY.md — implemented algorithms, CLI methodology, validation datasets, parity criteria, and the paper benchmark protocol. This is the authoritative technical reference.
• BENCHMARKS.md — measured results and their provenance.
• NOTICE — third-party attributions (MACS3, htslib).
• docs/archive/ and plans/ — superseded handoffs and development-history plans. Retained for provenance; not normative.

Status

Phase A (extraction) and Phase B streaming both landed. The MACS3 FRAG path (--macs3-frag-narrowpeak) now streams via RunMacs3FragPeakPipeline- FromSortedIterator and the sweep-line workspace by default — no fragment list materialization, ~tens of KB peak workspace memory. Falls back to the events-buffer workspace if the input violates chrom-grouped + start-sorted order.

Paired-end ChIP-seq/CUT&RUN support is also implemented. BAMPE parsing follows MACS3 3.0.3's proper-pair and primary-alignment semantics. BAMPE and BEDPE accept --keep-dup all or a positive coordinate cap and default to --keep-dup 1; BEDPE accepts MACS3's three-column fragment form and conventional six-column mate intervals. With -c, local lambda is the MACS3 maximum of the matched control's d/slocal/llocal endpoint windows and genome-wide background, including default scale-to-small and optional --scale-to large behavior.

MACS3-compatible --nolambda and --call-summits are opt-in on BED, FRAG, BAMPE, and BEDPE callpeak paths. --nolambda retains the effective-genome-size and treatment/control depth-scaling arithmetic but constructs only the constant genome-wide background, avoiding d/slocal/llocal window work. --call-summits reproduces MACS3's Savitzky-Golay derivative and peak-shape filter and emits MACS3-style 1a, 1b, ... subpeak names. Neither flag changes the default single-summit/local-lambda behavior.

See TECHNICAL-METHODOLOGY.md for the final complete-CLI ATAC, CUT&RUN, and ChIP-seq triplicates, including wall-time ranges, peak RSS, scaling, output-parity criteria, and provenance.

The legacy custom-caller path (--out-prefix → BuildBinnedCutSignal + CallPeaksOnBins) also streams (uint16 bins, sliding-window lambda, sparse BH).

Two test lanes guard regressions:

• make test-legacy-streaming-100k — --out-prefix exercises the streaming path (uint16 bins, sliding-window lambda, sparse BH). Currently runs the 100K fixture in ~155 MB RSS / 0.35 s and produces a byte-identical 121,592-peak narrowPeak vs upstream chromap_callpeaks.
• make test-macs3-frag-100k — --macs3-frag-narrowpeak + --macs3-frag-summits exercises the MACS3 FRAG path via the sweep-line workspace. Currently runs the 100K in ~55 MB RSS / 0.4 s and reproduces the frozen 2,047-peak narrowPeak + summits snapshot exactly.

FRAG retains the historical cutoff of 5 (= MACS3 -p 1e-5), and the 100K regression smoke is byte-identical to its MACS3 3.0.3 golden. Paired BAMPE/BEDPE input defaults to MACS3's -q 0.05 and dynamic mean-fragment peak geometry. Pass --preset shiftTAG to opt in to MACS3 -p 0.01 (bdgpeakcall_cutoff = 2). Place --preset shiftTAG before per-flag overrides if explicit flags should win.

The production Chromap/MorPHiC lane remains the MACS3 FRAG path, not the Signac BED profile. The current MorPHiC bulk-ATAC wrapper drives Chromap inline peak calling with shifted fragment rows, memory-source FRAG peaks, --macs3-frag-pvalue 0.01, --macs3-frag-min-length 200, --macs3-frag-max-gap 30, human effective genome size, and symmetric Tn5-shifted fragments. Those settings should remain the default suite-facing behavior unless a benchmark/runbook explicitly opts into another profile.

For paper-baseline compatibility, the BED callpeak surface exposes ResolveMacs3BedCallPeakProfile("signac-atac", ...) in include/rapidmacs/macs3_bed_callpeak.h. Downstream consumers should use that resolver instead of copying profile constants when they need to emulate Signac. The profile resolves to BED input, --nomodel, --extsize 200, --shift -100, effective genome size 2700000000, q-value 0.05, --llocal 10000, and --keep-dup 1. It is opt-in compatibility behavior, not the library or suite default.

Final full PBMC 3k ATAC complete-CLI medians (53,969,811 fragments; three recorded warm-cache runs after one warm-up, measured 2026-08-01):

Tool	Threads	Wall	Peak RSS	Speedup vs MACS3
MACS3 3.0.3	1	743.07 s	2.397 GiB	1.00x
RapidMACS	1	59.14 s	2.802 GiB	12.56x
RapidMACS	4	24.97 s	3.080 GiB	29.76x
RapidMACS	24	15.69 s	4.053 GiB	47.36x

Phase B substrate changes (legacy custom caller only — NOT yet applied to the MACS3 FRAG pipeline):

• obs[] is uint16_t (saturating at 65535/bin); 4× smaller than the prior int64_t. Saturation is reported in *.summary.tsv (saturated_bins, saturated_bin_fraction).
• Local lambda is computed as a sliding-window stream rather than materialized into a lambda[n] array.
• BH FDR runs over the sparse subset of bins with obs >= 1; we never build p_values[n] / q_values[n].
• FileFragmentIterator is strict: malformed start / end / count produce an error (path:line:reason) rather than silently defaulting.

Open follow-ups (in rough order):

• Re-link libchromap against RapidMACS — done. The in-tree src/peak_caller/ was deleted and Chromap-suite consumes third_party/rapidmacs/lib/librapidmacs.a via a git submodule, using the canonical name.

• libchromap concurrent-path adapter — landed in Chromap-suite branch relink-libmacs3 (commit 94cbba6). mapping_writer.cc now buffers fragments as macs3::FragmentRecord (~16 B each) instead of workspace events (~48 B each), and the driver consumes the buffer via macs3::WrapVectorFragmentIterator → RunMacs3FragPeakPipelineFromSortedIterator. Validated byte-identical on the 100K integration harness: chromap --call-macs3-frag-peaks narrowPeak (md5 e5c8835..., 2047 peaks) matches standalone chromap_callpeaks. Awaits merge to Chromap-suite master and full PBMC 3k re-validation.

• Tolerance-based parity harness vs upstream MACS3 (Python) and MACS2 is available at tests/run_parity_vs_macs.sh for FRAG/BAMPE/BEDPE and optional control input.

• Namespace rename chromap::peaks → macs3 (mechanical).

• Full PBMC 3k re-validation through the streamed MACS3 FRAG path via libchromap (now that the adapter on relink-libmacs3 has landed).

Output parity with MACS3

The validated callpeak CLI paths are byte-identical to MACS3 3.0.3. This includes every narrowPeak field and every summit BED field (names, coordinates, signal, p-score, and q-score) on both the 2,047-peak 100K ATAC fixture and the 50,274-peak full PBMC 3k dataset. Controlled CUT&RUN and ChIP-seq narrowPeak and summit files are likewise byte-identical.

The audited ArchR BED command (--nolambda --call-summits) and pycisTopic BEDPE command (also --nolambda --call-summits) are byte-identical to MACS3 3.0.3 on the complete PBMC 100K workflow input: 12,923 and 9,181 rows, respectively. Reproduction commands and checksums are recorded in benchmarks/WORKFLOW-COMPATIBILITY.md.

The former ATAC q-score discrepancy came from rounding the in-memory local lambda and using non-MACS3 score-table keys. The callpeak path now retains MACS3's float32 lambda for Poisson scoring and uses its float32 p/q-table arithmetic; five-decimal rounding occurs only when a diagnostic bedGraph is written. Summit BED scores follow MACS3's selected score column: p-score for -p, q-score for -q.

Pass --no-macs3-frag-summit-parity for a lower-memory shortcut that diverges on ~2% of summits at PBMC scale (max ~222 bp; BED3 / bedGraphs / score columns unaffected).

Diagnosed via an instrumented build of MACS3 v3.0.3: MACS3's internal pos_array (the array iterated by __close_peak_wo_subpeaks to find the median-tied summit) is per-event, not per-distinct-value. When a fragment endpoint lands inside a treat plateau without changing the running pileup (one fragment ends and another starts at the same position with net ΔX=0), MACS3 still emits a new pos_array entry. The bedGraph writer collapses those adjacent same-value runs, so the exported .bdg matches RapidMACS's collapsed treat track byte-for-byte, but the median over peak_content is computed against more sub-pieces than the default RapidMACS path sees.

--macs3-frag-summit-parity reproduces MACS3's pos_array directly (see src/macs3_pos_array.cc): per chrom, build treat_pv positions via the quick_pileup start/end parallel- walk (with start==end coincidence skip) on count-duplicated fragment endpoints, build ctrl_pv positions via the se_all_in_one_pileup walk at d=llocal (10000 default for FRAG-without-control), then take the sorted union. The summit argmax then splits each tied plateau at every interior pos_array position before the median tie-pick.

Worked example on the worst-offender (chr4:76382069-76382561, 222 bp Δ under default; 0 bp Δ under --macs3-frag-summit-parity):

• bedGraph has 3 plateaus of treat=30 → default ties=3, midindex=1 → summit=76382470.
• MACS3 pos_array has 6 sub-pieces of treat=30 (first plateau split at 76382245, 76382246 from llocal-window endpoints of fragments at 76377245 and 76387246; second plateau split at 76382470) → ties=6, midindex=2 → summit=76382248. --macs3-frag-summit-parity matches this exactly.

Layout

include/rapidmacs/   public headers (consumers include "rapidmacs/<name>.h")
src/                library implementation (.cc files)
cli/                standalone CLI source (cli/rapidmacs.cc)
bin/                build output: rapidmacs
lib/                build output: librapidmacs.a
tests/
  golden/100k/      committed baseline narrowPeak from chromap_callpeaks
  unit/             small C++ unit tests
  run_100k_smoke.sh parity smoke vs golden baseline
  run_unit_tests.sh compiles + runs all tests/unit/*.cc
third_party/htslib  git submodule pinned at v1.23 (kfunc + BAM I/O)
benchmarks/         paper benchmark drivers and the MACS3 reference builds
  run_10x_pbmc_example.sh  single-dataset reproduction, checksum-verified
  results/          every recorded measurement behind the paper
  provenance/       the exact harness used for each recorded run
docs/archive/       superseded handoff documents (history, not normative)
plans/              development-history runbooks (history, not normative)

Reproducing the published benchmark

benchmarks/run_10x_pbmc_example.sh reproduces the single-cell ATAC row of Table 1 from a public dataset. It downloads the 10x Cell Ranger ARC PBMC 3k ATAC fragments, verifies the file against the SHA-256 the paper used, calls peaks with released MACS3 3.0.3 and with RapidMACS at 1, 4, and 24 threads, and exits non-zero unless every output file carries the published narrowPeak and summits MD5 sums:

make -j
benchmarks/setup_macs3_reference.sh VARIANT=released
MACS3=.macs3-reference/macs3-3.0.3-venv/bin/macs3 \
  benchmarks/run_10x_pbmc_example.sh

Expect roughly 25 minutes, nearly all of it the one MACS3 run. On the paper's host this configuration is 747.55 s for MACS3 against 57.87, 24.13, and 16.37 s for RapidMACS, all four producing 50,003 identical peaks. Set SKIP_MACS3=1 to time RapidMACS alone and check it against the published checksums, or WARMUPS=1 REPEATS=3 for the paper's full protocol rather than the single timed run the example defaults to.

The complete six-configuration matrix, both MACS3 builds, the recorded results, and the per-run provenance are described in benchmarks/README.md and benchmarks/results/README.md.

Building

git submodule update --init --recursive   # first time only
make                                       # produces lib/librapidmacs.a + bin/rapidmacs

The deprecated lib/libmacs3.a / bin/macs3frag symlinks and the include/libmacs3/*.h forwarding headers have been removed. Chromap-suite consumes the canonical lib/librapidmacs.a and bin/rapidmacs names via the third_party/rapidmacs submodule; no consumer references the old names.

By default the build links against the system-installed htslib (the same library Chromap-suite/chromap_callpeaks resolves to, currently v1.13 on this workstation). The submodule is pinned at v1.23 (matching Chromap-suite's vendored copy) for self-contained future builds; to build it and link against it instead:

(cd third_party/htslib && autoreconf -i && ./configure && make)
make HTSLIB_DIR=third_party/htslib

Paired-end examples

-f FRAG accepts both five-column Cell Ranger/Chromap barcoded fragments (chrom,start,end,barcode,count) and seven-column Chromap bulk fragment BED (chrom,start,end,N,mapq,strand,count). The duplicate count is read from column 5 or column 7 according to the detected schema. As in MACS3, FRAG duplicate counts are already collapsed: an explicit non-all --keep-dup value warns and is overridden to all.

# ChIP-seq: control is required by the chip preset.
bin/rapidmacs --preset chip -t ctcf.bam -c input.bam \
  -n ctcf --out-prefix out/ctcf

# CUT&RUN: matched IgG is optional.
bin/rapidmacs --preset cutrun -t ctcf_cutrun.bam -c igg.bam \
  -n ctcf_cutrun --out-prefix out/ctcf_cutrun

# Explicit standard BEDPE input (six mate columns) is also accepted.
bin/rapidmacs -f BEDPE -t fragments.bedpe -q 0.05 \
  -n bedpe --out-prefix out/bedpe

Testing

make test                    # unit + ATAC + ChIP-seq + CUT&RUN smokes
make test-unit               # C++ unit tests
make test-legacy-streaming-100k
make test-macs3-frag-100k    # existing FRAG golden
make test-chipseq            # public ENCODE BAMPE treatment/input fixture
make test-cutrun             # public GEO BAMPE CTCF/IgG fixture

The 100K smoke baseline is a snapshot of chromap_callpeaks output on pbmc_unsorted_3k_100k/pbmc_unsorted_3k_100k_chromap/fragments.tsv, committed at tests/golden/100k/golden.narrowPeak (md5 87f274f23fe635092bf1635418b98bc9, 121,592 peaks).

Override the fragments path:

FIXTURE_FRAGS=/path/to/fragments.tsv.gz make test-macs3-frag-100k

Phase A → Phase B

Once Phase A is fully closed (libchromap relinked, full PBMC 3k re-validated), the next work is the streaming sidecar (StreamingPeakCaller) — see roadmap §2. The FragmentIterator interface in include/rapidmacs/fragments.h is the boundary that streaming will plug into; today it is in place but only exercised by unit tests, not by the production pipeline.

License

RapidMACS is distributed under the MIT license — see LICENSE.

RapidMACS is an independent implementation, not a port. It was written clean-room from MACS3's published algorithmic specification and validated by output comparison against MACS3 v3.0.3; no MACS3 code is included, in original or translated form. MACS3's source was consulted only for the handful of low-level details needed to reach byte-for-byte equivalence — chiefly the per-event pos_array summit tie-break in src/macs3_pos_array.cc. That independence is also why RapidMACS is faster: a specialized paired-end narrow-peak codebase can pick representations a general caller cannot, rather than inheriting MACS3's structure at a faster language's speed.

NOTICE records the acknowledgement of MACS3 and the htslib attribution. MACS was originally developed in X. Shirley Liu's group and has been maintained and extended since by Tao Liu and the MACS3 contributors; this project would not have been possible without that work or the permissive licence it carries. htslib is linked, not vendored, and carries its own MIT license.