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

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

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;
 
 //
@@ -163,8 +165,13 @@ 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
+    /// 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),

src/cli/call_peaks_opts.rs

@@ -7,6 +7,9 @@ pub struct CallPeaksOptions {
     #[clap(flatten)]
     pub input: InputBam<CallableFibers>,
 
+    #[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,35 @@ 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<String>,
+
+    /// Output the FDR table to this file
+    #[clap(long)]
+    pub fdr_table_out: Option<String>,
+
+    /// 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. 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<i32>,
@@ -46,10 +78,15 @@ 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,
+    /// 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,
@@ -69,24 +106,4 @@ pub struct CallPeaksOptions {
     /// Maximum number of grouping iterations for merging
     #[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<String>,
-
-    /// Output the FDR table to this file
-    #[clap(long)]
-    pub fdr_table_out: Option<String>,
-
-    /// 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,
 }

src/cli/union_peaks_opts.rs

@@ -0,0 +1,39 @@
+use crate::cli::{GlobalOpts, PeakMergeParams};
+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 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<String>,
+    /// 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<String>,
+    /// 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, default_value_t = 0.0, value_parser = frac_in_range)]
+    pub min_frac_support: f64,
+    #[clap(flatten)]
+    pub merge: PeakMergeParams,
+    #[clap(flatten)]
+    pub global: GlobalOpts,
+}
+
+fn frac_in_range(s: &str) -> Result<f64, String> {
+    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())
+    }
+}

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)?;
         }

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;

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

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};
@@ -16,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.merge.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()
@@ -157,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
     );