diff --git a/src/lib/commands/extract.rs b/src/lib/commands/extract.rs index 3c44c1089..47f0352f1 100644 --- a/src/lib/commands/extract.rs +++ b/src/lib/commands/extract.rs @@ -23,18 +23,13 @@ use crate::commands::common::{ }; use crate::fastq::FastqSegment; use crate::fastq::FastqSet; -use crate::fastq::ReadSetIterator; use crate::fastq_deinterleave::deinterleave; use crate::fastq_parse::strip_read_suffix; -use crate::logging::OperationTimer; use crate::sam::SamTag; -use crate::unified_pipeline::fastq_out_of_sync_error; use crate::validation::validate_input_exists; use anyhow::{Context, Result, bail, ensure}; use bstr::{BString, ByteSlice}; use clap::Parser; -use fgumi_bam_io::ProgressTracker; -use fgumi_bam_io::{RawBamWriter, create_raw_bam_writer}; use fgumi_raw_bam::UnmappedSamBuilder; use fgumi_raw_bam::fields::flags; use log::{debug, info}; @@ -45,17 +40,6 @@ use crate::read_structure::{ReadStructure, SegmentType}; use fgumi_bam_io::create_bam_reader; use fgumi_bam_io::{ChainedReader, InputFormat, TeeReader, is_stdin_path}; use fgumi_simd_fastq::SimdFastqReader; -use noodles::sam::header::Header; -use noodles::sam::header::record::value::Map; -use noodles::sam::header::record::value::map::ReadGroup; -use noodles::sam::header::record::value::map::builder::Builder; -use noodles::sam::header::record::value::map::header::group_order; -use noodles::sam::header::record::value::map::header::sort_order; -use noodles::sam::header::record::value::map::read_group::tag as rg_tag; -use noodles::sam::header::record::value::{ - Map as HeaderRecordMap, - map::{Header as HeaderRecord, Tag as HeaderTag}, -}; use std::fs::File; use std::io::{self, BufRead, BufReader}; use std::path::{Path, PathBuf}; @@ -684,9 +668,9 @@ impl Extract { /// and overrides the field before the extract step runs (see /// [`ExtractOptions::quality_encoding`]). Every other field maps directly /// from the identically-named CLI flag; `platform` (CLI default `"illumina"`) - /// becomes `Some(..)` so `build_fastq_header` always emits `@RG PL:`, matching - /// `Self::create_header`. `clipping_attribute` is intentionally dropped — it - /// does not apply to FASTQ input (there is no existing clipping to adjust). + /// becomes `Some(..)` so `build_fastq_header` always emits `@RG PL:`. + /// `clipping_attribute` is intentionally dropped — it does not apply to FASTQ + /// input (there is no existing clipping to adjust). #[must_use] pub fn to_extract_options(&self) -> ExtractOptions { ExtractOptions { @@ -715,16 +699,17 @@ impl Extract { } } - /// Run extract on the declarative chain builder (the `--threads` path). + /// Run extract on the declarative chain builder — the only execution path. /// - /// Hand-builds the [`ChainSpec`] (rather than using - /// [`ChainSpec::single_stage`], which is BAM-in/BAM-out): the source is a - /// FASTQ source — [`SourceSpec::InterleavedFastq`] for `--interleaved`, else + /// [`Command::execute`] always dispatches here, with or without `--threads` + /// (absent `--threads` runs the chain at a single worker). Hand-builds the + /// [`ChainSpec`] (rather than using [`ChainSpec::single_stage`], which is + /// BAM-in/BAM-out): the source is a FASTQ source — + /// [`SourceSpec::InterleavedFastq`] for `--interleaved`, else /// [`SourceSpec::Fastqs`] — and the sink is a BAM. The chain opens its own /// readers and detects the quality encoding in `ChainBuilder::open_source`; /// `read_streams`/`verify_crc` are BAM-reader knobs and inert here (the FASTQ - /// source carries its CRC policy in [`ExtractOptions`]). The no-`--threads` - /// serial loop in [`Command::execute`] is the in-process parity oracle. + /// source carries its CRC policy in [`ExtractOptions`]). /// /// [`ChainSpec`]: crate::pipeline::chains::ChainSpec /// [`ChainSpec::single_stage`]: crate::pipeline::chains::ChainSpec::single_stage @@ -863,70 +848,6 @@ impl Extract { Ok(()) } - /// Helper to conditionally add a tag/value pair to a read group - /// - /// If the value is Some, inserts the tag with the value into the read group builder. - /// If the value is None, returns the builder unchanged. - /// - /// # Arguments - /// * `rg` - The read group builder - /// * `tag` - The tag to insert - /// * `value` - Optional value to insert - /// - /// # Returns - /// The read group builder, potentially with the tag added - fn add_to_read_group( - rg: Builder, - tag: noodles::sam::header::record::value::map::tag::Other, - value: Option<&String>, - ) -> Builder { - if let Some(v) = value { rg.insert(tag, v.clone()) } else { rg } - } - - /// Create SAM header - fn create_header(&self, command_line: &str) -> Result
{ - let mut header = Header::builder(); - - // Sort and group order - let HeaderTag::Other(so_tag) = HeaderTag::from([b'S', b'O']) else { unreachable!() }; - let HeaderTag::Other(go_tag) = HeaderTag::from([b'G', b'O']) else { unreachable!() }; - let map = HeaderRecordMap::::builder() - .insert(so_tag, sort_order::UNSORTED) - .insert(go_tag, group_order::QUERY) - .build()?; - header = header.set_header(map); - - // Add comments - for comment in &self.comment { - header = header.add_comment(comment.clone()); - } - - // Create read group - let mut rg = Map::::builder(); - rg = Self::add_to_read_group(rg, rg_tag::SAMPLE, Some(&self.sample.clone())); - rg = Self::add_to_read_group(rg, rg_tag::LIBRARY, Some(&self.library.clone())); - rg = Self::add_to_read_group(rg, rg_tag::BARCODE, self.barcode.as_ref()); - rg = Self::add_to_read_group(rg, rg_tag::PLATFORM, Some(&self.platform)); - rg = Self::add_to_read_group(rg, rg_tag::PLATFORM_UNIT, self.platform_unit.as_ref()); - rg = Self::add_to_read_group(rg, rg_tag::PLATFORM_MODEL, self.platform_model.as_ref()); - rg = - Self::add_to_read_group(rg, rg_tag::SEQUENCING_CENTER, self.sequencing_center.as_ref()); - rg = Self::add_to_read_group( - rg, - rg_tag::PREDICTED_MEDIAN_INSERT_SIZE, - self.predicted_insert_size.map(|i| i.to_string()).as_ref(), - ); - rg = Self::add_to_read_group(rg, rg_tag::DESCRIPTION, self.description.as_ref()); - rg = Self::add_to_read_group(rg, rg_tag::PRODUCED_AT, self.run_date.as_ref()); - - header = header.add_read_group(self.read_group_id.clone(), rg.build()?); - - // Add @PG record - header = crate::commands::common::add_pg_to_builder(header, command_line)?; - - Ok(header.build()) - } - /// Joins byte slices with a separator, pre-allocating capacity. /// Returns empty `BString` if iterator is empty. fn join_bytes_with_separator<'a>( @@ -1024,43 +945,6 @@ impl Extract { Ok((name_part.to_vec(), None)) } - /// Validates that all read names match across the read sets - fn validate_read_names_match(read_sets: &[FastqSet]) -> Result<()> { - if read_sets.is_empty() { - return Ok(()); - } - - // Extract the read name from the first header (removing @ prefix if present) - let first_header = &read_sets[0].header; - let first_name = if first_header.starts_with(b"@") { - &first_header[1..] - } else { - first_header.as_slice() - }; - - // Strip space comments and /1, /2 suffixes for comparison - let first_name_part = strip_read_suffix(first_name); - - // Check that all other read sets have the same name - for (i, read_set) in read_sets.iter().enumerate().skip(1) { - let header = &read_set.header; - let name = if header.starts_with(b"@") { &header[1..] } else { header.as_slice() }; - - let name_part = strip_read_suffix(name); - - if name_part != first_name_part { - bail!( - "Read names do not match across FASTQs: '{}' vs '{}' (FASTQ index 0 vs {})", - String::from_utf8_lossy(first_name_part), - String::from_utf8_lossy(name_part), - i - ); - } - } - - Ok(()) - } - /// Context for a read name that BAM cannot represent. /// /// The builder's own error reports the length and the limit but not which @@ -1084,234 +968,6 @@ impl Extract { }; format!("could not write the record for read {shown}") } - - /// Builds one raw BAM record for a template segment, substituting a single `N` @ Q2 when - /// the segment has no sequence. - /// - /// The name comes from the input FASTQ header (plus the UMI when `--annotate-read-names` - /// is set), so an over-long name is bad input rather than a bug: fail cleanly via - /// `try_build_record` instead of panicking partway through writing the output BAM, and - /// attach [`Self::read_name_too_long_context`] so the message names the offending read. - /// - /// Used by the serial-oracle `make_raw_records` path; the chain path's - /// `make_raw_records_from_fastq_set` builds records through the same - /// `UnmappedSamBuilder` API, which the parity tests keep in step. - fn build_template_record( - builder: &mut UnmappedSamBuilder, - name: &[u8], - flag: u16, - seq: &[u8], - quals: &[u8], - encoding: QualityEncoding, - ) -> Result<()> { - if seq.is_empty() { - builder.try_build_record(name, flag, b"N", &[2u8]) - } else { - builder.try_build_record(name, flag, seq, &encoding.to_standard_numeric(quals)) - } - .with_context(|| Self::read_name_too_long_context(name)) - } - - /// Write raw BAM records from a read set directly to a writer. - /// - /// Uses `UnmappedSamBuilder` to construct records as raw bytes, - /// bypassing `RecordBuf` allocation and encoding overhead. - /// - /// Returns the number of records written. - #[allow(clippy::too_many_lines)] - fn make_raw_records( - &self, - read_set: &FastqSet, - encoding: QualityEncoding, - builder: &mut UnmappedSamBuilder, - writer: &mut RawBamWriter, - ) -> Result { - let templates: Vec<&FastqSegment> = read_set.template_segments().collect(); - - let read_name = String::from_utf8_lossy(&read_set.header); - ensure!(!templates.is_empty(), "No template segments found for read: {read_name}"); - - // Extract various barcode types as BString - use optimized join - let cell_barcode_bs = Self::join_bytes_with_separator( - read_set.cell_barcode_segments().map(|s| s.seq.as_slice()), - b'-', - ); - let cell_quals_bs = Self::join_bytes_with_separator( - read_set.cell_barcode_segments().map(|s| s.quals.as_slice()), - b' ', - ); - let sample_barcode_bs = Self::join_bytes_with_separator( - read_set.sample_barcode_segments().map(|s| s.seq.as_slice()), - b'-', - ); - let sample_quals_bs = Self::join_bytes_with_separator( - read_set.sample_barcode_segments().map(|s| s.quals.as_slice()), - b' ', - ); - let umi_bs = Self::join_bytes_with_separator( - read_set.molecular_barcode_segments().map(|s| s.seq.as_slice()), - b'-', - ); - let umi_qual_bs = Self::join_bytes_with_separator( - read_set.molecular_barcode_segments().map(|s| s.quals.as_slice()), - b' ', - ); - - // Extract UMI from read name if requested - let (read_name_bytes, umi_from_name) = - Self::extract_read_name_and_umi(&read_set.header, self.extract_umis_from_read_names)?; - - // Prepare final UMI as BString - avoid format! and unnecessary allocations - let final_umi_bs: BString = match (umi_bs.is_empty(), &umi_from_name) { - (true, Some(from_name)) => BString::from(from_name.as_slice()), - (true, None) => BString::default(), - (false, Some(from_name)) => { - let mut combined = Vec::with_capacity(from_name.len() + 1 + umi_bs.len()); - combined.extend_from_slice(from_name); - combined.push(b'-'); - combined.extend_from_slice(umi_bs.as_bytes()); - BString::from(combined) - } - (false, None) => umi_bs, - }; - - let num_templates = templates.len(); - - for (index, template) in templates.iter().enumerate() { - // Compute flags for unmapped reads - let mut flag = flags::UNMAPPED; - if num_templates == 2 { - flag |= flags::PAIRED | flags::MATE_UNMAPPED; - if index == 0 { - flag |= flags::FIRST_SEGMENT; - } else { - flag |= flags::LAST_SEGMENT; - } - } - - // Set read name (optionally with UMI annotation) - let annotated_name: Option> = if self.annotate_read_names - && !final_umi_bs.is_empty() - { - let mut name = Vec::with_capacity(read_name_bytes.len() + 1 + final_umi_bs.len()); - name.extend_from_slice(&read_name_bytes); - name.push(b'+'); - name.extend_from_slice(final_umi_bs.as_bytes()); - Some(name) - } else { - None - }; - let final_read_name: &[u8] = annotated_name.as_deref().unwrap_or(&read_name_bytes); - - Self::build_template_record( - builder, - final_read_name, - flag, - &template.seq, - &template.quals, - encoding, - )?; - - // Append tags - // Read group - builder.append_string_tag(SamTag::RG, self.read_group_id.as_bytes()); - - // Cell barcode - if !cell_barcode_bs.is_empty() { - builder.append_string_tag(SamTag::CB, cell_barcode_bs.as_bytes()); - } - - if !cell_quals_bs.is_empty() && self.store_cell_quals { - builder.append_string_tag(SamTag::CY, cell_quals_bs.as_bytes()); - } - - // Sample barcode - if !sample_barcode_bs.is_empty() { - builder.append_string_tag(SamTag::BC, sample_barcode_bs.as_bytes()); - } - - if self.store_sample_barcode_qualities && !sample_quals_bs.is_empty() { - builder.append_string_tag(SamTag::QT, sample_quals_bs.as_bytes()); - } - - // UMI - if !final_umi_bs.is_empty() { - builder.append_string_tag(SamTag::RX, final_umi_bs.as_bytes()); - - // Single tag for all concatenated UMIs (if specified) - if let Some(single_tag) = self.single_tag { - builder.append_string_tag(single_tag, final_umi_bs.as_bytes()); - } - - // Only add UMI qualities if not extracted from read names - if umi_from_name.is_none() && !umi_qual_bs.is_empty() && self.store_umi_quals { - builder.append_string_tag(SamTag::QX, umi_qual_bs.as_bytes()); - } - } - - // Write the record directly - writer.write_raw_record(builder.as_bytes())?; - builder.clear(); - } - - Ok(num_templates as u64) - } - - /// Process records in single-threaded mode - /// - /// Returns the number of records written. - fn process_singlethreaded( - &self, - fq_iterators: &mut [ReadSetIterator], - writer: &mut RawBamWriter, - encoding: QualityEncoding, - ) -> Result { - let progress = ProgressTracker::new("Processed records").with_interval(1_000_000); - let mut read_pair_count: u64 = 0; - let mut builder = UnmappedSamBuilder::new(); - - // One record per stream at this position; `1` = a record was read, `0` = - // the stream ended. Naming the streams (rather than a vague "some ended - // before others") matches the threaded pipeline and the help text's - // promise to name which input ended first. Hoisted out of the loop and - // reset each iteration to avoid a per-template allocation on this fast path. - let mut counts = vec![0usize; fq_iterators.len()]; - loop { - counts.fill(0); - let mut next_read_sets = Vec::with_capacity(fq_iterators.len()); - for (idx, iter) in fq_iterators.iter_mut().enumerate() { - match iter.next() { - Some(Ok(rec)) => { - counts[idx] = 1; - next_read_sets.push(rec); - } - Some(Err(e)) => return Err(e), - None => {} // stream ended; leave its count at 0 - } - } - - // Every stream ended together: clean EOF. - if counts.iter().all(|count| *count == 0) { - break; - } - // Some streams still had a record while others ended: out of sync. - if counts.iter().any(|count| *count != counts[0]) { - return Err(fastq_out_of_sync_error(&counts).into()); - } - - // Validate read names match across all FASTQs - Self::validate_read_names_match(&next_read_sets)?; - - let read_set = FastqSet::combine_readsets(next_read_sets); - let num_records = self.make_raw_records(&read_set, encoding, &mut builder, writer)?; - - read_pair_count += num_records; - progress.log_if_needed(num_records); - } - - progress.log_final(); - Ok(read_pair_count) - } } /// Pool quality strings for encoding detection across the first @@ -1402,11 +1058,9 @@ fn sample_detection_quals_from_stream( /// Open the decompressed FASTQ input readers for a run. /// /// De-interleaves a single interleaved stream into an R1/R2 pair and consumes a -/// pre-sampled stdin reader when present. Both the serial oracle and the chain -/// reach it through [`detect_encoding_and_open_fastq_readers`] — the serial -/// `Extract::execute` path directly, and the chain FASTQ source via -/// `ChainBuilder::open_source` — so both open readers identically, the parity -/// contract for the `--threads` cutover. +/// pre-sampled stdin reader when present. Reached through +/// [`detect_encoding_and_open_fastq_readers`] from the chain FASTQ source +/// (`ChainBuilder::open_source`), the only execution path for extract. fn open_fastq_input_readers( inputs: &[PathBuf], interleaved: bool, @@ -1447,9 +1101,9 @@ fn open_fastq_input_readers( // workers would spawn up to N×threads decode threads ON TOP OF the pipeline's // own worker pool — over-subscription the pre-cutover chain deliberately // avoided (it passed 1 here, "the pipeline framework provides the - // parallelism"). The serial oracle reaches this branch with `decomp_threads` - // == 1 (it is single-threaded), so hardcoding 1 leaves it unchanged; per-file - // BGZF-decode parallelism, if wanted, is a separate benchmarked change. + // parallelism"). A no-`--threads` run reaches this branch through the chain at + // a single worker, so hardcoding 1 leaves it unchanged; per-file BGZF-decode + // parallelism, if wanted, is a separate benchmarked change. inputs .iter() .map(|path| open_fastq_reader(path, 1, async_reader, check_crc, no_check_crc)) @@ -1457,9 +1111,8 @@ fn open_fastq_input_readers( } /// Detect the input FASTQ quality encoding and open the decompressed input -/// readers in one pass, mirroring `Extract::execute`'s detection + reader-open -/// exactly so the serial oracle and the `--threads` chain path see identical -/// readers and encoding. +/// readers in one pass. Called by the chain FASTQ source +/// (`ChainBuilder::open_source`), the only execution path for extract. /// /// stdin is sampled once and its bytes replayed into the returned reader (there /// is no second open); file inputs are sampled through separate readers so the @@ -1505,69 +1158,20 @@ pub(crate) fn detect_encoding_and_open_fastq_readers( impl Command for Extract { fn execute(&self, command_line: &str) -> Result<()> { - // Validate inputs + // Validate inputs. self.validate()?; - // `--threads N` routes onto the declarative chain builder — the chain - // opens its own FASTQ readers and detects the quality encoding inside - // `ChainBuilder::open_source`. The no-`--threads` serial loop below is the - // in-process parity oracle. `--threads 1` takes the chain like any other - // `Some(n)`: the old dispatch was `is_parallel()`, which was already true - // for every `Threads(n)` including `Threads(1)`, so this - // `threads.is_some()` predicate is behaviorally identical — the cutover - // swaps the parallel *implementation* (the now-deleted legacy - // `process_with_pipeline`) for the chain builder, it does not change which - // `--threads` values run in parallel. - if self.threading.threads.is_some() { - return self.execute_chain(command_line); - } - - let timer = OperationTimer::new("Extracting UMIs"); - let read_structures = self.get_read_structures()?; - - // Serial oracle: detect the quality encoding and open the input readers - // in one pass — the same helper the chain FASTQ source uses, so both - // paths open identical readers (de-interleaving / stdin-replay included) - // and apply the same detected encoding. - let (fq_readers, encoding) = detect_encoding_and_open_fastq_readers( - &self.inputs, - self.interleaved, - self.threading.num_threads(), - self.async_reader, - self.check_crc, - self.no_check_crc, - )?; - - // Create header with @PG record - let header = self.create_header(command_line)?; - - let fq_sources: Vec>> = fq_readers - .into_iter() - .map(|fq| SimdFastqReader::with_capacity(fq, BUFFER_SIZE)) - .collect(); - - let mut fq_iterators: Vec = fq_sources - .into_iter() - .zip(read_structures.iter()) - .map(|(source, rs)| ReadSetIterator::new(rs.clone(), source, Vec::new())) - .collect(); - - // Raw BAM writer (single-threaded oracle output). - let writer_threads = self.threading.num_threads(); - let mut writer = create_raw_bam_writer( - &self.output, - &header, - writer_threads, - self.compression.compression_level, - )?; - - let count = self.process_singlethreaded(&mut fq_iterators, &mut writer, encoding)?; - - // Flush and finish the writer so any write error surfaces here, not on drop. - writer.finish()?; - - timer.log_completion(count); - Ok(()) + // The declarative chain is the only execution path: `execute` always + // dispatches to `execute_chain`, with or without `--threads` (absent + // `--threads` runs the chain at a single worker). The chain opens its own + // FASTQ readers and detects the quality encoding inside + // `ChainBuilder::open_source`, and emits every user-facing diagnostic — + // the CRC-verify log, the `Starting Extract` banner + Input/Output lines, + // the `Extracting UMIs` `OperationTimer`, and the records-emitted summary + // — inside `ChainBuilder::add_extract` and its finalize hook. 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. + self.execute_chain(command_line) } } @@ -1580,11 +1184,12 @@ impl Command for Extract { /// /// - **Header-synthesis fields** (`sample`, `library`, `platform`, /// `platform_unit`, `read_group_id`, `comments`) map directly to the -/// corresponding `@RG` and `@CO` entries written by `Extract::create_header`. +/// corresponding `@RG` and `@CO` entries written by the chain's +/// `build_fastq_header`. /// - **Behavior options** control tag output and name annotation in the same /// way as the identically-named [`Extract`] CLI flags. /// -/// Constructed by `Extract::execute` (T5.7) from the parsed CLI struct and +/// Constructed by `Extract::to_extract_options` from the parsed CLI struct and /// placed into [`crate::pipeline::chains::StageOptionsBag::extract`]. #[derive(Debug, Clone)] #[allow(clippy::struct_excessive_bools)] @@ -1621,9 +1226,8 @@ pub struct ExtractOptions { /// On the chain path this is a placeholder at spec-construction time: the /// FASTQ source detects the encoding while opening its readers (in /// `ChainBuilder::open_source`) and overrides this field before the extract - /// step runs, so both the serial oracle and the chain apply the same - /// source-detected encoding. Tests that construct `ExtractOptions` directly - /// set the encoding they want to exercise. + /// step runs, so the chain applies the source-detected encoding. Tests that + /// construct `ExtractOptions` directly set the encoding they want to exercise. pub quality_encoding: QualityEncoding, /// Store UMI base qualities in the `QX` tag. pub store_umi_quals: bool, @@ -1642,7 +1246,7 @@ pub struct ExtractOptions { pub async_reader: bool, /// Force BGZF input CRC verification (`--check-crc`). Consumed when the chain /// FASTQ source opens its readers, so `--check-crc`/`--no-check-crc` reach - /// `open_fastq_reader` on the `--threads` path exactly as on the serial path. + /// `open_fastq_reader` on the chain path. pub check_crc: bool, /// Disable BGZF input CRC verification (`--no-check-crc`). See `check_crc`. pub no_check_crc: bool, @@ -1655,17 +1259,15 @@ pub struct ExtractOptions { /// tags, and produces one `RawRecord` per template segment. The caller wraps /// the result in a [`crate::template::Template`] for the typed-step pipeline. /// -/// KNOWN DUPLICATION: this is the chain path's copy of the FASTQ→`RawRecord` -/// per-read loop, near-identical to `Extract::make_raw_records` (the serial -/// oracle's, which writes to a `RawBamWriter`). The two differ only in output -/// target, so tag-extraction/flag logic can drift between them. The -/// `chain_matches_serial_oracle*` integration tests guard against drift by -/// comparing chain-vs-oracle output: the base case covers RX/RG, and -/// `chain_matches_serial_oracle_barcode_and_annotate_tags` covers the -/// CB/CY, BC/QT, QX, `--single-tag`, and `--annotate-read-names` branches. -/// Unifying the two onto one shared per-read loop is a worthwhile follow-up. -/// (The former third copy, `make_raw_records_static`, was deleted with the -/// legacy threaded pipeline when extract cut over to the chain builder.) +/// This is the single per-read FASTQ→`RawRecord` builder for extract: the chain +/// is the only execution path, so there is no second copy to keep in sync (the +/// legacy serial `make_raw_records` and the earlier `make_raw_records_static` +/// were both removed with the single-threaded path). Its output is pinned by the +/// `no_threads_matches_threaded*` integration tests, which assert byte-identical +/// output across worker counts and pin the emitted tag values: the base case +/// covers RX/RG, and `no_threads_matches_threaded_barcode_and_annotate_tags` +/// covers the CB/CY, BC/QT, QX, `--single-tag`, and `--annotate-read-names` +/// branches. /// /// # Errors /// @@ -1758,7 +1360,7 @@ pub(crate) fn make_raw_records_from_fastq_set( // derived from the input FASTQ header (plus an optional `+`), so an // over-long (>=255-byte) name is bad input, not a bug — fail cleanly with // context instead of panicking partway through the batch and truncating - // the output BAM. Mirrors the standalone path's `build_template_record`. + // the output BAM. if template.seq.is_empty() { builder.try_build_record(final_read_name, flag, b"N", &[2u8]) } else { @@ -3207,11 +2809,13 @@ mod tests { } #[test] - #[should_panic(expected = "Read names do not match")] + #[should_panic(expected = "FASTQ read name mismatch")] fn test_extract_fails_on_underscore_read_number_suffix() { // A `read_1`/`read_2` pair previously passed validation (`_1`/`_2` were stripped) - // yet was written with two different QNAMEs, violating the SAM spec. Validation - // now only strips `/`, so the mismatch is reported instead. + // yet was written with two different QNAMEs, violating the SAM spec. Only `/` + // is stripped, so the mismatch is reported instead. The rejection comes from the + // chain's FASTQ zip step (the only execution path), which names the two disagreeing + // read names. let tmp = TempDir::new().expect("failed to create temp dir"); let r1 = create_fastq(&tmp, "r1.fq", &[("read_1", "AAAAAAAAAA", "==========")]); let r2 = create_fastq(&tmp, "r2.fq", &[("read_2", "CCCCCCCCCC", "##########")]); @@ -3254,7 +2858,7 @@ mod tests { } #[rstest] - #[case::fast_path(ThreadingOptions::none())] + #[case::single_worker(ThreadingOptions::none())] #[case::threaded(ThreadingOptions::new(2))] fn test_extract_strips_read_number_suffix_from_written_qname( #[case] threading: ThreadingOptions, @@ -3262,10 +2866,10 @@ mod tests { // End-to-end: paired FASTQs with old-style `/1` and `/2` read-number // suffixes must produce a single shared QNAME (no suffix) for both mates, // as required by the SAM spec. Guards the wiring from - // `extract_read_name_and_umi` through to the written BAM QNAME in both the - // serial-oracle (`make_raw_records`) and the chain + // `extract_read_name_and_umi` through to the written BAM QNAME in the chain // (`make_raw_records_from_fastq_set` + `strip_read_suffix`) record-building - // paths — this test is parameterized over threading to exercise both. + // path — parameterized over worker count (no-`--threads` single worker and + // `--threads 2`) to exercise both concurrency configurations. let tmp = TempDir::new().expect("failed to create temp dir"); let r1 = create_fastq(&tmp, "r1.fq", &[("SRR001.1/1", "AAAAAAAAAA", "==========")]); let r2 = create_fastq(&tmp, "r2.fq", &[("SRR001.1/2", "CCCCCCCCCC", "##########")]); @@ -3320,9 +2924,9 @@ mod tests { /// Rust panic partway through writing the output BAM and leave a truncated /// file behind. /// - /// Both record-building paths are covered (the test is parameterized over - /// threading): `make_raw_records` (serial oracle) and - /// `make_raw_records_from_fastq_set` (chain, `--threads`). + /// The chain record-builder (`make_raw_records_from_fastq_set`) is exercised + /// across worker counts (the test is parameterized over threading: + /// no-`--threads` single worker and `--threads 2`). #[rstest] #[case::at_limit(254, true)] #[case::one_over(255, false)] @@ -3598,7 +3202,7 @@ mod tests { } #[test] - #[should_panic(expected = "Read names do not match")] + #[should_panic(expected = "FASTQ read name mismatch")] fn test_fail_mismatched_read_names() { let tmp = TempDir::new().expect("failed to create temp dir"); let r1 = create_fastq(&tmp, "r1.fq", &[("q1", "AAAAAAAAAA", "==========")]); @@ -4910,11 +4514,11 @@ mod tests { assert_eq!(format, CompressionFormat::Bgzf); } - /// Test that extraction works correctly across all threading modes. - /// This parameterized test ensures both the single-threaded fast path (None) - /// and the multi-threaded pipeline (Some(1), Some(2)) produce correct results. + /// Test that extraction works correctly across all worker-count modes. + /// This parameterized test ensures the chain produces correct results at a + /// single worker (no `--threads`) and at `--threads 1` / `--threads 2`. #[rstest] - #[case::fast_path(ThreadingOptions::none())] + #[case::single_worker(ThreadingOptions::none())] #[case::pipeline_1(ThreadingOptions::new(1))] #[case::pipeline_2(ThreadingOptions::new(2))] fn test_threading_modes(#[case] threading: ThreadingOptions) -> Result<()> { @@ -4977,11 +4581,11 @@ mod tests { Ok(()) } - /// Verifies that the chain (`--threads`) path emits the same tags (RX, QX, RG) as the - /// serial oracle, ensuring `make_raw_records_from_fastq_set` stays in sync with - /// `make_raw_records`. + /// Verifies that the chain emits the correct tags (RX, QX, RG) regardless of + /// worker count, exercising `make_raw_records_from_fastq_set` under both a + /// no-`--threads` single worker and `--threads 2`. #[rstest] - #[case::fast_path(ThreadingOptions::none())] + #[case::single_worker(ThreadingOptions::none())] #[case::threaded(ThreadingOptions::new(2))] fn test_threading_modes_emit_correct_tags(#[case] threading: ThreadingOptions) -> Result<()> { let tmp = TempDir::new()?; @@ -5088,12 +4692,12 @@ mod tests { } /// A single interleaved input (`--interleaved`) must produce exactly the same - /// extracted records as the equivalent pair of separate R1/R2 FASTQ files, on - /// both the single-threaded and the multi-threaded pipeline paths. Distinct R1 - /// and R2 read structures (a 2 bp inline UMI on R1, template-only R2) prove the - /// split routes R1 → structure[0] and R2 → structure[1]. + /// extracted records as the equivalent pair of separate R1/R2 FASTQ files, + /// across worker counts (no-`--threads` single worker and `--threads 2`). + /// Distinct R1 and R2 read structures (a 2 bp inline UMI on R1, template-only + /// R2) prove the split routes R1 → structure[0] and R2 → structure[1]. #[rstest] - #[case::single_threaded(ThreadingOptions::none())] + #[case::single_worker(ThreadingOptions::none())] #[case::multi_threaded(ThreadingOptions::new(2))] fn interleaved_input_matches_two_file_input(#[case] threading: ThreadingOptions) { let tmp = TempDir::new().expect("failed to create temp dir"); @@ -5348,9 +4952,10 @@ mod tests { /// `open_fastq_reader` and `deinterleave` splits the *decompressed* bytes into /// R1/R2, so BGZF is never decoded twice or line-split as raw compressed /// bytes. The BGZF and plaintext interleaved runs must produce the same - /// records, on both the serial and chain paths. + /// records, across worker counts (no-`--threads` single worker and + /// `--threads 2`). #[rstest] - #[case::single_threaded(ThreadingOptions::none())] + #[case::single_worker(ThreadingOptions::none())] #[case::multi_threaded(ThreadingOptions::new(2))] fn interleaved_bgzf_input_matches_plaintext(#[case] threading: ThreadingOptions) { let tmp = TempDir::new().expect("temp dir"); @@ -5472,9 +5077,9 @@ mod tests { } /// `Extract::to_extract_options` must project every CLI field onto the - /// matching `ExtractOptions` field so the chain path sees the same options - /// the serial oracle uses. `quality_encoding` is a placeholder (`Standard`) - /// because the chain FASTQ source overrides it with the detected encoding. + /// matching `ExtractOptions` field so the chain path sees every option the + /// user set. `quality_encoding` is a placeholder (`Standard`) because the + /// chain FASTQ source overrides it with the detected encoding. #[test] fn to_extract_options_maps_all_fields() { let extract = Extract { diff --git a/src/lib/pipeline/chains/commands/extract.rs b/src/lib/pipeline/chains/commands/extract.rs index b41638984..13cc4f816 100644 --- a/src/lib/pipeline/chains/commands/extract.rs +++ b/src/lib/pipeline/chains/commands/extract.rs @@ -4,12 +4,12 @@ //! helper that synthesizes an unmapped-BAM `@HD`/`@RG`/`@CO` header from //! [`ExtractOptions`], both consumed by `ChainBuilder::add_extract`. //! -//! The header builder reproduces the same `@HD`, `@RG`, and `@CO` records -//! that [`Extract::create_header`] emits, minus the `@PG` record — -//! [`ChainBuilder::new`] adds that uniformly for all chains. +//! `build_fastq_header` is the only unmapped-BAM header synthesizer for extract +//! (the chain is extract's only execution path). It emits the `@HD`, `@RG`, and +//! `@CO` records minus the `@PG` record — [`ChainBuilder::new`] adds that +//! uniformly for all chains. //! //! [`ExtractOptions`]: crate::commands::extract::ExtractOptions -//! [`Extract::create_header`]: crate::commands::extract::Extract //! [`ChainBuilder::new`]: crate::pipeline::chains::builder::ChainBuilder use std::sync::Arc; @@ -59,13 +59,13 @@ impl FinalizeHook for ExtractFinalizeHook { } // ───────────────────────────────────────────────────────────────────────────── -// build_fastq_header — reproduce Extract::create_header sans @PG +// build_fastq_header — the unmapped-BAM header synthesizer (sans @PG) // ───────────────────────────────────────────────────────────────────────────── /// Conditionally add a read-group tag value to a [`Builder`]. /// /// If `value` is `Some`, inserts the tag with the value. Otherwise returns -/// the builder unchanged. Mirrors `Extract::add_to_read_group`. +/// the builder unchanged. fn add_to_read_group( rg: Builder, tag: noodles::sam::header::record::value::map::tag::Other, @@ -84,11 +84,10 @@ fn add_to_read_group( /// The header does **not** include a `@PG` record — [`ChainBuilder::new`] /// injects that uniformly for every chain. /// -/// This reproduces the header structure of [`Extract::create_header`] so the -/// two code paths produce byte-identical headers (modulo `@PG`). +/// This is the sole header synthesizer for extract (the chain is extract's only +/// execution path). /// /// [`ChainBuilder::new`]: crate::pipeline::chains::builder::ChainBuilder -/// [`Extract::create_header`]: crate::commands::extract::Extract pub(crate) fn build_fastq_header(extract_opts: &ExtractOptions) -> Result
{ let mut header = Header::builder(); @@ -111,11 +110,11 @@ pub(crate) fn build_fastq_header(extract_opts: &ExtractOptions) -> Result
, out: &Path| { let mut args: Vec = vec![ @@ -770,10 +772,11 @@ fn chain_matches_serial_oracle_barcode_and_annotate_tags(#[case] threads: usize) assert_eq!(string_tag(rec, SamTag::BC), b"CCC", "sample barcode"); assert_eq!(string_tag(rec, SamTag::RX), b"GGG", "UMI"); assert_eq!(string_tag(rec, SamTag::MI), b"GGG", "--single-tag copy of the UMI"); - // The store_* flags are on, so the quality tags must be present. - assert!(rec.data().get(&Tag::from(SamTag::CY)).is_some(), "cell-barcode quals (CY)"); - assert!(rec.data().get(&Tag::from(SamTag::QT)).is_some(), "sample-barcode quals (QT)"); - assert!(rec.data().get(&Tag::from(SamTag::QX)).is_some(), "UMI quals (QX)"); + // The store_* flags are on, so the quality tags carry the per-segment quality + // slices verbatim (the input is unambiguously Phred+33, so no encoding shift). + assert_eq!(string_tag(rec, SamTag::CY), b"+,-", "cell-barcode quals (CY)"); + assert_eq!(string_tag(rec, SamTag::QT), b"./0", "sample-barcode quals (QT)"); + assert_eq!(string_tag(rec, SamTag::QX), b"123", "UMI quals (QX)"); // --annotate-read-names appends `+` to the read name. let name: &[u8] = rec.name().expect("read name must be present").as_ref(); assert_eq!(name, &b"readX+GGG"[..], "annotated read name"); @@ -847,11 +850,11 @@ fn test_parallel_parse_determinism() { // Unified Pipeline Path Tests // ============================================================================ -/// Test BGZF+sync: multithreaded output matches single-threaded content. +/// Test BGZF+sync: multi-worker output matches single-worker content. /// -/// This verifies the new BGZF+synchronized code path (which didn't exist before -/// the unified pipeline) produces correct output by comparing against the -/// single-threaded fast-path result. +/// This verifies the BGZF+synchronized code path produces correct output by +/// comparing a `--threads 1` (single-worker) chain run against a multi-worker +/// chain run over BGZF-compressed input. #[test] fn test_bgzf_sync_multithreaded_matches_single_threaded() { let tmp = TempDir::new().unwrap(); @@ -870,7 +873,7 @@ fn test_bgzf_sync_multithreaded_matches_single_threaded() { let r1 = create_bgzf_fastq(&tmp, "r1.fq.bgz", &records_r1); let r2 = create_bgzf_fastq(&tmp, "r2.fq.bgz", &records_r2); - // Run single-threaded (fast-path) + // Run at a single worker (--threads 1) let output_st = tmp.path().join("output_st.bam"); let cmd = Extract::try_parse_from([ "extract", @@ -892,7 +895,7 @@ fn test_bgzf_sync_multithreaded_matches_single_threaded() { "1", ]) .expect("failed to parse extract args"); - cmd.execute("fgumi extract").expect("Failed to execute single-threaded extract"); + cmd.execute("fgumi extract").expect("Failed to execute single-worker extract"); // Run multithreaded (BGZF+sync through unified pipeline) let output_threaded = tmp.path().join("output_mt.bam"); @@ -1600,9 +1603,8 @@ fn run_extract_pair(r1: &Path, r2: &Path, output: &Path, threads: usize) -> anyh cmd.execute("fgumi extract") } -/// Run `fgumi extract` over a paired FASTQ input on the single-threaded fast -/// path (no `--threads`), returning the command result. This path runs -/// `process_singlethreaded` rather than the 7-step pipeline. +/// Run `fgumi extract` over a paired FASTQ input with no `--threads` flag (the +/// chain at a single worker), returning the command result. fn run_extract_pair_single_threaded(r1: &Path, r2: &Path, output: &Path) -> anyhow::Result<()> { let cmd = Extract::try_parse_from([ "extract", @@ -1852,11 +1854,11 @@ fn test_extract_rejects_mismatched_fastq_pair(#[case] flavor: FastqFlavor, #[cas } } -/// The single-threaded fast path (no `--threads`) must reject a mismatched pair -/// with a message that names which stream ended first, matching the threaded -/// pipeline and the help text's promise (#773). `--threads N` (even `N == 1`) -/// runs the pipeline, so `process_singlethreaded` is reachable only with no -/// `--threads` flag at all — the other rejection test never exercises it. +/// A no-`--threads` run (the chain at a single worker) must reject a mismatched +/// pair with a message that names which stream ended first, matching the +/// multi-worker pipeline and the help text's promise (#773). This pins the +/// directional out-of-sync wording on the no-`--threads` invocation specifically +/// — the other rejection test always passes `--threads N`. #[rstest] #[case::plain(FastqFlavor::Plain)] #[case::gzip(FastqFlavor::Gzip)] diff --git a/tests/integration/test_extract_cutover_parity.rs b/tests/integration/test_extract_cutover_parity.rs new file mode 100644 index 000000000..82ced729e --- /dev/null +++ b/tests/integration/test_extract_cutover_parity.rs @@ -0,0 +1,412 @@ +//! Parity gate for the `extract` command's single-threaded-path retirement (C4): +//! `Extract::execute` no longer has a serial in-process 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** (`extract_no_threads_routes_through_chain`). +//! A no-`--threads` run now emits the chain's `"Starting Extract"` banner (and +//! the `Input:`/`Output:` lines) that `ChainBuilder::add_extract` logs and the +//! retired serial tail never did. This is the genuine RED/GREEN discriminator: +//! before the removal a no-`--threads` run took the serial path and logged only +//! the `Extracting UMIs` timer + `Processed records` progress; after it, the +//! chain logs the banner. +//! +//! 2. **Output parity with the pre-removal serial path** +//! (`cutover_matches_baseline`). The current build's `extract` output — records +//! (byte-identical, modulo the `@PG` line) — must match the frozen 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 (every record's template sequence and combined `RX`/`RG` tags +//! asserted directly) rather than skipping — the exact fallback discipline of +//! `test_sort_cutover_parity.rs`. Cases +//! cover the three extract input shapes the chain must support: two-file input, +//! interleaved input (`-p/--interleaved`), and BGZF input with `--check-crc`, +//! each over a read structure carrying a UMI segment (`5M+T`). +//! +//! **Not a RED/GREEN gate for the parity half.** Because the removed serial loop +//! and the chain were already output-equivalent, the baseline byte-parity check +//! passes on both sides of the change — like the sort/retag cutovers it guards +//! equivalence, it does not observe a regression the cutover introduces. The +//! banner check in (1) is the part that flips RED→GREEN across the removal. + +use std::ffi::OsStr; +use std::fs::File; +use std::io::Write; +use std::path::{Path, PathBuf}; +use std::process::Command; + +use noodles::sam::alignment::RecordBuf; +use noodles::sam::alignment::record::data::field::Tag; +use noodles::sam::alignment::record_buf::data::field::Value; +use noodles_bgzf::io::Writer as BgzfWriter; +use rstest::rstest; +use tempfile::TempDir; + +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 = 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 +/// (extract synthesizes a fresh header with exactly one `@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. extract's contract is the +/// exact UMI/quality/tag bytes it writes, so that masking risk is what must be +/// avoided. +fn decompressed_records_without_pg(path: &Path) -> Vec { + let file = 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 +} + +/// A pair of `(name, sequence, quality)` FASTQ records: R1 carries UMI `AAAAA` +/// (template `CCCCC`), R2 carries UMI `GGGGG` (template `TTTTT`), so a `5M+T` +/// read structure yields template `CCCCC`/`TTTTT` and a combined `RX` of +/// `AAAAA-GGGGG` on both mates. +/// +/// Every quality string is **position-distinct** and differs across the four +/// records, so the self-consistency oracle's per-base quality check +/// (`assert_self_consistent`) has teeth: a transposed base, a swapped mate, or a +/// dropped/duplicated quality byte lands on a different Phred value and fails. +/// All bytes are unambiguous Phred+33 printable ASCII (`0`=Q15 … `T`=Q51). +const R1_RECORDS: &[(&str, &str, &str)] = + &[("pair0", "AAAAACCCCC", "0123456789"), ("pair1", "TTTTTGGGGG", "KLMNOPQRST")]; +const R2_RECORDS: &[(&str, &str, &str)] = + &[("pair0", "GGGGGTTTTT", "9876543210"), ("pair1", "CCCCCAAAAA", "TSRQPONMLK")]; + +/// Write plain-text FASTQ records to `dir/name`, returning the path. +fn write_plain_fastq(dir: &Path, name: &str, records: &[(&str, &str, &str)]) -> PathBuf { + let path = dir.join(name); + let mut file = File::create(&path).unwrap_or_else(|e| panic!("create {name}: {e}")); + for (rname, seq, qual) in records { + writeln!(file, "@{rname}\n{seq}\n+\n{qual}").expect("write fastq record"); + } + path +} + +/// Write BGZF-compressed FASTQ records to `dir/name`, returning the path. Used by +/// the `--check-crc` case, where the input must be BGZF for the CRC policy to bite. +fn write_bgzf_fastq(dir: &Path, name: &str, records: &[(&str, &str, &str)]) -> PathBuf { + let path = dir.join(name); + let file = File::create(&path).unwrap_or_else(|e| panic!("create {name}: {e}")); + let mut writer = BgzfWriter::new(file); + for (rname, seq, qual) in records { + writeln!(writer, "@{rname}\n{seq}\n+\n{qual}").expect("write bgzf fastq record"); + } + writer.finish().expect("finish bgzf fastq"); + path +} + +/// The extract input shape a case exercises. Each is a distinct chain-source +/// topology (`SourceSpec::Fastqs` vs `SourceSpec::InterleavedFastq`) and/or a +/// distinct reader-open policy (`--check-crc`), so all three must retain parity. +#[derive(Clone, Copy, Debug)] +enum InputShape { + /// Two plain FASTQ files (R1 + R2). + TwoFile, + /// One interleaved plain FASTQ file (`-p/--interleaved`). + Interleaved, + /// Two BGZF FASTQ files run with `--check-crc`. + CheckCrcBgzf, +} + +/// Materializes the input FASTQ(s) for `shape` under `dir` and returns +/// `(input_paths, interleaved, check_crc)` ready to pass to `run_extract`. +fn build_inputs(shape: InputShape, dir: &Path) -> (Vec, bool, bool) { + match shape { + InputShape::TwoFile => { + let r1 = write_plain_fastq(dir, "r1.fq", R1_RECORDS); + let r2 = write_plain_fastq(dir, "r2.fq", R2_RECORDS); + (vec![r1, r2], false, false) + } + InputShape::Interleaved => { + // R1/R2 records alternate in one physical stream. + let interleaved: Vec<(&str, &str, &str)> = + R1_RECORDS.iter().zip(R2_RECORDS.iter()).flat_map(|(a, b)| [*a, *b]).collect(); + let path = write_plain_fastq(dir, "interleaved.fq", &interleaved); + (vec![path], true, false) + } + InputShape::CheckCrcBgzf => { + let r1 = write_bgzf_fastq(dir, "r1.fq.gz", R1_RECORDS); + let r2 = write_bgzf_fastq(dir, "r2.fq.gz", R2_RECORDS); + (vec![r1, r2], false, true) + } + } +} + +/// Runs ` extract --inputs --output --read-structures +/// 5M+T [5M+T] --sample s --library l --compression-level 1 [--interleaved] +/// [--check-crc]` (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_extract( + bin: &Path, + inputs: &[PathBuf], + output: &Path, + interleaved: bool, + check_crc: bool, +) -> std::process::Output { + let mut cmd = Command::new(bin); + cmd.env("RUST_LOG", "info"); + cmd.arg("extract").arg("--inputs"); + cmd.args(inputs.iter().map(|p| p.as_os_str())); + cmd.args([OsStr::new("--output"), output.as_os_str(), OsStr::new("--read-structures")]); + // One `5M+T` per read: two structures for a two-read template (two-file and + // interleaved both describe two reads), matching validation's requirement. + if interleaved || inputs.len() == 2 { + cmd.args([OsStr::new("5M+T"), OsStr::new("5M+T")]); + } else { + cmd.arg("5M+T"); + } + cmd.args([ + OsStr::new("--sample"), + OsStr::new("s"), + OsStr::new("--library"), + OsStr::new("l"), + OsStr::new("--compression-level"), + OsStr::new("1"), + ]); + if interleaved { + cmd.arg("--interleaved"); + } + if check_crc { + cmd.arg("--check-crc"); + } + cmd.output().unwrap_or_else(|e| panic!("failed to spawn `{}` extract: {e}", bin.display())) +} + +/// A no-`--threads` run now routes through the declarative chain, which logs the +/// `"Starting Extract"` banner (and `Input:`/`Output:` lines) from +/// `ChainBuilder::add_extract`. The retired serial tail logged no such line — it +/// emitted only the `Extracting UMIs` timer + `Processed records` progress — so +/// this is the RED (pre-removal) → GREEN (post-removal) discriminator. +#[test] +fn extract_no_threads_routes_through_chain() { + let dir = TempDir::new().expect("temp dir"); + let (inputs, interleaved, check_crc) = build_inputs(InputShape::TwoFile, dir.path()); + let output = dir.path().join("out.bam"); + + let current_bin = Path::new(env!("CARGO_BIN_EXE_fgumi")); + let out = run_extract(current_bin, &inputs, &output, interleaved, check_crc); + let stderr = String::from_utf8_lossy(&out.stderr); + assert!(out.status.success(), "no-threads extract must succeed; stderr:\n{stderr}"); + assert!( + stderr.contains("Starting Extract"), + "a no-`--threads` extract must route through the chain (which logs the \ + `Starting Extract` banner); the serial path is retired. stderr:\n{stderr}" + ); +} + +/// Output parity of the post-cutover chain against the pre-removal serial +/// baseline binary, across the three extract input shapes — plus the +/// always-available self-consistency oracle when no baseline is set. +#[rstest] +#[case::two_file(InputShape::TwoFile)] +#[case::interleaved(InputShape::Interleaved)] +#[case::check_crc_bgzf(InputShape::CheckCrcBgzf)] +fn cutover_matches_baseline(#[case] shape: InputShape) { + let dir = TempDir::new().expect("temp dir"); + let (inputs, interleaved, check_crc) = build_inputs(shape, dir.path()); + + let current_bin = Path::new(env!("CARGO_BIN_EXE_fgumi")); + let current_out = dir.path().join("current.bam"); + let current = run_extract(current_bin, &inputs, ¤t_out, interleaved, check_crc); + let current_stderr = String::from_utf8_lossy(¤t.stderr); + assert!( + current.status.success(), + "current extract ({shape:?}) must succeed; stderr:\n{current_stderr}" + ); + + if let Some(baseline) = baseline_bin() { + let baseline_out = dir.path().join("baseline.bam"); + let base = run_extract(&baseline, &inputs, &baseline_out, interleaved, check_crc); + assert!( + base.status.success(), + "baseline extract ({shape:?}) failed; stderr:\n{}", + String::from_utf8_lossy(&base.stderr) + ); + + assert_eq!( + decompressed_records_without_pg(¤t_out), + decompressed_records_without_pg(&baseline_out), + "chain extract output ({shape:?}) diverges from the pre-removal serial baseline \ + binary ({}) after stripping @PG — a real cutover parity bug, not something to relax", + baseline.display(), + ); + } else { + eprintln!( + "SKIP baseline half of cutover_matches_baseline[{shape:?}]: FGUMI_BASELINE_BIN is \ + unset or does not name an existing file — running self-consistency oracle instead" + ); + assert_self_consistent(¤t_out); + } +} + +/// Always-available oracle used when no baseline binary is set (the default in +/// CI): the chain's output must carry exactly the records the known input + +/// `5M+T` read structure prescribe — for *every* input template pair, two mates +/// carrying the template bases (and their qualities) after the 5 bp UMI, named +/// after the input pair, flagged paired/first/last, and both tagged with the +/// combined R1-R2 UMI in `RX` and the default `RG`. Every produced record is +/// checked (full record identity — QNAME, paired/first/last flags, template +/// sequence, per-base template qualities, `RX`, `RG` — not just count), so a +/// per-record divergence in any pair fails the case. The prior version asserted +/// only sequence/`RX`/`RG`, so a wrong QNAME, a swapped first/last flag, or a +/// corrupted quality value slipped through. +fn assert_self_consistent(output: &Path) { + let (_, recs) = read_bam_output(output); + assert_eq!(recs.len(), R1_RECORDS.len() * 2, "two BAM records per input template pair"); + + let string_tag = |rec: &RecordBuf, tag: SamTag| -> Vec { + match rec.data().get(&Tag::from(tag)) { + Some(Value::String(s)) => s.to_vec(), + other => panic!("expected a string {tag:?} tag, got {other:?}"), + } + }; + // The template's Phred qualities are the input quality string past the 5 bp + // UMI, decoded from Phred+33 to the numeric scores `quality_scores()` returns. + let expected_template_qual = + |qual: &str| qual.as_bytes()[5..].iter().map(|q| q - 33).collect::>(); + + // Derive the expected output from the known input `5M+T` structure so the + // oracle reads as a spec: each read splits into a 5 bp UMI and the template + // remainder; the two mates of a pair share `RX = -`, the + // QNAME of the input pair, and paired flags with R1 first / R2 last. + for (pair_idx, (r1, r2)) in R1_RECORDS.iter().zip(R2_RECORDS.iter()).enumerate() { + let (r1_umi, r1_template) = r1.1.split_at(5); + let (r2_umi, r2_template) = r2.1.split_at(5); + let expected_rx = format!("{r1_umi}-{r2_umi}").into_bytes(); + // Both mates of a pair carry the R1 record's name (paired FASTQ input). + let expected_name = r1.0.as_bytes(); + + let r1_out = &recs[pair_idx * 2]; + let r2_out = &recs[pair_idx * 2 + 1]; + assert_eq!( + r1_out.sequence().as_ref(), + r1_template.as_bytes(), + "pair {pair_idx}: R1 template after 5M UMI" + ); + assert_eq!( + r2_out.sequence().as_ref(), + r2_template.as_bytes(), + "pair {pair_idx}: R2 template after 5M UMI" + ); + assert_eq!( + r1_out.quality_scores().as_ref(), + expected_template_qual(r1.2).as_slice(), + "pair {pair_idx}: R1 template qualities after 5M UMI" + ); + assert_eq!( + r2_out.quality_scores().as_ref(), + expected_template_qual(r2.2).as_slice(), + "pair {pair_idx}: R2 template qualities after 5M UMI" + ); + + // Paired flags: both segmented, R1 first / R2 last (never the reverse). + let r1_flags = r1_out.flags(); + let r2_flags = r2_out.flags(); + assert!(r1_flags.is_segmented(), "pair {pair_idx}: R1 must be paired"); + assert!(r2_flags.is_segmented(), "pair {pair_idx}: R2 must be paired"); + assert!( + r1_flags.is_first_segment() && !r1_flags.is_last_segment(), + "pair {pair_idx}: R1 must be the first segment, got {r1_flags:?}" + ); + assert!( + r2_flags.is_last_segment() && !r2_flags.is_first_segment(), + "pair {pair_idx}: R2 must be the last segment, got {r2_flags:?}" + ); + + for (mate, rec) in [("R1", r1_out), ("R2", r2_out)] { + assert_eq!( + rec.name().map(|n| n.to_vec()), + Some(expected_name.to_vec()), + "pair {pair_idx}: {mate} QNAME must be the input pair name" + ); + assert_eq!( + string_tag(rec, SamTag::RX), + expected_rx, + "pair {pair_idx}: {mate} combined R1-R2 UMI in RX" + ); + assert_eq!( + string_tag(rec, SamTag::RG), + b"A", + "pair {pair_idx}: {mate} default read group id" + ); + } + } +}