From c5696c106659404b003c8c27dac59d5f1bc3aa4e Mon Sep 17 00:00:00 2001 From: "Mitchell R. Vollger" Date: Wed, 19 Aug 2026 13:25:32 -0600 Subject: [PATCH 01/11] feat: add `ft union-peaks` to merge peak calls from many BED files Running mock-fire and call-peaks by hand needs an intermediate coordinate sort and index, which mock-fire cannot produce (it writes in read-name order), so combine the two into one command that takes N BED files samtools-merge style. Each input BED is one sample: its intervals become FIRE elements on mock fibers and go straight into the peak caller in memory. Peaks are reported with the number and names of the input BEDs that support them. To share the peak caller, extract `call_peaks_for_chrom` out of `call_peaks` and move the CLI values it reads into `PeakCallingParams`, so a command that never opens a BAM can call it (`CallPeaksOptions` flattens `InputBam`, whose fields also carry env-var defaults that would leak in). `ft call-peaks` output is byte-identical before and after on ctcf.bam, all.bam, and NAPA.bam. The window is sized to islands of input intervals instead of whole chromosomes: the pileup track costs ~56 bytes a base, which is tens of GB genome-wide. --- .../docs/mm-ml-per-group-passthrough.md | 2 +- src/cli.rs | 7 + src/cli/union_peaks_opts.rs | 27 ++ src/main.rs | 3 + src/subcommands.rs | 2 + src/subcommands/call_peaks/mod.rs | 4 +- src/subcommands/call_peaks/peaks.rs | 199 ++++++--- src/subcommands/mock_fire.rs | 4 +- src/subcommands/union_peaks.rs | 383 ++++++++++++++++++ tests/regression.rs | 2 + tests/regression/union_peaks.rs | 136 +++++++ 11 files changed, 698 insertions(+), 71 deletions(-) create mode 100644 src/cli/union_peaks_opts.rs create mode 100644 src/subcommands/union_peaks.rs create mode 100644 tests/regression/union_peaks.rs diff --git a/molecular-annotation/docs/mm-ml-per-group-passthrough.md b/molecular-annotation/docs/mm-ml-per-group-passthrough.md index 785b5f719..20c50f94a 100644 --- a/molecular-annotation/docs/mm-ml-per-group-passthrough.md +++ b/molecular-annotation/docs/mm-ml-per-group-passthrough.md @@ -18,7 +18,7 @@ on intent: `FiberseqData::serialize_annotations`. - `write_record_with_basemods` / `write_mm_ml` — destructively re-encodes MM/ML canonically from the model. Used by base-mod producers (predict-m6a, - ddda-to-m6a, strip-basemods, mock-fire, fibertig synthesis). + ddda-to-m6a, strip-basemods, mock-fire, union-peaks, fibertig synthesis). The split is a proxy for one boolean: "did this code path change base mods?" It works, but the decision lives in the caller's choice of function — a footgun diff --git a/src/cli.rs b/src/cli.rs index 34644a493..0c2d0975e 100644 --- a/src/cli.rs +++ b/src/cli.rs @@ -23,6 +23,7 @@ mod pileup_opts; mod predict_opts; mod qc_opts; mod strip_basemods_opts; +mod union_peaks_opts; mod validate_opts; // include the subcommand modules as top level functions and structs in the cli module @@ -45,6 +46,7 @@ pub use pileup_opts::*; pub use predict_opts::*; pub use qc_opts::*; pub use strip_basemods_opts::*; +pub use union_peaks_opts::*; pub use validate_opts::ValidateOptions; // @@ -165,6 +167,11 @@ pub enum Commands { /// Each interval in the BED becomes a FIRE element. The 4th column groups intervals into the same mock read. #[clap(name = "mock-fire")] MockFire(MockFireOptions), + /// Combine peak calls from many BED files into one union peak set. + /// Each input BED is one sample, and peaks are reported with the number and names of + /// the input BEDs that support them. + #[clap(name = "union-peaks", visible_aliases = &["union"])] + UnionPeaks(UnionPeaksOptions), /// Benchmark fiber iterator performance (hidden command for testing) #[clap(hide = true)] Benchmark(BenchmarkOptions), diff --git a/src/cli/union_peaks_opts.rs b/src/cli/union_peaks_opts.rs new file mode 100644 index 000000000..be7a040e5 --- /dev/null +++ b/src/cli/union_peaks_opts.rs @@ -0,0 +1,27 @@ +use crate::cli::GlobalOpts; +use clap::Args; +use std::fmt::Debug; + +#[derive(Args, Debug)] +pub struct UnionPeaksOptions { + /// Input BED files, one per sample. + /// Every interval becomes a FIRE element on a mock fiber for that sample, and + /// peaks are called across all the samples at once. Overlapping intervals within + /// one file are merged first, so a file can add at most 1 to a peak's support. + #[clap(required = true, num_args = 1..)] + pub beds: Vec, + /// Output BED file with union peaks + #[clap(short, long, default_value = "-")] + pub out: String, + /// Sample names, comma separated, one per input BED [default: input file basenames] + #[clap(long, value_delimiter = ',')] + pub names: Vec, + /// Minimum number of input BEDs that must overlap a peak for it to be reported + #[clap(short = 'n', long, default_value_t = 1)] + pub min_support: usize, + /// Rolling window size for finding local maxima (in base pairs) + #[clap(long, default_value_t = 200)] + pub window_size: usize, + #[clap(flatten)] + pub global: GlobalOpts, +} diff --git a/src/main.rs b/src/main.rs index 652330526..c62ffa0d4 100644 --- a/src/main.rs +++ b/src/main.rs @@ -138,6 +138,9 @@ pub fn main() -> Result<(), Error> { Some(Commands::MockFire(mock_fire_opts)) => { subcommands::mock_fire::run_mock_fire(mock_fire_opts)?; } + Some(Commands::UnionPeaks(union_peaks_opts)) => { + subcommands::union_peaks::run_union_peaks(union_peaks_opts)?; + } Some(Commands::Benchmark(benchmark_opts)) => { subcommands::benchmark::run_benchmark(benchmark_opts)?; } diff --git a/src/subcommands.rs b/src/subcommands.rs index ec233ea21..5691a8a2b 100644 --- a/src/subcommands.rs +++ b/src/subcommands.rs @@ -26,6 +26,8 @@ pub mod predict_m6a; pub mod qc; /// Remove base modifications from a bam record pub mod strip_basemods; +/// Call union peaks across many BED files +pub mod union_peaks; /// Create mock BAM from reference FASTA pub mod pg_inject; diff --git a/src/subcommands/call_peaks/mod.rs b/src/subcommands/call_peaks/mod.rs index 6eb94d728..ee4616b81 100644 --- a/src/subcommands/call_peaks/mod.rs +++ b/src/subcommands/call_peaks/mod.rs @@ -5,7 +5,9 @@ pub use fdr::{ fdr_table, lookup_fdr, read_fdr_table, write_fdr_table, FdrEntry, IncrementalFdrBuilder, PileupRecord, }; -pub use peaks::{call_peaks, reciprocal_overlap_raw}; +pub use peaks::{ + call_peaks, call_peaks_for_chrom, reciprocal_overlap_raw, Peak, PeakCallingParams, +}; use crate::cli::CallPeaksOptions; use crate::subcommands::pileup::{FireTrack, FireTrackOptions}; diff --git a/src/subcommands/call_peaks/peaks.rs b/src/subcommands/call_peaks/peaks.rs index 3a9dbfc83..7d7281bbe 100644 --- a/src/subcommands/call_peaks/peaks.rs +++ b/src/subcommands/call_peaks/peaks.rs @@ -1,6 +1,7 @@ use super::chrom_names_and_lengths; use super::fdr::{lookup_fdr, FdrEntry}; use crate::cli::CallPeaksOptions; +use crate::fiber::FiberseqData; use crate::subcommands::pileup::{ FiberseqPileup, FiberseqPileupOptions, FireTrack, FireTrackOptions, }; @@ -34,6 +35,44 @@ pub fn reciprocal_overlap_raw( (overlap_len / a_len).min(overlap_len / b_len) } +/// Everything peak calling needs from the CLI, decoupled from `CallPeaksOptions` +/// (which flattens `InputBam`/`FiberFilters`, so commands that never read a BAM +/// cannot build one). See `ft union-peaks` for the other caller. +#[derive(Debug, Clone, Copy)] +pub struct PeakCallingParams { + pub window_size: usize, + pub min_fire_coverage: i32, + pub min_cov: Option, + pub max_cov: Option, + pub sd_cov: f64, + pub max_fdr: f64, + pub min_fire_frac: Option, + pub min_fire_frac_filter: f64, + pub min_frac_overlap: f64, + pub min_reciprocal_overlap: f64, + pub high_reciprocal_overlap: f64, + pub max_grouping_iterations: usize, +} + +impl From<&CallPeaksOptions> for PeakCallingParams { + fn from(o: &CallPeaksOptions) -> Self { + Self { + window_size: o.window_size, + min_fire_coverage: o.min_fire_coverage, + min_cov: o.min_cov, + max_cov: o.max_cov, + sd_cov: o.sd_cov, + max_fdr: o.max_fdr, + min_fire_frac: o.min_fire_frac, + min_fire_frac_filter: o.min_fire_frac_filter, + min_frac_overlap: o.min_frac_overlap, + min_reciprocal_overlap: o.min_reciprocal_overlap, + high_reciprocal_overlap: o.high_reciprocal_overlap, + max_grouping_iterations: o.max_grouping_iterations, + } + } +} + /// Filtering thresholds for peak calling #[derive(Debug, Clone, Copy)] struct PeakThresholds { @@ -518,17 +557,20 @@ fn merge_peaks_single_iteration<'a>( /// Phase 1: High reciprocal overlap (default 90%, configurable via --high-reciprocal-overlap) /// Phase 2: FIRE element overlap (default 50%, configurable via --min-frac-overlap) /// Phase 3: Reciprocal overlap (default 75%, configurable via --min-reciprocal-overlap) -fn merge_peaks_iterative<'a>(mut peaks: Vec>, opts: &CallPeaksOptions) -> Vec> { +fn merge_peaks_iterative<'a>( + mut peaks: Vec>, + params: &PeakCallingParams, +) -> Vec> { let initial_count = peaks.len(); // Phase 1: High reciprocal overlap log::debug!( " Phase 1: Merging peaks with reciprocal overlap >= {}", - opts.high_reciprocal_overlap + params.high_reciprocal_overlap ); - for iteration in 0..opts.max_grouping_iterations { + for iteration in 0..params.max_grouping_iterations { let prev_count = peaks.len(); - peaks = merge_peaks_single_iteration(peaks, 0.0, opts.high_reciprocal_overlap); + peaks = merge_peaks_single_iteration(peaks, 0.0, params.high_reciprocal_overlap); log::debug!( " Iteration {}: {} -> {} peaks", iteration + 1, @@ -544,11 +586,11 @@ fn merge_peaks_iterative<'a>(mut peaks: Vec>, opts: &CallPeaksOptions) // Phase 2: FIRE element overlap log::debug!( " Phase 2: Merging peaks with FIRE element overlap >= {}", - opts.min_frac_overlap + params.min_frac_overlap ); - for iteration in 0..opts.max_grouping_iterations { + for iteration in 0..params.max_grouping_iterations { let prev_count = peaks.len(); - peaks = merge_peaks_single_iteration(peaks, opts.min_frac_overlap, 0.0); + peaks = merge_peaks_single_iteration(peaks, params.min_frac_overlap, 0.0); log::debug!( " Iteration {}: {} -> {} peaks", iteration + 1, @@ -564,11 +606,11 @@ fn merge_peaks_iterative<'a>(mut peaks: Vec>, opts: &CallPeaksOptions) // Phase 3: High reciprocal overlap again log::debug!( " Phase 3: Merging peaks with reciprocal overlap >= {}", - opts.min_reciprocal_overlap + params.min_reciprocal_overlap ); - for iteration in 0..opts.max_grouping_iterations { + for iteration in 0..params.max_grouping_iterations { let prev_count = peaks.len(); - peaks = merge_peaks_single_iteration(peaks, 0.0, opts.min_reciprocal_overlap); + peaks = merge_peaks_single_iteration(peaks, 0.0, params.min_reciprocal_overlap); log::debug!( " Iteration {}: {} -> {} peaks", iteration + 1, @@ -591,6 +633,75 @@ fn merge_peaks_iterative<'a>(mut peaks: Vec>, opts: &CallPeaksOptions) peaks } +/// Call and merge peaks over one reference window using an in-memory fiber stream. +/// +/// `chrom_start`/`chrom_end` bound the pileup track, so a caller can size the track +/// to its data instead of to the whole chromosome (`ft call-peaks` passes the whole +/// chromosome, `ft union-peaks` passes one island of BED intervals). `emit` is called +/// once per merged peak in output order. Returns `(peaks_before_merge, peaks_after_merge)`. +pub fn call_peaks_for_chrom( + chrom: &str, + chrom_start: usize, + chrom_end: usize, + fibers: impl Iterator, + params: &PeakCallingParams, + fdr_table: &[FdrEntry], + mut emit: impl FnMut(&Peak) -> Result<()>, +) -> Result<(usize, usize)> { + // Create FiberseqPileupOptions for peak calling + let pileup_opts = FiberseqPileupOptions { + fire_track_opts: FireTrackOptions { + no_nuc: false, + no_msp: false, + m6a: false, + cpg: false, + callable_fibers: true, + shuffle: false, + random_shuffle: false, + shuffle_seed: None, + rolling_max: Some(params.window_size), + track_fire_elements: true, // Enable FIRE element tracking for peak calling + }, + rolling_max: Some(params.window_size), + haps: false, + per_base: false, + keep_zeros: false, + min_fire_coverage: Some(params.min_fire_coverage), + }; + + // Process fibers to build the track + let mut pileup = FiberseqPileup::new(chrom, chrom_start, chrom_end, pileup_opts, &None, None); + pileup.add_fibers(fibers); + + // Calculate coverage thresholds + let (median, std_dev, _) = pileup.all_data.median_and_std_coverage(); + let min_cov = params.min_cov.unwrap_or_else(|| { + let calculated_min = (median - params.sd_cov * std_dev).round() as i32; + calculated_min.max(4) // DEFAULT_MIN_COVERAGE = 4 + }); + let max_cov = params + .max_cov + .unwrap_or_else(|| (median + params.sd_cov * std_dev).round() as i32); + + // Find local maxima and filter by threshold (FDR or FIRE fraction) + let peaks = Peak::from_pileup( + &pileup, + fdr_table, + params.max_fdr, + params.min_fire_frac, + params.min_fire_frac_filter, + min_cov, + max_cov, + ); + + let peaks_before = peaks.len(); + let merged_peaks = merge_peaks_iterative(peaks, params); + for peak in &merged_peaks { + emit(peak)?; + } + Ok((peaks_before, merged_peaks.len())) +} + /// Call peaks using FDR table or FIRE fraction filtering /// /// This will: @@ -626,6 +737,7 @@ pub fn call_peaks( let mut writer = bio_io::writer(&opts.out)?; writeln!(writer, "{}", Peak::header())?; + let params = PeakCallingParams::from(&*opts); let mut total_peaks_before_merge = 0; let mut total_peaks_after_merge = 0; @@ -642,60 +754,18 @@ pub fn call_peaks( chrom_len ); - // Create FiberseqPileupOptions for peak calling - let pileup_opts = FiberseqPileupOptions { - fire_track_opts: FireTrackOptions { - no_nuc: false, - no_msp: false, - m6a: false, - cpg: false, - callable_fibers: true, - shuffle: false, - random_shuffle: false, - shuffle_seed: None, - rolling_max: Some(opts.window_size), - track_fire_elements: true, // Enable FIRE element tracking for peak calling - }, - rolling_max: Some(opts.window_size), - haps: false, - per_base: false, - keep_zeros: false, - min_fire_coverage: Some(opts.min_fire_coverage), - }; - - // Process fibers to build the track - let mut pileup = - FiberseqPileup::new(&chrom, 0, chrom_len as usize, pileup_opts, &None, None); let fibers = opts.input.fetch_fibers(bam, &chrom, None, None)?; - pileup.add_fibers(fibers); - - // Calculate coverage thresholds - let (median, std_dev, _) = pileup.all_data.median_and_std_coverage(); - let min_cov = opts.min_cov.unwrap_or_else(|| { - let calculated_min = (median - opts.sd_cov * std_dev).round() as i32; - calculated_min.max(4) // DEFAULT_MIN_COVERAGE = 4 - }); - let max_cov = opts - .max_cov - .unwrap_or_else(|| (median + opts.sd_cov * std_dev).round() as i32); - - // Find local maxima and filter by threshold (FDR or FIRE fraction) - let peaks = Peak::from_pileup( - &pileup, + let (peaks_before, peaks_after) = call_peaks_for_chrom( + &chrom, + 0, + chrom_len as usize, + fibers, + ¶ms, fdr_table, - opts.max_fdr, - opts.min_fire_frac, - opts.min_fire_frac_filter, - min_cov, - max_cov, - ); - - let peaks_before = peaks.len(); + |peak| Ok(writeln!(writer, "{}", peak)?), + )?; total_peaks_before_merge += peaks_before; - - // Merge peaks for this chromosome - let merged_peaks = merge_peaks_iterative(peaks, opts); - total_peaks_after_merge += merged_peaks.len(); + total_peaks_after_merge += peaks_after; // Summary info statement per chromosome log::info!( @@ -703,13 +773,8 @@ pub fn call_peaks( chrom, chrom_len / 1_000_000, peaks_before, - merged_peaks.len(), + peaks_after, ); - - // Write merged peaks for this chromosome - for peak in &merged_peaks { - writeln!(writer, "{}", peak)?; - } } log::info!("Total peaks before merging: {}", total_peaks_before_merge); diff --git a/src/subcommands/mock_fire.rs b/src/subcommands/mock_fire.rs index 4e998462d..4291dc52e 100644 --- a/src/subcommands/mock_fire.rs +++ b/src/subcommands/mock_fire.rs @@ -33,7 +33,7 @@ fn group_bed_by_name_and_chrom( } /// Create a BAM header from BED records -fn create_header_from_bed(bed_records: &[BedRecord]) -> Header { +pub(crate) fn create_header_from_bed(bed_records: &[BedRecord]) -> Header { let mut header = Header::new(); // Collect unique chromosomes and their max positions @@ -61,7 +61,7 @@ fn create_header_from_bed(bed_records: &[BedRecord]) -> Header { } /// Create a mock BAM record with FIRE elements -fn create_mock_fire_record( +pub(crate) fn create_mock_fire_record( read_name: &str, intervals: &[BedRecord], header_view: &HeaderView, diff --git a/src/subcommands/union_peaks.rs b/src/subcommands/union_peaks.rs new file mode 100644 index 000000000..79f46c4d9 --- /dev/null +++ b/src/subcommands/union_peaks.rs @@ -0,0 +1,383 @@ +use crate::cli::UnionPeaksOptions; +use crate::fiber::FiberseqData; +use crate::subcommands::call_peaks::{call_peaks_for_chrom, PeakCallingParams}; +use crate::subcommands::mock_fire::{create_header_from_bed, create_mock_fire_record}; +use crate::utils::bio_io::{self, read_bed_regions, BedRecord}; +use crate::utils::input_bam::FiberFilters; +use anyhow::{bail, Context, Result}; +use rust_htslib::bam::HeaderView; +use std::collections::{BTreeMap, HashSet}; +use std::io::Write; + +/// Quality given to every synthesized FIRE element. Must stay at or above the pileup's +/// `MIN_FIRE_QUAL` (229) or the elements are silently ignored. +const MOCK_FIRE_QUALITY: u8 = 255; + +/// Slack added to `--window-size` when grouping intervals into islands. The gap has to +/// exceed the rolling-max window, otherwise a window could span two islands and the +/// islanded result would differ from a whole-chromosome one. +const ISLAND_PAD: i64 = 1000; + +/// Longest island we will build a mock fiber for. `Cigar::Equal(len)` packs the length +/// into 28 bits and wraps silently past that. +const MAX_ISLAND_LEN: i64 = 1 << 28; + +/// One input BED, collapsed to non-overlapping intervals per chromosome. +struct Sample { + name: String, + by_chrom: BTreeMap>, +} + +/// Sample name for a BED path: the basename with compression and BED-ish suffixes removed. +fn sample_name_from_path(path: &str) -> String { + let mut name = std::path::Path::new(path) + .file_name() + .map_or_else(|| path.to_string(), |f| f.to_string_lossy().into_owned()); + for suffix in [".gz", ".bgz"] { + if let Some(stripped) = name.strip_suffix(suffix) { + name = stripped.to_string(); + break; + } + } + for suffix in [".bed", ".narrowPeak", ".broadPeak", ".tsv", ".txt"] { + if let Some(stripped) = name.strip_suffix(suffix) { + name = stripped.to_string(); + break; + } + } + name +} + +/// One sample name per input BED, from `--names` or the file basenames. +fn resolve_names(opts: &UnionPeaksOptions) -> Result> { + let names: Vec = if opts.names.is_empty() { + opts.beds.iter().map(|b| sample_name_from_path(b)).collect() + } else { + if opts.names.len() != opts.beds.len() { + bail!( + "--names has {} names but {} BED files were given", + opts.names.len(), + opts.beds.len() + ); + } + opts.names.clone() + }; + // Names go into the comma joined support column and into BAM read names, so a comma + // or whitespace in a name corrupts the output rather than just looking odd. + for name in &names { + if name.is_empty() || name.contains(',') || name.split_whitespace().count() != 1 { + bail!("sample name {name:?} is empty or has a comma or whitespace in it; pass clean names with --names"); + } + } + let mut seen = HashSet::new(); + for name in &names { + if !seen.insert(name) { + bail!("duplicate sample name {name:?}; pass unique names with --names"); + } + } + Ok(names) +} + +/// Read one BED into per-chromosome sorted, non-overlapping intervals. +fn load_sample(path: &str, name: String) -> Result { + let mut by_chrom: BTreeMap> = BTreeMap::new(); + + // A peak caller that found nothing leaves a zero byte file, which the decompression + // layer rejects with an opaque "file is too short" error. Take it as a sample with no + // peaks instead of failing the whole union. Only regular files, so a fifo or a + // process substitution still goes down the normal read path. + if std::fs::metadata(path).is_ok_and(|m| m.is_file() && m.len() == 0) { + log::warn!("{path} is empty, sample {name} will support no peaks"); + return Ok(Sample { name, by_chrom }); + } + + for rec in read_bed_regions(path).with_context(|| format!("failed to read BED file {path}"))? { + // read_bed_regions parses the coordinates as bare i64s, and a negative start makes + // create_mock_fire_record emit a record whose positions are all dropped later. + if rec.start < 0 || rec.end <= rec.start { + bail!( + "invalid interval in {}: {} {} {}", + path, + rec.chrom, + rec.start, + rec.end + ); + } + by_chrom + .entry(rec.chrom) + .or_default() + .push((rec.start, rec.end)); + } + + // Collapse overlaps within a file, so one sample can only ever add 1 to the support + // count of a peak. + for intervals in by_chrom.values_mut() { + intervals.sort_unstable(); + let mut merged: Vec<(i64, i64)> = Vec::with_capacity(intervals.len()); + for &(start, end) in intervals.iter() { + match merged.last_mut() { + Some(last) if start <= last.1 => last.1 = last.1.max(end), + _ => merged.push((start, end)), + } + } + *intervals = merged; + } + + Ok(Sample { name, by_chrom }) +} + +/// Group sorted intervals into maximal runs separated by no more than `gap` bases. +/// +/// Peaks are called one island at a time so the pileup track is sized to the data. A +/// whole-chromosome track costs ~56 bytes a base, which is tens of GB for a genome-wide +/// peak union. +fn islands(sorted: &[(i64, i64)], gap: i64) -> Vec<(i64, i64)> { + let mut out: Vec<(i64, i64)> = Vec::new(); + for &(start, end) in sorted { + match out.last_mut() { + Some(last) if start - last.1 <= gap => last.1 = last.1.max(end), + _ => out.push((start, end)), + } + } + out +} + +/// Samples with an interval overlapping `[start, end)`, in input order, plus the outer +/// span of those intervals. +/// +/// The pileup cannot answer this: its FIRE elements carry no fiber identity. +fn support_for<'a>( + samples: &'a [Sample], + chrom: &str, + start: i64, + end: i64, +) -> (Vec<&'a str>, i64, i64) { + let mut names = Vec::new(); + let (mut union_start, mut union_end) = (start, end); + for sample in samples { + let Some(intervals) = sample.by_chrom.get(chrom) else { + continue; + }; + let mut hit = false; + let first = intervals.partition_point(|iv| iv.1 <= start); + for &(iv_start, iv_end) in &intervals[first..] { + if iv_start >= end { + break; + } + hit = true; + union_start = union_start.min(iv_start); + union_end = union_end.max(iv_end); + } + if hit { + names.push(sample.name.as_str()); + } + } + (names, union_start, union_end) +} + +/// Merge peak calls from many BED files into one union peak set. +/// +/// Each input BED is one sample. Its intervals are synthesized into the same mock FIRE +/// records `ft mock-fire` writes, and those go straight into the `ft call-peaks` +/// machinery, so this is `ft mock-fire` piped into `ft call-peaks` without the +/// intermediate sort and index (which the mock BAM needs and which `ft mock-fire` alone +/// cannot satisfy, since it writes in read-name order). +/// +/// Reported `start`/`end` are therefore the peak caller's consensus (median) boundaries +/// of the overlapping input intervals, not their outer span; the outer span is reported +/// alongside as `union_start`/`union_end`. `--min-support` filters the output only, so +/// `-n 3` and `-n 1` plus a downstream filter agree. +pub fn run_union_peaks(opts: &UnionPeaksOptions) -> Result<()> { + if opts.window_size < 2 { + bail!("--window-size must be at least 2; a smaller window finds no local maxima"); + } + let names = resolve_names(opts)?; + + // Load every sample up front: islands are built across all of them at once. + let mut samples: Vec = Vec::with_capacity(opts.beds.len()); + let mut chrom_lengths: BTreeMap = BTreeMap::new(); + for (bed, name) in opts.beds.iter().zip(names) { + log::info!("Reading BED file: {bed}"); + let sample = load_sample(bed, name)?; + for (chrom, intervals) in &sample.by_chrom { + let max_end = intervals.last().map_or(0, |iv| iv.1); + let entry = chrom_lengths.entry(chrom.clone()).or_insert(0); + *entry = (*entry).max(max_end); + } + samples.push(sample); + } + if chrom_lengths.is_empty() { + bail!( + "no intervals found in any of the {} input BED files", + opts.beds.len() + ); + } + + // The header exists only so the mock records can resolve tids and name their targets. + let header_records: Vec = chrom_lengths + .iter() + .map(|(chrom, end)| BedRecord { + chrom: chrom.clone(), + start: 0, + end: *end, + name: None, + extra_fields: vec![], + }) + .collect(); + let header_view = HeaderView::from_header(&create_header_from_bed(&header_records)); + + // Mock fibers have no background to estimate an FDR from (a shuffled control would be + // built from the same synthetic fibers), so call in FIRE-fraction mode with the + // fraction threshold off and filter on sample support afterwards instead. + let params = PeakCallingParams { + window_size: opts.window_size, + min_fire_coverage: 1, // one sample is enough to score a position + min_cov: Some(1), // coverage bounds only set pass_coverage, which we do not emit + max_cov: Some(i32::MAX), + sd_cov: 5.0, + max_fdr: 1.0, + min_fire_frac: Some(0.0), + min_fire_frac_filter: 0.0, + min_frac_overlap: 0.5, + min_reciprocal_overlap: 0.75, + high_reciprocal_overlap: 0.90, + max_grouping_iterations: 10, + }; + + let mut writer = bio_io::writer(&opts.out)?; + writeln!( + writer, + "#chrom\tstart\tend\tname\tn_support\tfrac_support\tsupport\tunion_start\tunion_end\tpeak_summit" + )?; + + let n_inputs = samples.len(); + let min_support = opts.min_support.max(1); + let gap = opts.window_size as i64 + ISLAND_PAD; + let mut n_peaks = 0; + for chrom in chrom_lengths.keys() { + let mut all: Vec<(i64, i64)> = samples + .iter() + .filter_map(|s| s.by_chrom.get(chrom)) + .flatten() + .copied() + .collect(); + all.sort_unstable(); + + for (island_start, island_end) in islands(&all, gap) { + if island_end - island_start >= MAX_ISLAND_LEN { + bail!( + "intervals at {chrom}:{island_start}-{island_end} span more than {MAX_ISLAND_LEN} bp, which is too long for a mock fiber" + ); + } + // Every interval falls wholly inside one island, so this is a select, not a clip. + let records = samples + .iter() + .filter_map(|sample| { + let intervals: Vec = sample + .by_chrom + .get(chrom)? + .iter() + .filter(|(start, end)| *start < island_end && *end > island_start) + .map(|&(start, end)| BedRecord { + chrom: chrom.clone(), + start, + end, + name: None, + extra_fields: vec![], + }) + .collect(); + if intervals.is_empty() { + return None; + } + Some(create_mock_fire_record( + &sample.name, + &intervals, + &header_view, + MOCK_FIRE_QUALITY, + None, + )) + }) + .collect::>>()?; + if records.is_empty() { + continue; + } + + let fibers = + FiberseqData::from_records(records, &header_view, &FiberFilters::default()); + let mut peaks = Vec::new(); + call_peaks_for_chrom( + chrom, + island_start as usize, + island_end as usize, + fibers.into_iter(), + ¶ms, + &[], + |peak| { + peaks.push(( + peak.start, + peak.end, + peak.pileup.chrom_start + peak.peak_index, + )); + Ok(()) + }, + )?; + + // Merged peaks come back in whatever order the merge left them. + peaks.sort_unstable(); + for (start, end, summit) in peaks { + let (support, union_start, union_end) = + support_for(&samples, chrom, start as i64, end as i64); + if support.len() < min_support { + continue; + } + n_peaks += 1; + writeln!( + writer, + "{chrom}\t{start}\t{end}\tunion_peak_{n_peaks}\t{}\t{:.4}\t{}\t{union_start}\t{union_end}\t{summit}", + support.len(), + support.len() as f64 / n_inputs as f64, + support.join(","), + )?; + } + } + } + writer.flush()?; + + log::info!( + "{n_peaks} union peaks from {n_inputs} BED files written to {}", + opts.out + ); + Ok(()) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn islands_split_on_gaps_larger_than_the_window() { + let intervals = [(1000, 1200), (1500, 1700), (50000, 50200)]; + assert_eq!( + islands(&intervals, 200 + ISLAND_PAD), + vec![(1000, 1700), (50000, 50200)] + ); + } + + #[test] + fn support_counts_only_overlapping_samples() { + let sample = |name: &str, intervals: Vec<(i64, i64)>| Sample { + name: name.to_string(), + by_chrom: BTreeMap::from([("chr1".to_string(), intervals)]), + }; + let samples = [ + sample("s1", vec![(100, 200)]), + sample("s2", vec![(900, 1000)]), + sample("s3", vec![(150, 260)]), + ]; + assert_eq!( + support_for(&samples, "chr1", 140, 210), + (vec!["s1", "s3"], 100, 260) + ); + assert_eq!(support_for(&samples, "chr2", 140, 210), (vec![], 140, 210)); + } +} diff --git a/tests/regression.rs b/tests/regression.rs index 43eb9467a..be1821476 100644 --- a/tests/regression.rs +++ b/tests/regression.rs @@ -20,3 +20,5 @@ mod pileup; mod predict_m6a; #[path = "regression/qc.rs"] mod qc; +#[path = "regression/union_peaks.rs"] +mod union_peaks; diff --git a/tests/regression/union_peaks.rs b/tests/regression/union_peaks.rs new file mode 100644 index 000000000..cfaebdfea --- /dev/null +++ b/tests/regression/union_peaks.rs @@ -0,0 +1,136 @@ +use super::common::{ft, run}; +use tempfile::{NamedTempFile, TempPath}; + +/// s3 has two overlapping intervals (20000-20200, 20100-20400) so the tests can pin that +/// one file only ever contributes 1 to a peak's support count. +fn fixture() -> [TempPath; 3] { + let beds = [ + "chr1\t1000\t1200\nchr1\t5000\t5300\nchr1\t9000\t9100\nchr2\t100\t400\n", + "chr1\t1020\t1220\nchr1\t5050\t5350\nchr2\t150\t420\n", + "chr1\t1050\t1250\nchr1\t20000\t20200\nchr1\t20100\t20400\nchr2\t120\t380\n", + ]; + beds.map(|contents| { + let bed = NamedTempFile::with_suffix(".bed").unwrap(); + std::fs::write(bed.path(), contents).unwrap(); + bed.into_temp_path() + }) +} + +fn union_peaks(beds: &[TempPath], extra: &[&str]) -> String { + let mut args = vec!["union-peaks"]; + args.extend(beds.iter().map(|b| b.to_str().unwrap())); + args.extend(["--names", "s1,s2,s3"]); + args.extend(extra); + run(&args) +} + +// Peak boundaries are the peak caller's median boundaries, not the outer span of the +// supporting intervals, and support is counted per input file. Peak 4 is the collapse +// case: s3's two overlapping intervals must yield one peak with support 1, not two peaks +// or a support of 2. +#[test] +fn union_peaks_support_counts_and_boundaries() { + let beds = fixture(); + assert_eq!( + union_peaks(&beds, &[]), + "#chrom\tstart\tend\tname\tn_support\tfrac_support\tsupport\tunion_start\tunion_end\tpeak_summit\n\ + chr1\t1020\t1220\tunion_peak_1\t3\t1.0000\ts1,s2,s3\t1000\t1250\t1125\n\ + chr1\t5050\t5350\tunion_peak_2\t2\t0.6667\ts1,s2\t5000\t5350\t5175\n\ + chr1\t9000\t9100\tunion_peak_3\t1\t0.3333\ts1\t9000\t9100\t9050\n\ + chr1\t20000\t20400\tunion_peak_4\t1\t0.3333\ts3\t20000\t20400\t20200\n\ + chr2\t120\t400\tunion_peak_5\t3\t1.0000\ts1,s2,s3\t100\t420\t265\n" + ); +} + +// --min-support is an output filter, so it must drop exactly the low support rows and +// leave the surviving rows byte for byte the same as the unfiltered run. +#[test] +fn union_peaks_min_support_filters_only() { + let beds = fixture(); + let all = union_peaks(&beds, &[]); + let filtered = union_peaks(&beds, &["-n", "2"]); + let kept: Vec<&str> = filtered.lines().skip(1).collect(); + assert_eq!(kept.len(), 3, "got: {filtered}"); + for (line, expected) in kept.iter().zip(["1020\t1220", "5050\t5350", "120\t400"]) { + assert!(line.contains(expected), "{line} lacks {expected}"); + // same peak, only the sequential name changes + let coords = line.split('\t').take(3).collect::>().join("\t"); + assert!(all.contains(&coords), "{coords} not in unfiltered output"); + } +} + +// mock-fire writes in read-name order, so a hand rolled `mock-fire | samtools index | +// call-peaks` pipeline fails with "unsorted positions" whenever sample names sort against +// their positions. union-peaks never writes a BAM, so it has to work here. +#[test] +fn union_peaks_names_ordered_against_positions() { + let beds: [TempPath; 2] = [ + "chr1\t100\t200\nchr10\t120\t220\nchr2\t150\t250\n", + "chr1\t5000\t5100\nchr10\t100\t200\nchr2\t100\t200\n", + ] + .map(|contents| { + let bed = NamedTempFile::with_suffix(".bed").unwrap(); + std::fs::write(bed.path(), contents).unwrap(); + bed.into_temp_path() + }); + let out = run(&[ + "union-peaks", + beds[0].to_str().unwrap(), + beds[1].to_str().unwrap(), + "--names", + "zz,aa", + ]); + let chroms: Vec<&str> = out + .lines() + .skip(1) + .map(|l| l.split('\t').next().unwrap()) + .collect(); + assert_eq!(chroms, ["chr1", "chr1", "chr10", "chr2"], "got: {out}"); + assert_eq!(out.matches("zz,aa").count(), 2, "got: {out}"); +} + +// An empty BED is a peak caller that found nothing, not an error, but it still counts in +// the frac_support denominator. +#[test] +fn union_peaks_empty_input_is_a_sample_with_no_peaks() { + let beds = fixture(); + let empty = NamedTempFile::with_suffix(".bed").unwrap(); + let out = run(&[ + "union-peaks", + beds[0].to_str().unwrap(), + empty.path().to_str().unwrap(), + ]); + assert_eq!(out.lines().count(), 5, "got: {out}"); + assert!( + out.lines().skip(1).all(|l| l.contains("\t0.5000\t")), + "{out}" + ); +} + +#[test] +fn union_peaks_rejects_bad_input() { + let beds = fixture(); + let bad = NamedTempFile::with_suffix(".bed").unwrap(); + std::fs::write(bad.path(), "chr1\t500\t100\n").unwrap(); + let cases: Vec<(Vec<&str>, &str)> = vec![ + ( + vec![beds[0].to_str().unwrap(), "--names", "only-one,and,three"], + "--names has 3 names", + ), + (vec![bad.path().to_str().unwrap()], "invalid interval"), + ( + vec![beds[0].to_str().unwrap(), "--window-size", "1"], + "--window-size must be at least 2", + ), + ]; + for (args, expected) in cases { + let out = std::process::Command::new(ft()) + .arg("union-peaks") + .args(&args) + .output() + .unwrap(); + assert!(!out.status.success(), "{args:?} unexpectedly succeeded"); + let stderr = String::from_utf8_lossy(&out.stderr); + assert!(stderr.contains(expected), "{args:?} stderr: {stderr}"); + } +} From 96c1765d632bc7d1d32946a5f3a5ec055de8424e Mon Sep 17 00:00:00 2001 From: "Mitchell R. Vollger" Date: Fri, 21 Aug 2026 11:12:01 -0600 Subject: [PATCH 02/11] perf: binary search each island's intervals in union-peaks Rescanning every interval of every sample for every island cost islands x intervals per chromosome, which was ~77% of the run on dense input. The intervals are already merged and sorted, so their ends rise with their starts: seek to the first one reaching the island and stop at the first one past it, the way support_for already does. 20 samples x 100k intervals on one chromosome (100k islands): 151.8s -> 19.3s. Whole-genome (20 x 150k peaks, hg38): 269.9s -> 237.9s, RSS flat at ~340 MB. Output byte-identical on all three benchmark sets. Also correct the 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, only local maxima become peaks, and --min-support renumbers the sequential name column. --- src/cli/union_peaks_opts.rs | 3 ++- src/subcommands/union_peaks.rs | 27 ++++++++++++++++++--------- tests/regression/union_peaks.rs | 2 +- 3 files changed, 21 insertions(+), 11 deletions(-) diff --git a/src/cli/union_peaks_opts.rs b/src/cli/union_peaks_opts.rs index be7a040e5..8c771d882 100644 --- a/src/cli/union_peaks_opts.rs +++ b/src/cli/union_peaks_opts.rs @@ -19,7 +19,8 @@ pub struct UnionPeaksOptions { /// Minimum number of input BEDs that must overlap a peak for it to be reported #[clap(short = 'n', long, default_value_t = 1)] pub min_support: usize, - /// Rolling window size for finding local maxima (in base pairs) + /// Rolling window size for finding local maxima (in base pairs). + /// Only local maxima are kept, so at most one peak is reported per window. #[clap(long, default_value_t = 200)] pub window_size: usize, #[clap(flatten)] diff --git a/src/subcommands/union_peaks.rs b/src/subcommands/union_peaks.rs index 79f46c4d9..b5177e655 100644 --- a/src/subcommands/union_peaks.rs +++ b/src/subcommands/union_peaks.rs @@ -14,8 +14,10 @@ use std::io::Write; const MOCK_FIRE_QUALITY: u8 = 255; /// Slack added to `--window-size` when grouping intervals into islands. The gap has to -/// exceed the rolling-max window, otherwise a window could span two islands and the -/// islanded result would differ from a whole-chromosome one. +/// exceed the rolling-max window so that no window ever spans two islands. That bounds the +/// rolling-max failure mode; it does not make islanding identical to a whole-chromosome +/// run, since confining a sample's mock fiber to one island also drops the coverage it +/// would have contributed between islands, which can move a summit. const ISLAND_PAD: i64 = 1000; /// Longest island we will build a mock fiber for. `Cigar::Equal(len)` packs the length @@ -185,8 +187,11 @@ fn support_for<'a>( /// /// Reported `start`/`end` are therefore the peak caller's consensus (median) boundaries /// of the overlapping input intervals, not their outer span; the outer span is reported -/// alongside as `union_start`/`union_end`. `--min-support` filters the output only, so -/// `-n 3` and `-n 1` plus a downstream filter agree. +/// alongside as `union_start`/`union_end`, and `peak_summit` (the pileup's local maximum, +/// `peak_max` in `ft call-peaks` output) can sit outside `start`/`end` for the same reason. +/// Only local maxima become peaks, so at most one peak is reported per `--window-size` +/// bases. `--min-support` filters the output only, so `-n 3` and `-n 1` plus a downstream +/// filter agree apart from the sequential `name` column, which renumbers. pub fn run_union_peaks(opts: &UnionPeaksOptions) -> Result<()> { if opts.window_size < 2 { bail!("--window-size must be at least 2; a smaller window finds no local maxima"); @@ -269,15 +274,19 @@ pub fn run_union_peaks(opts: &UnionPeaksOptions) -> Result<()> { "intervals at {chrom}:{island_start}-{island_end} span more than {MAX_ISLAND_LEN} bp, which is too long for a mock fiber" ); } - // Every interval falls wholly inside one island, so this is a select, not a clip. + // Every interval falls wholly inside one island, so this is a select, not a + // clip. Intervals are merged and sorted, so their ends rise with their starts: + // binary search to the first one reaching this island and stop at the first one + // past it, the way support_for does. Scanning them all instead costs + // islands x intervals, which dominates the run on dense whole-genome input. let records = samples .iter() .filter_map(|sample| { - let intervals: Vec = sample - .by_chrom - .get(chrom)? + let sample_intervals = sample.by_chrom.get(chrom)?; + let first = sample_intervals.partition_point(|iv| iv.1 <= island_start); + let intervals: Vec = sample_intervals[first..] .iter() - .filter(|(start, end)| *start < island_end && *end > island_start) + .take_while(|(start, _)| *start < island_end) .map(|&(start, end)| BedRecord { chrom: chrom.clone(), start, diff --git a/tests/regression/union_peaks.rs b/tests/regression/union_peaks.rs index cfaebdfea..7b7ef8dec 100644 --- a/tests/regression/union_peaks.rs +++ b/tests/regression/union_peaks.rs @@ -43,7 +43,7 @@ fn union_peaks_support_counts_and_boundaries() { } // --min-support is an output filter, so it must drop exactly the low support rows and -// leave the surviving rows byte for byte the same as the unfiltered run. +// leave the surviving rows unchanged apart from the sequential name, which renumbers. #[test] fn union_peaks_min_support_filters_only() { let beds = fixture(); From ab58d997ed561a452b39675ac71b33ac456d6f57 Mon Sep 17 00:00:00 2001 From: "Mitchell R. Vollger" Date: Fri, 21 Aug 2026 12:28:25 -0600 Subject: [PATCH 03/11] docs: quantify what islanding costs in union-peaks Measured islanded output against a whole-chromosome mock-fire | sort | index | call-peaks run over 60 adversarial configurations: 11 summit shifts, 1 peak-set change. The shifts are tie-breaks inside a plateau of equal sample support, and the support columns never come from the pileup. --- src/subcommands/union_peaks.rs | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/subcommands/union_peaks.rs b/src/subcommands/union_peaks.rs index b5177e655..ecb779863 100644 --- a/src/subcommands/union_peaks.rs +++ b/src/subcommands/union_peaks.rs @@ -17,7 +17,11 @@ const MOCK_FIRE_QUALITY: u8 = 255; /// exceed the rolling-max window so that no window ever spans two islands. That bounds the /// rolling-max failure mode; it does not make islanding identical to a whole-chromosome /// run, since confining a sample's mock fiber to one island also drops the coverage it -/// would have contributed between islands, which can move a summit. +/// would have contributed between islands. Coverage is the denominator of the FIRE score, +/// so that can move a summit: measured at 11 summit shifts and 1 peak-set change in 60 +/// adversarial random configurations. The shifts are tie-breaks inside a plateau of equal +/// sample support (the summit lands on a maximum-support position either way, 99% of the +/// time), and n_support/support/union_start/union_end never come from the pileup at all. const ISLAND_PAD: i64 = 1000; /// Longest island we will build a mock fiber for. `Cigar::Equal(len)` packs the length From 07edebec2221d600e942eccfdac6d16d164f390c Mon Sep 17 00:00:00 2001 From: "Mitchell R. Vollger" Date: Fri, 21 Aug 2026 19:23:38 -0600 Subject: [PATCH 04/11] test: pin ft call-peaks output on the ctcf fixture call_peaks_for_chrom is now also union-peaks' engine; the snapshot makes the byte-identical refactor gate permanent. --- tests/regression.rs | 2 ++ tests/regression/call_peaks.rs | 14 ++++++++ ..._call_peaks__call_peaks_ctcf_snapshot.snap | 36 +++++++++++++++++++ 3 files changed, 52 insertions(+) create mode 100644 tests/regression/call_peaks.rs create mode 100644 tests/regression/snapshots/regression__call_peaks__call_peaks_ctcf_snapshot.snap diff --git a/tests/regression.rs b/tests/regression.rs index be1821476..eae6c2f36 100644 --- a/tests/regression.rs +++ b/tests/regression.rs @@ -1,3 +1,5 @@ +#[path = "regression/call_peaks.rs"] +mod call_peaks; #[path = "regression/center.rs"] mod center; #[path = "regression/common.rs"] diff --git a/tests/regression/call_peaks.rs b/tests/regression/call_peaks.rs new file mode 100644 index 000000000..39fa148ce --- /dev/null +++ b/tests/regression/call_peaks.rs @@ -0,0 +1,14 @@ +use super::common::{fixture, run}; + +/// Pin ft call-peaks output so refactors of the shared peak caller +/// (call_peaks_for_chrom is also union-peaks' engine) cannot drift it. +#[test] +fn call_peaks_ctcf_snapshot() { + let out = run(&[ + "call-peaks", + fixture("ctcf.bam").to_str().unwrap(), + "--min-fire-frac", + "0.5", + ]); + insta::assert_snapshot!(out); +} diff --git a/tests/regression/snapshots/regression__call_peaks__call_peaks_ctcf_snapshot.snap b/tests/regression/snapshots/regression__call_peaks__call_peaks_ctcf_snapshot.snap new file mode 100644 index 000000000..a8789765e --- /dev/null +++ b/tests/regression/snapshots/regression__call_peaks__call_peaks_ctcf_snapshot.snap @@ -0,0 +1,36 @@ +--- +source: tests/regression/call_peaks.rs +expression: out +--- +#chrom peak_start peak_end peak_max FDR coverage fire_coverage score nuc_coverage msp_coverage coverage_H1 fire_coverage_H1 score_H1 nuc_coverage_H1 msp_coverage_H1 coverage_H2 fire_coverage_H2 score_H2 nuc_coverage_H2 msp_coverage_H2 pass_coverage +chr4 3009358 3009585 3009467 1.0000000000 170 100 44.39899 30 140 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr4 3013814 3014143 3013962 1.0000000000 152 97 49.55543 35 117 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr4 3037395 3037600 3037451 1.0000000000 123 88 55.25256 11 112 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr4 3040474 3040684 3040577 1.0000000000 114 59 38.09894 16 98 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr4 3047477 3048067 3047700 1.0000000000 121 117 80.41655 4 117 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr4 3073949 3074404 3074219 1.0000000000 104 96 72.09167 1 103 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr4 3074493 3074780 3074587 1.0000000000 105 102 78.14516 0 105 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr4 3075480 3075719 3075628 1.0000000000 100 55 40.77604 18 82 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr4 3077847 3078142 3078032 1.0000000000 57 53 69.93767 3 54 0 0 -1.0 0 0 0 0 -1.0 0 0 true +chr4 3079695 3079972 3079819 1.0000000000 44 28 48.85315 5 39 0 0 -1.0 0 0 0 0 -1.0 0 0 true +chr4 3085370 3085641 3085552 1.0000000000 8 5 44.75336 0 8 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr4 3106679 3106892 3106783 1.0000000000 125 112 69.67896 6 119 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr11 5280593 5280862 5280676 1.0000000000 14 8 46.76774 4 10 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr11 5291318 5291539 5291448 1.0000000000 165 97 46.41454 31 134 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr19 45669116 45669382 45669237 1.0000000000 14 10 58.72945 1 13 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr19 45670933 45671222 45671096 1.0000000000 17 9 40.57255 4 13 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr19 45677410 45677730 45677567 1.0000000000 52 36 52.86568 4 48 0 0 -1.0 0 0 0 0 -1.0 0 0 true +chr19 45681546 45681799 45681658 1.0000000000 101 98 77.05413 2 99 0 0 -1.0 0 0 0 0 -1.0 0 0 true +chr19 45691846 45692083 45691970 1.0000000000 93 87 70.31076 3 90 0 0 -1.0 0 0 0 0 -1.0 0 0 true +chr19 45692182 45692751 45692606 1.0000000000 94 92 76.69306 0 94 0 0 -1.0 0 0 0 0 -1.0 0 0 true +chr19 45693394 45693604 45693510 1.0000000000 102 68 49.32525 4 98 0 0 -1.0 0 0 0 0 -1.0 0 0 true +chr19 45706447 45706762 45706633 1.0000000000 16 12 52.59830 3 13 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr19 45707151 45707519 45707383 1.0000000000 17 12 51.08542 2 15 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr19 45717239 45717778 45717439 1.0000000000 81 77 73.71977 1 80 0 0 -1.0 0 0 0 0 -1.0 0 0 true +chr19 45730798 45731169 45731005 1.0000000000 90 87 75.56393 2 88 0 0 -1.0 0 0 0 0 -1.0 0 0 true +chr19 45769332 45769709 45769497 1.0000000000 60 56 68.42448 0 60 0 0 -1.0 0 0 0 0 -1.0 0 0 true +chr19 45770241 45770482 45770388 1.0000000000 71 47 48.65365 5 66 0 0 -1.0 0 0 0 0 -1.0 0 0 true +chr19 45770577 45770938 45770673 1.0000000000 69 40 42.31902 25 44 0 0 -1.0 0 0 0 0 -1.0 0 0 true +chr19 45779528 45779801 45779699 1.0000000000 104 60 41.48181 18 86 0 0 -1.0 0 0 0 0 -1.0 0 0 true +chr19 45792728 45793091 45792945 1.0000000000 41 38 69.14936 0 41 0 0 -1.0 0 0 0 0 -1.0 0 0 true +chr19 45799119 45799332 45799284 1.0000000000 93 49 40.12040 15 78 0 0 -1.0 0 0 0 0 -1.0 0 0 true From 638e078c982c79039d450230d5c2a04833a86d0c Mon Sep 17 00:00:00 2001 From: "Mitchell R. Vollger" Date: Fri, 21 Aug 2026 19:23:38 -0600 Subject: [PATCH 05/11] chore: hide ft mock-fire and warn of its removal ft union-peaks replaces the mock-fire | sort | call-peaks pipeline. The command keeps working for existing scripts (StergachisLab fire_consensus_pipeline) but is hidden from help and warns at runtime; also updates a stale pre-callable-model comment. --- src/cli.rs | 2 +- src/subcommands/mock_fire.rs | 12 ++++++++---- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/src/cli.rs b/src/cli.rs index 0c2d0975e..aea581e9a 100644 --- a/src/cli.rs +++ b/src/cli.rs @@ -165,7 +165,7 @@ pub enum Commands { CallPeaks(CallPeaksOptions), /// Create a mock BAM file with FIRE elements from a BED file. /// Each interval in the BED becomes a FIRE element. The 4th column groups intervals into the same mock read. - #[clap(name = "mock-fire")] + #[clap(name = "mock-fire", hide = true)] MockFire(MockFireOptions), /// Combine peak calls from many BED files into one union peak set. /// Each input BED is one sample, and peaks are reported with the number and names of diff --git a/src/subcommands/mock_fire.rs b/src/subcommands/mock_fire.rs index 4291dc52e..152f877ab 100644 --- a/src/subcommands/mock_fire.rs +++ b/src/subcommands/mock_fire.rs @@ -141,10 +141,10 @@ pub(crate) fn create_mock_fire_record( let mut annot = MolecularAnnotations::from_record(&record); ma_io::add_msp_annotations(&mut annot, &starts, &lengths, None); ma_io::add_fire_annotations(&mut annot, &starts, &lengths, &quals); - // Mock records carry no m6A, so the read-time backfill would brand - // them NotCallable and --callable-fibers would drop every one. - // Write an explicit full-width Callable annotation instead: mock data is - // callable by construction. + // Sparse mock reads would fail the >=10-MSP minimum if the tag were + // derived, and --callable-fibers would drop them. Write an explicit + // full-width Callable annotation instead: mock data is callable by + // construction. { use crate::utils::input_bam::{FIRE_CALLABLE_MIN_AVE_MSP_SIZE, FIRE_CALLABLE_MIN_MSP}; ma_io::set_fiberseq_callable( @@ -159,6 +159,10 @@ pub(crate) fn create_mock_fire_record( } pub fn run_mock_fire(opts: &MockFireOptions) -> Result<()> { + log::warn!( + "ft mock-fire is deprecated and will be removed in a future release; \ + ft union-peaks replaces the mock-fire | sort | call-peaks pipeline" + ); log::info!("Reading BED file: {}", opts.bed); // Read BED file From 48e89ffc833eb549f477c1623780bb260af6cdc3 Mon Sep 17 00:00:00 2001 From: "Mitchell R. Vollger" Date: Fri, 21 Aug 2026 20:30:28 -0600 Subject: [PATCH 06/11] fix: harden union-peaks from the deep review - Book-ended intervals from different samples: the touching sample formed the peak but was silently dropped from support and the union span; support now counts exact abutment with the consensus core (whose bounds are medians of element edges). The consensus core itself is unchanged, so ft call-peaks output is untouched. - Coordinates past the mock-header/liftover range (int32 LN, u32 positions) error cleanly instead of vanishing with exit 0 or panicking. - Warn once when an island exceeds 10 Mb: memory scales with island size and dense input can grow islands to chromosome scale. - --min-support 0 errors instead of silently acting as 1; -n above the input count warns that no peak can be reported. - read_bed_regions skips track/browser lines, so MACS2/UCSC peak files work unedited. - Docs: within-file collapse also merges book-ended intervals; the 56 B/base figure gains its per-element term. Both fixes pinned by new regression tests. --- src/cli/union_peaks_opts.rs | 2 +- src/subcommands/union_peaks.rs | 42 ++++++++++++++++++++++++++------- src/utils/bio_io.rs | 6 ++++- tests/regression/union_peaks.rs | 41 ++++++++++++++++++++++++++++++++ 4 files changed, 81 insertions(+), 10 deletions(-) diff --git a/src/cli/union_peaks_opts.rs b/src/cli/union_peaks_opts.rs index 8c771d882..439bbe50c 100644 --- a/src/cli/union_peaks_opts.rs +++ b/src/cli/union_peaks_opts.rs @@ -6,7 +6,7 @@ use std::fmt::Debug; pub struct UnionPeaksOptions { /// Input BED files, one per sample. /// Every interval becomes a FIRE element on a mock fiber for that sample, and - /// peaks are called across all the samples at once. Overlapping intervals within + /// peaks are called across all the samples at once. Overlapping or book-ended intervals within /// one file are merged first, so a file can add at most 1 to a peak's support. #[clap(required = true, num_args = 1..)] pub beds: Vec, diff --git a/src/subcommands/union_peaks.rs b/src/subcommands/union_peaks.rs index ecb779863..2c3926fc6 100644 --- a/src/subcommands/union_peaks.rs +++ b/src/subcommands/union_peaks.rs @@ -23,6 +23,9 @@ const MOCK_FIRE_QUALITY: u8 = 255; /// sample support (the summit lands on a maximum-support position either way, 99% of the /// time), and n_support/support/union_start/union_end never come from the pileup at all. const ISLAND_PAD: i64 = 1000; +/// BAM reference lengths are int32 and the liftover stores positions as u32; cap +/// coordinates so the mock header (max end + 10000) and every position stay in range. +const MAX_COORD: i64 = i32::MAX as i64 - 10_000; /// Longest island we will build a mock fiber for. `Cigar::Equal(len)` packs the length /// into 28 bits and wraps silently past that. @@ -100,10 +103,11 @@ fn load_sample(path: &str, name: String) -> Result { for rec in read_bed_regions(path).with_context(|| format!("failed to read BED file {path}"))? { // read_bed_regions parses the coordinates as bare i64s, and a negative start makes // create_mock_fire_record emit a record whose positions are all dropped later. - if rec.start < 0 || rec.end <= rec.start { + if rec.start < 0 || rec.end <= rec.start || rec.end > MAX_COORD { bail!( - "invalid interval in {}: {} {} {}", + "invalid interval in {} (coordinates must be 0 <= start < end <= {}): {} {} {}", path, + MAX_COORD, rec.chrom, rec.start, rec.end @@ -135,7 +139,8 @@ fn load_sample(path: &str, name: String) -> Result { /// Group sorted intervals into maximal runs separated by no more than `gap` bases. /// /// Peaks are called one island at a time so the pileup track is sized to the data. A -/// whole-chromosome track costs ~56 bytes a base, which is tens of GB for a genome-wide +/// whole-chromosome track costs ~56 bytes a base (plus 24 bytes per overlapping sample +/// element on covered bases), which is tens of GB for a genome-wide /// peak union. fn islands(sorted: &[(i64, i64)], gap: i64) -> Vec<(i64, i64)> { let mut out: Vec<(i64, i64)> = Vec::new(); @@ -148,8 +153,10 @@ fn islands(sorted: &[(i64, i64)], gap: i64) -> Vec<(i64, i64)> { out } -/// Samples with an interval overlapping `[start, end)`, in input order, plus the outer -/// span of those intervals. +/// Samples with an interval overlapping or exactly abutting `[start, end)`, in input +/// order, plus the outer span of those intervals. Touch counts because the consensus +/// bounds are medians of element edges: an exactly book-ended element that formed the +/// peak can land flush against the reported boundary. /// /// The pileup cannot answer this: its FIRE elements carry no fiber identity. fn support_for<'a>( @@ -165,9 +172,9 @@ fn support_for<'a>( continue; }; let mut hit = false; - let first = intervals.partition_point(|iv| iv.1 <= start); + let first = intervals.partition_point(|iv| iv.1 < start); for &(iv_start, iv_end) in &intervals[first..] { - if iv_start >= end { + if iv_start > end { break; } hit = true; @@ -260,7 +267,17 @@ pub fn run_union_peaks(opts: &UnionPeaksOptions) -> Result<()> { )?; let n_inputs = samples.len(); - let min_support = opts.min_support.max(1); + if opts.min_support == 0 { + bail!("--min-support must be >= 1"); + } + if opts.min_support > n_inputs { + log::warn!( + "--min-support {} exceeds the {} input BEDs; no peak can be reported", + opts.min_support, + n_inputs + ); + } + let min_support = opts.min_support; let gap = opts.window_size as i64 + ISLAND_PAD; let mut n_peaks = 0; for chrom in chrom_lengths.keys() { @@ -273,6 +290,15 @@ pub fn run_union_peaks(opts: &UnionPeaksOptions) -> Result<()> { all.sort_unstable(); for (island_start, island_end) in islands(&all, gap) { + if island_end - island_start > 10_000_000 { + static BIG_ISLAND: std::sync::Once = std::sync::Once::new(); + BIG_ISLAND.call_once(|| { + log::warn!( + "an island spans {} Mb ({chrom}:{island_start}-{island_end}); memory scales with island size, and dense input with no gaps over {gap} bp can grow islands to chromosome scale", + (island_end - island_start) / 1_000_000 + ); + }); + } if island_end - island_start >= MAX_ISLAND_LEN { bail!( "intervals at {chrom}:{island_start}-{island_end} span more than {MAX_ISLAND_LEN} bp, which is too long for a mock fiber" diff --git a/src/utils/bio_io.rs b/src/utils/bio_io.rs index 427690b3b..875a36da2 100644 --- a/src/utils/bio_io.rs +++ b/src/utils/bio_io.rs @@ -596,7 +596,11 @@ pub fn read_bed_regions(bed_path: &str) -> Result> { for line in reader.lines() { let line = line?; - if line.starts_with('#') || line.trim().is_empty() { + if line.starts_with('#') + || line.starts_with("track") + || line.starts_with("browser") + || line.trim().is_empty() + { continue; } let tokens: Vec<&str> = line.split('\t').collect(); diff --git a/tests/regression/union_peaks.rs b/tests/regression/union_peaks.rs index 7b7ef8dec..cab89d16c 100644 --- a/tests/regression/union_peaks.rs +++ b/tests/regression/union_peaks.rs @@ -134,3 +134,44 @@ fn union_peaks_rejects_bad_input() { assert!(stderr.contains(expected), "{args:?} stderr: {stderr}"); } } + +// Exactly book-ended intervals from two samples form one peak whose consensus core is +// the upper-median element; the touching sample must still be counted (support counts +// exact abutment because consensus bounds are medians of element edges). Before this +// was fixed, bkA vanished from the output entirely. +#[test] +fn book_ended_samples_both_support_the_peak() { + let a = NamedTempFile::with_suffix(".bed").unwrap(); + std::fs::write(a.path(), "chr1\t100\t200\n").unwrap(); + let b = NamedTempFile::with_suffix(".bed").unwrap(); + std::fs::write(b.path(), "chr1\t200\t300\n").unwrap(); + let out = run(&[ + "union-peaks", + a.path().to_str().unwrap(), + b.path().to_str().unwrap(), + "--names", + "bkA,bkB", + ]); + let rows: Vec<&str> = out.lines().skip(1).collect(); + assert_eq!(rows.len(), 1, "one peak: {out}"); + let f: Vec<&str> = rows[0].split('\t').collect(); + assert_eq!( + (f[4], f[6], f[7], f[8]), + ("2", "bkA,bkB", "100", "300"), + "both samples support the peak and the union spans both: {out}" + ); +} + +// Coordinates past the mock-header/liftover range must fail loudly, not vanish. +#[test] +fn oversized_coordinates_error_cleanly() { + let bed = NamedTempFile::with_suffix(".bed").unwrap(); + std::fs::write(bed.path(), "chr1\t4294967296\t4294967496\n").unwrap(); + let out = std::process::Command::new(ft()) + .args(["union-peaks", bed.path().to_str().unwrap()]) + .output() + .unwrap(); + assert!(!out.status.success(), "must fail"); + let err = String::from_utf8_lossy(&out.stderr); + assert!(err.contains("invalid interval"), "clean error: {err}"); +} From 469ede860dce3468f7633126c56d49cf33e30b7b Mon Sep 17 00:00:00 2001 From: "Mitchell R. Vollger" Date: Fri, 21 Aug 2026 20:44:16 -0600 Subject: [PATCH 07/11] fix!: true median for peak consensus boundaries The consensus start/end were the upper median of the pooled FIRE element edges, which lands entirely on one element's bounds whenever contributors tie (the book-ended two-sample case sat flush on the second sample). Even counts now take the midpoint of the two middle values. On the ctcf fixture 6 of 31 call-peaks boundaries shift by 1-15 bp (summits and all statistics unchanged); the union-peaks two-sample core centers between its contributors. --- src/subcommands/call_peaks/peaks.rs | 19 ++++++++++++++---- ..._call_peaks__call_peaks_ctcf_snapshot.snap | 20 +++++++++---------- tests/regression/union_peaks.rs | 17 ++++++++++++++-- 3 files changed, 40 insertions(+), 16 deletions(-) diff --git a/src/subcommands/call_peaks/peaks.rs b/src/subcommands/call_peaks/peaks.rs index 7d7281bbe..db887668d 100644 --- a/src/subcommands/call_peaks/peaks.rs +++ b/src/subcommands/call_peaks/peaks.rs @@ -35,6 +35,17 @@ pub fn reciprocal_overlap_raw( (overlap_len / a_len).min(overlap_len / b_len) } +/// Median of a sorted slice: middle value for odd counts, midpoint of the two middle +/// values for even counts. +fn median(sorted: &[i64]) -> i64 { + let n = sorted.len(); + if n % 2 == 1 { + sorted[n / 2] + } else { + (sorted[n / 2 - 1] + sorted[n / 2]) / 2 + } +} + /// Everything peak calling needs from the CLI, decoupled from `CallPeaksOptions` /// (which flattens `InputBam`/`FiberFilters`, so commands that never read a BAM /// cannot build one). See `ft union-peaks` for the other caller. @@ -312,12 +323,12 @@ impl<'a> Peak<'a> { let end = pileup.chrom_start + positions[positions.len() - 1] + 1; (start, end) } else { - // Calculate median start and end from FIRE elements + // True median: for even counts, the midpoint of the two middle + // values. The upper median biased the consensus onto one + // element's bounds whenever contributors tied. starts.sort_unstable(); ends.sort_unstable(); - let median_start = starts[starts.len() / 2] as usize; - let median_end = ends[ends.len() / 2] as usize; - (median_start, median_end) + (median(&starts) as usize, median(&ends) as usize) } } else { // Fallback: FIRE element tracking not enabled diff --git a/tests/regression/snapshots/regression__call_peaks__call_peaks_ctcf_snapshot.snap b/tests/regression/snapshots/regression__call_peaks__call_peaks_ctcf_snapshot.snap index a8789765e..ce0a653b9 100644 --- a/tests/regression/snapshots/regression__call_peaks__call_peaks_ctcf_snapshot.snap +++ b/tests/regression/snapshots/regression__call_peaks__call_peaks_ctcf_snapshot.snap @@ -8,29 +8,29 @@ chr4 3013814 3014143 3013962 1.0000000000 152 97 49.55543 35 117 0 0 -1.0 0 0 0 chr4 3037395 3037600 3037451 1.0000000000 123 88 55.25256 11 112 0 0 -1.0 0 0 0 0 -1.0 0 0 false chr4 3040474 3040684 3040577 1.0000000000 114 59 38.09894 16 98 0 0 -1.0 0 0 0 0 -1.0 0 0 false chr4 3047477 3048067 3047700 1.0000000000 121 117 80.41655 4 117 0 0 -1.0 0 0 0 0 -1.0 0 0 false -chr4 3073949 3074404 3074219 1.0000000000 104 96 72.09167 1 103 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr4 3073949 3074403 3074219 1.0000000000 104 96 72.09167 1 103 0 0 -1.0 0 0 0 0 -1.0 0 0 false chr4 3074493 3074780 3074587 1.0000000000 105 102 78.14516 0 105 0 0 -1.0 0 0 0 0 -1.0 0 0 false chr4 3075480 3075719 3075628 1.0000000000 100 55 40.77604 18 82 0 0 -1.0 0 0 0 0 -1.0 0 0 false chr4 3077847 3078142 3078032 1.0000000000 57 53 69.93767 3 54 0 0 -1.0 0 0 0 0 -1.0 0 0 true -chr4 3079695 3079972 3079819 1.0000000000 44 28 48.85315 5 39 0 0 -1.0 0 0 0 0 -1.0 0 0 true +chr4 3079689 3079970 3079819 1.0000000000 44 28 48.85315 5 39 0 0 -1.0 0 0 0 0 -1.0 0 0 true chr4 3085370 3085641 3085552 1.0000000000 8 5 44.75336 0 8 0 0 -1.0 0 0 0 0 -1.0 0 0 false -chr4 3106679 3106892 3106783 1.0000000000 125 112 69.67896 6 119 0 0 -1.0 0 0 0 0 -1.0 0 0 false -chr11 5280593 5280862 5280676 1.0000000000 14 8 46.76774 4 10 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr4 3106678 3106892 3106783 1.0000000000 125 112 69.67896 6 119 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr11 5280590 5280847 5280676 1.0000000000 14 8 46.76774 4 10 0 0 -1.0 0 0 0 0 -1.0 0 0 false chr11 5291318 5291539 5291448 1.0000000000 165 97 46.41454 31 134 0 0 -1.0 0 0 0 0 -1.0 0 0 false -chr19 45669116 45669382 45669237 1.0000000000 14 10 58.72945 1 13 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr19 45669109 45669371 45669237 1.0000000000 14 10 58.72945 1 13 0 0 -1.0 0 0 0 0 -1.0 0 0 false chr19 45670933 45671222 45671096 1.0000000000 17 9 40.57255 4 13 0 0 -1.0 0 0 0 0 -1.0 0 0 false -chr19 45677410 45677730 45677567 1.0000000000 52 36 52.86568 4 48 0 0 -1.0 0 0 0 0 -1.0 0 0 true +chr19 45677398 45677730 45677567 1.0000000000 52 36 52.86568 4 48 0 0 -1.0 0 0 0 0 -1.0 0 0 true chr19 45681546 45681799 45681658 1.0000000000 101 98 77.05413 2 99 0 0 -1.0 0 0 0 0 -1.0 0 0 true chr19 45691846 45692083 45691970 1.0000000000 93 87 70.31076 3 90 0 0 -1.0 0 0 0 0 -1.0 0 0 true chr19 45692182 45692751 45692606 1.0000000000 94 92 76.69306 0 94 0 0 -1.0 0 0 0 0 -1.0 0 0 true chr19 45693394 45693604 45693510 1.0000000000 102 68 49.32525 4 98 0 0 -1.0 0 0 0 0 -1.0 0 0 true -chr19 45706447 45706762 45706633 1.0000000000 16 12 52.59830 3 13 0 0 -1.0 0 0 0 0 -1.0 0 0 false -chr19 45707151 45707519 45707383 1.0000000000 17 12 51.08542 2 15 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr19 45706445 45706757 45706633 1.0000000000 16 12 52.59830 3 13 0 0 -1.0 0 0 0 0 -1.0 0 0 false +chr19 45707150 45707486 45707383 1.0000000000 17 12 51.08542 2 15 0 0 -1.0 0 0 0 0 -1.0 0 0 false chr19 45717239 45717778 45717439 1.0000000000 81 77 73.71977 1 80 0 0 -1.0 0 0 0 0 -1.0 0 0 true chr19 45730798 45731169 45731005 1.0000000000 90 87 75.56393 2 88 0 0 -1.0 0 0 0 0 -1.0 0 0 true chr19 45769332 45769709 45769497 1.0000000000 60 56 68.42448 0 60 0 0 -1.0 0 0 0 0 -1.0 0 0 true chr19 45770241 45770482 45770388 1.0000000000 71 47 48.65365 5 66 0 0 -1.0 0 0 0 0 -1.0 0 0 true -chr19 45770577 45770938 45770673 1.0000000000 69 40 42.31902 25 44 0 0 -1.0 0 0 0 0 -1.0 0 0 true +chr19 45770571 45770931 45770673 1.0000000000 69 40 42.31902 25 44 0 0 -1.0 0 0 0 0 -1.0 0 0 true chr19 45779528 45779801 45779699 1.0000000000 104 60 41.48181 18 86 0 0 -1.0 0 0 0 0 -1.0 0 0 true -chr19 45792728 45793091 45792945 1.0000000000 41 38 69.14936 0 41 0 0 -1.0 0 0 0 0 -1.0 0 0 true +chr19 45792728 45793089 45792945 1.0000000000 41 38 69.14936 0 41 0 0 -1.0 0 0 0 0 -1.0 0 0 true chr19 45799119 45799332 45799284 1.0000000000 93 49 40.12040 15 78 0 0 -1.0 0 0 0 0 -1.0 0 0 true diff --git a/tests/regression/union_peaks.rs b/tests/regression/union_peaks.rs index cab89d16c..efcdb674f 100644 --- a/tests/regression/union_peaks.rs +++ b/tests/regression/union_peaks.rs @@ -35,7 +35,7 @@ fn union_peaks_support_counts_and_boundaries() { union_peaks(&beds, &[]), "#chrom\tstart\tend\tname\tn_support\tfrac_support\tsupport\tunion_start\tunion_end\tpeak_summit\n\ chr1\t1020\t1220\tunion_peak_1\t3\t1.0000\ts1,s2,s3\t1000\t1250\t1125\n\ - chr1\t5050\t5350\tunion_peak_2\t2\t0.6667\ts1,s2\t5000\t5350\t5175\n\ + chr1\t5025\t5325\tunion_peak_2\t2\t0.6667\ts1,s2\t5000\t5350\t5175\n\ chr1\t9000\t9100\tunion_peak_3\t1\t0.3333\ts1\t9000\t9100\t9050\n\ chr1\t20000\t20400\tunion_peak_4\t1\t0.3333\ts3\t20000\t20400\t20200\n\ chr2\t120\t400\tunion_peak_5\t3\t1.0000\ts1,s2,s3\t100\t420\t265\n" @@ -51,7 +51,7 @@ fn union_peaks_min_support_filters_only() { let filtered = union_peaks(&beds, &["-n", "2"]); let kept: Vec<&str> = filtered.lines().skip(1).collect(); assert_eq!(kept.len(), 3, "got: {filtered}"); - for (line, expected) in kept.iter().zip(["1020\t1220", "5050\t5350", "120\t400"]) { + for (line, expected) in kept.iter().zip(["1020\t1220", "5025\t5325", "120\t400"]) { assert!(line.contains(expected), "{line} lacks {expected}"); // same peak, only the sequential name changes let coords = line.split('\t').take(3).collect::>().join("\t"); @@ -175,3 +175,16 @@ fn oversized_coordinates_error_cleanly() { let err = String::from_utf8_lossy(&out.stderr); assert!(err.contains("invalid interval"), "clean error: {err}"); } + +// --min-frac-support is the fractional twin of -n: 0.5 of 3 inputs rounds up to 2. +#[test] +fn min_frac_support_filters_like_min_support() { + let beds = fixture(); + let out = union_peaks(&beds, &["--min-frac-support", "0.5"]); + let rows: Vec<&str> = out.lines().skip(1).collect(); + assert!(rows.iter().all(|r| { + let n: usize = r.split('\t').nth(4).unwrap().parse().unwrap(); + n >= 2 + })); + assert_eq!(rows.len(), 3, "peaks with support 1 are dropped: {out}"); +} From bb56c70e1b2be6127a9cf89f0001f30b8baae0a7 Mon Sep 17 00:00:00 2001 From: "Mitchell R. Vollger" Date: Fri, 21 Aug 2026 20:44:16 -0600 Subject: [PATCH 08/11] feat: --min-frac-support, the fractional twin of -n ceil(frac x n_inputs) folds into the same threshold, making the fire_consensus_pipeline's --min-fire-frac semantics expressible directly. --- src/cli/union_peaks_opts.rs | 13 +++++++++++++ src/subcommands/union_peaks.rs | 4 +++- 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/src/cli/union_peaks_opts.rs b/src/cli/union_peaks_opts.rs index 439bbe50c..470bf8074 100644 --- a/src/cli/union_peaks_opts.rs +++ b/src/cli/union_peaks_opts.rs @@ -19,6 +19,10 @@ pub struct UnionPeaksOptions { /// Minimum number of input BEDs that must overlap a peak for it to be reported #[clap(short = 'n', long, default_value_t = 1)] pub min_support: usize, + /// Minimum fraction of input BEDs that must overlap a peak for it to be + /// reported (0-1). Applied together with --min-support. + #[clap(long, value_parser = frac_in_range)] + pub min_frac_support: Option, /// Rolling window size for finding local maxima (in base pairs). /// Only local maxima are kept, so at most one peak is reported per window. #[clap(long, default_value_t = 200)] @@ -26,3 +30,12 @@ pub struct UnionPeaksOptions { #[clap(flatten)] pub global: GlobalOpts, } + +fn frac_in_range(s: &str) -> Result { + let v: f64 = s.parse().map_err(|e| format!("{e}"))?; + if (0.0..=1.0).contains(&v) { + Ok(v) + } else { + Err("must be between 0 and 1".to_string()) + } +} diff --git a/src/subcommands/union_peaks.rs b/src/subcommands/union_peaks.rs index 2c3926fc6..c8c65e09d 100644 --- a/src/subcommands/union_peaks.rs +++ b/src/subcommands/union_peaks.rs @@ -277,7 +277,9 @@ pub fn run_union_peaks(opts: &UnionPeaksOptions) -> Result<()> { n_inputs ); } - let min_support = opts.min_support; + let min_support = opts + .min_support + .max((opts.min_frac_support.unwrap_or(0.0) * n_inputs as f64).ceil() as usize); let gap = opts.window_size as i64 + ISLAND_PAD; let mut n_peaks = 0; for chrom in chrom_lengths.keys() { From cf6624151a66da6ea4d181d86959ecae5be58920 Mon Sep 17 00:00:00 2001 From: "Mitchell R. Vollger" Date: Fri, 21 Aug 2026 20:48:10 -0600 Subject: [PATCH 09/11] refactor: --min-frac-support defaults to 0 instead of Option --- src/cli/union_peaks_opts.rs | 4 ++-- src/subcommands/union_peaks.rs | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/cli/union_peaks_opts.rs b/src/cli/union_peaks_opts.rs index 470bf8074..a7b6a7819 100644 --- a/src/cli/union_peaks_opts.rs +++ b/src/cli/union_peaks_opts.rs @@ -21,8 +21,8 @@ pub struct UnionPeaksOptions { pub min_support: usize, /// Minimum fraction of input BEDs that must overlap a peak for it to be /// reported (0-1). Applied together with --min-support. - #[clap(long, value_parser = frac_in_range)] - pub min_frac_support: Option, + #[clap(long, default_value_t = 0.0, value_parser = frac_in_range)] + pub min_frac_support: f64, /// Rolling window size for finding local maxima (in base pairs). /// Only local maxima are kept, so at most one peak is reported per window. #[clap(long, default_value_t = 200)] diff --git a/src/subcommands/union_peaks.rs b/src/subcommands/union_peaks.rs index c8c65e09d..302655189 100644 --- a/src/subcommands/union_peaks.rs +++ b/src/subcommands/union_peaks.rs @@ -279,7 +279,7 @@ pub fn run_union_peaks(opts: &UnionPeaksOptions) -> Result<()> { } let min_support = opts .min_support - .max((opts.min_frac_support.unwrap_or(0.0) * n_inputs as f64).ceil() as usize); + .max((opts.min_frac_support * n_inputs as f64).ceil() as usize); let gap = opts.window_size as i64 + ISLAND_PAD; let mut n_peaks = 0; for chrom in chrom_lengths.keys() { From 8e678184366e108113cb6afbb460738fd49d5204 Mon Sep 17 00:00:00 2001 From: "Mitchell R. Vollger" Date: Fri, 21 Aug 2026 20:56:50 -0600 Subject: [PATCH 10/11] refactor: PeakCallingParams is the clap source of truth The NucleosomeParameters pattern: the peak-caller knobs live in one #[derive(Args)] struct flattened into CallPeaksOptions (flags unchanged), and union-peaks keeps constructing it directly. Deletes the 12-field From mapping that had to be kept in sync by hand. --- src/cli/call_peaks_opts.rs | 48 +++++++++++++++++------------ src/subcommands/call_peaks/fdr.rs | 2 +- src/subcommands/call_peaks/mod.rs | 17 +++++----- src/subcommands/call_peaks/peaks.rs | 43 ++------------------------ 4 files changed, 41 insertions(+), 69 deletions(-) diff --git a/src/cli/call_peaks_opts.rs b/src/cli/call_peaks_opts.rs index b1b0b9b98..37cd1f686 100644 --- a/src/cli/call_peaks_opts.rs +++ b/src/cli/call_peaks_opts.rs @@ -7,6 +7,9 @@ pub struct CallPeaksOptions { #[clap(flatten)] pub input: InputBam, + #[clap(flatten)] + pub peak_params: PeakCallingParams, + /// BED file with shuffled fiber positions (from bedtools shuffle) /// If not provided, will use all positions as real data (no FDR calculation) #[clap(short, long)] @@ -16,6 +19,31 @@ pub struct CallPeaksOptions { #[clap(short, long, default_value = "-")] pub out: String, + /// Minimum fraction of accessible bases in peak + #[clap(long, default_value = "0.0", hide = true)] + pub min_frac_accessible: f64, + + /// Skip the FDR table generation and use existing table + #[clap(long)] + pub fdr_table: Option, + + /// Output the FDR table to this file + #[clap(long)] + pub fdr_table_out: Option, + + /// Include nucleosome and MSP coverage in pileup (default: only FIRE coverage) + #[clap(long)] + pub include_nuc_msp: bool, + + /// Include haplotype-specific calls + #[clap(long)] + pub haps: bool, +} + +/// The knobs of the shared peak caller. Flattened into `CallPeaksOptions` for the +/// CLI; other callers (union-peaks) construct it directly with their own values. +#[derive(Args, Debug, Clone)] +pub struct PeakCallingParams { /// Maximum coverage threshold for filtering (optional) #[clap(long)] pub max_cov: Option, @@ -46,10 +74,6 @@ pub struct CallPeaksOptions { #[clap(long, default_value = "0.1")] pub min_fire_frac_filter: f64, - /// Minimum fraction of accessible bases in peak - #[clap(long, default_value = "0.0", hide = true)] - pub min_frac_accessible: f64, - /// Rolling window size for finding local maxima (in base pairs) #[clap(long, default_value = "200")] pub window_size: usize, @@ -70,22 +94,6 @@ pub struct CallPeaksOptions { #[clap(long, default_value = "10")] pub max_grouping_iterations: usize, - /// Skip the FDR table generation and use existing table - #[clap(long)] - pub fdr_table: Option, - - /// Output the FDR table to this file - #[clap(long)] - pub fdr_table_out: Option, - - /// Include nucleosome and MSP coverage in pileup (default: only FIRE coverage) - #[clap(long)] - pub include_nuc_msp: bool, - - /// Include haplotype-specific calls - #[clap(long)] - pub haps: bool, - /// Minimum FIRE coverage required to calculate a score (default: 4) #[clap(long, default_value = "4", hide = true)] pub min_fire_coverage: i32, diff --git a/src/subcommands/call_peaks/fdr.rs b/src/subcommands/call_peaks/fdr.rs index ba723c6ab..e3464c2c8 100644 --- a/src/subcommands/call_peaks/fdr.rs +++ b/src/subcommands/call_peaks/fdr.rs @@ -371,7 +371,7 @@ pub fn fdr_table( } // Build the final FDR table - fdr_builder.build(opts.max_fdr) + fdr_builder.build(opts.peak_params.max_fdr) } /// Write FDR table to TSV file diff --git a/src/subcommands/call_peaks/mod.rs b/src/subcommands/call_peaks/mod.rs index ee4616b81..14b37ab39 100644 --- a/src/subcommands/call_peaks/mod.rs +++ b/src/subcommands/call_peaks/mod.rs @@ -18,18 +18,18 @@ pub fn run_call_peaks(opts: &mut CallPeaksOptions) -> Result<()> { log::info!(" Input BAM: {}", opts.input.bam); log::info!(" Output: {}", opts.out); - if let Some(min_frac) = opts.min_fire_frac { + if let Some(min_frac) = opts.peak_params.min_fire_frac { log::info!(" Using FIRE fraction mode: min_fire_frac = {}", min_frac); } else { - log::info!(" Max FDR: {}", opts.max_fdr); + log::info!(" Max FDR: {}", opts.peak_params.max_fdr); } - log::info!(" Window size: {}", opts.window_size); + log::info!(" Window size: {}", opts.peak_params.window_size); let mut bam = opts.input.indexed_bam_reader(); let header = opts.input.header_view(); // Generate or load FDR table (skip if using FIRE fraction mode) - let fdr_table = if opts.min_fire_frac.is_some() { + let fdr_table = if opts.peak_params.min_fire_frac.is_some() { // FIRE fraction mode: use empty FDR table (won't be used for filtering) log::info!(" Skipping FDR calculation (using FIRE fraction threshold)"); Vec::new() @@ -159,19 +159,20 @@ fn process_chromosome_pileup_both( // Apply sd_cov thresholds if max_cov/min_cov are not explicitly set // Match Python behavior: minimum coverage defaults to 4 - let min_cov_threshold = opts.min_cov.unwrap_or_else(|| { - let calculated_min = (median - opts.sd_cov * std_dev).round() as i32; + let min_cov_threshold = opts.peak_params.min_cov.unwrap_or_else(|| { + let calculated_min = (median - opts.peak_params.sd_cov * std_dev).round() as i32; calculated_min.max(DEFAULT_MIN_COVERAGE) }); let max_cov_threshold = opts + .peak_params .max_cov - .unwrap_or_else(|| (median + opts.sd_cov * std_dev).round() as i32); + .unwrap_or_else(|| (median + opts.peak_params.sd_cov * std_dev).round() as i32); log::debug!( " Coverage: median={:.1}, std_dev={:.1} ({:.1} SDs), range=[{}, {}]", median, std_dev, - opts.sd_cov, + opts.peak_params.sd_cov, min_cov_threshold, max_cov_threshold ); diff --git a/src/subcommands/call_peaks/peaks.rs b/src/subcommands/call_peaks/peaks.rs index db887668d..337f09fc0 100644 --- a/src/subcommands/call_peaks/peaks.rs +++ b/src/subcommands/call_peaks/peaks.rs @@ -1,6 +1,7 @@ use super::chrom_names_and_lengths; use super::fdr::{lookup_fdr, FdrEntry}; use crate::cli::CallPeaksOptions; +pub use crate::cli::PeakCallingParams; use crate::fiber::FiberseqData; use crate::subcommands::pileup::{ FiberseqPileup, FiberseqPileupOptions, FireTrack, FireTrackOptions, @@ -46,44 +47,6 @@ fn median(sorted: &[i64]) -> i64 { } } -/// Everything peak calling needs from the CLI, decoupled from `CallPeaksOptions` -/// (which flattens `InputBam`/`FiberFilters`, so commands that never read a BAM -/// cannot build one). See `ft union-peaks` for the other caller. -#[derive(Debug, Clone, Copy)] -pub struct PeakCallingParams { - pub window_size: usize, - pub min_fire_coverage: i32, - pub min_cov: Option, - pub max_cov: Option, - pub sd_cov: f64, - pub max_fdr: f64, - pub min_fire_frac: Option, - pub min_fire_frac_filter: f64, - pub min_frac_overlap: f64, - pub min_reciprocal_overlap: f64, - pub high_reciprocal_overlap: f64, - pub max_grouping_iterations: usize, -} - -impl From<&CallPeaksOptions> for PeakCallingParams { - fn from(o: &CallPeaksOptions) -> Self { - Self { - window_size: o.window_size, - min_fire_coverage: o.min_fire_coverage, - min_cov: o.min_cov, - max_cov: o.max_cov, - sd_cov: o.sd_cov, - max_fdr: o.max_fdr, - min_fire_frac: o.min_fire_frac, - min_fire_frac_filter: o.min_fire_frac_filter, - min_frac_overlap: o.min_frac_overlap, - min_reciprocal_overlap: o.min_reciprocal_overlap, - high_reciprocal_overlap: o.high_reciprocal_overlap, - max_grouping_iterations: o.max_grouping_iterations, - } - } -} - /// Filtering thresholds for peak calling #[derive(Debug, Clone, Copy)] struct PeakThresholds { @@ -727,7 +690,7 @@ pub fn call_peaks( header: &rust_htslib::bam::HeaderView, fdr_table: &[FdrEntry], ) -> Result<()> { - if opts.min_fire_frac.is_some() { + if opts.peak_params.min_fire_frac.is_some() { log::info!("Calling peaks using FIRE fraction threshold"); } else { log::info!( @@ -748,7 +711,7 @@ pub fn call_peaks( let mut writer = bio_io::writer(&opts.out)?; writeln!(writer, "{}", Peak::header())?; - let params = PeakCallingParams::from(&*opts); + let params = opts.peak_params.clone(); let mut total_peaks_before_merge = 0; let mut total_peaks_after_merge = 0; From ec4442b73f1d1d7966d43c5082408178a73f5bcd Mon Sep 17 00:00:00 2001 From: "Mitchell R. Vollger" Date: Fri, 21 Aug 2026 21:00:45 -0600 Subject: [PATCH 11/11] feat: union-peaks exposes the shared merge knobs PeakMergeParams (window size + the three merge overlaps + iteration cap) splits out of PeakCallingParams and flattens into both CLIs; defaults are unchanged on both. The mode/coverage fields stay hardcoded in union-peaks, where FDR and real-fiber coverage semantics do not apply. --- src/cli/call_peaks_opts.rs | 19 ++++++++++++++----- src/cli/union_peaks_opts.rs | 8 +++----- src/subcommands/call_peaks/mod.rs | 2 +- src/subcommands/call_peaks/peaks.rs | 22 +++++++++++----------- src/subcommands/union_peaks.rs | 10 +++------- 5 files changed, 32 insertions(+), 29 deletions(-) diff --git a/src/cli/call_peaks_opts.rs b/src/cli/call_peaks_opts.rs index 37cd1f686..0ab1c4e27 100644 --- a/src/cli/call_peaks_opts.rs +++ b/src/cli/call_peaks_opts.rs @@ -41,9 +41,13 @@ pub struct CallPeaksOptions { } /// The knobs of the shared peak caller. Flattened into `CallPeaksOptions` for the -/// CLI; other callers (union-peaks) construct it directly with their own values. +/// CLI. union-peaks flattens only [`PeakMergeParams`] (shared semantics) and +/// hardcodes the mode/coverage fields, which assume real fibers. #[derive(Args, Debug, Clone)] pub struct PeakCallingParams { + #[clap(flatten)] + pub merge: PeakMergeParams, + /// Maximum coverage threshold for filtering (optional) #[clap(long)] pub max_cov: Option, @@ -74,6 +78,15 @@ pub struct PeakCallingParams { #[clap(long, default_value = "0.1")] pub min_fire_frac_filter: f64, + /// Minimum FIRE coverage required to calculate a score (default: 4) + #[clap(long, default_value = "4", hide = true)] + pub min_fire_coverage: i32, +} + +/// Local-max window and merge geometry: meaningful for any element source, so +/// union-peaks exposes these too. +#[derive(Args, Debug, Clone)] +pub struct PeakMergeParams { /// Rolling window size for finding local maxima (in base pairs) #[clap(long, default_value = "200")] pub window_size: usize, @@ -93,8 +106,4 @@ pub struct PeakCallingParams { /// Maximum number of grouping iterations for merging #[clap(long, default_value = "10")] pub max_grouping_iterations: usize, - - /// Minimum FIRE coverage required to calculate a score (default: 4) - #[clap(long, default_value = "4", hide = true)] - pub min_fire_coverage: i32, } diff --git a/src/cli/union_peaks_opts.rs b/src/cli/union_peaks_opts.rs index a7b6a7819..4bc7867a0 100644 --- a/src/cli/union_peaks_opts.rs +++ b/src/cli/union_peaks_opts.rs @@ -1,4 +1,4 @@ -use crate::cli::GlobalOpts; +use crate::cli::{GlobalOpts, PeakMergeParams}; use clap::Args; use std::fmt::Debug; @@ -23,10 +23,8 @@ pub struct UnionPeaksOptions { /// reported (0-1). Applied together with --min-support. #[clap(long, default_value_t = 0.0, value_parser = frac_in_range)] pub min_frac_support: f64, - /// Rolling window size for finding local maxima (in base pairs). - /// Only local maxima are kept, so at most one peak is reported per window. - #[clap(long, default_value_t = 200)] - pub window_size: usize, + #[clap(flatten)] + pub merge: PeakMergeParams, #[clap(flatten)] pub global: GlobalOpts, } diff --git a/src/subcommands/call_peaks/mod.rs b/src/subcommands/call_peaks/mod.rs index 14b37ab39..72224594c 100644 --- a/src/subcommands/call_peaks/mod.rs +++ b/src/subcommands/call_peaks/mod.rs @@ -23,7 +23,7 @@ pub fn run_call_peaks(opts: &mut CallPeaksOptions) -> Result<()> { } else { log::info!(" Max FDR: {}", opts.peak_params.max_fdr); } - log::info!(" Window size: {}", opts.peak_params.window_size); + log::info!(" Window size: {}", opts.peak_params.merge.window_size); let mut bam = opts.input.indexed_bam_reader(); let header = opts.input.header_view(); diff --git a/src/subcommands/call_peaks/peaks.rs b/src/subcommands/call_peaks/peaks.rs index 337f09fc0..8afcbd189 100644 --- a/src/subcommands/call_peaks/peaks.rs +++ b/src/subcommands/call_peaks/peaks.rs @@ -540,11 +540,11 @@ fn merge_peaks_iterative<'a>( // Phase 1: High reciprocal overlap log::debug!( " Phase 1: Merging peaks with reciprocal overlap >= {}", - params.high_reciprocal_overlap + params.merge.high_reciprocal_overlap ); - for iteration in 0..params.max_grouping_iterations { + for iteration in 0..params.merge.max_grouping_iterations { let prev_count = peaks.len(); - peaks = merge_peaks_single_iteration(peaks, 0.0, params.high_reciprocal_overlap); + peaks = merge_peaks_single_iteration(peaks, 0.0, params.merge.high_reciprocal_overlap); log::debug!( " Iteration {}: {} -> {} peaks", iteration + 1, @@ -560,11 +560,11 @@ fn merge_peaks_iterative<'a>( // Phase 2: FIRE element overlap log::debug!( " Phase 2: Merging peaks with FIRE element overlap >= {}", - params.min_frac_overlap + params.merge.min_frac_overlap ); - for iteration in 0..params.max_grouping_iterations { + for iteration in 0..params.merge.max_grouping_iterations { let prev_count = peaks.len(); - peaks = merge_peaks_single_iteration(peaks, params.min_frac_overlap, 0.0); + peaks = merge_peaks_single_iteration(peaks, params.merge.min_frac_overlap, 0.0); log::debug!( " Iteration {}: {} -> {} peaks", iteration + 1, @@ -580,11 +580,11 @@ fn merge_peaks_iterative<'a>( // Phase 3: High reciprocal overlap again log::debug!( " Phase 3: Merging peaks with reciprocal overlap >= {}", - params.min_reciprocal_overlap + params.merge.min_reciprocal_overlap ); - for iteration in 0..params.max_grouping_iterations { + for iteration in 0..params.merge.max_grouping_iterations { let prev_count = peaks.len(); - peaks = merge_peaks_single_iteration(peaks, 0.0, params.min_reciprocal_overlap); + peaks = merge_peaks_single_iteration(peaks, 0.0, params.merge.min_reciprocal_overlap); log::debug!( " Iteration {}: {} -> {} peaks", iteration + 1, @@ -633,10 +633,10 @@ pub fn call_peaks_for_chrom( shuffle: false, random_shuffle: false, shuffle_seed: None, - rolling_max: Some(params.window_size), + rolling_max: Some(params.merge.window_size), track_fire_elements: true, // Enable FIRE element tracking for peak calling }, - rolling_max: Some(params.window_size), + rolling_max: Some(params.merge.window_size), haps: false, per_base: false, keep_zeros: false, diff --git a/src/subcommands/union_peaks.rs b/src/subcommands/union_peaks.rs index 302655189..c4bfee9fb 100644 --- a/src/subcommands/union_peaks.rs +++ b/src/subcommands/union_peaks.rs @@ -204,7 +204,7 @@ fn support_for<'a>( /// bases. `--min-support` filters the output only, so `-n 3` and `-n 1` plus a downstream /// filter agree apart from the sequential `name` column, which renumbers. pub fn run_union_peaks(opts: &UnionPeaksOptions) -> Result<()> { - if opts.window_size < 2 { + if opts.merge.window_size < 2 { bail!("--window-size must be at least 2; a smaller window finds no local maxima"); } let names = resolve_names(opts)?; @@ -246,7 +246,7 @@ pub fn run_union_peaks(opts: &UnionPeaksOptions) -> Result<()> { // built from the same synthetic fibers), so call in FIRE-fraction mode with the // fraction threshold off and filter on sample support afterwards instead. let params = PeakCallingParams { - window_size: opts.window_size, + merge: opts.merge.clone(), min_fire_coverage: 1, // one sample is enough to score a position min_cov: Some(1), // coverage bounds only set pass_coverage, which we do not emit max_cov: Some(i32::MAX), @@ -254,10 +254,6 @@ pub fn run_union_peaks(opts: &UnionPeaksOptions) -> Result<()> { max_fdr: 1.0, min_fire_frac: Some(0.0), min_fire_frac_filter: 0.0, - min_frac_overlap: 0.5, - min_reciprocal_overlap: 0.75, - high_reciprocal_overlap: 0.90, - max_grouping_iterations: 10, }; let mut writer = bio_io::writer(&opts.out)?; @@ -280,7 +276,7 @@ pub fn run_union_peaks(opts: &UnionPeaksOptions) -> Result<()> { let min_support = opts .min_support .max((opts.min_frac_support * n_inputs as f64).ceil() as usize); - let gap = opts.window_size as i64 + ISLAND_PAD; + let gap = opts.merge.window_size as i64 + ISLAND_PAD; let mut n_peaks = 0; for chrom in chrom_lengths.keys() { let mut all: Vec<(i64, i64)> = samples