diff --git a/src/lib/commands/filter.rs b/src/lib/commands/filter.rs index ba3bef0fc..b0eb331d3 100644 --- a/src/lib/commands/filter.rs +++ b/src/lib/commands/filter.rs @@ -17,38 +17,26 @@ use crate::consensus_filter::{ filter_read, is_duplex_consensus, mask_bases, mask_duplex_bases, mask_methylation_depth_duplex_raw_with_tags, mask_methylation_depth_simplex_raw_with_tags, mask_strand_methylation_agreement_raw_with_ref_bases_and_tags, mean_base_quality_full_length, - retained_primary_masked_bases, template_passes, }; -use crate::grouper::{SingleRawRecordGrouper, TemplateGrouper}; -use crate::logging::OperationTimer; use crate::per_thread_accumulator::PerThreadAccumulator; -use crate::read_info::LibraryIndex; use crate::reference::ReferenceReader; use crate::tag_reversal::reverse_per_base_tags_raw; -use crate::template::TemplateBatch; -use crate::unified_pipeline::{ - BamPipelineConfig, BatchWeight, GroupKeyConfig, Grouper, MemoryEstimate, - run_bam_pipeline_from_reader, run_bam_pipeline_from_reader_with_secondary, -}; use crate::validation::validate_file_exists; -use ahash::AHashMap; use anyhow::{Result, bail}; use clap::Parser; -use fgumi_bam_io::create_bam_reader_for_pipeline_with_opts; use fgumi_raw_bam; use fgumi_raw_bam::{RawRecord, RawRecordView}; use log::info; use noodles::sam::Header; -use std::io; use std::path::{Path, PathBuf}; use std::sync::Arc; -use std::sync::atomic::{AtomicU64, Ordering}; +use std::sync::atomic::AtomicU64; use std::time::Instant; use crate::commands::command::Command; use crate::commands::common::{ BamIoOptions, CompressionOptions, QueueMemoryOptions, SchedulerOptions, ThreadingOptions, - build_pipeline_config, reject_output_collisions, serialize_raw_bam_records, + reject_output_collisions, }; /// Filters and masks consensus reads based on various quality metrics. @@ -356,34 +344,9 @@ impl Filter { } // ============================================================================ -// 7-Step Pipeline Types +// Pipeline Types // ============================================================================ -/// Result from processing a batch of records through raw-byte filtering. -struct FilterProcessedBatchRaw { - /// Records that passed filtering (raw bytes). - kept_records: Vec, - /// Records that failed filtering (raw bytes, if tracking rejects). - rejected_records: Vec, - /// Number of records processed. - records_count: u64, - /// Number of records that passed. - passed_count: u64, - /// Number of bases masked. - bases_masked: u64, -} - -impl MemoryEstimate for FilterProcessedBatchRaw { - fn estimate_heap_size(&self) -> usize { - let vec_overhead = std::mem::size_of::(); - let kept_outer = self.kept_records.capacity() * vec_overhead; - let kept_inner: usize = self.kept_records.iter().map(RawRecord::capacity).sum(); - let rejected_outer = self.rejected_records.capacity() * vec_overhead; - let rejected_inner: usize = self.rejected_records.iter().map(RawRecord::capacity).sum(); - kept_outer + kept_inner + rejected_outer + rejected_inner - } -} - /// Per-thread accumulator merged into final counts after pipeline completion. #[derive(Default)] pub(crate) struct CollectedFilterMetrics { @@ -444,100 +407,38 @@ impl Command for Filter { // Validate parameter counts (1-3 values for duplex support) self.validate_parameters()?; - if self.threading.threads.is_some() { - return self.execute_chain(command_line); - } - - let timer = OperationTimer::new("Filtering consensus reads"); - - info!("Starting Filter"); - info!("Input: {}", self.io.input.display()); - info!("Output: {}", self.io.output.display()); - match &self.reference { - Some(r) => info!("Reference: {}", r.display()), - None => info!("Reference: (tag regeneration disabled)"), - } - info!("Min reads: {:?}", self.min_reads); - info!("Max read error rate: {:?}", self.max_read_error_rate); - info!("Max base error rate: {:?}", self.max_base_error_rate); - if let Some(q) = self.min_base_quality { - info!("Min base quality: {q}"); - } - if let Some(q) = self.min_mean_base_quality { - info!("Min mean base quality: {q}"); - } - info!("Max no-call fraction: {}", self.max_no_call_fraction); - if !self.min_methylation_depth.is_empty() { - info!("Min methylation depth: {:?}", self.min_methylation_depth); - } - if self.require_strand_methylation_agreement { - info!("Require strand methylation agreement: true"); - } - if let Some(frac) = self.min_conversion_fraction { - info!("Min conversion fraction: {frac}"); - } - if let Some(mode) = &self.methylation_mode { - info!("Methylation mode: {mode:?}"); - } - - self.io.log_effective_check_crc(); - - // Open input using streaming-capable reader for pipeline use - let (reader, header) = create_bam_reader_for_pipeline_with_opts( - &self.io.input, - self.io.pipeline_reader_opts(), - )?; - - // Synthesize @HD VN:1.6 SO:unsorted when the input lacks one (match fgbio, - // which never passes a header-less BAM straight through). - let header = crate::commands::common::ensure_hd_record(header)?; - // FILT3-02: fgbio's FilterConsensusReads calls Bams.requireQueryGrouped. - // Filtering is template-based, so coordinate-sorted input silently scatters - // mates and corrupts the both-primaries-pass logic. Reject it like fgbio. - crate::commands::common::require_query_grouped( - &header, - &self.io.input.display().to_string(), - )?; - - // Add @PG record with PP chaining to input's last program - let header = crate::commands::common::add_pg_record(header, command_line)?; - - let track_rejects = self.rejects.is_some(); - let threads = self.threading.threads.unwrap_or(1); - - // Route to appropriate 7-step pipeline mode - let total_reads = if self.filter_by_template { - self.execute_threads_mode_template(threads, reader, header, track_rejects)? - } else { - self.execute_threads_mode_single_read(threads, reader, header, track_rejects)? - }; - - timer.log_completion(total_reads); - Ok(()) + // The declarative chain is the only execution path. `execute_chain` emits + // the CRC-verify log, and `ChainBuilder::add_filter` (plus its finalize + // hooks) emits the `Starting Filter` banner + `OperationTimer`, the + // query-grouped input check, `@PG` injection, the summary banner, and the + // stats-file/rejects hooks (running any of those here first would + // double-log and pre-consume stdin), so `execute` does only the pre-flight + // validation above and then dispatches. Absent `--threads` runs the chain + // at a single worker. + self.execute_chain(command_line) } } impl Filter { - /// Run the filter stage on the declarative chain builder (the `--threads N` - /// path). + /// Run the filter stage on the declarative chain builder. /// - /// Replaces the hand-rolled unified-pipeline construction in `execute` for - /// the threaded case. The chain opens its own source, validates the - /// query-grouped input ordering, injects `@PG`, runs the shared filter - /// process step(s), writes the kept records (and, when configured, the - /// rejects BAM and stats file) via its finalize hooks — all through the - /// same shared helpers as the non-chain path, so the two orchestrations - /// stay in parity. The no-`--threads` path keeps its own unified pipeline - /// in `execute`, which is the in-process parity oracle for this one (see - /// `test_filter_chain_matches_single_threaded`). + /// The chain is the only execution path: `execute` always dispatches here, + /// with or without `--threads` (absent `--threads` runs the chain at a single + /// worker). The chain opens its own source, validates the query-grouped input + /// ordering, injects `@PG`, runs the shared filter process step(s), and writes + /// the kept records (and, when configured, the rejects BAM and stats file) via + /// its finalize hooks. All user-facing diagnostics — the `Starting Filter` + /// banner, the `OperationTimer`, the parameter log lines, and the + /// summary/stats finalize hooks — are emitted inside + /// `ChainBuilder::add_filter`, so `execute_chain` itself only surfaces the CRC + /// policy and runs the pipeline (mirrors `dedup`/`retag` `execute_chain`). fn execute_chain(&self, command_line: &str) -> Result<()> { use crate::pipeline::chains::{ ChainSpec, SingleStageContext, Stage, StageOptionsBag, build_for, }; // `add_filter` re-emits the timer/banner/threading log lines but not the - // CRC-verify status; emit it here so the --threads path reports it once, - // matching the non-chain path. + // CRC-verify status; emit it here so the run reports it exactly once. self.io.log_effective_check_crc(); let stage_opts = @@ -557,12 +458,10 @@ impl Filter { impl FilterOptions { /// Build the shared filter config, reference, per-thread metrics - /// accumulators, and progress counter used by both the chain builder and - /// the non-chain unified pipeline. + /// accumulators, and progress counter used by the chain builder. /// - /// The [`BamPipelineConfig`] is *not* built here: the chain builder derives - /// it from its [`crate::pipeline::chains::ChainSpec`], and the non-chain run - /// builds it via [`Filter::build_filter_pipeline_config`]. + /// The `BamPipelineConfig` is *not* built here: the chain builder derives it + /// from its [`crate::pipeline::chains::ChainSpec`]. pub(crate) fn setup_pipeline(&self, num_threads: usize) -> Result { let config = Arc::new(FilterConfig::new( &self.min_reads, @@ -770,356 +669,6 @@ impl FilterOptions { } impl Filter { - // ======================================================================== - // 7-Step Unified Pipeline Implementation - // ======================================================================== - - /// Build the pipeline configuration for the non-chain (`unified_pipeline`) - /// filter run. The chain builder assembles its own [`BamPipelineConfig`] - /// from the [`crate::pipeline::chains::ChainSpec`], so this lives on - /// [`Filter`] (which carries the io/compression/scheduler/queue flags) - /// rather than on [`FilterOptions`] (tuning knobs only). - fn build_filter_pipeline_config( - &self, - num_threads: usize, - header: &Header, - ) -> Result { - let mut pipeline_config = build_pipeline_config( - &self.scheduler_opts, - &self.compression, - &self.queue_memory, - &self.io, - num_threads, - )?; - // Template mode groups by QNAME and reads the key's `name_hash` fast-path - // (see `TemplateGrouper`), so it needs the key. Single-read mode uses - // `SingleRawRecordGrouper`, which discards the decoded record entirely and - // never reads the key — so skip its per-record computation with `None`. - pipeline_config.group_key_config = if self.filter_by_template { - Some(GroupKeyConfig::new_raw_no_cell(LibraryIndex::from_header(header))) - } else { - None - }; - Ok(pipeline_config) - } - - /// Build the shared filter config, reference, metrics queue, and progress - /// counter. Delegates to [`FilterOptions::setup_pipeline`] so the chain - /// builder and the non-chain run share one implementation. - pub(crate) fn setup_pipeline(&self, num_threads: usize) -> Result { - self.to_filter_options().setup_pipeline(num_threads) - } - - /// Build the process closure captures. Delegates to - /// [`FilterOptions::process_captures`]. - pub(crate) fn process_captures( - &self, - setup: &FilterPipelineSetup, - header: &Header, - ) -> FilterProcessCaptures { - self.to_filter_options().process_captures(setup, header) - } - - /// Run the filter pipeline with the given grouper and process function. - /// - /// This is the common pipeline executor shared by single-read and template modes. - /// It builds the serialize function, runs the pipeline (with optional secondary - /// output for rejects), and aggregates metrics. - fn run_filter_pipeline( - &self, - pipeline_config: BamPipelineConfig, - setup: FilterPipelineSetup, - reader: Box, - header: Header, - grouper_fn: GrouperFn, - process_fn: ProcessFn, - ) -> Result - where - G: Send + BatchWeight + MemoryEstimate + 'static, - GrouperFn: FnOnce(&Header) -> Box + Send>, - ProcessFn: Fn(G) -> io::Result + Send + Sync + 'static, - { - let collected_for_serialize = Arc::clone(&setup.collected_metrics); - let progress = Arc::clone(&setup.progress_counter); - - // Primary serialize: write kept records - let serialize_fn = move |processed: FilterProcessedBatchRaw, - _header: &Header, - output: &mut Vec| - -> io::Result { - collected_for_serialize.with_slot(|m| { - m.total_records += processed.records_count; - m.passed_records += processed.passed_count; - m.failed_records += processed.records_count - processed.passed_count; - m.total_bases_masked += processed.bases_masked; - }); - - // Interim "Processed N records" heartbeat, advanced once per batch - // here rather than once per record inside the parallel `process_fn` — - // the per-record `fetch_add` was a contended shared atomic (cacheline - // ping-pong across pipeline workers). Log-only; the authoritative total - // is aggregated from `records_count` below. - advance_progress(&progress, processed.records_count); - - serialize_raw_bam_records(&processed.kept_records, output) - }; - - if let Some(rejects_path) = &self.rejects { - // Secondary serialize: write rejected records - let secondary_serialize_fn = - |batch: &FilterProcessedBatchRaw, buf: &mut Vec| -> io::Result { - serialize_raw_bam_records(&batch.rejected_records, buf) - }; - - run_bam_pipeline_from_reader_with_secondary( - pipeline_config, - reader, - header, - &self.io.output, - None, - rejects_path, - None, // secondary uses resolved primary output header - grouper_fn, - process_fn, - serialize_fn, - secondary_serialize_fn, - )?; - } else { - run_bam_pipeline_from_reader( - pipeline_config, - reader, - header, - &self.io.output, - None, - grouper_fn, - process_fn, - serialize_fn, - )?; - } - - // Aggregate metrics - let mut total_reads = 0u64; - let mut passed_reads = 0u64; - let mut failed_reads = 0u64; - let mut total_bases_masked = 0u64; - - for slot in setup.collected_metrics.slots() { - let m = slot.lock(); - total_reads += m.total_records; - passed_reads += m.passed_records; - failed_reads += m.failed_records; - total_bases_masked += m.total_bases_masked; - } - - if let Some(stats_path) = &self.stats { - self.write_filter_stats(stats_path, total_reads, passed_reads, failed_reads)?; - } - - info!("Processed {total_reads} reads; kept {passed_reads} and rejected {failed_reads}"); - if self.rejects.is_some() && failed_reads > 0 { - info!("Wrote {failed_reads} rejected records to rejects file"); - } - info!("Total bases masked: {total_bases_masked}"); - - Ok(total_reads) - } - - /// Execute using the 7-step unified pipeline (single-read mode, raw bytes). - /// - /// Each record is filtered independently without template awareness. - fn execute_threads_mode_single_read( - &self, - num_threads: usize, - reader: Box, - header: Header, - track_rejects: bool, - ) -> Result { - let setup = self.setup_pipeline(num_threads)?; - let pipeline_config = self.build_filter_pipeline_config(num_threads, &header)?; - let ctx = self.process_captures(&setup, &header); - - let grouper_fn = move |_header: &Header| { - Box::new(SingleRawRecordGrouper::new()) as Box + Send> - }; - - let process_fn = move |mut record: RawRecord| -> io::Result { - let mut kept_records: Vec = Vec::new(); - let mut rejected_records: Vec = Vec::new(); - let mut passed_count = 0u64; - - let (masked, pass) = Self::process_record_raw( - &mut record, - &ctx.config, - ctx.reference.as_deref(), - &ctx.header, - ctx.should_reverse_tags, - ctx.min_base_quality, - ctx.require_single_strand_agreement, - ctx.min_mean_base_quality, - ctx.max_no_call_fraction, - ctx.methylation_depth_thresholds.as_ref(), - ctx.require_strand_methylation_agreement, - ctx.min_conversion_fraction, - ctx.methylation_mode, - &ctx.ref_names, - ) - .map_err(io::Error::other)?; - - // Match fgbio's `maskedBases`: count a record's masked bases only when it is a - // retained primary read (FilterConsensusReads.scala:207-219). In this per-record - // streaming mode the record is its own single-read "template", so `pass` is the - // template result; secondary/supplementary reads still contribute nothing. - let bases_masked = - retained_primary_masked_bases(std::slice::from_ref(&record), &[masked], pass); - - if pass { - passed_count = 1; - kept_records.push(record); - } else if track_rejects { - rejected_records.push(record); - } - - Ok(FilterProcessedBatchRaw { - kept_records, - rejected_records, - records_count: 1, - passed_count, - bases_masked, - }) - }; - - self.run_filter_pipeline(pipeline_config, setup, reader, header, grouper_fn, process_fn) - } - - /// Execute using the 7-step unified pipeline (template-aware mode). - /// - /// All primary reads in a template must pass for the template to pass. - fn execute_threads_mode_template( - &self, - num_threads: usize, - reader: Box, - header: Header, - track_rejects: bool, - ) -> Result { - let setup = self.setup_pipeline(num_threads)?; - let pipeline_config = self.build_filter_pipeline_config(num_threads, &header)?; - let ctx = self.process_captures(&setup, &header); - - #[cfg(test)] - const BATCH_SIZE: usize = 50; - #[cfg(not(test))] - const BATCH_SIZE: usize = 1000; - - let grouper_fn = move |_header: &Header| { - Box::new(TemplateGrouper::new(BATCH_SIZE)) - as Box + Send> - }; - - let process_fn = move |batch: TemplateBatch| -> io::Result { - let mut kept_records: Vec = Vec::new(); - let mut rejected_records: Vec = Vec::new(); - let mut total_records = 0u64; - let mut passed_count = 0u64; - let mut bases_masked = 0u64; - - for template in batch { - let mut template_records: Vec = template.into_records(); - let mut pass_map: AHashMap = - AHashMap::with_hasher(crate::hashing::deterministic_state()); - let mut masked_by_record: Vec = Vec::with_capacity(template_records.len()); - - for (idx, record) in template_records.iter_mut().enumerate() { - total_records += 1; - - let (masked, pass) = Self::process_record_raw( - record, - &ctx.config, - ctx.reference.as_deref(), - &ctx.header, - ctx.should_reverse_tags, - ctx.min_base_quality, - ctx.require_single_strand_agreement, - ctx.min_mean_base_quality, - ctx.max_no_call_fraction, - ctx.methylation_depth_thresholds.as_ref(), - ctx.require_strand_methylation_agreement, - ctx.min_conversion_fraction, - ctx.methylation_mode, - &ctx.ref_names, - ) - .map_err(io::Error::other)?; - masked_by_record.push(masked); - pass_map.insert(idx, pass); - } - - let template_pass = template_passes(&template_records, &pass_map); - - // fgbio tallies masked bases only over the primary reads of retained - // templates (FilterConsensusReads.scala:207-219); a dropped template and - // any secondary/supplementary read contribute nothing. - bases_masked += retained_primary_masked_bases( - &template_records, - &masked_by_record, - template_pass, - ); - - for (idx, record) in template_records.into_iter().enumerate() { - let flags = RawRecordView::new(&record).flags(); - let is_primary = (flags & fgumi_raw_bam::flags::SECONDARY) == 0 - && (flags & fgumi_raw_bam::flags::SUPPLEMENTARY) == 0; - - if is_primary { - if template_pass { - passed_count += 1; - kept_records.push(record); - } else if track_rejects { - rejected_records.push(record); - } - } else { - let record_pass = pass_map.get(&idx).copied().unwrap_or(false); - if template_pass && record_pass { - passed_count += 1; - kept_records.push(record); - } else if track_rejects { - rejected_records.push(record); - } - } - } - } - - Ok(FilterProcessedBatchRaw { - kept_records, - rejected_records, - records_count: total_records, - passed_count, - bases_masked, - }) - }; - - self.run_filter_pipeline(pipeline_config, setup, reader, header, grouper_fn, process_fn) - } - - /// Write filtering statistics to a file. - fn write_filter_stats( - &self, - path: &std::path::Path, - total: u64, - passed: u64, - failed: u64, - ) -> Result<()> { - use std::fs::File; - use std::io::Write; - - let mut file = File::create(path)?; - writeln!(file, "total_reads\t{total}")?; - writeln!(file, "passed_reads\t{passed}")?; - writeln!(file, "failed_reads\t{failed}")?; - #[allow(clippy::cast_precision_loss)] - let pass_rate = if total > 0 { passed as f64 / total as f64 } else { 0.0 }; - writeln!(file, "pass_rate\t{pass_rate:.4}")?; - Ok(()) - } - /// Process a single raw BAM record: reverse tags, mask bases, regenerate alignment /// tags, and check filters. /// @@ -1381,40 +930,12 @@ impl Filter { /// Validates that parameter vectors have 1-3 values and are in valid ranges. /// - /// Delegates to [`FilterOptions::validate_parameters`] so the chain builder - /// and the non-chain run share one implementation. + /// Delegates to [`FilterOptions::validate_parameters`] so `execute`'s + /// pre-flight and the chain builder share one implementation. pub(crate) fn validate_parameters(&self) -> Result<()> { self.to_filter_options().validate_parameters() } } - -/// Advances the shared interim-progress `counter` by `records` and logs a -/// "Processed N records" heartbeat when the running total crosses a -/// 1,000,000-record boundary. -/// -/// This is log-only — the authoritative record total is aggregated from each -/// batch's `records_count` and logged at completion — so callers may batch the -/// update coarsely to keep the shared atomic off the per-record hot path. -fn advance_progress(counter: &AtomicU64, records: u64) { - let before = counter.fetch_add(records, Ordering::Relaxed); - if let Some(total) = progress_heartbeat_total(before, records) { - info!("Processed {total} records"); - } -} - -/// Returns the running total to announce in a "Processed N records" heartbeat -/// when advancing the interim-progress counter from `before` by `records`, or -/// `None` when the advance does not cross a 1,000,000-record boundary. -/// -/// This is the pure boundary-crossing decision behind [`advance_progress`]'s -/// heartbeat: a non-`None` result means a heartbeat should be emitted for that -/// exact running total. Extracting it keeps the emit-or-not behavior directly -/// testable without capturing log output. -fn progress_heartbeat_total(before: u64, records: u64) -> Option { - let after = before + records; - if after / 1_000_000 > before / 1_000_000 { Some(after) } else { None } -} - #[cfg(test)] #[allow(clippy::float_cmp)] mod tests { @@ -1635,44 +1156,6 @@ mod tests { use noodles::sam::alignment::record_buf::RecordBuf; use rstest::rstest; - /// `advance_progress` accumulates exactly across successive batched advances, - /// regardless of whether any given advance crosses a heartbeat boundary. The - /// counter is the authoritative side of the function; the heartbeat cadence - /// is verified separately by [`test_progress_heartbeat_total`]. - #[test] - fn test_advance_progress_accumulates_and_crosses_boundary() { - let counter = AtomicU64::new(0); - advance_progress(&counter, 1_000_000); // 0 -> 1_000_000: crosses, logs - assert_eq!(counter.load(Ordering::Relaxed), 1_000_000); - advance_progress(&counter, 5); // 1_000_000 -> 1_000_005: no crossing - assert_eq!(counter.load(Ordering::Relaxed), 1_000_005); - advance_progress(&counter, 1_000_000); // crosses again - assert_eq!(counter.load(Ordering::Relaxed), 2_000_005); - } - - /// The heartbeat fires exactly when an advance crosses a 1,000,000-record - /// boundary, and reports the *running total* (not the batch size). Each case - /// mirrors a step of the accumulation walk plus the boundary edge cases: - /// a first crossing, a within-window increment that must stay silent, a - /// repeated crossing at the next boundary, a no-op advance, an advance that - /// lands exactly on a boundary, and a single advance that spans several - /// boundaries (still one heartbeat, at the final total). - #[rstest] - #[case::first_crossing(0, 1_000_000, Some(1_000_000))] - #[case::within_window_is_silent(1_000_000, 5, None)] - #[case::repeated_crossing(1_000_005, 999_995, Some(2_000_000))] - #[case::zero_advance_is_silent(1_000_005, 0, None)] - #[case::partial_within_window_is_silent(0, 999_999, None)] - #[case::lands_exactly_on_boundary(999_999, 1, Some(1_000_000))] - #[case::spans_multiple_boundaries(500_000, 2_500_000, Some(3_000_000))] - fn test_progress_heartbeat_total( - #[case] before: u64, - #[case] records: u64, - #[case] expected: Option, - ) { - assert_eq!(progress_heartbeat_total(before, records), expected); - } - /// Helper function to create a Filter command with commonly used test defaults. fn create_filter_with_paths(input: PathBuf, output: PathBuf, reference: PathBuf) -> Filter { Filter { @@ -3792,13 +3275,11 @@ mod tests { Ok(()) } - /// Non-template filtering through the *multi-threaded* pipeline. With - /// `filter_by_template: false`, `build_filter_pipeline_config` leaves - /// `group_key_config = None`, so the threaded run must route records through - /// `SingleRawRecordGrouper` — which discards the decoded key entirely — and - /// filter each read independently. Every other non-template execution test - /// runs single-threaded, so this pins the `> 1`-thread path over that `None` - /// branch. + /// Non-template filtering through the *multi-threaded* chain. With + /// `filter_by_template: false`, the chain's single-read filter step filters + /// each record independently (no queryname grouping). Every other + /// non-template execution test runs the chain at a single worker, so this + /// pins the `> 1`-thread configuration of that single-read step. #[test] fn test_filter_execute_non_template_mode_multithreaded() -> Result<()> { let dir = TempDir::new()?; diff --git a/src/lib/pipeline/chains/builder.rs b/src/lib/pipeline/chains/builder.rs index fe22a5ad4..4f0d5a123 100644 --- a/src/lib/pipeline/chains/builder.rs +++ b/src/lib/pipeline/chains/builder.rs @@ -1140,10 +1140,10 @@ impl<'a> ChainBuilder<'a> { /// pass entirely; for a SAM-first sort that pass runs once per record in /// `ParseSamChunk` and would otherwise be pure waste. This changes only the /// discarded key, never the record bytes, so output is unchanged. - /// (The legacy single-threaded path uses `new_raw_no_cell`, which still - /// pays the combined aux pass; name-hash-only is strictly less work.) All - /// other first stages (group/dedup/consensus) need the full position/cell - /// key, so they fall through to [`Self::bam_group_key_config`]. + /// (The full-key config `new_raw_no_cell` still pays the combined aux pass; + /// name-hash-only is strictly less work.) All other first stages + /// (group/dedup/consensus) need the full position/cell key, so they fall + /// through to [`Self::bam_group_key_config`]. /// /// Every arm below uses a DEFAULT `LibraryIndex`, never `from_header`: /// `name_hash_only` never reads `library_index` (see `name_hash_key` in @@ -1152,11 +1152,7 @@ impl<'a> ChainBuilder<'a> { /// `library_index` at all), so resolving it from the header is pure waste /// for these stages — and `LibraryIndex::from_header` panics on a header /// with more than 65,535 distinct `@RG` libraries, a needless crash risk - /// this avoids. (Some of these stages' serial oracles independently build - /// a full-header group key of their own — e.g. filter-by-template's - /// `build_filter_pipeline_config` calls `LibraryIndex::from_header` — so - /// this isn't "a panic the oracle never has"; it is simply dead work this - /// arm has no reason to repeat.) + /// this avoids. fn source_group_key_config(&self) -> fgumi_bam_io::GroupKeyConfig { match self.spec.stages.first() { Some( @@ -4437,11 +4433,10 @@ impl<'a> ChainBuilder<'a> { // FILT3-02: filtering is template-based, so coordinate-sorted input // silently scatters mates and corrupts the both-primaries-pass logic. - // Reject it here exactly as `Filter::execute`'s legacy path does, so - // the two orchestrations of the filter stage cannot drift on accepted - // orders, error text, or logging (mirrors `add_group`'s - // `require_group_input_ordering` call, shared verbatim with - // `Group::execute`). + // Reject it on the chain — the only filter execution path — so filter + // rejects mis-ordered input before any record is processed (mirrors + // `add_group`'s `require_group_input_ordering` call, shared verbatim + // with `Group::execute`). crate::commands::common::require_query_grouped( &self.header, &input_path.display().to_string(), diff --git a/src/lib/pipeline/chains/commands/filter.rs b/src/lib/pipeline/chains/commands/filter.rs index 78acc95f4..70687684e 100644 --- a/src/lib/pipeline/chains/commands/filter.rs +++ b/src/lib/pipeline/chains/commands/filter.rs @@ -225,7 +225,7 @@ pub(crate) fn build_filter_step_single_no_rejects( let mut record = decoded.into_raw_bytes(); let (bases_masked, pass) = process_record_raw_call(&mut record, &captures).map_err(io::Error::other)?; - // Match the legacy path's fgbio-parity "Total bases masked" tally: + // Match fgbio's "Total bases masked" tally: // count masked bases only in a retained primary read (0 for a // rejected read / secondary / supplementary), not the raw count. bases_masked_total += retained_primary_masked_bases( @@ -302,7 +302,7 @@ pub(crate) fn build_filter_step_single_with_rejects( let mut record = decoded.into_raw_bytes(); let (bases_masked, pass) = process_record_raw_call(&mut record, &captures) .map_err(io::Error::other)?; - // Match the legacy path's fgbio-parity "Total bases masked" tally: + // Match fgbio's "Total bases masked" tally: // count masked bases only in a retained primary read (0 for a // rejected read / secondary / supplementary), not the raw count. bases_masked_total += @@ -375,7 +375,7 @@ pub(crate) fn build_filter_step_template_no_rejects( } let template_pass = template_passes(&template_records, &pass_map); - // Match the legacy path's fgbio-parity "Total bases masked" tally: + // Match fgbio's "Total bases masked" tally: // count masked bases only in retained primary reads of a retained // template (0 for a dropped template), not the raw per-record sum. bases_masked_total += retained_primary_masked_bases( @@ -466,7 +466,7 @@ pub(crate) fn build_filter_step_template_with_rejects( } let template_pass = template_passes(&template_records, &pass_map); - // Match the legacy path's fgbio-parity "Total bases masked" tally: + // Match fgbio's "Total bases masked" tally: // count masked bases only in retained primary reads of a retained // template (0 for a dropped template), not the raw per-record sum. bases_masked_total += retained_primary_masked_bases( diff --git a/src/lib/pipeline/steps/group/queryname.rs b/src/lib/pipeline/steps/group/queryname.rs index f52e21ac5..160d050d5 100644 --- a/src/lib/pipeline/steps/group/queryname.rs +++ b/src/lib/pipeline/steps/group/queryname.rs @@ -25,9 +25,8 @@ use crate::template::Template; /// Serial mutex acquisition; mirrors `GroupByMi::MAX_BATCHES_PER_LOCK`. const MAX_BATCHES_PER_LOCK: usize = 8; -/// Default target batch count. Mirrors legacy's -/// `TemplateGrouper::new(1000)` (the production batch size in -/// `Filter::execute_threads_mode_template`). +/// Default target batch count: the production batch size used when filter +/// groups templates on the chain (`fgumi filter --filter-by-template=true`). pub const DEFAULT_TARGET_BATCH_COUNT: usize = 1000; /// `Serial + ByItemOrdinal` queryname grouper. Records arriving in diff --git a/src/lib/pipeline/steps/parse/decode.rs b/src/lib/pipeline/steps/parse/decode.rs index 2d3c34d2f..a50ccd6bb 100644 --- a/src/lib/pipeline/steps/parse/decode.rs +++ b/src/lib/pipeline/steps/parse/decode.rs @@ -441,6 +441,103 @@ mod tests { } } + /// Guards the `name_hash_only` decode skip that the chain filter path relies + /// on (`source_group_key_config` routes `Stage::Filter` — and `Correct` — to + /// `GroupKeyConfig::name_hash_only`). filter groups templates through + /// `GroupByQueryname`, which reads only `key.name_hash`, and its process step + /// operates on the raw record bytes (untouched by the key config), so the skip + /// is output-safe **iff** `name_hash_only` produces the same `name_hash` the + /// full key would for every record — even when the fields the full key also + /// computes (library index from `RG`, unclipped position) vary across records. + /// + /// This is the in-repo replacement for the parity that + /// `test_filter_chain_matches_single_threaded_with_rg_and_cb_variation` used + /// to provide before the cutover made both of its runs take the same + /// `name_hash_only` chain: it pins the invariant directly against the two + /// decode-consumer key branches (`compute_group_key_from_raw` vs + /// `name_hash_key`) without needing `FGUMI_BASELINE_BIN`. + #[test] + fn name_hash_only_matches_full_key_grouping_when_rg_and_position_vary() { + use crate::sam::SamTag; + use fgumi_bam_io::{GroupKey, LibraryIndex}; + use fgumi_raw_bam::SamBuilder; + use fgumi_raw_bam::flags::{FIRST_SEGMENT, LAST_SEGMENT, PAIRED}; + use noodles::sam::alignment::record::data::field::Tag; + use noodles::sam::header::record::value::Map; + use noodles::sam::header::record::value::map::ReadGroup; + use noodles::sam::header::record::value::map::read_group::tag as rg_tag; + + // Two read groups in distinct libraries, so the full key's `library_idx` + // genuinely differs by RG — the field name_hash_only skips computing. + let mut header = noodles::sam::Header::builder(); + for (id, library) in [("RG1", "libA"), ("RG2", "libB")] { + let rg = Map::::builder() + .insert(rg_tag::LIBRARY, String::from(library)) + .build() + .expect("read group builds"); + header = header.add_read_group(bstr::BString::from(id), rg); + } + let lib = LibraryIndex::from_header(&header.build()); + let cb = Some(Tag::from([b'C', b'B'])); + + // Two paired templates, each R1+R2 sharing a name; the templates carry + // DIFFERENT read groups and DIFFERENT mapped positions, so the full key's + // library_idx and pos1 differ across them (the test is non-vacuous only if + // the skipped fields actually vary). + let mate = |name: &[u8], rg: &[u8], pos: i32, first: bool| -> fgumi_raw_bam::RawRecord { + let mut b = SamBuilder::new(); + b.read_name(name) + .flags(PAIRED | if first { FIRST_SEGMENT } else { LAST_SEGMENT }) + .ref_id(0) + .pos(pos) + .mapq(60) + .cigar_ops(&[4u32 << 4]) // 4M + .sequence(b"ACGT") + .qualities(&[30u8; 4]); + b.add_string_tag(SamTag::RG, rg); + b.build() + }; + let records = [ + mate(b"tmpl-A", b"RG1", 100, true), + mate(b"tmpl-A", b"RG1", 200, false), + mate(b"tmpl-B", b"RG2", 300, true), + mate(b"tmpl-B", b"RG2", 400, false), + ]; + + let full: Vec<_> = + records.iter().map(|r| compute_group_key_from_raw(r.as_ref(), &lib, cb)).collect(); + let name_only: Vec<_> = records.iter().map(|r| name_hash_key(r.as_ref())).collect(); + + for (i, (full_key, name_key)) in full.iter().zip(&name_only).enumerate() { + assert_eq!( + name_key.name_hash, full_key.name_hash, + "record {i}: name_hash_only must reproduce the full key's name_hash despite \ + differing RG/position, or filter's GroupByQueryname would group differently" + ); + // Everything except name_hash is left at the config-independent default, + // so the RG/position the full key computes cannot leak into the grouping + // key the chain actually uses. + assert_eq!( + *name_key, + GroupKey { name_hash: full_key.name_hash, ..GroupKey::default() }, + "record {i}: name_hash_only must leave all non-name fields at default" + ); + } + + // Mates of a template share a name_hash under BOTH configs, so template + // membership (and thus filter's both-primaries aggregation) is identical. + assert_eq!(name_only[0].name_hash, name_only[1].name_hash, "tmpl-A mates group together"); + assert_eq!(name_only[2].name_hash, name_only[3].name_hash, "tmpl-B mates group together"); + assert_ne!(name_only[0].name_hash, name_only[2].name_hash, "distinct templates stay apart"); + + // Non-vacuous: the full key really does populate the fields name_hash_only + // drops, and they differ across the two read groups / positions — so the + // name_hash parity above is a genuine skip guard, not a comparison of two + // all-default keys. + assert_ne!(full[0].library_idx, full[2].library_idx, "full key varies library_idx by RG"); + assert_ne!(full[0].pos1, full[2].pos1, "full key varies pos1 by mapped position"); + } + // ======================================================================== // Fail-closed validation of undersized / malformed records // ======================================================================== diff --git a/tests/integration/main.rs b/tests/integration/main.rs index 6963b56e0..05f91996c 100644 --- a/tests/integration/main.rs +++ b/tests/integration/main.rs @@ -37,6 +37,7 @@ mod test_extract_command; mod test_fastq_command; mod test_fastq_pipeline_memory_backpressure; mod test_filter_command; +mod test_filter_cutover_parity; mod test_group_command; mod test_group_determinism; mod test_input_source_matrix; diff --git a/tests/integration/test_filter_command.rs b/tests/integration/test_filter_command.rs index f7f038adf..b3b400855 100644 --- a/tests/integration/test_filter_command.rs +++ b/tests/integration/test_filter_command.rs @@ -527,8 +527,11 @@ fn test_filter_command_no_ref_mapped_reads_fails() { } ////////////////////////////////////////////////////////////////////////////// -// R3.1 chain-cutover parity tests: `--threads N` (declarative chain builder) -// vs. no-`--threads` (legacy unified-pipeline oracle). +// Worker-count-independence tests: filter always runs on the declarative chain +// builder (the legacy single-threaded path is retired), so a no-`--threads` run +// (the chain at a single worker) must produce output identical to a `--threads N` +// run. Byte-parity against the pre-removal binary lives in +// `test_filter_cutover_parity.rs`. ////////////////////////////////////////////////////////////////////////////// /// Run `filter` on `input` writing `output`, with `extra` args appended @@ -612,8 +615,8 @@ fn build_mixed_depth_templates(n: usize) -> Vec { records } -/// The chain (`--threads N`) path produces output identical to the non-chain -/// (no-`--threads`) path, for both `--threads 1` (the minimal chain engine) +/// A `--threads N` run produces output identical to a no-`--threads` run (the +/// chain at a single worker), for both `--threads 1` (the minimal chain engine) /// and `--threads 4` (genuinely parallel), and for both `filter-by-template` /// step factories (template mode groups via `GroupByQueryname` before /// filtering; single-read mode filters each record independently -- distinct @@ -636,8 +639,8 @@ fn test_filter_chain_matches_single_threaded( let template_flag = if filter_by_template { "true" } else { "false" }; let threads_str = threads.to_string(); - let oracle_out = temp_dir.path().join("oracle.bam"); - filter_run(&input_bam, &oracle_out, &ref_path, &["--filter-by-template", template_flag]); + let single_worker_out = temp_dir.path().join("single_worker.bam"); + filter_run(&input_bam, &single_worker_out, &ref_path, &["--filter-by-template", template_flag]); let chain_out = temp_dir.path().join("chain.bam"); filter_run( @@ -651,7 +654,7 @@ fn test_filter_chain_matches_single_threaded( // normalizes the `@PG` command-line field that legitimately differs by // `--threads`, so a chain regression that changed `@RG`/`@SQ`/`@CO` or any // other header record — not only sort order — is caught too. - let (oracle_header, expected) = crate::helpers::read_bam_output(&oracle_out); + let (single_worker_header, expected) = crate::helpers::read_bam_output(&single_worker_out); let (chain_header, actual) = crate::helpers::read_bam_output(&chain_out); // Non-vacuous: filtering must actually drop records, or a pass-through // regression on BOTH paths would satisfy `actual == expected` silently. @@ -662,39 +665,34 @@ fn test_filter_chain_matches_single_threaded( assert_eq!( expected.len(), expected_kept, - "oracle must keep exactly {expected_kept} of 16 records \ + "single-worker run must keep exactly {expected_kept} of 16 records \ (filter_by_template={filter_by_template})" ); assert_eq!( actual, expected, "chain (threads={threads}, filter_by_template={filter_by_template}) output must match \ - the non-chain path record-for-record" + the single-worker chain record-for-record" ); // Guard against a vacuous header comparison: if a regression dropped `@HD` on // both paths, the header equality could still pass on two equally-broken - // headers, so assert the oracle actually declares `@HD` first. - assert!(oracle_header.header().is_some(), "oracle output must declare @HD"); + // headers, so assert the single-worker run actually declares `@HD` first. + assert!(single_worker_header.header().is_some(), "single-worker output must declare @HD"); assert_eq!( - chain_header, oracle_header, + chain_header, single_worker_header, "chain (threads={threads}, filter_by_template={filter_by_template}) output header must \ - match the non-chain path (complete normalized header, not just @HD)" + match the single-worker chain (complete normalized header, not just @HD)" ); } /// Header carrying two `@RG` lines with distinct `LB` values, for /// `test_filter_chain_matches_single_threaded_with_rg_and_cb_variation`. The -/// `source_group_key_config` perf fix routes the chain filter path's first -/// stage to a `name_hash_only` `GroupKeyConfig`, which never resolves a -/// per-record library index or walks the CIGAR for position during decode — -/// unlike the filter-by-template oracle's own decode config -/// (`Filter::build_filter_pipeline_config`'s `new_raw_no_cell`), which -/// resolves the library index and the position but, like `name_hash_only`, -/// never extracts `CB` either way. So RG (library index) and position are the -/// fields that genuinely differ between the two decode configs, and are what -/// this parity test needs to vary to be non-vacuous; the `CB` variation below -/// is incidental fixture realism, not a discriminating input — the oracle -/// never reads `CB` regardless of this change, so it cannot expose a -/// regression here. +/// `source_group_key_config` perf fix routes the chain filter's first stage to a +/// `name_hash_only` `GroupKeyConfig`, which never resolves a per-record library +/// index or walks the CIGAR for position during decode. This fixture varies the +/// `@RG`/`LB` (library index) and position fields precisely so a regression that +/// made those fields load-bearing for filtering — despite the skip — would change +/// the output and be caught. The per-template `CB` tag is incidental fixture +/// realism, not a discriminating input. fn create_consensus_header_with_read_groups( ref_name: &str, ref_len: usize, @@ -784,19 +782,22 @@ fn build_mixed_depth_templates_with_rg_and_cb(n: usize) -> Vec { records } -/// The chain (`--threads N`) path matches the non-chain oracle record-for-record -/// even when the discarded position/RG portion of the group key genuinely -/// varies per record (distinct `@RG`/`LB` per template and distinct positions -/// — the fields the oracle's decode config resolves but `name_hash_only` -/// does not, per `create_consensus_header_with_read_groups`'s doc comment). -/// `test_filter_chain_matches_single_threaded` already covers the -/// `threads`/`filter_by_template` matrix on RG-free, single-position records; -/// this test's unique contribution is exercising the exact fields -/// `source_group_key_config`'s `name_hash_only` routing stops computing for -/// the chain filter path, confirming the skip is genuinely invisible in the -/// output. Records also carry a per-template `CB` tag for fixture realism, -/// but it is not discriminating: the oracle's own decode config never -/// extracts `CB` either, so `CB` variation cannot expose a regression here. +/// A `--threads 4` run matches the single-worker run record-for-record when the +/// input carries per-template `@RG`/`LB` and position variation (both runs take +/// the chain post-cutover, so this is a worker-count-invariance check, not a +/// chain-vs-legacy one). `test_filter_chain_matches_single_threaded` already +/// covers the `threads`/`filter_by_template` matrix on RG-free, single-position +/// records; this test's unique contribution is holding that invariance with the +/// RG/position fields present in the input. +/// +/// It does NOT by itself guard `source_group_key_config`'s `name_hash_only` skip +/// (the perf routing that stops computing library index/position during decode +/// for `Stage::Filter`): post-cutover both runs here take that same skip, so a bug +/// in it would change both sides equally. That invariant is pinned directly at the +/// decode consumer by +/// `decode::tests::name_hash_only_matches_full_key_grouping_when_rg_and_position_vary`. +/// The per-template `CB` tag below is incidental fixture realism; the decode config +/// never extracts `CB`. #[test] fn test_filter_chain_matches_single_threaded_with_rg_and_cb_variation() { let temp_dir = TempDir::new().unwrap(); @@ -809,8 +810,8 @@ fn test_filter_chain_matches_single_threaded_with_rg_and_cb_variation() { build_mixed_depth_templates_with_rg_and_cb(8), ); - let oracle_out = temp_dir.path().join("oracle.bam"); - filter_run(&input_bam, &oracle_out, &ref_path, &["--filter-by-template", "true"]); + let single_worker_out = temp_dir.path().join("single_worker.bam"); + filter_run(&input_bam, &single_worker_out, &ref_path, &["--filter-by-template", "true"]); let chain_out = temp_dir.path().join("chain.bam"); filter_run( @@ -820,46 +821,46 @@ fn test_filter_chain_matches_single_threaded_with_rg_and_cb_variation() { &["--filter-by-template", "true", "--threads", "4"], ); - let (oracle_header, expected) = crate::helpers::read_bam_output(&oracle_out); + let (single_worker_header, expected) = crate::helpers::read_bam_output(&single_worker_out); let (chain_header, actual) = crate::helpers::read_bam_output(&chain_out); // Non-vacuous: filter-by-template drops the 4 odd templates whole (keeps 8 // of 16), same as the RG/CB-free matrix test. assert_eq!( expected.len(), 8, - "oracle must keep exactly 8 of 16 records with RG/CB variation present" + "single-worker run must keep exactly 8 of 16 records with RG/CB variation present" ); assert_eq!( actual, expected, - "chain output must match the non-chain path record-for-record with RG/CB variation \ + "chain output must match the single-worker chain record-for-record with RG/CB variation \ present, proving source_group_key_config's name_hash_only skip for Stage::Filter \ changes no output" ); assert_eq!( - chain_header, oracle_header, - "chain output header must match the non-chain path with RG/CB variation present" + chain_header, single_worker_header, + "chain output header must match the single-worker chain with RG/CB variation present" ); } /// The `--rejects` output BAM matches record-for-record between the chain and -/// oracle paths, not just the kept output. Template mode over +/// single-worker runs, not just the kept output. Template mode over /// `build_mixed_depth_templates` rejects both mates of every odd-indexed /// template (8 records across 4 templates), so the rejects comparison is /// non-vacuous. #[test] -fn test_filter_chain_rejects_bam_matches_single_threaded() { +fn test_filter_chain_rejects_bam_matches_single_worker() { let temp_dir = TempDir::new().unwrap(); let input_bam = temp_dir.path().join("input.bam"); let ref_path = create_test_reference(temp_dir.path()); create_consensus_bam(&input_bam, build_mixed_depth_templates(8)); - let oracle_out = temp_dir.path().join("oracle.bam"); - let oracle_rejects = temp_dir.path().join("oracle.rejects.bam"); + let single_worker_out = temp_dir.path().join("single_worker.bam"); + let single_worker_rejects = temp_dir.path().join("single_worker.rejects.bam"); filter_run( &input_bam, - &oracle_out, + &single_worker_out, &ref_path, - &["--rejects", oracle_rejects.to_str().unwrap()], + &["--rejects", single_worker_rejects.to_str().unwrap()], ); let chain_out = temp_dir.path().join("chain.bam"); @@ -871,34 +872,39 @@ fn test_filter_chain_rejects_bam_matches_single_threaded() { &["--rejects", chain_rejects.to_str().unwrap(), "--threads", "4"], ); - let expected_kept = read_filter_records(&oracle_out); + let expected_kept = read_filter_records(&single_worker_out); let actual_kept = read_filter_records(&chain_out); - assert!(!expected_kept.is_empty(), "oracle kept output must be non-empty"); + assert!(!expected_kept.is_empty(), "single-worker kept output must be non-empty"); assert_eq!(actual_kept, expected_kept, "kept output must match record-for-record"); - let expected_rejects = read_filter_records(&oracle_rejects); + let expected_rejects = read_filter_records(&single_worker_rejects); let actual_rejects = read_filter_records(&chain_rejects); assert!( !expected_rejects.is_empty(), - "oracle rejects output must be non-empty (guard against a vacuous pass)" + "single-worker rejects output must be non-empty (guard against a vacuous pass)" ); assert_eq!( actual_rejects, expected_rejects, - "rejects output must match record-for-record between chain and oracle" + "rejects output must match record-for-record between chain and single-worker runs" ); } -/// The `--stats` file is byte-identical between the chain and oracle paths. +/// The `--stats` file is byte-identical between the chain and single-worker runs. #[test] -fn test_filter_chain_stats_file_matches_single_threaded() { +fn test_filter_chain_stats_file_matches_single_worker() { let temp_dir = TempDir::new().unwrap(); let input_bam = temp_dir.path().join("input.bam"); let ref_path = create_test_reference(temp_dir.path()); create_consensus_bam(&input_bam, build_mixed_depth_templates(8)); - let oracle_out = temp_dir.path().join("oracle.bam"); - let oracle_stats = temp_dir.path().join("oracle.stats.txt"); - filter_run(&input_bam, &oracle_out, &ref_path, &["--stats", oracle_stats.to_str().unwrap()]); + let single_worker_out = temp_dir.path().join("single_worker.bam"); + let single_worker_stats = temp_dir.path().join("single_worker.stats.txt"); + filter_run( + &input_bam, + &single_worker_out, + &ref_path, + &["--stats", single_worker_stats.to_str().unwrap()], + ); let chain_out = temp_dir.path().join("chain.bam"); let chain_stats = temp_dir.path().join("chain.stats.txt"); @@ -909,21 +915,25 @@ fn test_filter_chain_stats_file_matches_single_threaded() { &["--stats", chain_stats.to_str().unwrap(), "--threads", "4"], ); - let oracle_content = fs::read_to_string(&oracle_stats).expect("read oracle stats"); + let single_worker_content = + fs::read_to_string(&single_worker_stats).expect("read single-worker stats"); let chain_content = fs::read_to_string(&chain_stats).expect("read chain stats"); - assert!(!oracle_content.trim().is_empty(), "oracle stats file should not be empty"); - assert!(oracle_content.contains("total_reads"), "stats should contain total_reads"); + assert!( + !single_worker_content.trim().is_empty(), + "single-worker stats file should not be empty" + ); + assert!(single_worker_content.contains("total_reads"), "stats should contain total_reads"); // Non-vacuous: the stats must record that filtering actually rejected reads, // or a pass-through regression on both paths would still byte-match. - let failed = oracle_content + let failed = single_worker_content .lines() .find_map(|l| l.strip_prefix("failed_reads\t")) .and_then(|v| v.trim().parse::().ok()) .expect("stats file must report failed_reads"); - assert!(failed > 0, "filtering must reject at least one read; stats:\n{oracle_content}"); + assert!(failed > 0, "filtering must reject at least one read; stats:\n{single_worker_content}"); assert_eq!( - oracle_content, chain_content, - "stats file must be byte-identical across chain and oracle modes" + single_worker_content, chain_content, + "stats file must be byte-identical across chain and single-worker modes" ); } @@ -966,14 +976,14 @@ fn create_multi_batch_filter_input(path: &Path, n: usize) { create_consensus_bam(path, records); } -/// The chain's multi-worker path (`--threads 4`) must match the single-threaded -/// oracle across `GroupByQueryname`'s batch machinery: both its cross-input-batch +/// The chain's multi-worker run (`--threads 4`) must match the single-worker +/// run across `GroupByQueryname`'s batch machinery: both its cross-input-batch /// `current_template` carry-over (a paired template whose two mates land in /// consecutive input record-batches) AND its 1000-template output-emit boundary. /// `create_multi_batch_filter_input(1200)` emits 1200 paired templates (2400 /// records) to exercise both -- the small-fixture parity tests above (well under /// one batch, single-record templates) exercise neither. Compare `--threads 4` -/// output record-for-record against the oracle. +/// output record-for-record against the single-worker run. #[test] fn test_filter_chain_threads4_matches_single_threaded_multi_batch() { let temp_dir = TempDir::new().unwrap(); @@ -981,31 +991,31 @@ fn test_filter_chain_threads4_matches_single_threaded_multi_batch() { let ref_path = create_test_reference(temp_dir.path()); create_multi_batch_filter_input(&input_bam, 1200); - let oracle_out = temp_dir.path().join("oracle.bam"); - filter_run(&input_bam, &oracle_out, &ref_path, &[]); + let single_worker_out = temp_dir.path().join("single_worker.bam"); + filter_run(&input_bam, &single_worker_out, &ref_path, &[]); let chain_out = temp_dir.path().join("chain.bam"); filter_run(&input_bam, &chain_out, &ref_path, &["--threads", "4"]); - let expected = read_filter_records(&oracle_out); + let expected = read_filter_records(&single_worker_out); let actual = read_filter_records(&chain_out); // 1200 paired templates, all passing -> all 2400 records kept on both paths. assert_eq!( expected.len(), 2400, - "expected all 2400 records (1200 paired passing templates) in the oracle output" + "expected all 2400 records (1200 paired passing templates) in the single-worker output" ); assert_eq!( actual, expected, - "chain (--threads 4) output must match the single-threaded oracle record-for-record \ + "chain (--threads 4) output must match the single-worker chain record-for-record \ across GroupByQueryname input-batch carry-over and output-emit boundaries" ); } -/// FILT3-02 on the chain path: mirrors `test_filter_rejects_coordinate_sorted_input` -/// but with `--threads 4`, proving the `require_query_grouped` guard added to -/// `add_filter` (builder.rs) rejects coordinate-sorted input on the chain path -/// exactly as the legacy path does -- not just when the legacy oracle is used. +/// FILT3-02 with `--threads 4`: mirrors `test_filter_rejects_coordinate_sorted_input` +/// (a no-`--threads` run) but at four workers, proving the `require_query_grouped` +/// guard in `add_filter` (builder.rs) rejects coordinate-sorted input regardless of +/// worker count (the chain is the only filter execution path). #[test] fn test_filter_chain_rejects_coordinate_sorted_input() { let temp_dir = TempDir::new().unwrap(); @@ -1042,14 +1052,13 @@ fn test_filter_chain_rejects_coordinate_sorted_input() { assert!(msg.contains("queryname sorted or query grouped"), "unexpected error message: {msg}"); } -/// FILT3-02 on the chain path, header-less variant: mirrors -/// `test_filter_rejects_headerless_input` but with `--threads 4`. Both paths now -/// synthesize `@HD VN:1.6 SO:unsorted` for header-less input — the legacy path -/// via `ensure_hd_record` in `execute()`, the chain path via `ensure_hd_record` -/// in `ChainBuilder::new` (before `add_pg_record`) — and `require_query_grouped` -/// rejects `SO:unsorted`, so the chain path rejects header-less input with the -/// same message. Pins chain-path header-less rejection so a future -/// header-handling refactor cannot silently break it. +/// FILT3-02 header-less variant with `--threads 4`: mirrors +/// `test_filter_rejects_headerless_input` (a no-`--threads` run) but at four +/// workers. The chain synthesizes `@HD VN:1.6 SO:unsorted` for header-less input +/// via `ensure_hd_record` in `ChainBuilder::new` (before `add_pg_record`), and +/// `require_query_grouped` rejects `SO:unsorted`, so header-less input is rejected +/// with the same message regardless of worker count. Pins header-less rejection so +/// a future header-handling refactor cannot silently break it. #[test] fn test_filter_chain_rejects_headerless_input() { let temp_dir = TempDir::new().unwrap(); diff --git a/tests/integration/test_filter_cutover_parity.rs b/tests/integration/test_filter_cutover_parity.rs new file mode 100644 index 000000000..6e1a36db9 --- /dev/null +++ b/tests/integration/test_filter_cutover_parity.rs @@ -0,0 +1,508 @@ +//! Parity gate for the `filter` command's single-threaded-path retirement (C4): +//! `Filter::execute` no longer has a serial in-process unified-pipeline loop +//! reached when `--threads` is absent — it *always* routes through the +//! declarative chain builder. This test proves that cutover lost nothing +//! user-observable. +//! +//! Two independent things are checked here: +//! +//! 1. **The cutover actually happened** (`filter_no_threads_routes_through_chain`). +//! A no-`--threads` run now emits the chain-only `"Using pipeline with N +//! threads"` banner that `ChainBuilder::add_filter` logs and the retired +//! unified-pipeline tail never did. This is the genuine RED/GREEN +//! discriminator: before the removal a no-`--threads` run took the serial +//! path and printed no such line; after it, the chain does. +//! +//! 2. **Output parity with the pre-removal serial path** +//! (`cutover_matches_baseline`). The current build's `filter` output — kept +//! records and (when written) the rejects BAM (both byte-identical modulo the +//! `@PG` line) and the `--stats` TSV — must match the frozen owned-serial +//! baseline binary. The baseline path comes from `FGUMI_BASELINE_BIN`; when it +//! is unset (or names a missing file) the case degrades to a self-consistency +//! oracle (min-reads filtering, base masking, and NM regeneration asserted +//! directly) rather than skipping — the exact fallback discipline of +//! `test_sort_cutover_parity.rs`. Cases cover `--min-reads` filtering (records +//! kept vs rejected), `--ref`-based NM/UQ/MD regeneration, base masking, and +//! both `filter-by-template` modes, plus `--rejects` and `--stats`. +//! +//! **Not a RED/GREEN gate for the parity half.** Because the removed serial loop +//! and the chain were already output-equivalent (see the in-tree +//! worker-count-independence tests in `test_filter_command.rs`), the baseline +//! byte-parity check passes on both sides of the change — like the sort cutover +//! it guards equivalence, it does not observe a regression the cutover +//! introduces. The chain-banner check in (1) is the part that flips RED→GREEN +//! across the removal. + +use std::ffi::OsStr; +use std::path::Path; +use std::process::Command; + +use fgumi_raw_bam::{RawRecord, SamBuilder, flags}; +use rstest::rstest; +use tempfile::TempDir; + +use crate::helpers::bam_generator::{create_minimal_header, create_test_reference, write_bam}; +use crate::helpers::read_bam_output; +use fgumi_lib::sam::SamTag; + +/// Resolves the saved pre-removal serial baseline binary to compare against. +/// +/// The path comes solely from `FGUMI_BASELINE_BIN`; when unset (or naming a +/// missing file) this returns `None` — "no baseline oracle available", never a +/// silent pass. Callers layer the baseline byte-parity check on top of the +/// always-available self-consistency oracle; a missing baseline drops only the +/// byte-parity half, it never skips the case. No hardcoded fallback: a baseline +/// binary is host-specific and must never be a path committed into the repo. +fn baseline_bin() -> Option { + let path = std::path::PathBuf::from(std::env::var_os("FGUMI_BASELINE_BIN")?); + if path.is_file() { + return Some(path); + } + eprintln!( + "FGUMI_BASELINE_BIN={} does not name an existing file; baseline oracle unavailable", + path.display() + ); + None +} + +/// Removes every `@PG` line from a SAM header text blob. +/// +/// Both the current build and the baseline binary stamp a single `@PG` line +/// whose `VN` (git-describe version) and `CL` (command line, naming argv[0] — a +/// different binary path for each side — and the per-run output path) necessarily +/// differ between two independent invocations. Stripping the whole line is +/// correct here because neither side emits more than one `@PG` for these inputs +/// (the test BAMs carry no pre-existing `@PG`). +fn strip_pg_lines(text: &str) -> String { + if text.is_empty() { + return String::new(); + } + let mut lines: Vec<&str> = text.split('\n').collect(); + let had_trailing_newline = lines.last() == Some(&""); + if had_trailing_newline { + lines.pop(); + } + lines.retain(|line| !line.starts_with("@PG")); + let mut out = lines.join("\n"); + if had_trailing_newline { + out.push('\n'); + } + out +} + +/// Reads a BAM's raw BGZF stream, decompresses it, strips the `@PG` line(s) from +/// the embedded SAM header text, and returns the resulting bytes (new header + +/// unmodified `n_ref`/reference-list/record bytes). +/// +/// Comparing the *decompressed* BAM binary — rather than `RecordBuf`-parsed +/// records or raw file bytes — keeps everything except the `@PG` line an exact, +/// uninterpreted byte comparison: BGZF block boundaries differ between the two +/// writers even for identical logical content (so raw file bytes never match), +/// while parsing into `RecordBuf` and re-encoding could mask a real tag-order or +/// binary-layout regression by normalizing it away. filter regenerates NM/UQ/MD +/// and masks bases in place, so its output record bytes must match the pre-removal +/// path exactly, which is what this comparison pins. +fn decompressed_records_without_pg(path: &Path) -> Vec { + let file = std::fs::File::open(path).unwrap_or_else(|e| panic!("open {}: {e}", path.display())); + let mut reader = noodles::bgzf::io::Reader::new(std::io::BufReader::new(file)); + let mut raw = Vec::new(); + std::io::Read::read_to_end(&mut reader, &mut raw) + .unwrap_or_else(|e| panic!("decompress BGZF stream for {}: {e}", path.display())); + + assert!( + raw.len() >= 8 && &raw[0..4] == b"BAM\x01", + "{} does not decompress to a BAM binary stream (missing magic)", + path.display() + ); + let l_text = i32::from_le_bytes(raw[4..8].try_into().expect("4 bytes")); + let l_text = usize::try_from(l_text).expect("l_text is non-negative"); + let text_start = 8; + let text_end = text_start + l_text; + assert!(raw.len() >= text_end, "{} header text runs past end of stream", path.display()); + + let text = String::from_utf8_lossy(&raw[text_start..text_end]); + let stripped_text = strip_pg_lines(&text); + + let mut out = Vec::with_capacity(raw.len()); + out.extend_from_slice(b"BAM\x01"); + let new_l_text = i32::try_from(stripped_text.len()).expect("stripped header text fits i32"); + out.extend_from_slice(&new_l_text.to_le_bytes()); + out.extend_from_slice(stripped_text.as_bytes()); + out.extend_from_slice(&raw[text_end..]); // n_ref, reference list, and all records, untouched + out +} + +/// One mapped consensus read: `name`, `flags`, position, per-consensus depth (`cD`) +/// and error (`cE`), and per-base depth (`cd`)/error (`ce`) arrays. The sequence +/// is 8 bases; a mismatch to the reference or a masked base changes the +/// regenerated `NM`, and a low per-base depth drives base masking. +fn consensus_read( + name: &str, + flags: u16, + pos: i32, + depth: i32, + seq: &[u8], + per_base_depth: &[u16], +) -> RawRecord { + let mut b = SamBuilder::new(); + b.read_name(name.as_bytes()) + .flags(flags) + .ref_id(0) + .pos(pos) + .mapq(60) + .cigar_ops(&[u32::try_from(seq.len()).expect("seq len fits u32") << 4]) // M + .sequence(seq) + .qualities(&vec![35u8; seq.len()]); + b.add_int_tag(SamTag::CD, depth).add_float_tag(SamTag::CE, 0.0_f32); + b.add_array_u16(SamTag::CD_BASES, per_base_depth) + .add_array_u16(SamTag::CE_BASES, &vec![0u16; seq.len()]); + b.build() +} + +/// The parity corpus, at reference-aligned positions (`create_test_reference` +/// writes `ACGTACGT` repeated, so an 8-base window at any multiple of 8 is +/// `ACGTACGT`). Exercises every behavior the brief names in one input. +/// +/// Four single-end reads (each its own template, so they behave identically in +/// both `--filter-by-template` modes): +/// - `pass`: depth 10, exact match, no masking — kept, `NM` 0. +/// - `low_depth`: depth 1 (< `--min-reads 3`) — rejected by the read-level filter. +/// - `masked`: depth 10 but one per-base depth 1 (< 3) — that base masked to `N`, +/// read still kept (1 `N` of 8 = 0.125 < `--max-no-call-fraction 0.2`), `NM` +/// regenerated over the masked base. +/// - `mismatch`: depth 10, first base `T` vs reference `A` — kept, `NM` 1 after +/// regeneration. +/// +/// Plus one paired template (`pair` R1 + R2, adjacent as query-grouped input +/// requires) whose R1 passes (depth 10) but R2 fails (depth 1). This is the +/// record that makes `--filter-by-template` genuinely different from single-read +/// mode: template mode drops the *whole* template (both mates rejected, R1 too, +/// under fgbio's "all primaries must pass" rule), while single-read mode keeps R1 +/// and drops only R2. Without it the `template` and `single_read` cases would +/// exercise identical behavior. +fn build_parity_records() -> Vec { + let paired_r1 = flags::PAIRED | flags::FIRST_SEGMENT; + let paired_r2 = flags::PAIRED | flags::LAST_SEGMENT; + vec![ + consensus_read("pass", 0, 0, 10, b"ACGTACGT", &[10; 8]), + consensus_read("low_depth", 0, 8, 1, b"ACGTACGT", &[1; 8]), + consensus_read("masked", 0, 16, 10, b"ACGTACGT", &[10, 10, 10, 1, 10, 10, 10, 10]), + consensus_read("mismatch", 0, 24, 10, b"TCGTACGT", &[10; 8]), + consensus_read("pair", paired_r1, 32, 10, b"ACGTACGT", &[10; 8]), + consensus_read("pair", paired_r2, 40, 1, b"ACGTACGT", &[1; 8]), + ] +} + +/// Writes the parity corpus to `dir/in.bam` (query-grouped header, which filter +/// requires) and returns its path. +fn write_parity_input(dir: &Path) -> std::path::PathBuf { + let input = dir.join("in.bam"); + write_bam(&input, &create_minimal_header("chr1", 10_000), &build_parity_records()); + input +} + +/// Runs ` filter -i -o --ref --min-reads 3 +/// --max-no-call-fraction 0.2 --compression-level 1 --filter-by-template +/// [--rejects ] [--stats ]` (no `--threads`, so the current build takes the +/// post-cutover chain path and the baseline takes its serial path) with +/// `RUST_LOG=info`, and returns the process output for stderr assertions. +fn run_filter( + bin: &Path, + input: &Path, + output: &Path, + ref_path: &Path, + filter_by_template: bool, + rejects: Option<&Path>, + stats: Option<&Path>, +) -> std::process::Output { + let mut cmd = Command::new(bin); + cmd.env("RUST_LOG", "info").args([ + OsStr::new("filter"), + OsStr::new("-i"), + input.as_os_str(), + OsStr::new("-o"), + output.as_os_str(), + OsStr::new("--ref"), + ref_path.as_os_str(), + OsStr::new("--min-reads"), + OsStr::new("3"), + OsStr::new("--max-no-call-fraction"), + OsStr::new("0.2"), + OsStr::new("--compression-level"), + OsStr::new("1"), + OsStr::new("--filter-by-template"), + OsStr::new(if filter_by_template { "true" } else { "false" }), + ]); + if let Some(r) = rejects { + cmd.args([OsStr::new("--rejects"), r.as_os_str()]); + } + if let Some(s) = stats { + cmd.args([OsStr::new("--stats"), s.as_os_str()]); + } + cmd.output().unwrap_or_else(|e| panic!("failed to spawn `{}` filter: {e}", bin.display())) +} + +/// A no-`--threads` run now routes through the declarative chain, which logs the +/// `"Using pipeline with N threads"` banner from `ChainBuilder::add_filter`. The +/// retired unified-pipeline tail logged no such line, so this is the RED +/// (pre-removal) → GREEN (post-removal) discriminator for the cutover. +#[test] +fn filter_no_threads_routes_through_chain() { + let dir = TempDir::new().expect("temp dir"); + let input = write_parity_input(dir.path()); + let ref_path = create_test_reference(dir.path()); + let output = dir.path().join("out.bam"); + + let current_bin = Path::new(env!("CARGO_BIN_EXE_fgumi")); + let out = run_filter(current_bin, &input, &output, &ref_path, true, None, None); + let stderr = String::from_utf8_lossy(&out.stderr); + assert!(out.status.success(), "no-threads filter must succeed; stderr:\n{stderr}"); + assert!( + stderr.contains("Using pipeline with 1 threads"), + "a no-`--threads` filter must route through the chain (which logs the pipeline banner); \ + the serial path is retired. stderr:\n{stderr}" + ); + assert!( + stderr.contains("Starting Filter"), + "the chain must still emit the `Starting Filter` banner; stderr:\n{stderr}" + ); +} + +/// Output parity of the post-cutover chain against the pre-removal serial +/// baseline binary, across `filter-by-template` modes and the `--rejects` / +/// `--stats` outputs — plus the always-available self-consistency oracle when no +/// baseline is set. +/// +/// `#[case]` args: a label, `filter_by_template`, `with_rejects`, `with_stats`. +#[rstest] +#[case::template(true, false, false)] +#[case::single_read(false, false, false)] +#[case::template_rejects(true, true, false)] +#[case::template_stats(true, false, true)] +#[case::single_read_rejects(false, true, false)] +#[case::single_read_stats(false, false, true)] +fn cutover_matches_baseline( + #[case] filter_by_template: bool, + #[case] with_rejects: bool, + #[case] with_stats: bool, +) { + let dir = TempDir::new().expect("temp dir"); + let input = write_parity_input(dir.path()); + let ref_path = create_test_reference(dir.path()); + + let current_bin = Path::new(env!("CARGO_BIN_EXE_fgumi")); + let current_out = dir.path().join("current.bam"); + let current_rejects = dir.path().join("current.rejects.bam"); + let current_stats = dir.path().join("current.stats.txt"); + let cur_rej = with_rejects.then(|| current_rejects.clone()); + let cur_stats = with_stats.then(|| current_stats.clone()); + let current = run_filter( + current_bin, + &input, + ¤t_out, + &ref_path, + filter_by_template, + cur_rej.as_deref(), + cur_stats.as_deref(), + ); + let current_stderr = String::from_utf8_lossy(¤t.stderr); + assert!(current.status.success(), "current filter must succeed; stderr:\n{current_stderr}"); + + if let Some(baseline) = baseline_bin() { + let baseline_out = dir.path().join("baseline.bam"); + let baseline_rejects = dir.path().join("baseline.rejects.bam"); + let baseline_stats = dir.path().join("baseline.stats.txt"); + let base_rej = with_rejects.then(|| baseline_rejects.clone()); + let base_stats = with_stats.then(|| baseline_stats.clone()); + let base = run_filter( + &baseline, + &input, + &baseline_out, + &ref_path, + filter_by_template, + base_rej.as_deref(), + base_stats.as_deref(), + ); + assert!( + base.status.success(), + "baseline filter failed; stderr:\n{}", + String::from_utf8_lossy(&base.stderr) + ); + + assert_eq!( + decompressed_records_without_pg(¤t_out), + decompressed_records_without_pg(&baseline_out), + "chain filter output diverges from the pre-removal serial baseline binary ({}) \ + after stripping @PG — a real cutover parity bug, not something to relax", + baseline.display(), + ); + if with_rejects { + assert_eq!( + decompressed_records_without_pg(¤t_rejects), + decompressed_records_without_pg(&baseline_rejects), + "chain --rejects output diverges from the serial baseline binary" + ); + } + if with_stats { + assert_eq!( + std::fs::read_to_string(¤t_stats).expect("current stats"), + std::fs::read_to_string(&baseline_stats).expect("baseline stats"), + "chain --stats TSV diverges from the serial baseline binary" + ); + } + } else { + eprintln!( + "SKIP baseline half of cutover_matches_baseline[template={filter_by_template}, \ + rejects={with_rejects}, stats={with_stats}]: FGUMI_BASELINE_BIN is unset or does \ + not name an existing file — running self-consistency oracle instead" + ); + assert_self_consistent( + filter_by_template, + ¤t_out, + cur_rej.as_deref(), + cur_stats.as_deref(), + ); + } +} + +/// Always-available oracle used when no baseline binary is set: the chain's +/// output must show exactly the `--min-reads` filtering (including the +/// template-drop that distinguishes the two modes), base masking, and NM +/// regeneration the parity corpus prescribes. +/// +/// Kept/rejected sets differ by mode because of the `pair` template (R1 passes, +/// R2 fails at depth 1): +/// - template mode: keeps `pass`, `masked`, `mismatch` (3); rejects `low_depth` +/// and both `pair` mates (the whole template drops even though R1 passed). +/// - single-read mode: additionally keeps `pair` R1 (4 kept); rejects `low_depth` +/// and `pair` R2 (2). +fn assert_self_consistent( + filter_by_template: bool, + output: &Path, + rejects: Option<&Path>, + stats: Option<&Path>, +) { + let (total, passed, failed): (usize, usize, usize) = + if filter_by_template { (6, 3, 3) } else { (6, 4, 2) }; + + let paired_r1 = flags::PAIRED | flags::FIRST_SEGMENT; + let paired_r2 = flags::PAIRED | flags::LAST_SEGMENT; + + // Assert the *exact* kept set as a (read name, flag bits) multiset, not just + // names and a count. This pins the identity of every kept record: in + // single-read mode `pair` R1 is kept with FIRST_SEGMENT set (never R2 nor a + // flag-stripped duplicate), while template mode drops the whole `pair` + // template. The three single-end passers (`pass`, `masked`, `mismatch`) keep + // their flags of 0; `low_depth` is absent from either mode. + let (_, kept) = read_bam_output(output); + let expected_kept: Vec<(&str, u16)> = if filter_by_template { + vec![("pass", 0), ("masked", 0), ("mismatch", 0)] + } else { + vec![("pass", 0), ("masked", 0), ("mismatch", 0), ("pair", paired_r1)] + }; + assert_eq!( + sorted_name_flags(&kept), + sorted_expected(&expected_kept), + "filter_by_template={filter_by_template}: kept (read name, flags) multiset mismatch" + ); + + // Base masking: the `masked` read carries exactly one N, at the base whose + // per-base depth (index 3) fell below --min-reads — every other base intact. + let masked = kept.iter().find(|r| record_name(r) == "masked").expect("masked read present"); + let seq: Vec = masked.sequence().as_ref().to_vec(); + let n_positions: Vec = + seq.iter().enumerate().filter(|&(_, &b)| b == b'N').map(|(i, _)| i).collect(); + assert_eq!(n_positions, vec![3], "exactly base index 3 must be masked to N, nothing else"); + + // --ref NM/UQ/MD regeneration: exact-match read is NM 0; the single-base + // mismatch read is NM 1; masking base 3 of the `masked` read (an exact match + // pre-masking) turns that base into an N/ref difference, so its NM regenerates + // to 1 as well. + let pass = kept.iter().find(|r| record_name(r) == "pass").expect("pass read present"); + assert_eq!(nm_tag(pass), Some(0), "exact-match read NM must be regenerated to 0"); + let mismatch = + kept.iter().find(|r| record_name(r) == "mismatch").expect("mismatch read present"); + assert_eq!(nm_tag(mismatch), Some(1), "mismatch read NM must be regenerated to 1"); + assert_eq!(nm_tag(masked), Some(1), "masked read NM must regenerate over the masked base"); + + if let Some(rejects_path) = rejects { + // Assert the *exact* rejected set as a (read name, flag bits) multiset, + // so a duplicated or flag-substituted `low_depth` record cannot pass. In + // template mode the whole `pair` template is rejected (R1 too, with + // FIRST_SEGMENT); in single-read mode only `pair` R2 (LAST_SEGMENT) is. + let (_, rejected) = read_bam_output(rejects_path); + let expected_rejected: Vec<(&str, u16)> = if filter_by_template { + vec![("low_depth", 0), ("pair", paired_r1), ("pair", paired_r2)] + } else { + vec![("low_depth", 0), ("pair", paired_r2)] + }; + assert_eq!( + sorted_name_flags(&rejected), + sorted_expected(&expected_rejected), + "filter_by_template={filter_by_template}: rejected (read name, flags) multiset mismatch" + ); + } + + if let Some(stats_path) = stats { + let tsv = std::fs::read_to_string(stats_path).expect("read stats"); + // Parse the `keyvalue` rows and compare each field's value *exactly*. + // A substring match (`contains("total_reads\t6")`) would also accept a + // wrong `total_reads\t60`, so the count must be pinned to the precise + // value, not merely a prefix of it. + let fields = parse_stats_fields(&tsv); + for (name, expected) in + [("total_reads", total), ("passed_reads", passed), ("failed_reads", failed)] + { + let value = fields + .get(name) + .unwrap_or_else(|| panic!("stats missing `{name}` row; got:\n{tsv}")); + assert_eq!(value, &expected.to_string(), "stats `{name}` value mismatch; got:\n{tsv}"); + } + } +} + +/// Parse a `keyvalue` stats TSV (as written by `filter --stats`) into a map +/// of field name to its raw value string, so each field can be compared exactly +/// rather than by substring. A non-empty row without a tab is malformed and fails +/// loudly with the offending line. +fn parse_stats_fields(tsv: &str) -> std::collections::HashMap { + tsv.lines() + .filter(|line| !line.is_empty()) + .map(|line| { + line.split_once('\t').map_or_else( + || panic!("malformed stats row (no tab): {line:?}"), + |(key, value)| (key.to_string(), value.to_string()), + ) + }) + .collect() +} + +/// The read name of a parsed `RecordBuf`. +fn record_name(record: &noodles::sam::alignment::RecordBuf) -> String { + record.name().map(|n| String::from_utf8_lossy(n.as_ref()).into_owned()).unwrap_or_default() +} + +/// The sorted `(read name, flag bits)` multiset of a record set, so kept/rejected +/// identities can be compared exactly and order-independently. +fn sorted_name_flags(records: &[noodles::sam::alignment::RecordBuf]) -> Vec<(String, u16)> { + let mut pairs: Vec<(String, u16)> = + records.iter().map(|r| (record_name(r), r.flags().bits())).collect(); + pairs.sort(); + pairs +} + +/// The same sorted multiset built from a `(name, flags)` expectation table. +fn sorted_expected(expected: &[(&str, u16)]) -> Vec<(String, u16)> { + let mut pairs: Vec<(String, u16)> = + expected.iter().map(|&(name, flags)| (name.to_string(), flags)).collect(); + pairs.sort(); + pairs +} + +/// Reads the integer `NM` aux tag off a parsed `RecordBuf`, or `None` when absent +/// or not an integer. +fn nm_tag(record: &noodles::sam::alignment::RecordBuf) -> Option { + use noodles::sam::alignment::record::data::field::Tag; + let nm: SamTag = "NM".parse().expect("valid NM tag"); + record.data().get(&Tag::from(nm))?.as_int() +}