feat: add ft union-peaks to merge peak calls from many BED files

[mrvollger]

Combines ft mock-fire and ft call-peaks into one command that takes any number of BED files and returns a single union peak set.

ft union-peaks sample1.bed sample2.bed sample3.bed -o union.bed

Each input BED is one sample: its intervals become FIRE elements on mock fibers, which go straight into the existing peak caller in memory. Doing this by hand needs an intermediate coordinate sort and index, which ft mock-fire cannot produce on its own (it writes in read-name order, so samtools index fails with Unsorted positions whenever read names sort against positions).

Output

#chrom  start  end    name          n_support  frac_support  support   union_start  union_end  peak_summit
chr1    1020   1220   union_peak_1  3          1.0000        s1,s2,s3  1000         1250       1125
chr1    5050   5350   union_peak_2  2          0.6667        s1,s2     5000         5350       5175

start/end are the peak caller's consensus (median) boundaries of the overlapping input intervals; union_start/union_end are their outer span. Overlapping intervals within one file are collapsed first, so a file can add at most 1 to n_support. --min-support filters the output only, so -n 3 and -n 1 plus a downstream filter agree.

Flags: -o/--out, --names (default: file basenames), -n/--min-support, --min-frac-support, and the shared peak-caller knobs (--window-size plus the merge thresholds), exposed through the same PeakMergeParams struct call-peaks uses, with identical defaults.

Refactor

To share the peak caller, call_peaks_for_chrom is extracted out of call_peaks and the CLI values it reads move into a PeakCallingParams struct. CallPeaksOptions flattens InputBam, which cannot be built by a command that never opens a BAM, and whose FiberFilters fields carry env= defaults (MIN_MSP, MIN_AVE_MSP_SIZE, FT_MIN_ML_SCORE) that would otherwise leak into union peak calls.

Peaks are emitted through a callback, so Peak and its Display impl are untouched and ft call-peaks output cannot drift.

Gate: ft call-peaks --min-fire-frac 0.5 output is byte-identical before and after the refactor on tests/data/{ctcf,all,NAPA}.bam.

The pileup window is sized to islands of input intervals rather than whole chromosomes; the track costs ~56 bytes a base, which is tens of GB genome-wide. On the test fixture, islanded output matches a hand-run mock-fire | samtools sort | samtools index | call-peaks pipeline exactly, except where that pipeline double-counts one file's overlapping intervals.

Tests

5 regression tests in tests/regression/union_peaks.rs (support counts and boundaries, within-file collapse, --min-support, names ordered against positions, empty input, error paths) plus 2 unit tests for islands/support_for. cargo fmt, cargo clippy -p fibertools-rs, and cargo test -p fibertools-rs are clean.

mock-fire deprecation

ft mock-fire is now hidden from help and warns at runtime that it will be removed: union-peaks replaces the mock-fire | samtools sort | call-peaks pipeline (the one user found in the wild, StergachisLab/fire_consensus_pipeline, keeps working meanwhile — its --min-fire-coverage 1 --min-fire-frac 0.01 maps to union-peaks' built-ins plus --min-support). The module code stays as union-peaks' engine.

Refactor gate, made permanent

ft call-peaks output on the ctcf fixture is now snapshot-pinned (tests/regression/call_peaks.rs), so the byte-identical guarantee for the shared call_peaks_for_chrom path is enforced by CI rather than a one-off manual check.

Deep-review outcome

A multi-agent review compared this against the legacy pipeline empirically: on ctcf.bed.gz, mock-fire | sort | call-peaks --min-fire-coverage 1 --min-fire-frac 0.01 and ft union-peaks produce identical peak sets (boundaries and summits byte-equal, legacy fire_coverage == n_support in the multi-sample case), deterministically, at ~0.01 s vs ~10 s. Every divergence is intended and documented except two that this PR now fixes: exactly book-ended intervals from different samples no longer drop the touching sample from support (fixed at the root too: consensus boundaries now use a true median — even element counts take the midpoint of the two middle values instead of the upper median, shifting 10 of 31 ctcf call-peaks boundaries by 1-33 bp with summits and all statistics unchanged), and out-of-range coordinates (past int32 header / u32 liftover space) error cleanly instead of vanishing. Also added: a warning when an island exceeds 10 Mb (memory scales with island size), --min-support misuse guards, and track/browser line tolerance in the BED reader.

Migration note for fire_consensus_pipeline: union-peaks derives the sample from the input file and ignores the BED 4th column — passing the pipeline's single merged.4col.bed.gz gives one sample (n_support ≡ 1). Split per-sample BEDs to migrate. Its --min-fire-frac 0.01 maps directly to the new --min-frac-support (the fractional twin of -n; both fold into one threshold); residual differences appear only past ~100 samples (the 1% floor in the peak caller itself).

[mrvollger]

Adversarial review + scale test

15 agents across 4 lenses (correctness, refactor safety, CLI/docs, genome-scale perf), each finding independently refuted before being kept. 10 candidate findings, 9 refuted, 1 survived.

Fixed in 8980f7f

Per-island interval selection rescanned every interval of every sample for every island — islands x intervals per chromosome, ~77% of the run on dense input (profiler: 79.8% of main-thread stacks in that block vs 10.3% in call_peaks_for_chrom). Now binary-searches, mirroring support_for.

Input	Intervals	Islands	Before	After
20 samples, spread, 1 chrom	2,000,000	100,000	151.8 s	19.3 s
20 samples, whole genome (hg38)	3,000,000	—	269.9 s	237.9 s

Output byte-identical before/after on all three benchmark sets (genome: 2,071,199 peak rows). Peak RSS flat at ~340 MB genome-wide — islanding is a ~1,780x memory reduction over the whole-chromosome pileup it wraps.

Also corrected four docs the review found overpromising: islanding bounds the rolling-max failure mode but is not identical to a whole-chromosome run; peak_summit can fall outside start/end (same median-boundary reason as peak_max in call-peaks); only local maxima become peaks, so at most one per --window-size; --min-support renumbers the sequential name column.

Refuted (checked, not real)

• Islanding not result-preserving — divergence traces to pre-existing ft mock-fire fibers spanning min-start..max-end; 0/60 fuzz differences once that artifact is removed.
• High-support peaks lost to a lone sample's nearby interval — mechanism false (isolated all-FIRE intervals score 99.0, the lowest all-FIRE score); 0/291 fuzz cases, 0/60 realistic fixtures.
• peak_summit outside [start, end) — pre-existing call-peaks behavior, reproducible with the pre-PR binary; the value is the true argmax, and clamping would center metaplots on a boundary pixel.
• Lexicographic chromosome order breaks bedtools -sorted — the repro used a hand-written karyotypic -g; it reproduces with bedtools sort's own output.
• MAX_ISLAND_LEN unreachable on human data — intentional; it guards the 28-bit Cigar::Equal wraparound that ft mock-fire itself lacks.
• Plus: --min-support renumbering, --help not documenting medians, and the deliberate simplifications (hardcoded merge thresholds, no --min-frac-support, callback instead of an owned PeakRecord).

Refactor safety was re-checked on the paths the original gate did not cover (FDR mode, --shuffled, --fdr-table-out, --haps, explicit coverage bounds) against a git worktree build of the parent commit.

[mrvollger]

Depends on #130. Rebase this branch after #130 merges.

Why it depends. Both branches change the same three files: src/subcommands/mock_fire.rs, src/subcommands/call_peaks/mod.rs, src/subcommands/call_peaks/peaks.rs. The conflict is one line: call_peaks_for_chrom sets fiber_coverage: true (peaks.rs:658), and #130 renames that field to callable_coverage.

Two things to check after the rebase.

1. feat!: fiberseq_callable MA tag + filtered/unfiltered ft qc columns #130 drops NotCallable fibers in call-peaks. Mock fibers have no m6A and often fewer than 10 MSPs in one island, so they would all be dropped and union-peaks would emit nothing. It works only because feat!: fiberseq_callable MA tag + filtered/unfiltered ft qc columns #130 stamps mock records full-width Callable in create_mock_fire_record. Add a union-peaks regression test so a later change to that stamp does not silently empty the output.

2. feat!: fiberseq_callable MA tag + filtered/unfiltered ft qc columns #130 makes the summit fix possible. Coverage is the denominator of the FIRE score, and today it is the fiber MSP extent, which for a mock fiber is the first to last peak in its island. That makes the summit depend on how we chunked the genome (measured: 11 summit shifts and 1 peak-set change in 60 random configs). After feat!: fiberseq_callable MA tag + filtered/unfiltered ft qc columns #130 coverage is the alignment span instead, so pinning each island mock record to the island bounds gives a constant denominator of N samples. The score then becomes proportional to sample support, and the summit is a maximum-support position by definition. Do this as a follow-up on top of feat!: fiberseq_callable MA tag + filtered/unfiltered ft qc columns #130.