From 1d6d56373a8553aa819c3d9e9a09175f2be721e9 Mon Sep 17 00:00:00 2001 From: "Mitchell R. Vollger" Date: Wed, 16 Sep 2026 11:02:55 -0600 Subject: [PATCH 1/2] fix: ft call-peaks --haps fills the H1/H2 columns The flag was parsed on CallPeaksOptions but call_peaks_for_chrom hardcoded haps: false, so the pileup never built haplotype tracks and every H1/H2 column was zero. Move haps into PeakCallingParams so it reaches the pileup. union-peaks sets it false since BED intervals carry no HP tag. Closes #140 --- src/cli/call_peaks_opts.rs | 8 ++++---- src/subcommands/call_peaks/peaks.rs | 2 +- src/subcommands/union_peaks.rs | 1 + tests/regression/call_peaks.rs | 19 +++++++++++++++++++ 4 files changed, 25 insertions(+), 5 deletions(-) diff --git a/src/cli/call_peaks_opts.rs b/src/cli/call_peaks_opts.rs index 0ab1c4e2..28ffb6e7 100644 --- a/src/cli/call_peaks_opts.rs +++ b/src/cli/call_peaks_opts.rs @@ -34,10 +34,6 @@ pub struct CallPeaksOptions { /// 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 @@ -81,6 +77,10 @@ pub struct PeakCallingParams { /// Minimum FIRE coverage required to calculate a score (default: 4) #[clap(long, default_value = "4", hide = true)] pub min_fire_coverage: i32, + + /// Include haplotype-specific calls + #[clap(long)] + pub haps: bool, } /// Local-max window and merge geometry: meaningful for any element source, so diff --git a/src/subcommands/call_peaks/peaks.rs b/src/subcommands/call_peaks/peaks.rs index 4bb5be8f..c40e4457 100644 --- a/src/subcommands/call_peaks/peaks.rs +++ b/src/subcommands/call_peaks/peaks.rs @@ -637,7 +637,7 @@ pub fn call_peaks_for_chrom( track_fire_elements: true, // Enable FIRE element tracking for peak calling }, rolling_max: Some(params.merge.window_size), - haps: false, + haps: params.haps, per_base: false, keep_zeros: false, min_fire_coverage: Some(params.min_fire_coverage), diff --git a/src/subcommands/union_peaks.rs b/src/subcommands/union_peaks.rs index c4bfee9f..15be4eae 100644 --- a/src/subcommands/union_peaks.rs +++ b/src/subcommands/union_peaks.rs @@ -254,6 +254,7 @@ pub fn run_union_peaks(opts: &UnionPeaksOptions) -> Result<()> { max_fdr: 1.0, min_fire_frac: Some(0.0), min_fire_frac_filter: 0.0, + haps: false, // BED intervals carry no HP tag }; let mut writer = bio_io::writer(&opts.out)?; diff --git a/tests/regression/call_peaks.rs b/tests/regression/call_peaks.rs index 39fa148c..018fd747 100644 --- a/tests/regression/call_peaks.rs +++ b/tests/regression/call_peaks.rs @@ -12,3 +12,22 @@ fn call_peaks_ctcf_snapshot() { ]); insta::assert_snapshot!(out); } + +/// `--haps` must fill the H1/H2 columns from HP tags (issue #140: the flag was +/// parsed but never reached the pileup, so every haplotype column was zero). +#[test] +fn call_peaks_haps_fills_haplotype_columns() { + let out = run(&[ + "call-peaks", + fixture("NAPA.bam").to_str().unwrap(), + "--haps", + "--min-fire-frac", + "0.5", + ]); + let peak = out.lines().find(|l| !l.starts_with('#')).expect("one peak"); + let cols: Vec<&str> = peak.split('\t').collect(); + // coverage, coverage_H1, coverage_H2 + assert_eq!(cols[5], "95"); + assert_eq!(cols[10], "45"); + assert_eq!(cols[15], "14"); +} From 6957ae33aa7686800b97cb69760660dc6788fc12 Mon Sep 17 00:00:00 2001 From: "Mitchell R. Vollger" Date: Fri, 18 Sep 2026 09:36:09 -0600 Subject: [PATCH 2/2] perf: skip FIRE element tracking on the haplotype pileup tracks Only all_data's fire_elements are read (peak boundaries in call-peaks). With --haps now live, the H1/H2 tracks were allocating a Vec per base that nothing used, about a third of the tripled call-peaks memory. --- src/subcommands/pileup.rs | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/subcommands/pileup.rs b/src/subcommands/pileup.rs index c8457845..9c6cc0d2 100644 --- a/src/subcommands/pileup.rs +++ b/src/subcommands/pileup.rs @@ -699,19 +699,25 @@ impl<'a> FiberseqPileup<'a> { &None, ); let (hap1_data, hap2_data) = if pileup_opts.haps { + // Only all_data's FIRE elements are ever read (peak boundaries in + // call-peaks), so do not pay for them on the haplotype tracks. + let hap_opts = FireTrackOptions { + track_fire_elements: false, + ..fire_track_opts.clone() + }; ( Some(FireTrack::new( chrom.to_string(), chrom_start, chrom_end, - fire_track_opts.clone(), + hap_opts.clone(), &None, )), Some(FireTrack::new( chrom.to_string(), chrom_start, chrom_end, - fire_track_opts.clone(), + hap_opts, &None, )), )