From 2ad088341114ce3ad04b4132f0df79e054876c10 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 09:39:43 -0700 Subject: [PATCH 01/10] Supply CRAM references to legacy sampling --- modkit-core/src/command_utils.rs | 43 +- modkit-core/src/entropy/mod.rs | 12 +- modkit-core/src/entropy/subcommand.rs | 26 +- modkit-core/src/modbam_util/subcommands.rs | 53 ++- modkit-core/src/pileup/duplex.rs | 14 +- modkit-core/src/pileup/subcommand.rs | 83 ++-- modkit-core/src/reads_sampler/mod.rs | 130 +++++- .../src/reads_sampler/sampling_schedule.rs | 82 +++- modkit-core/src/thresholds.rs | 81 +++- modkit-core/src/util.rs | 57 ++- modkit/tests/test_cram_reference_consumers.rs | 373 ++++++++++++++++++ 11 files changed, 857 insertions(+), 97 deletions(-) create mode 100644 modkit/tests/test_cram_reference_consumers.rs diff --git a/modkit-core/src/command_utils.rs b/modkit-core/src/command_utils.rs index 29ee4689..19736d0f 100644 --- a/modkit-core/src/command_utils.rs +++ b/modkit-core/src/command_utils.rs @@ -13,7 +13,7 @@ use crate::mod_base_code::{DnaBase, ModCodeRepr}; use crate::motifs::motif_bed::RegexMotif; use crate::position_filter::StrandedPositionFilter; use crate::threshold_mod_caller::MultipleThresholdModCaller; -use crate::thresholds::calc_threshold_from_bam; +use crate::thresholds::calc_threshold_from_bam_with_reference; use crate::util::{create_out_directory, Region}; pub fn parse_per_mod_thresholds( @@ -141,6 +141,44 @@ pub fn get_threshold_from_options( position_filter: Option<&StrandedPositionFilter<()>>, only_mapped: bool, suppress_progress: bool, +) -> anyhow::Result { + get_threshold_from_options_with_reference( + in_bam, + None, + threads, + interval_size, + sample_frac, + num_reads, + no_filtering, + filter_percentile, + seed, + region, + per_mod_thresholds, + edge_filter, + collapse_method, + position_filter, + only_mapped, + suppress_progress, + ) +} + +pub(crate) fn get_threshold_from_options_with_reference( + in_bam: &PathBuf, + reference_fasta: Option<&PathBuf>, + threads: usize, + interval_size: u32, + sample_frac: Option, + num_reads: usize, + no_filtering: bool, + filter_percentile: f32, + seed: Option, + region: Option<&Region>, + per_mod_thresholds: Option>, + edge_filter: Option<&EdgeFilter>, + collapse_method: Option<&CollapseMethod>, + position_filter: Option<&StrandedPositionFilter<()>>, + only_mapped: bool, + suppress_progress: bool, ) -> anyhow::Result { if no_filtering { info!("not performing filtering"); @@ -157,8 +195,9 @@ pub fn get_threshold_from_options( (None, Some(num_reads)) } }; - let per_base_thresholds = calc_threshold_from_bam( + let per_base_thresholds = calc_threshold_from_bam_with_reference( in_bam, + reference_fasta, threads, interval_size, sample_frac, diff --git a/modkit-core/src/entropy/mod.rs b/modkit-core/src/entropy/mod.rs index 5de1ee4b..3626e262 100644 --- a/modkit-core/src/entropy/mod.rs +++ b/modkit-core/src/entropy/mod.rs @@ -25,7 +25,10 @@ use crate::read_ids_to_base_mod_probs::{PositionModCalls, ReadBaseModProfile}; use crate::reads_sampler::sampling_schedule::ReferenceSequencesLookup; use crate::threshold_mod_caller::MultipleThresholdModCaller; use crate::thresholds::percentile_linear_interp; -use crate::util::{record_is_not_primary, ReferenceRecord, Strand}; +use crate::util::{ + record_is_not_primary, set_reference_for_cram_indexed_reader, + ReferenceRecord, Strand, +}; mod methylation_entropy; pub mod subcommand; @@ -1486,11 +1489,13 @@ struct Message { fn process_bam_fp( bam_fp: &PathBuf, + reference_fasta: &PathBuf, fetch_definition: FetchDefinition, caller: Arc, io_threads: usize, ) -> anyhow::Result> { let mut reader = bam::IndexedReader::from_path(bam_fp)?; + set_reference_for_cram_indexed_reader(&mut reader, Some(reference_fasta))?; reader.set_threads(io_threads)?; reader.fetch(fetch_definition)?; @@ -1583,9 +1588,11 @@ pub(super) fn process_entropy_window( io_threads: usize, caller: Arc, bam_fps: &[PathBuf], + reference_fasta: &PathBuf, ) -> anyhow::Result { let bam_fp = &bam_fps[0]; - let reader = bam::IndexedReader::from_path(bam_fp)?; + let mut reader = bam::IndexedReader::from_path(bam_fp)?; + set_reference_for_cram_indexed_reader(&mut reader, Some(reference_fasta))?; let chrom_id = entropy_windows.chrom_id; drop(reader); @@ -1594,6 +1601,7 @@ pub(super) fn process_entropy_window( .map(|fp| { process_bam_fp( fp, + reference_fasta, entropy_windows.get_fetch_definition(), caller.clone(), io_threads, diff --git a/modkit-core/src/entropy/subcommand.rs b/modkit-core/src/entropy/subcommand.rs index cbfe1257..c16a35c7 100644 --- a/modkit-core/src/entropy/subcommand.rs +++ b/modkit-core/src/entropy/subcommand.rs @@ -22,8 +22,8 @@ use crate::reads_sampler::sampling_schedule::{ }; use crate::threshold_mod_caller::MultipleThresholdModCaller; use crate::thresholds::{ - calculate_threshold_with_fallback, get_modbase_probs_from_bam, - log_calculated_thresholds, + calculate_threshold_with_fallback, + get_modbase_probs_from_bam_with_reference, log_calculated_thresholds, }; use crate::util::{ format_errors_table, get_master_progress_bar, get_ticker, MutOpMax, @@ -194,13 +194,18 @@ impl MethylationEntropy { bail!("min-valid-coverage must be at least 1") } for bam_fp in self.in_bams.iter() { - IdxStats::check_any_mapped_reads(&bam_fp, None, None) - .with_context(|| { - format!( - "did not find any mapped reads in {bam_fp:?}, perform \ + IdxStats::check_any_mapped_reads_with_reference( + &bam_fp, + Some(&self.reference_fasta), + None, + None, + ) + .with_context(|| { + format!( + "did not find any mapped reads in {bam_fp:?}, perform \ alignment first" - ) - })?; + ) + })?; } let mut writer: Box = @@ -341,6 +346,7 @@ impl MethylationEntropy { let (snd, rcv) = crossbeam::channel::bounded(10_000); let bam_fps = self.in_bams.clone(); + let reference_fasta = self.reference_fasta.clone(); let min_coverage = self.min_valid_coverage; let threads = self.threads; let io_threads = self.io_threads.unwrap_or(threads); @@ -387,6 +393,7 @@ impl MethylationEntropy { io_threads, threshold_caller.clone(), &bam_fps, + &reference_fasta, ) }) .collect::>(); @@ -482,8 +489,9 @@ impl MethylationEntropy { HashMap::::new(); for in_bam in self.in_bams.iter() { let (per_base_thresholds, explicit_canonical_probs) = - get_modbase_probs_from_bam( + get_modbase_probs_from_bam_with_reference( in_bam, + Some(&self.reference_fasta), self.threads, 1_000_000, None, diff --git a/modkit-core/src/modbam_util/subcommands.rs b/modkit-core/src/modbam_util/subcommands.rs index 509b270b..feda6a8a 100644 --- a/modkit-core/src/modbam_util/subcommands.rs +++ b/modkit-core/src/modbam_util/subcommands.rs @@ -20,8 +20,8 @@ use rustc_hash::FxHashMap; use crate::adjust::adjust_modbam; use crate::command_utils::{ get_bam_writer, get_motif_lookup_from_parts, get_serial_reader, - get_threshold_from_options, parse_edge_filter_input, parse_forward_motifs, - parse_per_mod_thresholds, parse_raw_motifs, + get_threshold_from_options_with_reference, parse_edge_filter_input, + parse_forward_motifs, parse_per_mod_thresholds, parse_raw_motifs, parse_raw_thresholds_string_with_default, parse_thresholds, using_stream, }; use crate::errs::{MkError, MkResult}; @@ -54,7 +54,8 @@ use crate::summarize::ModSummary; use crate::util::{ add_modkit_pg_records, filter_reference_records, format_errors_table, get_master_progress_bar, get_subroutine_progress_bar, get_targets, - get_ticker, ReferenceRecord, Region, DEFAULT_NUM_READS, + get_ticker, preflight_cram_input, reader_is_cram, + set_reference_for_cram_reader, ReferenceRecord, Region, DEFAULT_NUM_READS, }; use crate::writers::{ MultiTableWriter, OutWriter, SampledProbs, TableWriter, TsvWriter, @@ -529,6 +530,9 @@ pub struct Adjust { /// File path to new BAM file to be created. Can be a path to a file or one /// of `-` or `stdin` to specify a stream from standard output. out_bam: String, + /// Reference sequence in FASTA format for CRAM decoding. (alias: 'ref') + #[arg(long = "reference", alias = "ref", short = 'r')] + reference_fasta: Option, /// Output debug logs to file at this path. #[clap(help_heading = "Output Options")] #[arg(long, alias = "log")] @@ -705,12 +709,19 @@ impl Adjust { let _handle = init_logging(self.log_filepath.as_ref()); let io_threadpool = tpool::ThreadPool::new(self.threads as u32)?; let mut reader = get_serial_reader(self.in_bam.as_str())?; + set_reference_for_cram_reader( + &mut reader, + self.reference_fasta.as_ref(), + )?; + if !using_stream(&self.in_bam) && reader_is_cram(&reader) { + preflight_cram_input( + Path::new(&self.in_bam), + self.reference_fasta.as_ref(), + )?; + } reader.set_thread_pool(&io_threadpool)?; let mut header = bam::Header::from_template(reader.header()); add_modkit_pg_records(&mut header); - let mut bam_writer = - get_bam_writer(&self.out_bam, &header, self.output_sam)?; - bam_writer.set_thread_pool(&io_threadpool)?; let methods = if let Some(convert) = &self.convert { let convert = convert @@ -829,8 +840,9 @@ impl Adjust { .build() .with_context(|| "failed to make threadpool")?; pool.install(|| { - get_threshold_from_options( + get_threshold_from_options_with_reference( &Path::new(&self.in_bam).to_path_buf(), + self.reference_fasta.as_ref(), self.threads, self.sampling_interval_size, None, @@ -853,6 +865,10 @@ impl Adjust { None }; + let mut bam_writer = + get_bam_writer(&self.out_bam, &header, self.output_sam)?; + bam_writer.set_thread_pool(&io_threadpool)?; + adjust_modbam( &mut reader, &mut bam_writer, @@ -2205,6 +2221,9 @@ pub struct CallMods { /// Output BAM, can be a path to a file or one of `-` or /// `stdin` to specify a stream from standard input. out_bam: String, + /// Reference sequence in FASTA format for CRAM decoding. (alias: 'ref') + #[arg(long = "reference", alias = "ref", short = 'r')] + reference_fasta: Option, /// Specify a file for debug logs to be written to, otherwise ignore them. /// Setting a file is recommended. #[arg(long, alias = "log")] @@ -2372,12 +2391,19 @@ impl CallMods { let _handle = init_logging(self.log_filepath.as_ref()); let io_threadpool = tpool::ThreadPool::new(self.threads as u32)?; let mut reader = get_serial_reader(&self.in_bam)?; + set_reference_for_cram_reader( + &mut reader, + self.reference_fasta.as_ref(), + )?; + if !using_stream(&self.in_bam) && reader_is_cram(&reader) { + preflight_cram_input( + Path::new(&self.in_bam), + self.reference_fasta.as_ref(), + )?; + } reader.set_thread_pool(&io_threadpool)?; let mut header = bam::Header::from_template(reader.header()); add_modkit_pg_records(&mut header); - let mut bam_writer = - get_bam_writer(&self.out_bam, &header, self.output_sam)?; - bam_writer.set_thread_pool(&io_threadpool)?; let edge_filter = self .edge_filter @@ -2419,8 +2445,9 @@ impl CallMods { .build() .with_context(|| "failed to make threadpool")?; pool.install(|| { - get_threshold_from_options( + get_threshold_from_options_with_reference( &Path::new(&self.in_bam).to_path_buf(), + self.reference_fasta.as_ref(), self.threads, self.sampling_interval_size, self.sampling_frac, @@ -2439,6 +2466,10 @@ impl CallMods { })? }; + let mut bam_writer = + get_bam_writer(&self.out_bam, &header, self.output_sam)?; + bam_writer.set_thread_pool(&io_threadpool)?; + adjust_modbam( &mut reader, &mut bam_writer, diff --git a/modkit-core/src/pileup/duplex.rs b/modkit-core/src/pileup/duplex.rs index d86ba76b..3eb40fc6 100644 --- a/modkit-core/src/pileup/duplex.rs +++ b/modkit-core/src/pileup/duplex.rs @@ -1,6 +1,6 @@ use std::cmp::Ordering; use std::collections::HashMap; -use std::path::Path; +use std::path::{Path, PathBuf}; use anyhow::bail; use derive_new::new; @@ -13,7 +13,9 @@ use crate::motifs::motif_bed::MotifInfo; use crate::pileup::{get_forward_read_base, PileupIter, PileupNumericOptions}; use crate::read_cache::DuplexReadCache; use crate::threshold_mod_caller::MultipleThresholdModCaller; -use crate::util::record_is_not_primary; +use crate::util::{ + record_is_not_primary, set_reference_for_cram_indexed_reader, +}; /// Summarizes the duplex (hemi) methylation patterns for /// a genomic interval @@ -206,9 +208,12 @@ impl DuplexFeatureVector { // todo this function should be removed in favor of a more // generic version in pileup/mod.rs -pub fn process_region_duplex_batch + Copy>( +pub(crate) fn process_region_duplex_batch_with_reference< + T: AsRef + Copy, +>( chromosome_coordintes: &MultiChromCoordinates, bam_fp: T, + reference_fasta: Option<&PathBuf>, caller: &MultipleThresholdModCaller, pileup_numeric_options: &PileupNumericOptions, force_allow: bool, @@ -222,6 +227,7 @@ pub fn process_region_duplex_batch + Copy>( .map(|chrom_coords| { process_region_duplex( bam_fp, + reference_fasta, chrom_coords.chrom_tid, chrom_coords.start_pos, chrom_coords.end_pos, @@ -239,6 +245,7 @@ pub fn process_region_duplex_batch + Copy>( fn process_region_duplex>( bam_fp: T, + reference_fasta: Option<&PathBuf>, chrom_tid: u32, start_pos: u32, end_pos: u32, @@ -261,6 +268,7 @@ fn process_region_duplex>( }; let mut bam_reader = bam::IndexedReader::from_path(bam_fp)?; + set_reference_for_cram_indexed_reader(&mut bam_reader, reference_fasta)?; let chrom_name = String::from_utf8_lossy(bam_reader.header().tid2name(chrom_tid)) .to_string(); diff --git a/modkit-core/src/pileup/subcommand.rs b/modkit-core/src/pileup/subcommand.rs index 4503f99e..01d5961a 100644 --- a/modkit-core/src/pileup/subcommand.rs +++ b/modkit-core/src/pileup/subcommand.rs @@ -16,7 +16,7 @@ use rust_htslib::bam::{self, HeaderView, Read}; use modkit_logging::init_logging; use crate::command_utils::{ - get_motif_lookup_from_parts, get_threshold_from_options, + get_motif_lookup_from_parts, get_threshold_from_options_with_reference, parse_edge_filter_input, parse_per_base_thresholds, parse_per_mod_thresholds, parse_raw_motifs, parse_raw_thresholds_string_with_default, parse_thresholds, @@ -33,7 +33,9 @@ use crate::mod_base_code::{ }; use crate::motifs::motif_bed::{MotifInfo, RegexMotif}; use crate::pileup::bedrmod::BedRModArgs; -use crate::pileup::duplex::{process_region_duplex_batch, DuplexModBasePileup}; +use crate::pileup::duplex::{ + process_region_duplex_batch_with_reference, DuplexModBasePileup, +}; use crate::pileup::pileup_processor::{ CountsMatrix, DnaAllContext, DnaCpGCombineStrands, DnaCytosineCombine, DnaModOption, DnaPileupWorker, Dynamic, GenericPileupWorker, PileupWorker, @@ -50,7 +52,8 @@ use crate::sample_probs::{ use crate::util::{ create_out_directory, filter_reference_records, get_master_progress_bar, get_master_progress_bar_fancy, get_subroutine_progress_bar, get_targets, - get_ticker, reader_is_bam, reader_is_cram, Region, + get_ticker, reader_is_bam, reader_is_cram, + set_reference_for_cram_indexed_reader, Region, }; use crate::writers::{ BedMethylWriter, BedMethylWriter2, MultipleMotifBedmethylWriter, @@ -1991,17 +1994,19 @@ impl DuplexModBamPileup { ); } // do this first so we fail when the file isn't readable - let header = - bam::IndexedReader::from_path(&self.in_bam).map(|reader| { - if !reader_is_bam(&reader) { - info!( - "\ - detected non-BAM input format, please consider using BAM, \ - CRAM may be unstable" - ); - } - reader.header().to_owned() - })?; + let mut reader = bam::IndexedReader::from_path(&self.in_bam)?; + set_reference_for_cram_indexed_reader( + &mut reader, + Some(&self.reference_fasta), + )?; + if !reader_is_bam(&reader) { + info!( + "\ + detected non-BAM input format, please consider using BAM, \ + CRAM may be unstable" + ); + } + let header = reader.header().to_owned(); // options parsing below let region = self @@ -2052,8 +2057,9 @@ impl DuplexModBamPileup { .transpose()?; // use the path here instead of passing the reader directly to avoid // potentially changing mutable internal state of the reader. - IdxStats::check_any_mapped_reads( + IdxStats::check_any_mapped_reads_with_reference( &self.in_bam, + Some(&self.reference_fasta), region.as_ref(), position_filter.as_ref(), ) @@ -2132,26 +2138,6 @@ impl DuplexModBamPileup { bail!("motif must be palindromic for pileup-hemi") } - let mut writer: Box> = - if let Some(out_fp) = self.out_bed.as_ref() { - create_out_directory(out_fp)?; - let fh = std::fs::File::create(out_fp) - .context("failed to make output file")?; - let writer = BufWriter::new(fh); - Box::new(BedMethylWriter::new( - writer, - self.mixed_delimiters, - false, - )?) - } else { - let writer = BufWriter::new(std::io::stdout()); - Box::new(BedMethylWriter::new( - writer, - self.mixed_delimiters, - false, - )?) - }; - let pool = rayon::ThreadPoolBuilder::new() .num_threads(self.threads) .build() @@ -2173,8 +2159,9 @@ impl DuplexModBamPileup { parse_thresholds(raw_threshold, per_mod_thresholds)? } else { pool.install(|| { - get_threshold_from_options( + get_threshold_from_options_with_reference( &self.in_bam, + Some(&self.reference_fasta), self.threads, self.sampling_interval_size, self.sampling_frac, @@ -2236,6 +2223,26 @@ impl DuplexModBamPileup { } } + let mut writer: Box> = + if let Some(out_fp) = self.out_bed.as_ref() { + create_out_directory(out_fp)?; + let fh = std::fs::File::create(out_fp) + .context("failed to make output file")?; + let writer = BufWriter::new(fh); + Box::new(BedMethylWriter::new( + writer, + self.mixed_delimiters, + false, + )?) + } else { + let writer = BufWriter::new(std::io::stdout()); + Box::new(BedMethylWriter::new( + writer, + self.mixed_delimiters, + false, + )?) + }; + let (snd, rx) = bounded(self.queue_size); let reference_records = if let Some(pf) = position_filter.as_ref() { pf.optimize_reference_records(reference_records, self.interval_size) @@ -2252,6 +2259,7 @@ impl DuplexModBamPileup { )?; let in_bam_fp = self.in_bam.clone(); + let reference_fasta = self.reference_fasta.clone(); let master_progress = MultiProgress::new(); if self.suppress_progress { @@ -2299,9 +2307,10 @@ impl DuplexModBamPileup { .into_par_iter() .progress_with(chunk_progress) .map(|multi_chrom_coords| { - process_region_duplex_batch( + process_region_duplex_batch_with_reference( multi_chrom_coords, &in_bam_fp, + Some(&reference_fasta), &threshold_caller, &pileup_options, force_allow, diff --git a/modkit-core/src/reads_sampler/mod.rs b/modkit-core/src/reads_sampler/mod.rs index 7dfc6b4b..191e1e35 100644 --- a/modkit-core/src/reads_sampler/mod.rs +++ b/modkit-core/src/reads_sampler/mod.rs @@ -20,7 +20,9 @@ use crate::reads_sampler::sampling_schedule::{ }; use crate::record_processor::{RecordProcessor, WithRecords}; use crate::util::{ - get_master_progress_bar, get_targets, get_ticker, ReferenceRecord, Region, + get_master_progress_bar, get_targets, get_ticker, + set_reference_for_cram_indexed_reader, set_reference_for_cram_reader, + ReferenceRecord, Region, }; use record_sampler::RecordSampler; @@ -41,6 +43,43 @@ pub fn get_sampled_read_ids_to_base_mod_probs( only_mapped: bool, suppress_progress: bool, ) -> anyhow::Result +where + P::Output: Moniod + WithRecords, +{ + get_sampled_read_ids_to_base_mod_probs_with_reference::

( + bam_fp, + None, + reader_threads, + interval_size, + sample_frac, + num_reads, + seed, + region, + collapse_method, + edge_filter, + position_filter, + only_mapped, + suppress_progress, + ) +} + +pub(crate) fn get_sampled_read_ids_to_base_mod_probs_with_reference< + P: RecordProcessor, +>( + bam_fp: &PathBuf, + reference_fasta: Option<&PathBuf>, + reader_threads: usize, + interval_size: u32, + sample_frac: Option, + num_reads: Option, + seed: Option, + region: Option<&Region>, + collapse_method: Option<&CollapseMethod>, + edge_filter: Option<&EdgeFilter>, + position_filter: Option<&StrandedPositionFilter<()>>, + only_mapped: bool, + suppress_progress: bool, +) -> anyhow::Result where P::Output: Moniod + WithRecords, { @@ -51,22 +90,29 @@ where chunks" ); let schedule = match (sample_frac, num_reads) { - (_, Some(num_reads)) => SamplingSchedule::from_num_reads( - bam_fp, - num_reads, - region, - position_filter, - !only_mapped, - ), - (Some(frac), _) => SamplingSchedule::from_sample_frac( - bam_fp, - frac as f32, - region, - position_filter, - !only_mapped, - ), - (None, None) => SamplingSchedule::from_sample_frac( + (_, Some(num_reads)) => { + SamplingSchedule::from_num_reads_with_reference( + bam_fp, + reference_fasta, + num_reads, + region, + position_filter, + !only_mapped, + ) + } + (Some(frac), _) => { + SamplingSchedule::from_sample_frac_with_reference( + bam_fp, + reference_fasta, + frac as f32, + region, + position_filter, + !only_mapped, + ) + } + (None, None) => SamplingSchedule::from_sample_frac_with_reference( bam_fp, + reference_fasta, 1.0, region, position_filter, @@ -76,6 +122,7 @@ where let mut read_ids_to_base_mod_calls = sample_reads_base_mod_calls_over_regions::

( bam_fp, + reference_fasta, interval_size, (reader_threads as f32 * 1.5).floor() as usize, region, @@ -94,6 +141,10 @@ where read_ids_to_base_mod_calls.len() ); let mut reader = bam::IndexedReader::from_path(bam_fp)?; + set_reference_for_cram_indexed_reader( + &mut reader, + reference_fasta, + )?; reader.set_threads(reader_threads)?; reader.fetch(bam::FetchDefinition::Unmapped)?; let num_reads_unmapped = num_reads.map(|nr| { @@ -138,6 +189,7 @@ where ); } let mut reader = bam::Reader::from_path(bam_fp)?; + set_reference_for_cram_reader(&mut reader, reference_fasta)?; reader.set_threads(reader_threads)?; let record_sampler = RecordSampler::new_from_options(sample_frac, num_reads, seed); @@ -162,6 +214,7 @@ where /// an entire sorted, aligned BAM. Only uses primary alignments fn sample_reads_base_mod_calls_over_regions( bam_fp: &PathBuf, + reference_fasta: Option<&PathBuf>, interval_size: u32, batch_size: usize, region: Option<&Region>, @@ -175,7 +228,8 @@ fn sample_reads_base_mod_calls_over_regions( where P::Output: Moniod + WithRecords, { - let reader = bam::IndexedReader::from_path(bam_fp)?; + let mut reader = bam::IndexedReader::from_path(bam_fp)?; + set_reference_for_cram_indexed_reader(&mut reader, reference_fasta)?; let header = reader.header(); let contigs = get_targets(header, region) @@ -232,6 +286,7 @@ where .map(|multi_coords| { run_batch::

( bam_fp, + reference_fasta, multi_coords, sampling_schedule, collapse_method, @@ -258,6 +313,7 @@ where fn run_batch( bam_fp: &PathBuf, + reference_fasta: Option<&PathBuf>, batch: Vec<(ChromCoordinates, CountOrSample)>, sampling_schedule: &SamplingSchedule, collapse_method: Option<&CollapseMethod>, @@ -293,9 +349,9 @@ where } CountOrSample::All => RecordSampler::new_passthrough(), }; - - match sample_reads_from_interval::

( + match sample_reads_from_interval_with_reference::

( bam_fp, + reference_fasta, cc.chrom_tid, cc.start_pos, cc.end_pos, @@ -351,10 +407,46 @@ pub(crate) fn sample_reads_from_interval( allow_non_primary: bool, kmer_size: Option, ) -> anyhow::Result +where + P::Output: Moniod, +{ + sample_reads_from_interval_with_reference::

( + bam_fp, + None, + chrom_tid, + start, + end, + prev_end, + record_sampler, + collapse_method, + edge_filter, + position_filter, + only_mapped, + allow_non_primary, + kmer_size, + ) +} + +pub(crate) fn sample_reads_from_interval_with_reference( + bam_fp: &PathBuf, + reference_fasta: Option<&PathBuf>, + chrom_tid: u32, + start: u32, + end: u32, + prev_end: Option, + record_sampler: RecordSampler, + collapse_method: Option<&CollapseMethod>, + edge_filter: Option<&EdgeFilter>, + position_filter: Option<&StrandedPositionFilter<()>>, + only_mapped: bool, + allow_non_primary: bool, + kmer_size: Option, +) -> anyhow::Result where P::Output: Moniod, { let mut bam_reader = bam::IndexedReader::from_path(bam_fp)?; + set_reference_for_cram_indexed_reader(&mut bam_reader, reference_fasta)?; bam_reader.fetch(bam::FetchDefinition::Region( chrom_tid as i32, start as i64, diff --git a/modkit-core/src/reads_sampler/sampling_schedule.rs b/modkit-core/src/reads_sampler/sampling_schedule.rs index 21a99d29..72847a54 100644 --- a/modkit-core/src/reads_sampler/sampling_schedule.rs +++ b/modkit-core/src/reads_sampler/sampling_schedule.rs @@ -19,7 +19,10 @@ use crate::interval_chunks::{ChromCoordinates, MultiChromCoordinates}; use crate::monoid::Moniod; use crate::position_filter::StrandedPositionFilter; use crate::reads_sampler::record_sampler::RecordSampler; -use crate::util::{get_ticker, reader_is_bam, ReferenceRecord, Region}; +use crate::util::{ + get_ticker, reader_is_bam, set_reference_for_cram_indexed_reader, + ReferenceRecord, Region, +}; /// Count is an exact count, Sample is a fraction to sample #[derive(Debug, PartialEq, Copy, Clone)] @@ -174,8 +177,27 @@ impl SamplingSchedule { region: Option<&Region>, position_filter: Option<&StrandedPositionFilter<()>>, include_unmapped: bool, + ) -> anyhow::Result { + Self::from_num_reads_with_reference( + bam_fp, + None, + num_reads, + region, + position_filter, + include_unmapped, + ) + } + + pub(crate) fn from_num_reads_with_reference>( + bam_fp: T, + reference_fasta: Option<&PathBuf>, + num_reads: usize, + region: Option<&Region>, + position_filter: Option<&StrandedPositionFilter<()>>, + include_unmapped: bool, ) -> anyhow::Result { let mut reader = bam::IndexedReader::from_path(bam_fp)?; + set_reference_for_cram_indexed_reader(&mut reader, reference_fasta)?; let header = reader.header().to_owned(); let index_stats = IdxStats::new_from_reader(&mut reader, region, position_filter)?; @@ -318,17 +340,37 @@ impl SamplingSchedule { } } + #[allow(dead_code)] pub fn from_sample_frac>( bam_fp: T, sample_frac: f32, region: Option<&Region>, position_filter: Option<&StrandedPositionFilter<()>>, include_unmapped: bool, + ) -> anyhow::Result { + Self::from_sample_frac_with_reference( + bam_fp, + None, + sample_frac, + region, + position_filter, + include_unmapped, + ) + } + + pub(crate) fn from_sample_frac_with_reference>( + bam_fp: T, + reference_fasta: Option<&PathBuf>, + sample_frac: f32, + region: Option<&Region>, + position_filter: Option<&StrandedPositionFilter<()>>, + include_unmapped: bool, ) -> anyhow::Result { if sample_frac > 1.0 { bail!("sample fraction must be <= 1") } let mut reader = bam::IndexedReader::from_path(bam_fp)?; + set_reference_for_cram_indexed_reader(&mut reader, reference_fasta)?; let index_stats = IdxStats::new_from_reader(&mut reader, region, position_filter)?; drop(reader); @@ -627,22 +669,43 @@ pub(crate) struct IdxStats { } impl IdxStats { - pub(crate) fn check_any_mapped_reads( + pub(crate) fn check_any_mapped_reads_with_reference( bam_fp: &PathBuf, + reference_fasta: Option<&PathBuf>, region: Option<&Region>, position_filter: Option<&StrandedPositionFilter<()>>, ) -> anyhow::Result { - Self::new_from_path(bam_fp, region, position_filter) - .map(|idx_stats| idx_stats.mapped_read_count > 0) + Self::new_from_path_with_reference( + bam_fp, + reference_fasta, + region, + position_filter, + ) + .map(|idx_stats| idx_stats.mapped_read_count > 0) } pub(crate) fn new_from_path( bam_fp: &PathBuf, region: Option<&Region>, position_filter: Option<&StrandedPositionFilter<()>>, + ) -> anyhow::Result { + Self::new_from_path_with_reference( + bam_fp, + None, + region, + position_filter, + ) + } + + pub(crate) fn new_from_path_with_reference( + bam_fp: &PathBuf, + reference_fasta: Option<&PathBuf>, + region: Option<&Region>, + position_filter: Option<&StrandedPositionFilter<()>>, ) -> anyhow::Result { let mut reader = bam::IndexedReader::from_path(bam_fp) .context("could not create reader for getting mapping stats")?; + set_reference_for_cram_indexed_reader(&mut reader, reference_fasta)?; Self::new_from_reader(&mut reader, region, position_filter) } @@ -679,7 +742,7 @@ impl IdxStats { }) .transpose()?; - let is_bam = reader_is_bam(&reader); + let is_bam = reader_is_bam(reader); if is_bam { let idx_stats = reader.index_stats().context("failed to get index stats")?; @@ -842,7 +905,14 @@ impl ReferenceSequencesLookup { let mut tid_to_id = HashMap::new(); let idxs = bam_fps .iter() - .map(|fp| IdxStats::new_from_path(fp, None, None)) + .map(|fp| { + IdxStats::new_from_path_with_reference( + fp, + Some(reference_fasta_fp), + None, + None, + ) + }) .collect::>>()?; let reader = bam::IndexedReader::from_path(&bam_fps[0])?; let header = reader.header(); diff --git a/modkit-core/src/thresholds.rs b/modkit-core/src/thresholds.rs index 84cbfe0e..9beb0d8f 100644 --- a/modkit-core/src/thresholds.rs +++ b/modkit-core/src/thresholds.rs @@ -6,7 +6,7 @@ use crate::mod_bam::{CollapseMethod, EdgeFilter}; use crate::mod_base_code::{DnaBase, ModCodeRepr}; use crate::position_filter::StrandedPositionFilter; use crate::read_ids_to_base_mod_probs::ReadIdsToBaseModProbs; -use crate::reads_sampler::get_sampled_read_ids_to_base_mod_probs; +use crate::reads_sampler::get_sampled_read_ids_to_base_mod_probs_with_reference; use crate::threshold_mod_caller::MultipleThresholdModCaller; use crate::util::Region; use anyhow::{Context, Result as AnyhowResult}; @@ -122,20 +122,56 @@ pub fn calc_threshold_from_bam( only_mapped: bool, suppress_progress: bool, ) -> AnyhowResult> { - let (can_base_probs, explicit_can_probs) = get_modbase_probs_from_bam( + calc_threshold_from_bam_with_reference( bam_fp, + None, threads, interval_size, sample_frac, num_reads, + filter_percentile, seed, region, - collapse_method, edge_filter, + collapse_method, position_filter, only_mapped, suppress_progress, - )?; + ) +} + +pub(crate) fn calc_threshold_from_bam_with_reference( + bam_fp: &PathBuf, + reference_fasta: Option<&PathBuf>, + threads: usize, + interval_size: u32, + sample_frac: Option, + num_reads: Option, + filter_percentile: f32, + seed: Option, + region: Option<&Region>, + edge_filter: Option<&EdgeFilter>, + collapse_method: Option<&CollapseMethod>, + position_filter: Option<&StrandedPositionFilter<()>>, + only_mapped: bool, + suppress_progress: bool, +) -> AnyhowResult> { + let (can_base_probs, explicit_can_probs) = + get_modbase_probs_from_bam_with_reference( + bam_fp, + reference_fasta, + threads, + interval_size, + sample_frac, + num_reads, + seed, + region, + collapse_method, + edge_filter, + position_filter, + only_mapped, + suppress_progress, + )?; calculate_threshold_with_fallback( can_base_probs, &explicit_can_probs, @@ -157,8 +193,43 @@ pub fn get_modbase_probs_from_bam( only_mapped: bool, suppress_progress: bool, ) -> AnyhowResult<(HashMap>, HashMap)> { - get_sampled_read_ids_to_base_mod_probs::( + get_modbase_probs_from_bam_with_reference( + bam_fp, + None, + threads, + interval_size, + sample_frac, + num_reads, + seed, + region, + collapse_method, + edge_filter, + position_filter, + only_mapped, + suppress_progress, + ) +} + +pub(crate) fn get_modbase_probs_from_bam_with_reference( + bam_fp: &PathBuf, + reference_fasta: Option<&PathBuf>, + threads: usize, + interval_size: u32, + sample_frac: Option, + num_reads: Option, + seed: Option, + region: Option<&Region>, + collapse_method: Option<&CollapseMethod>, + edge_filter: Option<&EdgeFilter>, + position_filter: Option<&StrandedPositionFilter<()>>, + only_mapped: bool, + suppress_progress: bool, +) -> AnyhowResult<(HashMap>, HashMap)> { + get_sampled_read_ids_to_base_mod_probs_with_reference::< + ReadIdsToBaseModProbs, + >( bam_fp, + reference_fasta, threads, interval_size, sample_frac, diff --git a/modkit-core/src/util.rs b/modkit-core/src/util.rs index 7323a1a7..925de514 100644 --- a/modkit-core/src/util.rs +++ b/modkit-core/src/util.rs @@ -26,7 +26,7 @@ use prettytable::row; use regex::Regex; use rust_htslib::bam::{ self, ext::BamRecordExtensions, header::HeaderRecord, record::Aux, - HeaderView, Read, + HeaderView, Read as BamRead, }; use rustc_hash::FxHashMap; use substring::Substring; @@ -852,20 +852,71 @@ pub fn get_reference_mod_strand( } #[inline] -pub(crate) fn reader_is_bam(reader: &bam::IndexedReader) -> bool { +pub(crate) fn reader_is_bam(reader: &R) -> bool { unsafe { (*reader.htsfile()).format.format == rust_htslib::htslib::htsExactFormat_bam } } #[inline] -pub(crate) fn reader_is_cram(reader: &bam::IndexedReader) -> bool { +pub(crate) fn reader_is_cram(reader: &R) -> bool { unsafe { (*reader.htsfile()).format.format == rust_htslib::htslib::htsExactFormat_cram } } +pub(crate) fn set_reference_for_cram_reader( + reader: &mut bam::Reader, + reference_fasta: Option<&PathBuf>, +) -> anyhow::Result<()> { + if reader_is_cram(reader) { + if let Some(reference_fasta) = reference_fasta { + reader + .set_reference(reference_fasta) + .context("failed to set CRAM reference")?; + } + } + Ok(()) +} + +pub(crate) fn set_reference_for_cram_indexed_reader( + reader: &mut bam::IndexedReader, + reference_fasta: Option<&PathBuf>, +) -> anyhow::Result<()> { + if reader_is_cram(reader) { + if let Some(reference_fasta) = reference_fasta { + reader + .set_reference(reference_fasta) + .context("failed to set CRAM reference")?; + } + } + Ok(()) +} + +/// Check that the first record of a CRAM file can be decoded before creating +/// an output file. When no explicit reference is supplied, htslib may still +/// resolve an embedded, cached, or remotely available reference. +pub(crate) fn preflight_cram_input( + input: &Path, + reference_fasta: Option<&PathBuf>, +) -> anyhow::Result<()> { + let mut reader = bam::Reader::from_path(input) + .with_context(|| format!("failed to open input {input:?}"))?; + if !reader_is_cram(&reader) { + return Ok(()); + } + set_reference_for_cram_reader(&mut reader, reference_fasta)?; + let mut record = bam::Record::new(); + if let Some(result) = reader.read(&mut record) { + result.with_context(|| { + "failed to decode CRAM input; provide --reference when the CRAM \ + depends on an external reference" + })?; + } + Ok(()) +} + pub(crate) const KMER_SIZE: usize = 50; #[derive(Copy, Clone)] diff --git a/modkit/tests/test_cram_reference_consumers.rs b/modkit/tests/test_cram_reference_consumers.rs new file mode 100644 index 00000000..e03dc7ef --- /dev/null +++ b/modkit/tests/test_cram_reference_consumers.rs @@ -0,0 +1,373 @@ +use std::fs; +use std::path::Path; +use std::process::{Command, Output}; + +const BAM: &str = "../tests/resources/bc_anchored_10_reads.sorted.bam"; +const CRAM: &str = "../tests/resources/bc_anchored_10_reads.sorted.cram"; +const REFERENCE: &str = "../tests/resources/CGI_ladder_3.6kb_ref.fa"; +const REFERENCE_INDEPENDENT_CRAM: &str = + "../tests/resources/bc_anchored_10_reads_unmapped.cram"; + +fn run_modkit(args: &[&str]) -> Output { + let exe = Path::new(env!("CARGO_BIN_EXE_modkit")); + assert!(exe.exists()); + Command::new(exe).args(args).output().unwrap() +} + +fn run_modkit_without_reference_resolution( + args: &[&str], + empty_reference_cache: &Path, +) -> Output { + let exe = Path::new(env!("CARGO_BIN_EXE_modkit")); + assert!(exe.exists()); + let reference_pattern = empty_reference_cache.join("%2s/%2s/%s"); + Command::new(exe) + .args(args) + .env("REF_PATH", &reference_pattern) + .env("REF_CACHE", &reference_pattern) + .output() + .unwrap() +} + +fn assert_success(output: &Output) { + assert!( + output.status.success(), + "command failed:\n{}", + String::from_utf8_lossy(&output.stderr) + ); +} + +fn normalized_sam_records(path: &Path) -> Vec { + fs::read_to_string(path) + .unwrap() + .lines() + .filter(|line| !line.starts_with('@')) + .map(|line| { + let mut fields = line.split('\t').collect::>(); + let mut aux = fields.split_off(11); + aux.sort_unstable(); + fields.extend(aux); + fields.join("\t") + }) + .collect() +} + +fn threshold_messages(output: &Output) -> Vec<&str> { + std::str::from_utf8(&output.stderr) + .unwrap() + .lines() + .filter(|line| { + line.contains("Threshold of") + || line.contains("calculated thresholds:") + }) + .collect() +} + +fn assert_entropy_parity(bam_output: &Path, cram_output: &Path) { + let bam = fs::read_to_string(bam_output).unwrap(); + let cram = fs::read_to_string(cram_output).unwrap(); + let bam_lines = bam.lines().collect::>(); + let cram_lines = cram.lines().collect::>(); + assert_eq!(bam_lines.len(), cram_lines.len()); + for (bam_line, cram_line) in bam_lines.iter().zip(cram_lines) { + let bam_fields = bam_line.split('\t').collect::>(); + let cram_fields = cram_line.split('\t').collect::>(); + assert_eq!(bam_fields.len(), 6); + assert_eq!(cram_fields.len(), 6); + assert_eq!(&bam_fields[..3], &cram_fields[..3]); + assert_eq!(&bam_fields[4..], &cram_fields[4..]); + let bam_entropy = bam_fields[3].parse::().unwrap(); + let cram_entropy = cram_fields[3].parse::().unwrap(); + assert!( + (bam_entropy - cram_entropy).abs() <= 1e-6, + "entropy differed: {bam_entropy} versus {cram_entropy}" + ); + } +} + +#[test] +fn call_and_adjust_mods_match_bam_and_cram() { + let temp_dir = tempfile::tempdir().unwrap(); + for command in ["call-mods", "adjust-mods"] { + let help = run_modkit(&[command, "--help"]); + assert_success(&help); + assert!(String::from_utf8_lossy(&help.stdout) + .contains("-r, --reference ")); + } + + let call_bam = temp_dir.path().join("call.bam.sam"); + let call_cram = temp_dir.path().join("call.cram.sam"); + + let bam_output = run_modkit(&[ + "call-mods", + BAM, + call_bam.to_str().unwrap(), + "--output-sam", + "--sampling-frac", + "1", + "--seed", + "7", + "--threads", + "1", + "--sampling-interval-size", + "20", + "--reference", + REFERENCE, + "--suppress-progress", + ]); + let cram_output = run_modkit(&[ + "call-mods", + CRAM, + call_cram.to_str().unwrap(), + "--output-sam", + "--sampling-frac", + "1", + "--seed", + "7", + "--threads", + "1", + "--sampling-interval-size", + "20", + "--ref", + REFERENCE, + "--suppress-progress", + ]); + assert_success(&bam_output); + assert_success(&cram_output); + let bam_records = normalized_sam_records(&call_bam); + assert_eq!(bam_records.len(), 10); + assert_eq!(bam_records, normalized_sam_records(&call_cram)); + + let adjust_bam = temp_dir.path().join("adjust.bam.sam"); + let adjust_cram = temp_dir.path().join("adjust.cram.sam"); + let bam_output = run_modkit(&[ + "adjust-mods", + BAM, + adjust_bam.to_str().unwrap(), + "--output-sam", + "--filter-probs", + "--filter-threshold", + "0.55", + "--threads", + "1", + "-r", + REFERENCE, + "--suppress-progress", + ]); + let cram_output = run_modkit(&[ + "adjust-mods", + CRAM, + adjust_cram.to_str().unwrap(), + "--output-sam", + "--filter-probs", + "--filter-threshold", + "0.55", + "--threads", + "1", + "--reference", + REFERENCE, + "--suppress-progress", + ]); + assert_success(&bam_output); + assert_success(&cram_output); + let bam_records = normalized_sam_records(&adjust_bam); + assert_eq!(bam_records.len(), 10); + assert_eq!(bam_records, normalized_sam_records(&adjust_cram)); + + let sampled_adjust_cram = temp_dir.path().join("adjust.sampled.cram.sam"); + let output = run_modkit(&[ + "adjust-mods", + CRAM, + sampled_adjust_cram.to_str().unwrap(), + "--output-sam", + "--filter-probs", + "--num-reads", + "10042", + "--threads", + "1", + "--sampling-interval-size", + "20", + "--ref", + REFERENCE, + "--suppress-progress", + ]); + assert_success(&output); + assert_eq!(normalized_sam_records(&sampled_adjust_cram).len(), 10); +} + +#[test] +fn pileup_hemi_and_entropy_match_bam_and_cram() { + let temp_dir = tempfile::tempdir().unwrap(); + let hemi_bam = temp_dir.path().join("hemi.bam.bed"); + let hemi_cram = temp_dir.path().join("hemi.cram.bed"); + let include_bed = "../tests/resources/include-pos-1-site.bed"; + + let run_hemi = |input: &str, output: &Path| { + run_modkit(&[ + "pileup-hemi", + input, + "--out-bed", + output.to_str().unwrap(), + "--ref", + REFERENCE, + "--cpg", + "--region", + "oligo_1512_adapters", + "--include-bed", + include_bed, + "--sampling-frac", + "1", + "--seed", + "7", + "--threads", + "1", + "--interval-size", + "20", + "--sampling-interval-size", + "20", + "--suppress-progress", + ]) + }; + let bam_output = run_hemi(BAM, &hemi_bam); + let cram_output = run_hemi(CRAM, &hemi_cram); + assert_success(&bam_output); + assert_success(&cram_output); + let bam_thresholds = threshold_messages(&bam_output); + let cram_thresholds = threshold_messages(&cram_output); + assert!(!bam_thresholds.is_empty()); + assert_eq!(bam_thresholds, cram_thresholds); + assert!(String::from_utf8_lossy(&bam_output.stderr) + .contains("Processed ~10 reads")); + assert!(String::from_utf8_lossy(&cram_output.stderr) + .contains("Processed ~10 reads")); + assert_eq!(fs::read(&hemi_bam).unwrap(), fs::read(&hemi_cram).unwrap()); + + let entropy_bam = temp_dir.path().join("entropy.bam.bed"); + let entropy_cram = temp_dir.path().join("entropy.cram.bed"); + let run_entropy = |input: &str, output: &Path| { + run_modkit(&[ + "entropy", + "--in-bam", + input, + "--out-bed", + output.to_str().unwrap(), + "--ref", + REFERENCE, + "--cpg", + "--min-coverage", + "1", + "--num-reads", + "10042", + "--threads", + "1", + "--io-threads", + "1", + "--suppress-progress", + "--force", + ]) + }; + let bam_output = run_entropy(BAM, &entropy_bam); + let cram_output = run_entropy(CRAM, &entropy_cram); + assert_success(&bam_output); + assert_success(&cram_output); + let bam_thresholds = threshold_messages(&bam_output); + let cram_thresholds = threshold_messages(&cram_output); + assert!(!bam_thresholds.is_empty()); + assert_eq!(bam_thresholds, cram_thresholds); + assert!(String::from_utf8_lossy(&bam_output.stderr) + .contains("6 windows processed successfully")); + assert!(String::from_utf8_lossy(&cram_output.stderr) + .contains("6 windows processed successfully")); + assert_entropy_parity(&entropy_bam, &entropy_cram); +} + +#[test] +fn missing_external_reference_fails_before_output_without_rejecting_all_cram() { + let temp_dir = tempfile::tempdir().unwrap(); + let cache_dir = temp_dir.path().join("empty-reference-cache"); + fs::create_dir(&cache_dir).unwrap(); + let disguised_cram = temp_dir.path().join("reference-dependent.bam"); + fs::copy(CRAM, &disguised_cram).unwrap(); + let disguised_cram = disguised_cram.to_str().unwrap(); + + for (command, transform_args) in [ + ("call-mods", vec!["--filter-threshold", "0.5"]), + ("adjust-mods", vec!["--ignore", "m"]), + ] { + for preexisting_output in [false, true] { + let output_path = temp_dir.path().join(format!( + "{command}-{}.sam", + if preexisting_output { "existing" } else { "absent" } + )); + if preexisting_output { + fs::write(&output_path, b"sentinel\n").unwrap(); + } + let mut args = vec![ + command, + disguised_cram, + output_path.to_str().unwrap(), + "--output-sam", + "--threads", + "1", + "--suppress-progress", + ]; + args.extend(transform_args.iter().copied()); + let output = + run_modkit_without_reference_resolution(&args, &cache_dir); + assert!(!output.status.success()); + assert!(String::from_utf8_lossy(&output.stderr) + .contains("failed to decode CRAM input")); + if preexisting_output { + assert_eq!(fs::read(&output_path).unwrap(), b"sentinel\n"); + } else { + assert!(!output_path.exists()); + } + } + } + + for (command, transform_args) in [ + ("call-mods", vec!["--filter-threshold", "0.5"]), + ("adjust-mods", vec!["--ignore", "m"]), + ] { + let output_path = temp_dir + .path() + .join(format!("reference-independent-{command}.sam")); + let mut args = vec![ + command, + REFERENCE_INDEPENDENT_CRAM, + output_path.to_str().unwrap(), + "--output-sam", + "--threads", + "1", + "--suppress-progress", + ]; + args.extend(transform_args); + let output = run_modkit_without_reference_resolution(&args, &cache_dir); + assert_success(&output); + assert!(output_path.exists()); + assert_eq!(normalized_sam_records(&output_path).len(), 10); + } + + let hemi_output = temp_dir.path().join("missing-ref-hemi.bed"); + let output = run_modkit(&[ + "pileup-hemi", + CRAM, + "--out-bed", + hemi_output.to_str().unwrap(), + "--cpg", + ]); + assert!(!output.status.success()); + assert!(!hemi_output.exists()); + + let entropy_output = temp_dir.path().join("missing-ref-entropy.bed"); + let output = run_modkit(&[ + "entropy", + "--in-bam", + CRAM, + "--out-bed", + entropy_output.to_str().unwrap(), + "--cpg", + ]); + assert!(!output.status.success()); + assert!(!entropy_output.exists()); +} From ac084130c12dbba96a2f3f271ddd935c3b4bfddc Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 09:57:34 -0700 Subject: [PATCH 02/10] Restore duplex batch compatibility API --- modkit-core/src/pileup/duplex.rs | 56 ++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/modkit-core/src/pileup/duplex.rs b/modkit-core/src/pileup/duplex.rs index 3eb40fc6..af1e2d93 100644 --- a/modkit-core/src/pileup/duplex.rs +++ b/modkit-core/src/pileup/duplex.rs @@ -208,6 +208,30 @@ impl DuplexFeatureVector { // todo this function should be removed in favor of a more // generic version in pileup/mod.rs +#[allow(dead_code)] +pub fn process_region_duplex_batch + Copy>( + chromosome_coordintes: &MultiChromCoordinates, + bam_fp: T, + caller: &MultipleThresholdModCaller, + pileup_numeric_options: &PileupNumericOptions, + force_allow: bool, + max_depth: u32, + motif: MotifInfo, + edge_filter: Option<&EdgeFilter>, +) -> Vec> { + process_region_duplex_batch_with_reference( + chromosome_coordintes, + bam_fp, + None, + caller, + pileup_numeric_options, + force_allow, + max_depth, + motif, + edge_filter, + ) +} + pub(crate) fn process_region_duplex_batch_with_reference< T: AsRef + Copy, >( @@ -350,3 +374,35 @@ fn process_region_duplex>( skipped_records, }) } + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn legacy_duplex_batch_signature_compiles() { + fn call_legacy_signature( + chromosome_coordinates: &MultiChromCoordinates, + bam_fp: &Path, + caller: &MultipleThresholdModCaller, + pileup_numeric_options: &PileupNumericOptions, + force_allow: bool, + max_depth: u32, + motif: MotifInfo, + edge_filter: Option<&EdgeFilter>, + ) -> Vec> { + process_region_duplex_batch( + chromosome_coordinates, + bam_fp, + caller, + pileup_numeric_options, + force_allow, + max_depth, + motif, + edge_filter, + ) + } + + let _ = call_legacy_signature; + } +} From 0d828f1efb7ee862cd5cbb7e3895c0444c6f99e9 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 10:52:32 -0700 Subject: [PATCH 03/10] Honor entropy no-filtering mode --- modkit-core/src/entropy/subcommand.rs | 12 ++- modkit/tests/test_entropy.rs | 121 ++++++++++++++++++++++++++ 2 files changed, 132 insertions(+), 1 deletion(-) diff --git a/modkit-core/src/entropy/subcommand.rs b/modkit-core/src/entropy/subcommand.rs index c16a35c7..5adfd36b 100644 --- a/modkit-core/src/entropy/subcommand.rs +++ b/modkit-core/src/entropy/subcommand.rs @@ -60,7 +60,12 @@ pub struct MethylationEntropy { window_size: usize, /// Do not perform any filtering, include all mod base calls in output. #[clap(help_heading = "Filtering Options")] - #[arg(group = "thresholds", long, default_value_t = false)] + #[arg( + group = "thresholds", + long, + conflicts_with = "mod_thresholds", + default_value_t = false + )] no_filtering: bool, /// Sample this many reads when estimating the filtering threshold. Reads /// will be sampled evenly across aligned genome. If a region is @@ -462,6 +467,11 @@ impl MethylationEntropy { &self, pool: &rayon::ThreadPool, ) -> anyhow::Result { + if self.no_filtering { + info!("not performing filtering"); + return Ok(MultipleThresholdModCaller::new_passthrough()); + } + let per_mod_thresholds = self .mod_thresholds .as_ref() diff --git a/modkit/tests/test_entropy.rs b/modkit/tests/test_entropy.rs index cf67ab32..190460c4 100644 --- a/modkit/tests/test_entropy.rs +++ b/modkit/tests/test_entropy.rs @@ -1,3 +1,5 @@ +use std::{fs, path::Path, process::Command}; + use crate::common::run_modkit; mod common; @@ -39,3 +41,122 @@ fn test_entropy_regression() { // check_against_expected_text_file(windows.to_str().unwrap(), // "../tests/resources/expected_entropy_windows.bed"); } + +#[test] +fn no_filtering_bypasses_sampling_and_matches_zero_threshold() { + let temp_dir = tempfile::tempdir().unwrap(); + let no_filter_output = temp_dir.path().join("no_filter.bed"); + let zero_threshold_output = temp_dir.path().join("zero_threshold.bed"); + let executable = Path::new(env!("CARGO_BIN_EXE_modkit")); + + let run = |output: &Path, threshold_args: &[&str]| { + Command::new(executable) + .args([ + "entropy", + "--in-bam", + "../tests/resources/bc_anchored_10_reads.sorted.bam", + "--out-bed", + output.to_str().unwrap(), + "--ref", + "../tests/resources/CGI_ladder_3.6kb_ref.fa", + "--cpg", + "--min-coverage", + "1", + "--threads", + "1", + "--io-threads", + "1", + "--num-reads", + "0", + "--suppress-progress", + ]) + .args(threshold_args) + .output() + .unwrap() + }; + + let no_filter = run(&no_filter_output, &["--no-filtering"]); + assert!( + no_filter.status.success(), + "no-filtering failed:\n{}", + String::from_utf8_lossy(&no_filter.stderr) + ); + let no_filter_stderr = String::from_utf8_lossy(&no_filter.stderr); + assert!(no_filter_stderr.contains("not performing filtering")); + assert!(!no_filter_stderr.contains("calculated thresholds:")); + + let zero_threshold = + run(&zero_threshold_output, &["--filter-threshold", "0"]); + assert!( + zero_threshold.status.success(), + "zero threshold failed:\n{}", + String::from_utf8_lossy(&zero_threshold.stderr) + ); + + let no_filter_text = fs::read_to_string(no_filter_output).unwrap(); + let zero_threshold_text = + fs::read_to_string(zero_threshold_output).unwrap(); + let no_filter_rows = no_filter_text.lines().collect::>(); + let zero_threshold_rows = zero_threshold_text.lines().collect::>(); + assert_eq!(no_filter_rows.len(), zero_threshold_rows.len()); + for (no_filter_row, zero_threshold_row) in + no_filter_rows.into_iter().zip(zero_threshold_rows) + { + let no_filter_fields = no_filter_row.split('\t').collect::>(); + let zero_threshold_fields = + zero_threshold_row.split('\t').collect::>(); + assert_eq!(&no_filter_fields[..3], &zero_threshold_fields[..3]); + assert_eq!(&no_filter_fields[4..], &zero_threshold_fields[4..]); + let no_filter_entropy = no_filter_fields[3].parse::().unwrap(); + let zero_threshold_entropy = + zero_threshold_fields[3].parse::().unwrap(); + assert!( + (no_filter_entropy - zero_threshold_entropy).abs() <= 1e-6, + "entropy differed: {no_filter_entropy} versus \ + {zero_threshold_entropy}" + ); + } +} + +#[test] +fn no_filtering_conflicts_with_per_mod_thresholds_before_output() { + let temp_dir = tempfile::tempdir().unwrap(); + let absent_output = temp_dir.path().join("absent.bed"); + let sentinel_output = temp_dir.path().join("sentinel.bed"); + fs::write(&sentinel_output, "keep-sentinel\n").unwrap(); + let executable = Path::new(env!("CARGO_BIN_EXE_modkit")); + + for (threshold_option, threshold_value, output) in [ + ("--mod-threshold", "m:0.5", &absent_output), + ("--mod-thresholds", "not-valid", &sentinel_output), + ] { + let result = Command::new(executable) + .args([ + "entropy", + "--in-bam", + "../tests/resources/bc_anchored_10_reads.sorted.bam", + "--out-bed", + output.to_str().unwrap(), + "--ref", + "../tests/resources/CGI_ladder_3.6kb_ref.fa", + "--cpg", + "--no-filtering", + threshold_option, + threshold_value, + "--suppress-progress", + ]) + .output() + .unwrap(); + + assert_eq!(result.status.code(), Some(2)); + assert!( + String::from_utf8_lossy(&result.stderr) + .contains("cannot be used with"), + "unexpected stderr:\n{}", + String::from_utf8_lossy(&result.stderr) + ); + } + + assert!(!absent_output.exists()); + assert_eq!(fs::read_to_string(sentinel_output).unwrap(), "keep-sentinel\n"); +} From ed9a6a205b82f2157952ba008f973683dcf7b2ea Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 09:53:55 -0700 Subject: [PATCH 04/10] Reject empty automatic threshold samples --- modkit-core/src/thresholds.rs | 21 +++++++++ modkit/tests/test_empty_threshold_sampling.rs | 43 +++++++++++++++++++ 2 files changed, 64 insertions(+) create mode 100644 modkit/tests/test_empty_threshold_sampling.rs diff --git a/modkit-core/src/thresholds.rs b/modkit-core/src/thresholds.rs index 9beb0d8f..b875ce49 100644 --- a/modkit-core/src/thresholds.rs +++ b/modkit-core/src/thresholds.rs @@ -250,6 +250,11 @@ pub(crate) fn calculate_threshold_with_fallback( max_explicit_canonical_prob: &HashMap, filter_percentile: f32, ) -> anyhow::Result> { + anyhow::ensure!( + !probs_per_base.is_empty(), + "cannot calculate automatic thresholds because no modification \ + probabilities were sampled" + ); probs_per_base .iter_mut() .map(|(dna_base, mod_base_probs)| { @@ -328,4 +333,20 @@ mod thresolds_tests { assert_eq!(t, 0.95); } } + + #[test] + fn empty_sample_does_not_produce_an_empty_threshold_map() { + let error = calculate_threshold_with_fallback( + Default::default(), + &Default::default(), + 0.1, + ) + .expect_err("an empty sample must not disable filtering"); + + assert_eq!( + error.to_string(), + "cannot calculate automatic thresholds because no modification \ + probabilities were sampled" + ); + } } diff --git a/modkit/tests/test_empty_threshold_sampling.rs b/modkit/tests/test_empty_threshold_sampling.rs new file mode 100644 index 00000000..e0df9851 --- /dev/null +++ b/modkit/tests/test_empty_threshold_sampling.rs @@ -0,0 +1,43 @@ +use std::path::Path; +use std::process::Command; + +#[test] +fn pileup_hemi_rejects_an_empty_automatic_threshold_sample_before_output() { + let temp_dir = tempfile::tempdir().unwrap(); + let output_path = temp_dir.path().join("pileup.bed"); + let executable = Path::new(env!("CARGO_BIN_EXE_modkit")); + + let output = Command::new(executable) + .args([ + "pileup-hemi", + "../tests/resources/duplex_modcalls_sort.bam", + "--out-bed", + output_path.to_str().unwrap(), + "--ref", + "../tests/resources/GRCh38_chr20.fa", + "--motif", + "CG", + "0", + "--region", + "chr20:22,613,835-22,640,468", + "--sampling-frac", + "0", + "--suppress-progress", + ]) + .output() + .expect("failed to run modkit pileup-hemi"); + + assert!(!output.status.success(), "an empty sample must fail"); + assert!( + String::from_utf8_lossy(&output.stderr).contains( + "cannot calculate automatic thresholds because no modification \ + probabilities were sampled" + ), + "unexpected stderr: {}", + String::from_utf8_lossy(&output.stderr) + ); + assert!( + !output_path.exists(), + "failure must occur before creating the scientific output" + ); +} From 06cf545ecc78f8c4b5598517b5f11431222a3bd8 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 11:03:27 -0700 Subject: [PATCH 05/10] Protect entropy outputs until setup completes --- modkit-core/Cargo.toml | 1 + modkit-core/src/entropy/subcommand.rs | 60 ++-- modkit-core/src/entropy/writers.rs | 128 ++++++-- modkit/tests/test_entropy.rs | 410 ++++++++++++++++++++++++++ 4 files changed, 555 insertions(+), 44 deletions(-) diff --git a/modkit-core/Cargo.toml b/modkit-core/Cargo.toml index 881c7d87..57c1e22e 100644 --- a/modkit-core/Cargo.toml +++ b/modkit-core/Cargo.toml @@ -29,6 +29,7 @@ indexmap = "2.2.6" indicatif = { version = "0.17.1", features = ["rayon"] } itertools = "0.12.1" lazy_static = "1.4" +libc = "0.2" linear-map = "1.2.0" log = "0.4.0" log-once = "0.4.0" diff --git a/modkit-core/src/entropy/subcommand.rs b/modkit-core/src/entropy/subcommand.rs index 5adfd36b..3a7f5cbc 100644 --- a/modkit-core/src/entropy/subcommand.rs +++ b/modkit-core/src/entropy/subcommand.rs @@ -213,33 +213,6 @@ impl MethylationEntropy { })?; } - let mut writer: Box = - match (self.out_bed.as_ref(), self.regions_fp.is_some()) { - (Some(out_fp), false) => Box::new( - WindowsWriter::new_file(out_fp, self.header, self.verbose) - .context("failed to make writer to file")?, - ), - (Some(out_dir), true) => Box::new( - RegionsWriter::new( - out_dir, - self.prefix.as_ref(), - self.header, - self.verbose, - ) - .context( - "failed to make regions writer, output must be a \ - directory", - )?, - ), - (None, false) => Box::new( - WindowsWriter::new_stdout(self.header, self.verbose) - .context("failed to make writer to stdout")?, - ), - (None, true) => { - bail!("must provide output directory with regions") - } - }; - let pool = rayon::ThreadPoolBuilder::new() .num_threads(self.threads) .build()?; @@ -348,6 +321,39 @@ impl MethylationEntropy { let threshold_caller = self.get_threshold_caller(&pool).map(|c| Arc::new(c))?; + let mut writer: Box = + match (self.out_bed.as_ref(), self.regions_fp.is_some()) { + (Some(out_fp), false) => Box::new( + WindowsWriter::new_file( + out_fp, + self.header, + self.verbose, + self.force, + ) + .context("failed to make writer to file")?, + ), + (Some(out_dir), true) => Box::new( + RegionsWriter::new( + out_dir, + self.prefix.as_ref(), + self.header, + self.verbose, + self.force, + ) + .context( + "failed to make regions writer, output must be a \ + directory", + )?, + ), + (None, false) => Box::new( + WindowsWriter::new_stdout(self.header, self.verbose) + .context("failed to make writer to stdout")?, + ), + (None, true) => { + bail!("must provide output directory with regions") + } + }; + let (snd, rcv) = crossbeam::channel::bounded(10_000); let bam_fps = self.in_bams.clone(); diff --git a/modkit-core/src/entropy/writers.rs b/modkit-core/src/entropy/writers.rs index 586a563b..76746794 100644 --- a/modkit-core/src/entropy/writers.rs +++ b/modkit-core/src/entropy/writers.rs @@ -6,10 +6,10 @@ use indicatif::ProgressBar; use log::debug; use rustc_hash::FxHashMap; use std::collections::HashMap; -use std::fs::File; -use std::io::{stdout, BufWriter, Write}; +use std::fs::{File, OpenOptions}; +use std::io::{stdout, BufWriter, ErrorKind, Write}; use std::ops::AddAssign; -use std::path::PathBuf; +use std::path::{Component, Path, PathBuf}; #[inline(always)] fn write_entropy_windows( @@ -142,6 +142,72 @@ pub(super) trait EntropyWriter { const WINDOWS_HEADER: &'static str = "\ #chrom\tstart\tend\tentropy\tstrand\tnum_reads\n"; +fn preflight_output_file(out_fp: &Path) -> anyhow::Result { + match std::fs::symlink_metadata(out_fp) { + Ok(metadata) if metadata.file_type().is_file() => Ok(true), + Ok(_) => bail!( + "entropy output target must be a regular file and may not be a \ + symbolic link: {}", + out_fp.display() + ), + Err(error) if error.kind() == ErrorKind::NotFound => Ok(false), + Err(error) => Err(error.into()), + } +} + +fn open_output_file_untruncated( + out_fp: &Path, + existed: bool, +) -> std::io::Result { + let mut options = OpenOptions::new(); + options.write(true); + if !existed { + options.create_new(true); + } + #[cfg(unix)] + { + use std::os::unix::fs::OpenOptionsExt; + options.custom_flags(libc::O_NOFOLLOW); + } + options.open(out_fp) +} + +fn open_output_file(out_fp: &Path, force: bool) -> anyhow::Result { + let existed = preflight_output_file(out_fp)?; + if existed && !force { + bail!( + "entropy output already exists, use --force to overwrite it: {}", + out_fp.display() + ) + } + let file = open_output_file_untruncated(out_fp, existed)?; + if force { + file.set_len(0)?; + } + Ok(file) +} + +fn validate_regions_prefix( + prefix: Option<&String>, +) -> anyhow::Result> { + let Some(prefix) = prefix else { + return Ok(None); + }; + let path = Path::new(prefix); + let mut components = path.components(); + match (components.next(), components.next()) { + (Some(Component::Normal(component)), None) + if !component.is_empty() && path.as_os_str() == component => + { + Ok(Some(prefix)) + } + _ => bail!( + "entropy regions prefix must be exactly one non-empty filename \ + component" + ), + } +} + pub(super) struct WindowsWriter { output: BufWriter, verbose: bool, @@ -152,8 +218,9 @@ impl WindowsWriter { out_fp: &PathBuf, header: bool, verbose: bool, + force: bool, ) -> anyhow::Result { - let mut output = BufWriter::new(File::create(out_fp)?); + let mut output = BufWriter::new(open_output_file(out_fp, force)?); if header { output.write(WINDOWS_HEADER.as_bytes())?; } @@ -186,27 +253,54 @@ impl RegionsWriter { prefix: Option<&String>, header: bool, verbose: bool, + force: bool, ) -> anyhow::Result { - if out_dir.is_file() { + let prefix = validate_regions_prefix(prefix)?; + if out_dir.exists() && !out_dir.is_dir() { bail!("regions output location must be a directory") } - std::fs::create_dir_all(out_dir)?; + if !out_dir.exists() { + std::fs::create_dir_all(out_dir)?; + } debug_assert!(out_dir.exists(), "out_dir should exist now"); - let mut regions_bed_out = if let Some(p) = prefix { - let fp = out_dir.join(format!("{p}_regions.bed")); - BufWriter::new(File::create(fp)?) + let regions_fp = if let Some(p) = prefix { + out_dir.join(format!("{p}_regions.bed")) } else { - let fp = out_dir.join("regions.bed"); - BufWriter::new(File::create(fp)?) + out_dir.join("regions.bed") }; - - let mut windows_bed_out = if let Some(p) = prefix { - let fp = out_dir.join(format!("{p}_windows.bedgraph")); - BufWriter::new(File::create(fp)?) + let windows_fp = if let Some(p) = prefix { + out_dir.join(format!("{p}_windows.bedgraph")) } else { - let fp = out_dir.join("windows.bedgraph"); - BufWriter::new(File::create(fp)?) + out_dir.join("windows.bedgraph") }; + let regions_existed = preflight_output_file(®ions_fp)?; + let windows_existed = preflight_output_file(&windows_fp)?; + if !force && (regions_existed || windows_existed) { + bail!( + "entropy region output already exists, use --force to \ + overwrite it" + ) + } + let regions_file = + open_output_file_untruncated(®ions_fp, regions_existed)?; + let windows_file = + match open_output_file_untruncated(&windows_fp, windows_existed) { + Ok(file) => file, + Err(error) => { + drop(regions_file); + if !regions_existed { + std::fs::remove_file(®ions_fp)?; + } + return Err(error.into()); + } + }; + if force { + regions_file.set_len(0)?; + windows_file.set_len(0)?; + } + + let mut regions_bed_out = BufWriter::new(regions_file); + let mut windows_bed_out = BufWriter::new(windows_file); if header { windows_bed_out.write(WINDOWS_HEADER.as_bytes())?; diff --git a/modkit/tests/test_entropy.rs b/modkit/tests/test_entropy.rs index 190460c4..e635623d 100644 --- a/modkit/tests/test_entropy.rs +++ b/modkit/tests/test_entropy.rs @@ -4,6 +4,26 @@ use crate::common::run_modkit; mod common; +fn entropy_command() -> Command { + let mut command = Command::new(env!("CARGO_BIN_EXE_modkit")); + command.args([ + "entropy", + "--in-bam", + "../tests/resources/bc_anchored_10_reads.sorted.bam", + "--ref", + "../tests/resources/CGI_ladder_3.6kb_ref.fa", + "--cpg", + "--min-coverage", + "1", + "--threads", + "1", + "--io-threads", + "1", + "--suppress-progress", + ]); + command +} + #[test] fn test_entropy_help() { run_modkit(&["entropy", "--help"]).expect("entropy help"); @@ -160,3 +180,393 @@ fn no_filtering_conflicts_with_per_mod_thresholds_before_output() { assert!(!absent_output.exists()); assert_eq!(fs::read_to_string(sentinel_output).unwrap(), "keep-sentinel\n"); } + +#[test] +fn threshold_failure_does_not_create_or_truncate_entropy_output() { + let temp_dir = tempfile::tempdir().unwrap(); + let absent_output = temp_dir.path().join("absent.bed"); + let sentinel_output = temp_dir.path().join("sentinel.bed"); + let sentinel = b"existing scientific output\n"; + fs::write(&sentinel_output, sentinel).unwrap(); + + for output_path in [&absent_output, &sentinel_output] { + let output = entropy_command() + .args([ + "--out-bed", + output_path.to_str().unwrap(), + "--num-reads", + "0", + "--header", + "--force", + ]) + .output() + .unwrap(); + assert!( + !output.status.success(), + "empty automatic-threshold setup unexpectedly succeeded" + ); + } + + assert!(!absent_output.exists()); + assert_eq!(fs::read(sentinel_output).unwrap(), sentinel); +} + +#[test] +fn entropy_file_output_requires_force() { + let temp_dir = tempfile::tempdir().unwrap(); + let output_path = temp_dir.path().join("entropy.bed"); + let sentinel = b"existing scientific output\n"; + fs::write(&output_path, sentinel).unwrap(); + + let rejected = entropy_command() + .args([ + "--out-bed", + output_path.to_str().unwrap(), + "--filter-threshold", + "0", + ]) + .output() + .unwrap(); + assert!(!rejected.status.success()); + assert_eq!(fs::read(&output_path).unwrap(), sentinel); + + let overwritten = entropy_command() + .args([ + "--out-bed", + output_path.to_str().unwrap(), + "--filter-threshold", + "0", + "--force", + ]) + .output() + .unwrap(); + assert!( + overwritten.status.success(), + "forced output failed:\n{}", + String::from_utf8_lossy(&overwritten.stderr) + ); + assert_ne!(fs::read(output_path).unwrap(), sentinel); +} + +#[test] +fn entropy_regions_preflight_both_prefixed_targets() { + let temp_dir = tempfile::tempdir().unwrap(); + let output_dir = temp_dir.path().join("entropy"); + fs::create_dir(&output_dir).unwrap(); + let unrelated = output_dir.join("keep.txt"); + let unrelated_contents = b"unrelated contents\n"; + fs::write(&unrelated, unrelated_contents).unwrap(); + + let regions_output = output_dir.join("sample_regions.bed"); + let windows_output = output_dir.join("sample_windows.bedgraph"); + let sentinel = b"existing windows output\n"; + fs::write(&windows_output, sentinel).unwrap(); + + let run_regions = |prefix: &str, force: bool| { + let mut command = entropy_command(); + command.args([ + "--out-bed", + output_dir.to_str().unwrap(), + "--regions", + "../tests/resources/entropy_test_regions.bed", + "--prefix", + prefix, + "--filter-threshold", + "0", + ]); + if force { + command.arg("--force"); + } + command.output().unwrap() + }; + + let rejected = run_regions("sample", false); + assert!(!rejected.status.success()); + assert!(!regions_output.exists()); + assert_eq!(fs::read(&windows_output).unwrap(), sentinel); + assert_eq!(fs::read(&unrelated).unwrap(), unrelated_contents); + + let fresh = run_regions("fresh", false); + assert!( + fresh.status.success(), + "fresh prefixed output failed:\n{}", + String::from_utf8_lossy(&fresh.stderr) + ); + assert!(output_dir.join("fresh_regions.bed").exists()); + assert!(output_dir.join("fresh_windows.bedgraph").exists()); + assert_eq!(fs::read(&unrelated).unwrap(), unrelated_contents); + + let forced = run_regions("sample", true); + assert!( + forced.status.success(), + "forced prefixed output failed:\n{}", + String::from_utf8_lossy(&forced.stderr) + ); + assert!(regions_output.exists()); + assert_ne!(fs::read(windows_output).unwrap(), sentinel); + assert_eq!(fs::read(unrelated).unwrap(), unrelated_contents); +} + +#[test] +fn entropy_regions_rejects_escaping_prefixes_before_forced_mutation() { + let temp_dir = tempfile::tempdir().unwrap(); + let output_dir = temp_dir.path().join("entropy"); + fs::create_dir(&output_dir).unwrap(); + fs::create_dir(output_dir.join("a")).unwrap(); + let unrelated = output_dir.join("keep.txt"); + let unrelated_contents = b"unrelated contents\n"; + fs::write(&unrelated, unrelated_contents).unwrap(); + + let absolute_prefix = temp_dir.path().join("absolute-x"); + let absolute_prefix_string = absolute_prefix.to_str().unwrap().to_string(); + let sentinel = b"outside scientific output\n"; + let invalid_cases = [ + ( + "../x".to_string(), + [ + temp_dir.path().join("x_regions.bed"), + temp_dir.path().join("x_windows.bedgraph"), + ], + ), + ( + "a/b".to_string(), + [ + output_dir.join("a/b_regions.bed"), + output_dir.join("a/b_windows.bedgraph"), + ], + ), + ( + absolute_prefix_string.clone(), + [ + format!("{absolute_prefix_string}_regions.bed").into(), + format!("{absolute_prefix_string}_windows.bedgraph").into(), + ], + ), + ]; + + for (prefix, targets) in invalid_cases { + for target in &targets { + fs::write(target, sentinel).unwrap(); + } + let rejected = entropy_command() + .args([ + "--out-bed", + output_dir.to_str().unwrap(), + "--regions", + "../tests/resources/entropy_test_regions.bed", + "--prefix", + &prefix, + "--filter-threshold", + "0", + "--force", + ]) + .output() + .unwrap(); + assert!(!rejected.status.success(), "prefix {prefix:?} was accepted"); + assert!( + String::from_utf8_lossy(&rejected.stderr) + .contains("exactly one non-empty filename component"), + "unexpected stderr for {prefix:?}:\n{}", + String::from_utf8_lossy(&rejected.stderr) + ); + for target in targets { + assert_eq!(fs::read(target).unwrap(), sentinel); + } + assert_eq!(fs::read(&unrelated).unwrap(), unrelated_contents); + } + + let ordinary = entropy_command() + .args([ + "--out-bed", + output_dir.to_str().unwrap(), + "--regions", + "../tests/resources/entropy_test_regions.bed", + "--prefix", + "ordinary", + "--filter-threshold", + "0", + "--force", + ]) + .output() + .unwrap(); + assert!( + ordinary.status.success(), + "ordinary prefix failed:\n{}", + String::from_utf8_lossy(&ordinary.stderr) + ); + assert!(output_dir.join("ordinary_regions.bed").is_file()); + assert!(output_dir.join("ordinary_windows.bedgraph").is_file()); + assert_eq!(fs::read(unrelated).unwrap(), unrelated_contents); +} + +#[test] +fn entropy_regions_second_open_failure_preserves_first_target() { + let temp_dir = tempfile::tempdir().unwrap(); + let output_dir = temp_dir.path().join("entropy"); + fs::create_dir(&output_dir).unwrap(); + + // NAME_MAX is 255 bytes on the supported Unix platforms. The first + // suffix brings this component to exactly 255 bytes, while the second + // makes it 260 bytes and must fail before the first target is truncated. + let prefix = "x".repeat(243); + let regions_output = output_dir.join(format!("{prefix}_regions.bed")); + let sentinel = b"existing scientific output\n"; + fs::write(®ions_output, sentinel).unwrap(); + + let rejected = entropy_command() + .args([ + "--out-bed", + output_dir.to_str().unwrap(), + "--regions", + "../tests/resources/entropy_test_regions.bed", + "--prefix", + &prefix, + "--filter-threshold", + "0", + "--force", + ]) + .output() + .unwrap(); + + assert!(!rejected.status.success()); + assert_eq!(fs::read(regions_output).unwrap(), sentinel); +} + +#[cfg(unix)] +#[test] +fn entropy_file_force_rejects_existing_and_dangling_symlinks() { + use std::os::unix::fs::symlink; + + let temp_dir = tempfile::tempdir().unwrap(); + let sentinel = b"outside scientific output\n"; + let outside_output = temp_dir.path().join("outside.bed"); + fs::write(&outside_output, sentinel).unwrap(); + let linked_output = temp_dir.path().join("linked.bed"); + symlink(&outside_output, &linked_output).unwrap(); + + let run = |output_path: &Path| { + entropy_command() + .args([ + "--out-bed", + output_path.to_str().unwrap(), + "--filter-threshold", + "0", + "--force", + ]) + .output() + .unwrap() + }; + + let linked = run(&linked_output); + assert!(!linked.status.success()); + assert!( + String::from_utf8_lossy(&linked.stderr).contains("regular file"), + "unexpected stderr:\n{}", + String::from_utf8_lossy(&linked.stderr) + ); + assert_eq!(fs::read(&outside_output).unwrap(), sentinel); + assert!(fs::symlink_metadata(&linked_output) + .unwrap() + .file_type() + .is_symlink()); + + let missing_output = temp_dir.path().join("missing.bed"); + let dangling_output = temp_dir.path().join("dangling.bed"); + symlink(&missing_output, &dangling_output).unwrap(); + let dangling = run(&dangling_output); + assert!(!dangling.status.success()); + assert!( + String::from_utf8_lossy(&dangling.stderr).contains("regular file"), + "unexpected stderr:\n{}", + String::from_utf8_lossy(&dangling.stderr) + ); + assert!(!missing_output.exists()); + assert!(fs::symlink_metadata(dangling_output) + .unwrap() + .file_type() + .is_symlink()); +} + +#[cfg(unix)] +#[test] +fn entropy_regions_force_rejects_symlinks_before_pair_mutation() { + use std::os::unix::fs::symlink; + + let temp_dir = tempfile::tempdir().unwrap(); + let output_dir = temp_dir.path().join("entropy"); + fs::create_dir(&output_dir).unwrap(); + let sentinel = b"outside scientific output\n"; + + let run = |prefix: &str| { + entropy_command() + .args([ + "--out-bed", + output_dir.to_str().unwrap(), + "--regions", + "../tests/resources/entropy_test_regions.bed", + "--prefix", + prefix, + "--filter-threshold", + "0", + "--force", + ]) + .output() + .unwrap() + }; + + let outside_regions = temp_dir.path().join("outside-regions.bed"); + let outside_windows = temp_dir.path().join("outside-windows.bedgraph"); + fs::write(&outside_regions, sentinel).unwrap(); + fs::write(&outside_windows, sentinel).unwrap(); + let linked_regions = output_dir.join("linked_regions.bed"); + let linked_windows = output_dir.join("linked_windows.bedgraph"); + symlink(&outside_regions, &linked_regions).unwrap(); + symlink(&outside_windows, &linked_windows).unwrap(); + + let linked = run("linked"); + assert!(!linked.status.success()); + assert!(String::from_utf8_lossy(&linked.stderr).contains("regular file")); + assert_eq!(fs::read(&outside_regions).unwrap(), sentinel); + assert_eq!(fs::read(&outside_windows).unwrap(), sentinel); + assert!(fs::symlink_metadata(&linked_regions) + .unwrap() + .file_type() + .is_symlink()); + assert!(fs::symlink_metadata(&linked_windows) + .unwrap() + .file_type() + .is_symlink()); + + let first_output = output_dir.join("second_regions.bed"); + let outside_second = temp_dir.path().join("outside-second.bedgraph"); + fs::write(&first_output, sentinel).unwrap(); + fs::write(&outside_second, sentinel).unwrap(); + symlink(&outside_second, output_dir.join("second_windows.bedgraph")) + .unwrap(); + let second = run("second"); + assert!(!second.status.success()); + assert_eq!(fs::read(&first_output).unwrap(), sentinel); + assert_eq!(fs::read(&outside_second).unwrap(), sentinel); + + let dangling_first = output_dir.join("dangling_regions.bed"); + let missing_second = temp_dir.path().join("missing-second.bedgraph"); + fs::write(&dangling_first, sentinel).unwrap(); + symlink(&missing_second, output_dir.join("dangling_windows.bedgraph")) + .unwrap(); + let dangling = run("dangling"); + assert!(!dangling.status.success()); + assert_eq!(fs::read(dangling_first).unwrap(), sentinel); + assert!(!missing_second.exists()); +} + +#[test] +fn entropy_stdout_output_remains_available() { + let output = + entropy_command().args(["--filter-threshold", "0"]).output().unwrap(); + assert!( + output.status.success(), + "stdout output failed:\n{}", + String::from_utf8_lossy(&output.stderr) + ); + assert!(!output.stdout.is_empty()); +} From b63f53af779710436d46a82807aaa3b2a13bc0cf Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 10:31:57 -0700 Subject: [PATCH 06/10] Reject empty probability threshold samples --- modkit-core/src/sample_probs/mod.rs | 74 +++++++- modkit-core/src/summarize.rs | 7 + .../tests/test_empty_qual_hist_thresholds.rs | 161 ++++++++++++++++++ 3 files changed, 235 insertions(+), 7 deletions(-) create mode 100644 modkit/tests/test_empty_qual_hist_thresholds.rs diff --git a/modkit-core/src/sample_probs/mod.rs b/modkit-core/src/sample_probs/mod.rs index a0080af9..3d255266 100644 --- a/modkit-core/src/sample_probs/mod.rs +++ b/modkit-core/src/sample_probs/mod.rs @@ -110,6 +110,11 @@ pub(crate) struct QualHist { } impl QualHist { + pub(crate) fn has_probability_observations(&self) -> bool { + self.base_totals.iter().any(|count| *count > 0) + || self.mods_hists.iter().any(|hist| hist.total > 0) + } + pub(crate) fn clear(&mut self) { for ar in self.hist.iter_mut() { ar.iter_mut().for_each(|x| *x = 0u64); @@ -569,8 +574,11 @@ impl QualHist { max_thresholds_per_base: Option<[f32; 4]>, multi_progress: &MultiProgress, ) -> anyhow::Result<[f32; 4]> { - if self.ok_records == 0 { - bail!("Failed to sample any records to estimate threshold.") + if !self.has_probability_observations() { + bail!( + "cannot calculate automatic thresholds because no \ + modification probabilities were sampled" + ) } let mut base_thresholds = [0f32; 4]; let filter_percentile = if filter_percentile >= 1.0f32 { @@ -580,10 +588,10 @@ impl QualHist { } else { filter_percentile * 100f32 }; - for (base, vals) in QualHist::percentiles( - &self.hist, - &self.base_totals, - &vec![filter_percentile], + for (base, vals) in self.base_level_percentiles( + &[filter_percentile], + multi_progress, + false, ) { assert_eq!(vals.len(), 1); let (t, q) = (vals[0].threshold, vals[0].qual); @@ -623,11 +631,18 @@ impl QualHist { base_thresholds[base as usize] = t_fb; } } else { + let modified_probability_count = self + .mods_hists + .iter() + .filter(|hist| hist.dna_base == base) + .fold(0u64, |total, hist| total.saturating_add(hist.total)); + let probability_count = self.base_totals[base as usize] + .saturating_add(modified_probability_count); multi_progress.suspend(|| { info!( "setting threshold {t} (qual: {q}) for base {base}, \ percentile {}, {} total probabilites", - filter_percentile, self.base_totals[base as usize] + filter_percentile, probability_count ); }); base_thresholds[base as usize] = t; @@ -1759,3 +1774,48 @@ fn byte_to_bool_positions(b: u8, agg: &mut [u32; SIZE]) { agg[i] = count.saturating_add(1u32); } } + +#[cfg(test)] +mod empty_threshold_tests { + use indicatif::MultiProgress; + + use super::{ModHist, QualHist}; + use crate::mod_base_code::{DnaBase, METHYL_CYTOSINE}; + + fn all_modified_qual_hist() -> QualHist { + let mut hist = [0u64; 256]; + hist[128] = 1; + let mut qual_hist = QualHist::default(); + qual_hist.mods_hists.push(ModHist { + total: 1, + mod_code: METHYL_CYTOSINE, + dna_base: DnaBase::C, + hist, + }); + qual_hist.num_records_with_base_mods[DnaBase::C as usize] = 1; + qual_hist.ok_records = 1; + qual_hist + } + + #[test] + fn probability_observations_are_not_inferred_from_record_count() { + let mut qual_hist = QualHist::default(); + qual_hist.ok_records = 1; + + assert!(!qual_hist.has_probability_observations()); + let error = qual_hist + .get_base_thresholds(0.1, None, &MultiProgress::new()) + .expect_err("record count alone must not enable auto-thresholding"); + assert!(error.to_string().contains("no modification probabilities")); + } + + #[test] + fn all_modified_observations_support_automatic_thresholds() { + let qual_hist = all_modified_qual_hist(); + assert!(qual_hist.has_probability_observations()); + let thresholds = qual_hist + .get_base_thresholds(0.1, None, &MultiProgress::new()) + .expect("all-modified observations must produce a threshold"); + assert!(thresholds[DnaBase::C as usize] > 0.0); + } +} diff --git a/modkit-core/src/summarize.rs b/modkit-core/src/summarize.rs index 04e90d7b..6033a6f0 100644 --- a/modkit-core/src/summarize.rs +++ b/modkit-core/src/summarize.rs @@ -1,6 +1,7 @@ use std::collections::{HashMap, HashSet}; use std::path::PathBuf; +use anyhow::bail; use common_macros::hash_map; use derive_new::new; use indicatif::{MultiProgress, ParallelProgressIterator}; @@ -154,6 +155,12 @@ impl<'a> ModSummary<'a> { let base_thresholds = if let Some(x) = base_thresholds { x } else { + if !qual_hist.has_probability_observations() { + bail!( + "cannot calculate automatic thresholds because no \ + modification probabilities were sampled" + ) + } let mut agg = [0f32; 4]; let base_level_threshs = qual_hist.base_level_percentiles( &[filter_percentile], diff --git a/modkit/tests/test_empty_qual_hist_thresholds.rs b/modkit/tests/test_empty_qual_hist_thresholds.rs new file mode 100644 index 00000000..dac1ed5b --- /dev/null +++ b/modkit/tests/test_empty_qual_hist_thresholds.rs @@ -0,0 +1,161 @@ +use std::fs; +use std::path::Path; +use std::process::{Command, Output}; + +const MOD_BAM: &str = "../tests/resources/bc_anchored_10_reads.sorted.bam"; +const EMPTY_TAGS_BAM: &str = "../tests/resources/empty-tags.sorted.bam"; +const REFERENCE: &str = "../tests/resources/CGI_ladder_3.6kb_ref.fa"; +const EMPTY_THRESHOLD_ERROR: &str = + "cannot calculate automatic thresholds because no modification \ + probabilities were sampled"; + +fn run_modkit(args: &[&str]) -> Output { + Command::new(Path::new(env!("CARGO_BIN_EXE_modkit"))) + .args(args) + .output() + .unwrap() +} + +fn assert_success(output: &Output) { + assert!( + output.status.success(), + "command failed:\n{}", + String::from_utf8_lossy(&output.stderr) + ); +} + +fn assert_empty_auto_failure(output: &Output) { + assert!(!output.status.success()); + assert!( + String::from_utf8_lossy(&output.stderr).contains(EMPTY_THRESHOLD_ERROR), + "unexpected stderr:\n{}", + String::from_utf8_lossy(&output.stderr) + ); +} + +#[test] +fn pileup_rejects_empty_automatic_observations_but_controls_succeed() { + let temp_dir = tempfile::tempdir().unwrap(); + let include_bed = temp_dir.path().join("empty_site.bed"); + fs::write(&include_bed, "oligo_1512_adapters\t0\t1\n").unwrap(); + + let run = |output_path: &Path, threshold_args: &[&str]| { + Command::new(Path::new(env!("CARGO_BIN_EXE_modkit"))) + .args([ + "pileup", + MOD_BAM, + output_path.to_str().unwrap(), + "--include-bed", + include_bed.to_str().unwrap(), + "--modified-bases", + "C:m", + "--reference", + REFERENCE, + "--threads", + "1", + "--sampling-threads", + "1", + "--io-threads", + "1", + "--suppress-progress", + ]) + .args(threshold_args) + .output() + .unwrap() + }; + + let automatic_path = temp_dir.path().join("automatic.bed"); + let automatic = run(&automatic_path, &[]); + assert_empty_auto_failure(&automatic); + + for (name, threshold_args) in [ + ("explicit.bed", vec!["--filter-threshold", "0"]), + ("no_filter.bed", vec!["--no-filtering"]), + ] { + let output_path = temp_dir.path().join(name); + let output = run(&output_path, &threshold_args); + assert_success(&output); + assert_eq!(fs::metadata(output_path).unwrap().len(), 0); + } +} + +#[test] +fn summary_rejects_empty_automatic_observations_but_explicit_succeeds() { + let automatic = run_modkit(&[ + "summary", + EMPTY_TAGS_BAM, + "--num-reads", + "10", + "--threads", + "1", + "--io-threads", + "1", + "--tsv", + "--suppress-progress", + ]); + assert_empty_auto_failure(&automatic); + assert!(automatic.stdout.is_empty()); + + let explicit = run_modkit(&[ + "summary", + EMPTY_TAGS_BAM, + "--num-reads", + "10", + "--threads", + "1", + "--io-threads", + "1", + "--filter-threshold", + "0", + "--tsv", + "--suppress-progress", + ]); + assert_success(&explicit); + assert!(String::from_utf8_lossy(&explicit.stdout) + .contains("total_reads_used\t0")); +} + +#[test] +fn extract_rejects_empty_automatic_observations_but_controls_succeed() { + let temp_dir = tempfile::tempdir().unwrap(); + let include_bed = temp_dir.path().join("empty_site.bed"); + fs::write(&include_bed, "oligo_1512_adapters\t0\t1\n").unwrap(); + + let run = |output_path: &Path, threshold_args: &[&str]| { + Command::new(Path::new(env!("CARGO_BIN_EXE_modkit"))) + .args([ + "extract", + "calls", + MOD_BAM, + output_path.to_str().unwrap(), + "--include-bed", + include_bed.to_str().unwrap(), + "--threads", + "1", + "--io-threads", + "1", + "--sample-num-reads", + "10", + "--suppress-progress", + "--no-headers", + ]) + .args(threshold_args) + .output() + .unwrap() + }; + + let automatic_path = temp_dir.path().join("automatic.tsv"); + let automatic = run(&automatic_path, &[]); + assert_empty_auto_failure(&automatic); + assert!(!automatic_path.exists()); + + for (name, threshold_args) in [ + ("explicit.tsv", vec!["--filter-threshold", "0"]), + ("no_filter.tsv", vec!["--no-filtering"]), + ] { + let output_path = temp_dir.path().join(name); + let output = run(&output_path, &threshold_args); + assert_success(&output); + assert_eq!(fs::metadata(output_path).unwrap().len(), 0); + } +} From 456d543b870128d8ef3bca8771fd9384680c2f92 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 12:55:26 -0700 Subject: [PATCH 07/10] Fix extract CRAM reference plumbing --- modkit-core/src/extract/subcommand.rs | 58 +++-- modkit-core/src/extract/util.rs | 86 +++--- modkit/tests/test_extract_cram_reference.rs | 274 ++++++++++++++++++++ 3 files changed, 365 insertions(+), 53 deletions(-) create mode 100644 modkit/tests/test_extract_cram_reference.rs diff --git a/modkit-core/src/extract/subcommand.rs b/modkit-core/src/extract/subcommand.rs index 42b4d6ed..3dc3cf50 100644 --- a/modkit-core/src/extract/subcommand.rs +++ b/modkit-core/src/extract/subcommand.rs @@ -41,7 +41,10 @@ use crate::reads_sampler::sampling_schedule::SamplingSchedule; use crate::record_processor::WithRecords; use crate::sample_probs::calc_per_base_thresholds_from_indexed_hts_file; use crate::threshold_mod_caller::MultipleThresholdModCaller; -use crate::util::{format_errors_table, get_ticker, Region, KMER_SIZE}; +use crate::util::{ + format_errors_table, get_ticker, preflight_cram_input, reader_is_cram, + set_reference_for_cram_reader, Region, KMER_SIZE, +}; use crate::writers::TsvWriter; #[derive(Subcommand)] @@ -150,6 +153,13 @@ impl EntryExtractFull { .transpose()?; let mut reader = get_serial_reader(&self.input_args.in_bam)?; + set_reference_for_cram_reader(&mut reader, self.reference.as_ref())?; + if !self.using_stdin() && reader_is_cram(&reader) { + preflight_cram_input( + Path::new(&self.input_args.in_bam), + self.reference.as_ref(), + )?; + } let header = reader.header().to_owned(); let queue_size = self.input_args.queue_size; @@ -219,13 +229,16 @@ impl EntryExtractFull { (_, true) | (None, false) => None, (Some(num_reads), false) => { match bam::IndexedReader::from_path(&self.input_args.in_bam) { - Ok(_) => Some(SamplingSchedule::from_num_reads( - &self.input_args.in_bam, - num_reads, - region.as_ref(), - reference_position_filter.include_pos.as_ref(), - reference_position_filter.include_unmapped_reads, - )?), + Ok(_) => { + Some(SamplingSchedule::from_num_reads_with_reference( + &self.input_args.in_bam, + self.reference.as_ref(), + num_reads, + region.as_ref(), + reference_position_filter.include_pos.as_ref(), + reference_position_filter.include_unmapped_reads, + )?) + } Err(_) => { debug!( "cannot use sampling schedule without index, \ @@ -251,6 +264,7 @@ impl EntryExtractFull { let threads = self.input_args.threads; let mapped_only = self.input_args.mapped_only; let in_bam = self.input_args.in_bam.clone(); + let reference = self.reference.clone(); let kmer_size = self.input_args.kmer_size; let allow_non_primary = self.input_args.allow_non_primary; let remove_inferred = self.input_args.ignore_implicit; @@ -259,6 +273,7 @@ impl EntryExtractFull { super::util::run_extract_reads( reader, in_bam, + reference, references_and_intervals, schedule, collapse_method, @@ -522,6 +537,7 @@ impl EntryExtractCalls { { calc_per_base_thresholds_from_stream( &Path::new(&self.input_args.in_bam).to_path_buf(), + self.reference.as_ref(), self.sample_num_reads, false, position_filter, @@ -613,6 +629,13 @@ impl EntryExtractCalls { .transpose()?; let mut reader = get_serial_reader(&self.input_args.in_bam)?; + set_reference_for_cram_reader(&mut reader, self.reference.as_ref())?; + if !self.using_stdin() && reader_is_cram(&reader) { + preflight_cram_input( + Path::new(&self.input_args.in_bam), + self.reference.as_ref(), + )?; + } let header = reader.header().to_owned(); let tid_to_name = (0..header.target_count()) @@ -776,13 +799,16 @@ impl EntryExtractCalls { (_, true) | (None, false) => None, (Some(num_reads), false) => { match bam::IndexedReader::from_path(&self.input_args.in_bam) { - Ok(_) => Some(SamplingSchedule::from_num_reads( - &self.input_args.in_bam, - num_reads, - region.as_ref(), - reference_position_filter.include_pos.as_ref(), - reference_position_filter.include_unmapped_reads, - )?), + Ok(_) => { + Some(SamplingSchedule::from_num_reads_with_reference( + &self.input_args.in_bam, + self.reference.as_ref(), + num_reads, + region.as_ref(), + reference_position_filter.include_pos.as_ref(), + reference_position_filter.include_unmapped_reads, + )?) + } Err(_) => { debug!( "cannot use sampling schedule without index, \ @@ -811,6 +837,7 @@ impl EntryExtractCalls { let threads = self.input_args.threads; let mapped_only = self.input_args.mapped_only; let in_bam = self.input_args.in_bam.clone(); + let reference = self.reference.clone(); let kmer_size = self.input_args.kmer_size; let allow_non_primary = self.input_args.allow_non_primary; let remove_inferred = self.input_args.ignore_implicit; @@ -819,6 +846,7 @@ impl EntryExtractCalls { super::util::run_extract_reads( reader, in_bam, + reference, references_and_intervals, schedule, collapse_method, diff --git a/modkit-core/src/extract/util.rs b/modkit-core/src/extract/util.rs index e628a9b5..1b2ed661 100644 --- a/modkit-core/src/extract/util.rs +++ b/modkit-core/src/extract/util.rs @@ -15,14 +15,15 @@ use crate::read_ids_to_base_mod_probs::{ ModProfile, ReadBaseModProfile, ReadsBaseModProfile, }; use crate::reads_sampler::record_sampler::RecordSampler; -use crate::reads_sampler::sample_reads_from_interval; +use crate::reads_sampler::sample_reads_from_interval_with_reference; use crate::reads_sampler::sampling_schedule::SamplingSchedule; use crate::record_processor::WithRecords; use crate::sample_probs::QualHist; use crate::util::{ get_guage, get_master_progress_bar, get_query_name_string, get_reference_mod_strand, get_subroutine_progress_bar, get_targets, - get_ticker, record_is_primary, Region, Strand, + get_ticker, record_is_primary, set_reference_for_cram_indexed_reader, + set_reference_for_cram_reader, Region, Strand, }; use derive_new::new; use indicatif::{MultiProgress, ParallelProgressIterator, ProgressBar}; @@ -371,6 +372,7 @@ pub(super) fn load_regions( pub(super) fn run_extract_reads( mut reader: bam::Reader, in_bam: String, + reference_fasta: Option, references_and_intervals: Option, schedule: Option, collapse_method: Option, @@ -437,35 +439,38 @@ pub(super) fn run_extract_reads( .unwrap_or_else(|| { RecordSampler::new_passthrough() }); - let batch_result = sample_reads_from_interval::< - ReadsBaseModProfile, - >( - &bam_fp, - cc.chrom_tid(), - cc.start_pos(), - cc.end_pos(), - cc.prev_end(), - record_sampler, - collapse_method.as_ref(), - edge_filter.as_ref(), - None, - false, - allow_non_primary, - Some(kmer_size), - ) - .map(|reads_base_mod_profile| { - if remove_inferred { - reads_base_mod_profile.remove_inferred() - } else { - reads_base_mod_profile - } - }) - .map(|reads_base_mod_profile| { - reference_position_filter - .filter_read_base_mod_probs( - reads_base_mod_profile, - ) - }); + let batch_result = + sample_reads_from_interval_with_reference::< + ReadsBaseModProfile, + >( + &bam_fp, + reference_fasta.as_ref(), + cc.chrom_tid(), + cc.start_pos(), + cc.end_pos(), + cc.prev_end(), + record_sampler, + collapse_method.as_ref(), + edge_filter.as_ref(), + None, + false, + allow_non_primary, + Some(kmer_size), + ); + let batch_result = batch_result + .map(|reads_base_mod_profile| { + if remove_inferred { + reads_base_mod_profile.remove_inferred() + } else { + reads_base_mod_profile + } + }) + .map(|reads_base_mod_profile| { + reference_position_filter + .filter_read_base_mod_probs( + reads_base_mod_profile, + ) + }); let num_reads_success = batch_result .as_ref() @@ -503,13 +508,16 @@ pub(super) fn run_extract_reads( } else { debug!("processing unmapped reads"); } - let reader = bam::IndexedReader::from_path(&bam_fp) - .and_then(|mut reader| { - reader.fetch(FetchDefinition::Unmapped).map(|_| reader) - }) - .and_then(|mut reader| { - reader.set_threads(threads).map(|_| reader) - }); + let reader = (|| -> anyhow::Result { + let mut reader = bam::IndexedReader::from_path(&bam_fp)?; + set_reference_for_cram_indexed_reader( + &mut reader, + reference_fasta.as_ref(), + )?; + reader.fetch(FetchDefinition::Unmapped)?; + reader.set_threads(threads)?; + Ok(reader) + })(); match reader { Ok(mut reader) => { let (skip, fail) = process_records_to_chan( @@ -831,6 +839,7 @@ impl ReadModStatsProcessor { pub(super) fn calc_per_base_thresholds_from_stream( bam_fp: &PathBuf, + reference_fasta: Option<&PathBuf>, num_reads: usize, allow_non_primary: bool, stranded_position_filter: Option>, @@ -847,6 +856,7 @@ pub(super) fn calc_per_base_thresholds_from_stream( ); }); let mut records = bam::Reader::from_path(bam_fp)?; + set_reference_for_cram_reader(&mut records, reference_fasta)?; records.set_threads(io_threads)?; QualHist::from_records( records.records(), diff --git a/modkit/tests/test_extract_cram_reference.rs b/modkit/tests/test_extract_cram_reference.rs new file mode 100644 index 00000000..7af06bb5 --- /dev/null +++ b/modkit/tests/test_extract_cram_reference.rs @@ -0,0 +1,274 @@ +use std::fs; +use std::path::{Path, PathBuf}; +use std::process::{Command, Output}; + +const BAM: &str = "../tests/resources/bc_anchored_10_reads.sorted.bam"; +const CRAM: &str = "../tests/resources/bc_anchored_10_reads.sorted.cram"; +const REFERENCE: &str = "../tests/resources/CGI_ladder_3.6kb_ref.fa"; +const WRONG_REFERENCE: &str = "../tests/resources/genome_for_phased_test.fasta"; +const UNMAPPED_BAM: &str = + "../tests/resources/bc_anchored_10_reads.unmapped.bam"; +const UNMAPPED_CRAM: &str = + "../tests/resources/bc_anchored_10_reads_unmapped.cram"; + +fn run_modkit(args: &[String], empty_reference_cache: Option<&Path>) -> Output { + let exe = Path::new(env!("CARGO_BIN_EXE_modkit")); + assert!(exe.exists()); + let mut command = Command::new(exe); + command.args(args); + if let Some(empty_reference_cache) = empty_reference_cache { + let reference_pattern = empty_reference_cache.join("%2s/%2s/%s"); + command + .env("REF_PATH", &reference_pattern) + .env("REF_CACHE", &reference_pattern); + } + command.output().unwrap() +} + +fn run_extract( + command: &str, + input: &Path, + output: &Path, + reference: Option<&Path>, + extra_args: &[&str], +) -> Output { + let mut args = vec![ + "extract".to_string(), + command.to_string(), + input.to_str().unwrap().to_string(), + output.to_str().unwrap().to_string(), + "--threads".to_string(), + "1".to_string(), + "--io-threads".to_string(), + "1".to_string(), + "--no-headers".to_string(), + "--suppress-progress".to_string(), + "--force".to_string(), + ]; + if let Some(reference) = reference { + args.extend([ + "--reference".to_string(), + reference.to_str().unwrap().to_string(), + ]); + } + args.extend(extra_args.iter().map(|arg| arg.to_string())); + run_modkit(&args, None) +} + +fn assert_success(output: &Output) { + assert!( + output.status.success(), + "command failed:\n{}", + String::from_utf8_lossy(&output.stderr) + ); +} + +fn assert_parity(expected: &Path, observed: &Path, expected_rows: usize) { + let expected = fs::read(expected).unwrap(); + let observed = fs::read(observed).unwrap(); + assert_eq!( + row_count(&expected), + expected_rows, + "unexpected control row count" + ); + assert_eq!( + row_count(&observed), + expected_rows, + "unexpected CRAM row count" + ); + assert_eq!(expected, observed); +} + +fn row_count(output: &[u8]) -> usize { + std::str::from_utf8(output).unwrap().lines().count() +} + +fn unindexed_copies(temp_dir: &Path) -> (PathBuf, PathBuf) { + let bam = temp_dir.join("unindexed.bam"); + let cram = temp_dir.join("unindexed.cram"); + fs::copy(BAM, &bam).unwrap(); + fs::copy(CRAM, &cram).unwrap(); + assert!(!bam.with_extension("bam.bai").exists()); + assert!(!cram.with_extension("cram.crai").exists()); + (bam, cram) +} + +#[test] +fn mapped_cram_matches_bam_for_full_and_calls() { + let temp_dir = tempfile::tempdir().unwrap(); + let (unindexed_bam, unindexed_cram) = unindexed_copies(temp_dir.path()); + let reference = Path::new(REFERENCE); + + let full_bam_indexed = temp_dir.path().join("full.bam.indexed.tsv"); + let full_cram_indexed = temp_dir.path().join("full.cram.indexed.tsv"); + let full_bam_serial = temp_dir.path().join("full.bam.serial.tsv"); + let full_cram_serial = temp_dir.path().join("full.cram.serial.tsv"); + for (input, output, extra_args) in [ + (Path::new(BAM), full_bam_indexed.as_path(), &[][..]), + (Path::new(CRAM), full_cram_indexed.as_path(), &[][..]), + ( + unindexed_bam.as_path(), + full_bam_serial.as_path(), + &["--ignore-index"][..], + ), + ( + unindexed_cram.as_path(), + full_cram_serial.as_path(), + &["--ignore-index"][..], + ), + ] { + assert_success(&run_extract( + "full", + input, + output, + Some(reference), + extra_args, + )); + } + assert_parity(&full_bam_indexed, &full_cram_indexed, 218); + assert_parity(&full_bam_serial, &full_cram_serial, 218); + assert_parity(&full_bam_indexed, &full_bam_serial, 218); + + for (mode, mode_args) in [ + ("automatic", &[][..]), + ("explicit", &["--filter-threshold", "0"][..]), + ("no-filter", &["--no-filtering"][..]), + ] { + let bam_indexed = + temp_dir.path().join(format!("calls.{mode}.bam.indexed.tsv")); + let cram_indexed = + temp_dir.path().join(format!("calls.{mode}.cram.indexed.tsv")); + let bam_serial = + temp_dir.path().join(format!("calls.{mode}.bam.serial.tsv")); + let cram_serial = + temp_dir.path().join(format!("calls.{mode}.cram.serial.tsv")); + + for (input, output, extra_args) in [ + (Path::new(BAM), bam_indexed.as_path(), mode_args), + (Path::new(CRAM), cram_indexed.as_path(), mode_args), + (unindexed_bam.as_path(), bam_serial.as_path(), mode_args), + (unindexed_cram.as_path(), cram_serial.as_path(), mode_args), + ] { + assert_success(&run_extract( + "calls", + input, + output, + Some(reference), + extra_args, + )); + } + assert_parity(&bam_indexed, &cram_indexed, 109); + assert_parity(&bam_serial, &cram_serial, 109); + assert_parity(&bam_indexed, &bam_serial, 109); + } +} + +#[test] +fn indexed_cram_num_reads_schedule_uses_reference() { + let temp_dir = tempfile::tempdir().unwrap(); + for (command, mode_args, expected_rows) in [ + ("full", &[][..], 218), + ("calls", &["--filter-threshold", "0"][..], 109), + ] { + let output_path = temp_dir.path().join(format!("{command}.tsv")); + let mut args = mode_args.to_vec(); + args.extend(["--num-reads", "10042"]); + assert_success(&run_extract( + command, + Path::new(CRAM), + &output_path, + Some(Path::new(REFERENCE)), + &args, + )); + let output = fs::read(&output_path).unwrap(); + assert_eq!(row_count(&output), expected_rows); + } +} + +#[test] +fn unmapped_cram_without_reference_matches_bam() { + let temp_dir = tempfile::tempdir().unwrap(); + for (command, mode_args, expected_rows) in + [("full", &[][..], 218), ("calls", &["--no-filtering"][..], 109)] + { + let bam_output = temp_dir.path().join(format!("{command}.bam.tsv")); + let cram_output = temp_dir.path().join(format!("{command}.cram.tsv")); + let mut args = mode_args.to_vec(); + args.push("--ignore-index"); + assert_success(&run_extract( + command, + Path::new(UNMAPPED_BAM), + &bam_output, + None, + &args, + )); + assert_success(&run_extract( + command, + Path::new(UNMAPPED_CRAM), + &cram_output, + None, + &args, + )); + assert_parity(&bam_output, &cram_output, expected_rows); + } +} + +#[test] +fn reference_errors_do_not_mutate_extract_outputs() { + let temp_dir = tempfile::tempdir().unwrap(); + let empty_reference_cache = temp_dir.path().join("empty-reference-cache"); + fs::create_dir(&empty_reference_cache).unwrap(); + let disguised_cram = temp_dir.path().join("reference-dependent.bam"); + fs::copy(CRAM, &disguised_cram).unwrap(); + let missing_reference = temp_dir.path().join("missing-reference.fa"); + + for (reference_name, reference) in [ + ("none", None), + ("missing", Some(missing_reference.as_path())), + ("wrong", Some(Path::new(WRONG_REFERENCE))), + ] { + for command in ["full", "calls"] { + for preexisting in [false, true] { + let output_path = temp_dir.path().join(format!( + "{command}-{reference_name}-{}.tsv", + if preexisting { "existing" } else { "absent" } + )); + if preexisting { + fs::write(&output_path, b"sentinel\n").unwrap(); + } + let mut args = vec![ + "extract".to_string(), + command.to_string(), + disguised_cram.to_str().unwrap().to_string(), + output_path.to_str().unwrap().to_string(), + "--threads".to_string(), + "1".to_string(), + "--io-threads".to_string(), + "1".to_string(), + "--suppress-progress".to_string(), + "--force".to_string(), + ]; + if command == "calls" { + args.push("--no-filtering".to_string()); + } + if let Some(reference) = reference { + args.extend([ + "--reference".to_string(), + reference.to_str().unwrap().to_string(), + ]); + } + let output = + run_modkit(&args, Some(empty_reference_cache.as_path())); + assert!( + !output.status.success(), + "{command} unexpectedly accepted {reference_name} reference" + ); + if preexisting { + assert_eq!(fs::read(&output_path).unwrap(), b"sentinel\n"); + } else { + assert!(!output_path.exists()); + } + } + } + } +} From 07035502328047fda1615b3dd340acf7e4f88943 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 13:55:44 -0700 Subject: [PATCH 08/10] Handle indexed zero-reference extract inputs --- modkit-core/src/extract/util.rs | 32 +- modkit-core/src/interval_chunks.rs | 275 +++++++++++++++--- modkit-core/src/modbam_util/subcommands.rs | 3 +- modkit-core/src/sample_probs/mod.rs | 38 ++- .../test_extract_zero_reference_targets.rs | 265 +++++++++++++++++ 5 files changed, 551 insertions(+), 62 deletions(-) create mode 100644 modkit/tests/test_extract_zero_reference_targets.rs diff --git a/modkit-core/src/extract/util.rs b/modkit-core/src/extract/util.rs index 1b2ed661..a13baa49 100644 --- a/modkit-core/src/extract/util.rs +++ b/modkit-core/src/extract/util.rs @@ -321,6 +321,7 @@ pub(super) fn load_regions( input_args.interval_size ); let reference_records = get_targets(reader.header(), region); + let zero_reference_targets = reference_records.is_empty(); let reference_records = if let Some(pf) = include_positions.as_ref() { pf.optimize_reference_records( @@ -331,14 +332,29 @@ pub(super) fn load_regions( reference_records }; - let feeder = ReferenceIntervalBatchesFeeder::new( - reference_records, - (input_args.threads as f32 * 1.5f32).floor() as usize, - input_args.interval_size, - false, - None, - None, - )?; + let batch_size = + (input_args.threads as f32 * 1.5f32).floor() as usize; + let feeder = if include_unmapped_reads && zero_reference_targets + { + ReferenceIntervalBatchesFeeder:: + new_allowing_zero_reference_targets( + reference_records, + batch_size, + input_args.interval_size, + false, + None, + None, + )? + } else { + ReferenceIntervalBatchesFeeder::new( + reference_records, + batch_size, + input_args.interval_size, + false, + None, + None, + )? + }; Some(feeder) } Err(_) => { diff --git a/modkit-core/src/interval_chunks.rs b/modkit-core/src/interval_chunks.rs index a6aa5541..f2fb5557 100644 --- a/modkit-core/src/interval_chunks.rs +++ b/modkit-core/src/interval_chunks.rs @@ -179,9 +179,13 @@ impl TotalLength for Vec { impl TotalLength for ReferenceIntervalBatchesFeeder { fn total_length(&self) -> u64 { - self.contigs.iter().fold(self.curr_contig.length as u64, |agg, next| { - agg.saturating_add(next.length as u64) - }) + self.contigs.iter().fold( + self.curr_contig + .as_ref() + .map(|record| record.length as u64) + .unwrap_or(0), + |agg, next| agg.saturating_add(next.length as u64), + ) } } @@ -192,9 +196,9 @@ pub(crate) struct ReferenceIntervalBatchesFeeder { motifs: Option, position_filter: Option>, combine_strands: bool, - curr_contig: ReferenceRecord, + // None denotes exhaustion, including an explicit zero-target terminal. + curr_contig: Option, curr_position: u32, - done: bool, } impl ReferenceIntervalBatchesFeeder { @@ -205,10 +209,53 @@ impl ReferenceIntervalBatchesFeeder { combine_strands: bool, multi_motif_locations: Option, position_filter: Option>, + ) -> anyhow::Result { + Self::new_inner( + reference_records, + batch_size, + interval_size, + combine_strands, + multi_motif_locations, + position_filter, + false, + ) + } + + /// Allow a zero-reference input to produce an already exhausted feeder. + /// This is only appropriate when the caller separately processes unmapped + /// records after the mapped interval phase. + pub(crate) fn new_allowing_zero_reference_targets( + reference_records: Vec, + batch_size: usize, + interval_size: u32, + combine_strands: bool, + multi_motif_locations: Option, + position_filter: Option>, + ) -> anyhow::Result { + Self::new_inner( + reference_records, + batch_size, + interval_size, + combine_strands, + multi_motif_locations, + position_filter, + true, + ) + } + + fn new_inner( + reference_records: Vec, + batch_size: usize, + interval_size: u32, + combine_strands: bool, + multi_motif_locations: Option, + position_filter: Option>, + allow_zero_reference_targets: bool, ) -> anyhow::Result { if combine_strands & !multi_motif_locations.is_some() { bail!("cannot combine strands without a motif") } + let zero_reference_targets = reference_records.is_empty(); let mut contigs = reference_records.into_iter().collect::>(); let n_contigs = contigs.iter().map(|r| r.tid).unique().count(); @@ -224,10 +271,14 @@ impl ReferenceIntervalBatchesFeeder { contigs.len() ); } - let curr_contig = contigs - .pop_front() - .ok_or(anyhow!("should be at least 1 contig"))?; - let curr_position = curr_contig.start; + let curr_contig = contigs.pop_front(); + if curr_contig.is_none() + && !(allow_zero_reference_targets && zero_reference_targets) + { + return Err(anyhow!("should be at least 1 contig")); + } + let curr_position = + curr_contig.as_ref().map(|record| record.start).unwrap_or(0); Ok(Self { contigs, batch_size, @@ -237,17 +288,16 @@ impl ReferenceIntervalBatchesFeeder { position_filter, curr_contig, curr_position, - done: false, }) } fn update_current(&mut self) { if let Some(reference_record) = self.contigs.pop_front() { self.curr_position = reference_record.start; - self.curr_contig = reference_record; + self.curr_contig = Some(reference_record); } else { debug!("no more records to process"); - self.done = true; + self.curr_contig = None; } } @@ -260,29 +310,30 @@ impl ReferenceIntervalBatchesFeeder { let mut batch_length = 0u32; loop { - if self.done { + if self.curr_contig.is_none() { break; } else if ret.len() >= self.batch_size { break; } - debug_assert!(self.curr_position < self.curr_contig.end()); + let curr_contig = self.curr_contig.as_ref().expect( + "non-terminal feeder must have a current reference record", + ); + debug_assert!(self.curr_position < curr_contig.end()); let start = self.curr_position; - let tid = self.curr_contig.tid; + let tid = curr_contig.tid; // in the case where we're on a large chrom end will be < length, // but batch length will be equal to interval size - let end = std::cmp::min( - start + self.interval_size, - self.curr_contig.end(), - ); + let end = + std::cmp::min(start + self.interval_size, curr_contig.end()); // get the sequence here. let (focus_positions, end) = if let Some(lookup) = self.motifs.as_mut() { // todo change everything to u64 let range = (start as u64)..(end as u64); let (fps, end) = lookup.get_motif_positions( - &self.curr_contig.name, + &curr_contig.name, tid, - self.curr_contig.end(), + curr_contig.end(), range, self.position_filter.as_ref(), self.combine_strands, @@ -296,7 +347,7 @@ impl ReferenceIntervalBatchesFeeder { } else { (FocusPositions2::AllPositions, end) }; - let end = std::cmp::min(end, self.curr_contig.end()); + let end = std::cmp::min(end, curr_contig.end()); // in the "short contig" case, chrom_coords.len() will be less than // interval size so batch length will be less than // interval size for a few rounds @@ -305,7 +356,7 @@ impl ReferenceIntervalBatchesFeeder { start, end, focus_positions, - end == self.curr_contig.end(), + end == curr_contig.end(), ); batch_length += chrom_coords.len(); batch.push(chrom_coords); @@ -319,7 +370,7 @@ impl ReferenceIntervalBatchesFeeder { } // might need to update the pointers, check if we're at the end of // this contig - if end >= self.curr_contig.end() { + if end >= curr_contig.end() { self.update_current(); } else { self.curr_position = end; @@ -348,9 +399,13 @@ impl Iterator for ReferenceIntervalBatchesFeeder { impl TotalLength for ChromCoordinatesFeeder { fn total_length(&self) -> u64 { - self.contigs.iter().fold(self.curr_contig.length as u64, |agg, next| { - agg.saturating_add(next.length as u64) - }) + self.contigs.iter().fold( + self.curr_contig + .as_ref() + .map(|record| record.length as u64) + .unwrap_or(0), + |agg, next| agg.saturating_add(next.length as u64), + ) } } @@ -362,9 +417,10 @@ pub(crate) struct ChromCoordinatesFeeder { // TODO: make this a borrow position_filter: Option>, combine_strands: bool, - curr_contig: ReferenceRecord, + // None denotes exhaustion, including an explicit zero-target terminal. + curr_contig: Option, curr_position: u32, - done: bool, + zero_reference_target_terminal: bool, } impl ChromCoordinatesFeeder { @@ -374,11 +430,49 @@ impl ChromCoordinatesFeeder { motifs: Option, combine_strands: bool, position_filter: Option>, + ) -> anyhow::Result { + Self::new_inner( + reference_records, + interval_size, + motifs, + combine_strands, + position_filter, + false, + ) + } + + /// Allow a zero-reference input to produce an already exhausted feeder. + /// Inputs whose targets are removed by motif filtering remain errors. + pub(crate) fn new_allowing_zero_reference_targets( + reference_records: &[ReferenceRecord], + interval_size: u32, + motifs: Option, + combine_strands: bool, + position_filter: Option>, + ) -> anyhow::Result { + Self::new_inner( + reference_records, + interval_size, + motifs, + combine_strands, + position_filter, + true, + ) + } + + fn new_inner( + reference_records: &[ReferenceRecord], + interval_size: u32, + motifs: Option, + combine_strands: bool, + position_filter: Option>, + allow_zero_reference_targets: bool, ) -> anyhow::Result { debug!("there is {} reference record(s)", reference_records.len()); if combine_strands & !motifs.is_some() { bail!("cannot combine strands without a motif") } + let zero_reference_targets = reference_records.is_empty(); let mut contigs = reference_records .into_iter() .filter(|rr| { @@ -403,10 +497,18 @@ impl ChromCoordinatesFeeder { contigs.len() ); } - let curr_contig = contigs.pop_front().ok_or(anyhow!( - "should be at least 1 contig, are all of the reads unmapped?" - ))?; - let curr_position = curr_contig.start; + let curr_contig = contigs.pop_front(); + if curr_contig.is_none() + && !(allow_zero_reference_targets && zero_reference_targets) + { + return Err(anyhow!( + "should be at least 1 contig, are all of the reads unmapped?" + )); + } + let curr_position = + curr_contig.as_ref().map(|record| record.start).unwrap_or(0); + let zero_reference_target_terminal = + allow_zero_reference_targets && zero_reference_targets; Ok(Self { contigs, interval_size, @@ -415,18 +517,23 @@ impl ChromCoordinatesFeeder { combine_strands, curr_contig, curr_position, - done: false, + zero_reference_target_terminal, }) } fn update_current(&mut self, end: u32) { - if end >= self.curr_contig.end() { + let curr_contig_end = self + .curr_contig + .as_ref() + .expect("non-terminal feeder must have a current reference record") + .end(); + if end >= curr_contig_end { if let Some(rr) = self.contigs.pop_front() { self.curr_position = rr.start; - self.curr_contig = rr; + self.curr_contig = Some(rr); } else { debug!("feeder is done"); - self.done = true; + self.curr_contig = None; } } else { self.curr_position = end @@ -434,14 +541,14 @@ impl ChromCoordinatesFeeder { } fn get_next(&mut self) -> MkResult> { - if self.done { + let Some(curr_contig) = self.curr_contig.as_ref() else { return Ok(None); - } + }; let start = self.curr_position; - let tid = self.curr_contig.tid; + let tid = curr_contig.tid; let end = std::cmp::min( start.saturating_add(self.interval_size), - self.curr_contig.end(), + curr_contig.end(), ); let (focus_positions, end) = if let Some(lookup) = self.motifs.as_mut() @@ -449,9 +556,9 @@ impl ChromCoordinatesFeeder { // todo change everything to u64 let range = (start as u64)..(end as u64); let (fps, end) = lookup.get_motif_positions( - &self.curr_contig.name, + &curr_contig.name, tid, - self.curr_contig.end(), + curr_contig.end(), range, self.position_filter.as_ref(), self.combine_strands, @@ -504,13 +611,13 @@ impl ChromCoordinatesFeeder { (FocusPositions2::AllPositions, end) }; - let end = std::cmp::min(end, self.curr_contig.end()); + let end = std::cmp::min(end, curr_contig.end()); let chrom_coords = ChromCoordinates::new( tid, start, end, focus_positions, - end == self.curr_contig.end(), + end == curr_contig.end(), ); self.update_current(end); @@ -525,6 +632,10 @@ impl ChromCoordinatesFeeder { self.position_filter.is_some() } + pub(crate) fn is_zero_reference_target_terminal(&self) -> bool { + self.zero_reference_target_terminal + } + pub(crate) fn get_motif_bases(&self) -> Option<[DnaBase; 4]> { if let Some(motif_lookup) = self.ref_motifs() { let motif_primary_bases = motif_lookup @@ -660,8 +771,12 @@ where mod interval_chunks_tests { use rust_htslib::faidx; - use crate::interval_chunks::slice_dna_sequence; + use crate::interval_chunks::{ + slice_dna_sequence, ChromCoordinatesFeeder, + ReferenceIntervalBatchesFeeder, TotalLength, + }; use crate::test_utils::load_test_sequence; + use crate::util::ReferenceRecord; #[test] fn test_check_sequence_slicing_is_same_as_fetch() { @@ -676,6 +791,74 @@ mod interval_chunks_tests { assert_eq!(slice_a, slice_b); } + #[test] + fn strict_feeders_reject_zero_reference_targets() { + assert!(ReferenceIntervalBatchesFeeder::new( + Vec::new(), + 1, + 100, + false, + None, + None, + ) + .is_err()); + assert!( + ChromCoordinatesFeeder::new(&[], 100, None, false, None).is_err() + ); + } + + #[test] + fn interval_batches_zero_reference_terminal_is_exhausted() { + let mut feeder = + ReferenceIntervalBatchesFeeder::new_allowing_zero_reference_targets( + Vec::new(), + 1, + 100, + false, + None, + None, + ) + .unwrap(); + assert_eq!(feeder.total_length(), 0); + assert!(feeder.next().is_none()); + assert!(feeder.next().is_none()); + } + + #[test] + fn chrom_coordinates_zero_reference_terminal_and_clone_are_exhausted() { + let reference_record = + ReferenceRecord::new(0, 0, 1, "chr1".to_string()); + let active = + ChromCoordinatesFeeder::new_allowing_zero_reference_targets( + &[reference_record], + 100, + None, + false, + None, + ) + .unwrap(); + assert!(!active.is_zero_reference_target_terminal()); + + let mut feeder = + ChromCoordinatesFeeder::new_allowing_zero_reference_targets( + &[], + 100, + None, + false, + None, + ) + .unwrap(); + let mut cloned = feeder.clone(); + assert!(feeder.is_zero_reference_target_terminal()); + assert!(cloned.is_zero_reference_target_terminal()); + assert_eq!(feeder.total_length(), 0); + assert_eq!(cloned.total_length(), 0); + assert!(feeder.next().is_none()); + assert!(cloned.next().is_none()); + assert!(feeder.next().is_none()); + assert!(cloned.next().is_none()); + } + #[test] fn test_interval_chunks() { diff --git a/modkit-core/src/modbam_util/subcommands.rs b/modkit-core/src/modbam_util/subcommands.rs index feda6a8a..13cf9f0f 100644 --- a/modkit-core/src/modbam_util/subcommands.rs +++ b/modkit-core/src/modbam_util/subcommands.rs @@ -48,7 +48,7 @@ use crate::sample_probs::{ get_base_mods_quals_from_indexed_hts_file, run_extract_probs_workers, AlignedBaseAndModArgmaxProbs, AlignedBaseArgmaxProbs, BaseAndModArgmaxProbs, BaseArgmaxProbs, ExtractProbsWorker, ProbsExtractor, - QualHist, RegionMleProbs, + QualHist, RegionMleProbs, ZeroReferenceTargetMode, }; use crate::summarize::ModSummary; use crate::util::{ @@ -1945,6 +1945,7 @@ impl ModSummarize { region.as_ref(), edge_filter.as_ref(), &io_threadpool, + ZeroReferenceTargetMode::Reject, multi_progress.clone(), )?; diff --git a/modkit-core/src/sample_probs/mod.rs b/modkit-core/src/sample_probs/mod.rs index 3d255266..726dc100 100644 --- a/modkit-core/src/sample_probs/mod.rs +++ b/modkit-core/src/sample_probs/mod.rs @@ -1511,6 +1511,7 @@ pub(crate) fn calc_per_base_thresholds_from_indexed_hts_file( sampling_region, edge_filter, io_threadpool, + ZeroReferenceTargetMode::FinishMappedPhaseForUnmappedFallback, multi_progress.clone(), )?; if qual_hist.ok_records < 100 { @@ -1540,6 +1541,14 @@ pub(crate) fn calc_per_base_thresholds_from_indexed_hts_file( ) } +#[derive(Clone, Copy)] +pub(crate) enum ZeroReferenceTargetMode { + Reject, + /// Finish the mapped phase immediately so this consumer can run its + /// existing unmapped-record fallback. + FinishMappedPhaseForUnmappedFallback, +} + pub(crate) fn get_base_mods_quals_from_indexed_hts_file( bam_fp: &PathBuf, reference_fasta: Option<&PathBuf>, @@ -1557,6 +1566,7 @@ pub(crate) fn get_base_mods_quals_from_indexed_hts_file( sampling_region: Option<&Region>, edge_filter: Option<&EdgeFilter>, io_threadpool: &rust_htslib::tpool::ThreadPool, + zero_reference_target_mode: ZeroReferenceTargetMode, multi_progress: MultiProgress, ) -> anyhow::Result { let bam_reader = bam::IndexedReader::from_path(bam_fp)?; @@ -1579,13 +1589,27 @@ pub(crate) fn get_base_mods_quals_from_indexed_hts_file( mask, preload_references, )?; - let feeder = ChromCoordinatesFeeder::new( - &reference_records, - interval_size, - motif_lookup, - false, - stranded_position_filter.clone(), - )?; + let feeder = match zero_reference_target_mode { + ZeroReferenceTargetMode::Reject => ChromCoordinatesFeeder::new( + &reference_records, + interval_size, + motif_lookup, + false, + stranded_position_filter.clone(), + ), + ZeroReferenceTargetMode::FinishMappedPhaseForUnmappedFallback => { + ChromCoordinatesFeeder::new_allowing_zero_reference_targets( + &reference_records, + interval_size, + motif_lookup, + false, + stranded_position_filter.clone(), + ) + } + }?; + if feeder.is_zero_reference_target_terminal() { + return Ok(QualHist::default()); + } if let Some(motif_bases) = feeder.get_motif_bases() { for i in 0..n_workers { let worker: Box = if collect_mod_histograms diff --git a/modkit/tests/test_extract_zero_reference_targets.rs b/modkit/tests/test_extract_zero_reference_targets.rs new file mode 100644 index 00000000..21973562 --- /dev/null +++ b/modkit/tests/test_extract_zero_reference_targets.rs @@ -0,0 +1,265 @@ +use rust_htslib::bam::{self, Read}; +use std::collections::HashSet; +use std::fs; +use std::path::{Path, PathBuf}; +use std::process::{Command, Output}; + +const UNMAPPED_BAM: &str = + "../tests/resources/bc_anchored_10_reads.unmapped.bam"; +const UNMAPPED_CRAM: &str = + "../tests/resources/bc_anchored_10_reads_unmapped.cram"; +const UNMAPPED_CRAM_INDEX: &str = + "../tests/resources/bc_anchored_10_reads_unmapped.cram.crai"; +const EMPTY_THRESHOLD_ERROR: &str = + "cannot calculate automatic thresholds because no modification \ + probabilities were sampled"; + +struct InputCopies { + serial_bam: PathBuf, + indexed_bam: PathBuf, + serial_cram: PathBuf, + indexed_cram: PathBuf, + protected_bytes: Vec<(PathBuf, Vec)>, +} + +impl InputCopies { + fn new(temp_dir: &Path) -> Self { + let serial_bam = temp_dir.join("serial.bam"); + let indexed_bam = temp_dir.join("indexed.bam"); + let serial_cram = temp_dir.join("serial.cram"); + let indexed_cram = temp_dir.join("indexed.cram"); + let indexed_cram_index = temp_dir.join("indexed.cram.crai"); + + fs::copy(UNMAPPED_BAM, &serial_bam).unwrap(); + fs::copy(UNMAPPED_BAM, &indexed_bam).unwrap(); + bam::index::build(&indexed_bam, None, bam::index::Type::Bai, 1) + .unwrap(); + let indexed_bam_index = temp_dir.join("indexed.bam.bai"); + assert!(indexed_bam_index.is_file()); + + fs::copy(UNMAPPED_CRAM, &serial_cram).unwrap(); + fs::copy(UNMAPPED_CRAM, &indexed_cram).unwrap(); + fs::copy(UNMAPPED_CRAM_INDEX, &indexed_cram_index).unwrap(); + + assert!(!temp_dir.join("serial.bam.bai").exists()); + assert!(!temp_dir.join("serial.cram.crai").exists()); + + let protected_bytes = [ + &serial_bam, + &indexed_bam, + &indexed_bam_index, + &serial_cram, + &indexed_cram, + &indexed_cram_index, + ] + .into_iter() + .map(|path| (path.clone(), fs::read(path).unwrap())) + .collect(); + + Self { + serial_bam, + indexed_bam, + serial_cram, + indexed_cram, + protected_bytes, + } + } + + fn assert_unchanged(&self) { + for (path, expected) in &self.protected_bytes { + assert_eq!( + fs::read(path).unwrap(), + *expected, + "input or index changed: {}", + path.display() + ); + } + } +} + +fn run_extract( + command: &str, + input: &Path, + output: &Path, + serial: bool, + mode_args: &[&str], +) -> Output { + let mut args = vec![ + "extract".to_string(), + command.to_string(), + input.to_str().unwrap().to_string(), + output.to_str().unwrap().to_string(), + "--threads".to_string(), + "1".to_string(), + "--io-threads".to_string(), + "1".to_string(), + "--no-headers".to_string(), + "--suppress-progress".to_string(), + "--force".to_string(), + ]; + if command == "calls" { + args.extend(["--sample-num-reads".to_string(), "10".to_string()]); + } + if serial { + args.push("--ignore-index".to_string()); + } + args.extend(mode_args.iter().map(|arg| arg.to_string())); + Command::new(env!("CARGO_BIN_EXE_modkit")).args(args).output().unwrap() +} + +fn assert_success(output: &Output) { + assert!( + output.status.success(), + "command failed:\n{}", + String::from_utf8_lossy(&output.stderr) + ); +} + +fn assert_exact_output(output: &Output, path: &Path, expected_rows: usize) { + assert_success(output); + assert!( + String::from_utf8_lossy(&output.stderr).contains("processed 10 reads"), + "unmapped reads were not processed exactly once:\n{}", + String::from_utf8_lossy(&output.stderr) + ); + let bytes = fs::read(path).unwrap(); + assert_eq!( + std::str::from_utf8(&bytes).unwrap().lines().count(), + expected_rows + ); + let read_ids = std::str::from_utf8(&bytes) + .unwrap() + .lines() + .map(|line| line.split('\t').next().unwrap()) + .collect::>(); + assert_eq!(read_ids.len(), 10); +} + +#[test] +fn indexed_zero_target_bam_and_cram_match_serial_extract_modes() { + let temp_dir = tempfile::tempdir().unwrap(); + let inputs = InputCopies::new(temp_dir.path()); + + for (command, mode, mode_args, expected_rows) in [ + ("full", "full", &[][..], 218), + ("calls", "automatic", &[][..], 109), + ("calls", "explicit", &["--filter-threshold", "0"][..], 109), + ("calls", "no-filter", &["--no-filtering"][..], 109), + ] { + let mut expected_bytes = None; + for (format, serial_input, indexed_input) in [ + ("bam", &inputs.serial_bam, &inputs.indexed_bam), + ("cram", &inputs.serial_cram, &inputs.indexed_cram), + ] { + for (route, input, serial) in [ + ("serial", serial_input, true), + ("indexed", indexed_input, false), + ] { + let output_path = temp_dir + .path() + .join(format!("{format}-{mode}-{route}.tsv")); + let output = run_extract( + command, + input, + &output_path, + serial, + mode_args, + ); + assert_exact_output(&output, &output_path, expected_rows); + let observed = fs::read(&output_path).unwrap(); + if let Some(expected) = expected_bytes.as_ref() { + assert_eq!( + &observed, expected, + "{format} {mode} {route} output differed" + ); + } else { + expected_bytes = Some(observed); + } + } + } + } + + inputs.assert_unchanged(); +} + +fn write_without_modification_tags(input: &Path, output: &Path) { + let mut reader = bam::Reader::from_path(input).unwrap(); + let header = bam::Header::from_template(reader.header()); + let mut writer = + bam::Writer::from_path(output, &header, bam::Format::Bam).unwrap(); + for record in reader.records() { + let mut record = record.unwrap(); + for tag in [b"MM", b"ML", b"Mm", b"Ml"] { + if record.aux(tag).is_ok() { + record.remove_aux(tag).unwrap(); + } + } + writer.write(&record).unwrap(); + } +} + +#[test] +fn indexed_zero_target_automatic_threshold_preserves_empty_sample_contract() { + let temp_dir = tempfile::tempdir().unwrap(); + let serial_bam = temp_dir.path().join("no-observations-serial.bam"); + let indexed_bam = temp_dir.path().join("no-observations-indexed.bam"); + write_without_modification_tags(Path::new(UNMAPPED_BAM), &serial_bam); + fs::copy(&serial_bam, &indexed_bam).unwrap(); + bam::index::build(&indexed_bam, None, bam::index::Type::Bai, 1).unwrap(); + let indexed_bam_index = + temp_dir.path().join("no-observations-indexed.bam.bai"); + let protected_bytes = [ + (&serial_bam, fs::read(&serial_bam).unwrap()), + (&indexed_bam, fs::read(&indexed_bam).unwrap()), + (&indexed_bam_index, fs::read(&indexed_bam_index).unwrap()), + ]; + + for (route, input, serial) in [ + ("serial", serial_bam.as_path(), true), + ("indexed", indexed_bam.as_path(), false), + ] { + let output_path = + temp_dir.path().join(format!("automatic-{route}.tsv")); + let output = run_extract("calls", input, &output_path, serial, &[]); + assert!(!output.status.success()); + assert!( + String::from_utf8_lossy(&output.stderr) + .contains(EMPTY_THRESHOLD_ERROR), + "unexpected {route} automatic-threshold error:\n{}", + String::from_utf8_lossy(&output.stderr) + ); + assert!(!output_path.exists()); + } + + for mode_args in [&["--filter-threshold", "0"][..], &["--no-filtering"][..]] + { + let serial_output = temp_dir.path().join(format!( + "control-serial-{}.tsv", + mode_args[0].trim_start_matches('-') + )); + let indexed_output = temp_dir.path().join(format!( + "control-indexed-{}.tsv", + mode_args[0].trim_start_matches('-') + )); + assert_success(&run_extract( + "calls", + &serial_bam, + &serial_output, + true, + mode_args, + )); + assert_success(&run_extract( + "calls", + &indexed_bam, + &indexed_output, + false, + mode_args, + )); + assert_eq!(fs::read(&serial_output).unwrap(), Vec::::new()); + assert_eq!(fs::read(&indexed_output).unwrap(), Vec::::new()); + } + + for (path, expected) in protected_bytes { + assert_eq!(fs::read(path).unwrap(), expected); + } +} From b836f1bd7a04395e0d6b511b6d9a8d82349194f3 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Sun, 9 Aug 2026 20:15:24 -0700 Subject: [PATCH 09/10] Complete CRAM threshold reader coverage --- modkit-core/src/sample_probs/mod.rs | 46 ++++++++++++++++++- modkit/tests/test_cram_reference_consumers.rs | 41 ++++++++++++----- 2 files changed, 74 insertions(+), 13 deletions(-) diff --git a/modkit-core/src/sample_probs/mod.rs b/modkit-core/src/sample_probs/mod.rs index 726dc100..ad35f2e3 100644 --- a/modkit-core/src/sample_probs/mod.rs +++ b/modkit-core/src/sample_probs/mod.rs @@ -39,7 +39,9 @@ use crate::{ }; use crate::{ mod_base_code::DnaBase, - util::{qual_to_prob, reader_is_cram}, + util::{ + qual_to_prob, reader_is_cram, set_reference_for_cram_indexed_reader, + }, }; #[derive(new, Debug)] @@ -1518,6 +1520,7 @@ pub(crate) fn calc_per_base_thresholds_from_indexed_hts_file( multi_progress .suspend(|| info!("collecting probabilities from unmapped reads")); let mut records = bam::IndexedReader::from_path(bam_fp)?; + set_reference_for_cram_indexed_reader(&mut records, reference_fasta)?; records.fetch(FetchDefinition::Unmapped)?; let unmapped_qual_hist = QualHist::from_records( records.records(), @@ -1801,9 +1804,13 @@ fn byte_to_bool_positions(b: u8, agg: &mut [u32; SIZE]) { #[cfg(test)] mod empty_threshold_tests { + use std::path::PathBuf; + use indicatif::MultiProgress; - use super::{ModHist, QualHist}; + use super::{ + calc_per_base_thresholds_from_indexed_hts_file, ModHist, QualHist, + }; use crate::mod_base_code::{DnaBase, METHYL_CYTOSINE}; fn all_modified_qual_hist() -> QualHist { @@ -1842,4 +1849,39 @@ mod empty_threshold_tests { .expect("all-modified observations must produce a threshold"); assert!(thresholds[DnaBase::C as usize] > 0.0); } + + #[test] + fn indexed_unmapped_threshold_reader_applies_supplied_cram_reference() { + let bam_fp = PathBuf::from( + "../tests/resources/bc_anchored_10_reads_unmapped.cram", + ); + let missing_reference = + PathBuf::from("../tests/resources/missing-reference.fa"); + let io_threadpool = rust_htslib::tpool::ThreadPool::new(1).unwrap(); + + let error = calc_per_base_thresholds_from_indexed_hts_file( + &bam_fp, + Some(&missing_reference), + 0.1, + false, + false, + false, + None, + 10, + None, + None, + None, + false, + 1, + 1_000, + Some(7), + None, + None, + &io_threadpool, + MultiProgress::new(), + ) + .expect_err("a supplied CRAM reference must reach the fallback reader"); + + assert!(error.to_string().contains("failed to set CRAM reference")); + } } diff --git a/modkit/tests/test_cram_reference_consumers.rs b/modkit/tests/test_cram_reference_consumers.rs index e03dc7ef..e0901a13 100644 --- a/modkit/tests/test_cram_reference_consumers.rs +++ b/modkit/tests/test_cram_reference_consumers.rs @@ -103,14 +103,10 @@ fn call_and_adjust_mods_match_bam_and_cram() { BAM, call_bam.to_str().unwrap(), "--output-sam", - "--sampling-frac", - "1", - "--seed", - "7", + "--filter-threshold", + "0.5", "--threads", "1", - "--sampling-interval-size", - "20", "--reference", REFERENCE, "--suppress-progress", @@ -120,6 +116,27 @@ fn call_and_adjust_mods_match_bam_and_cram() { CRAM, call_cram.to_str().unwrap(), "--output-sam", + "--filter-threshold", + "0.5", + "--threads", + "1", + "--ref", + REFERENCE, + "--suppress-progress", + ]); + assert_success(&bam_output); + assert_success(&cram_output); + let bam_records = normalized_sam_records(&call_bam); + assert_eq!(bam_records.len(), 10); + assert_eq!(bam_records, normalized_sam_records(&call_cram)); + + let automatic_cram = temp_dir.path().join("call.automatic.cram.sam"); + let automatic_log = temp_dir.path().join("call.automatic.cram.log"); + let automatic_output = run_modkit(&[ + "call-mods", + CRAM, + automatic_cram.to_str().unwrap(), + "--output-sam", "--sampling-frac", "1", "--seed", @@ -130,13 +147,15 @@ fn call_and_adjust_mods_match_bam_and_cram() { "20", "--ref", REFERENCE, + "--log", + automatic_log.to_str().unwrap(), "--suppress-progress", ]); - assert_success(&bam_output); - assert_success(&cram_output); - let bam_records = normalized_sam_records(&call_bam); - assert_eq!(bam_records.len(), 10); - assert_eq!(bam_records, normalized_sam_records(&call_cram)); + assert_success(&automatic_output); + assert_eq!(normalized_sam_records(&automatic_cram).len(), 10); + let automatic_log = fs::read_to_string(automatic_log).unwrap(); + assert!(automatic_log.contains("sampling 100% of reads")); + assert!(automatic_log.contains("estimated pass threshold")); let adjust_bam = temp_dir.path().join("adjust.bam.sam"); let adjust_cram = temp_dir.path().join("adjust.cram.sam"); From 87645355bdf76111381bb0bbee9f3c89b7864996 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Sun, 9 Aug 2026 20:21:25 -0700 Subject: [PATCH 10/10] Isolate automatic CRAM reference test --- modkit/tests/test_cram_reference_consumers.rs | 43 +++++++++++-------- 1 file changed, 24 insertions(+), 19 deletions(-) diff --git a/modkit/tests/test_cram_reference_consumers.rs b/modkit/tests/test_cram_reference_consumers.rs index e0901a13..d2fed607 100644 --- a/modkit/tests/test_cram_reference_consumers.rs +++ b/modkit/tests/test_cram_reference_consumers.rs @@ -132,25 +132,30 @@ fn call_and_adjust_mods_match_bam_and_cram() { let automatic_cram = temp_dir.path().join("call.automatic.cram.sam"); let automatic_log = temp_dir.path().join("call.automatic.cram.log"); - let automatic_output = run_modkit(&[ - "call-mods", - CRAM, - automatic_cram.to_str().unwrap(), - "--output-sam", - "--sampling-frac", - "1", - "--seed", - "7", - "--threads", - "1", - "--sampling-interval-size", - "20", - "--ref", - REFERENCE, - "--log", - automatic_log.to_str().unwrap(), - "--suppress-progress", - ]); + let empty_reference_cache = temp_dir.path().join("empty-reference-cache"); + fs::create_dir(&empty_reference_cache).unwrap(); + let automatic_output = run_modkit_without_reference_resolution( + &[ + "call-mods", + CRAM, + automatic_cram.to_str().unwrap(), + "--output-sam", + "--sampling-frac", + "1", + "--seed", + "7", + "--threads", + "1", + "--sampling-interval-size", + "20", + "--ref", + REFERENCE, + "--log", + automatic_log.to_str().unwrap(), + "--suppress-progress", + ], + &empty_reference_cache, + ); assert_success(&automatic_output); assert_eq!(normalized_sam_records(&automatic_cram).len(), 10); let automatic_log = fs::read_to_string(automatic_log).unwrap();