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..aea581e9a 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; // @@ -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), diff --git a/src/cli/call_peaks_opts.rs b/src/cli/call_peaks_opts.rs index b1b0b9b98..0ab1c4e27 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,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, + + /// 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. 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, @@ -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, - - /// 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/cli/union_peaks_opts.rs b/src/cli/union_peaks_opts.rs new file mode 100644 index 000000000..4bc7867a0 --- /dev/null +++ b/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, + /// 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, + /// 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 { + 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/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/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 6eb94d728..72224594c 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}; @@ -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 ); diff --git a/src/subcommands/call_peaks/peaks.rs b/src/subcommands/call_peaks/peaks.rs index 3a9dbfc83..8afcbd189 100644 --- a/src/subcommands/call_peaks/peaks.rs +++ b/src/subcommands/call_peaks/peaks.rs @@ -1,6 +1,8 @@ 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, }; @@ -34,6 +36,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 + } +} + /// Filtering thresholds for peak calling #[derive(Debug, Clone, Copy)] struct PeakThresholds { @@ -273,12 +286,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 @@ -518,17 +531,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.merge.high_reciprocal_overlap ); - for iteration in 0..opts.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, opts.high_reciprocal_overlap); + peaks = merge_peaks_single_iteration(peaks, 0.0, params.merge.high_reciprocal_overlap); log::debug!( " Iteration {}: {} -> {} peaks", iteration + 1, @@ -544,11 +560,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.merge.min_frac_overlap ); - for iteration in 0..opts.max_grouping_iterations { + for iteration in 0..params.merge.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.merge.min_frac_overlap, 0.0); log::debug!( " Iteration {}: {} -> {} peaks", iteration + 1, @@ -564,11 +580,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.merge.min_reciprocal_overlap ); - for iteration in 0..opts.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, opts.min_reciprocal_overlap); + peaks = merge_peaks_single_iteration(peaks, 0.0, params.merge.min_reciprocal_overlap); log::debug!( " Iteration {}: {} -> {} peaks", iteration + 1, @@ -591,6 +607,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.merge.window_size), + track_fire_elements: true, // Enable FIRE element tracking for peak calling + }, + rolling_max: Some(params.merge.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: @@ -605,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!( @@ -626,6 +711,7 @@ pub fn call_peaks( let mut writer = bio_io::writer(&opts.out)?; writeln!(writer, "{}", Peak::header())?; + let params = opts.peak_params.clone(); let mut total_peaks_before_merge = 0; let mut total_peaks_after_merge = 0; @@ -642,60 +728,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 +747,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..152f877ab 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, @@ -141,10 +141,10 @@ 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 @@ 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 diff --git a/src/subcommands/union_peaks.rs b/src/subcommands/union_peaks.rs new file mode 100644 index 000000000..c4bfee9fb --- /dev/null +++ b/src/subcommands/union_peaks.rs @@ -0,0 +1,420 @@ +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 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. 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; +/// 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. +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 || rec.end > MAX_COORD { + bail!( + "invalid interval in {} (coordinates must be 0 <= start < end <= {}): {} {} {}", + path, + MAX_COORD, + 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 (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(); + 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 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>( + 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`, 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.merge.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 { + 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), + sd_cov: 5.0, + max_fdr: 1.0, + min_fire_frac: Some(0.0), + min_fire_frac_filter: 0.0, + }; + + 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(); + 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 + .max((opts.min_frac_support * n_inputs as f64).ceil() as usize); + 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 + .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 > 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" + ); + } + // 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 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() + .take_while(|(start, _)| *start < island_end) + .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/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.rs b/tests/regression.rs index 43eb9467a..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"] @@ -20,3 +22,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/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..ce0a653b9 --- /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 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 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 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 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 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 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 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 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 new file mode 100644 index 000000000..efcdb674f --- /dev/null +++ b/tests/regression/union_peaks.rs @@ -0,0 +1,190 @@ +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\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" + ); +} + +// --min-support is an output filter, so it must drop exactly the low support rows and +// leave the surviving rows unchanged apart from the sequential name, which renumbers. +#[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", "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"); + 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}"); + } +} + +// 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}"); +} + +// --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}"); +}