diff --git a/Cargo.lock b/Cargo.lock index 758960ba0..a5dff182d 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -681,6 +681,7 @@ dependencies = [ "fgumi-simd-fastq", "fgumi-sort", "fgumi-sort-cli", + "fgumi-tag", "fgumi-umi", "flate2", "fs4", diff --git a/Cargo.toml b/Cargo.toml index 4424db348..cf25ebc4a 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -148,6 +148,8 @@ rstest = "0" proptest = "1.10" criterion = { version = "0.8", features = ["html_reports"] } fgumi-raw-bam = { workspace = true, features = ["test-utils", "noodles"] } +fgumi-tag = { workspace = true } +flate2 = { version = "1.1", features = ["zlib-rs"] } noodles = { version = "0.111.0", features = ["bam", "sam"] } trybuild = "1.0" diff --git a/crates/fgumi-bam-io/src/grouping.rs b/crates/fgumi-bam-io/src/grouping.rs index 548af4454..f116688c2 100644 --- a/crates/fgumi-bam-io/src/grouping.rs +++ b/crates/fgumi-bam-io/src/grouping.rs @@ -282,6 +282,14 @@ impl DecodedRecord { self.data.as_ref() } + /// Immutable access to the underlying [`RawRecord`], for read-only + /// consumers (e.g. the FASTQ-encode step) that need the typed accessors + /// (`flags()`, sequence, qualities, tags) rather than the raw byte slice. + #[must_use] + pub fn record(&self) -> &RawRecord { + &self.data + } + /// Mutable access to the underlying [`RawRecord`]. Used by mid-chain /// `Parallel` Process steps that need to mutate record bytes (e.g., /// MQ bumping, tag rewriting) without rebuilding the `DecodedRecord` diff --git a/crates/fgumi-bam-io/src/lib.rs b/crates/fgumi-bam-io/src/lib.rs index 80e5fa2ab..5940d645c 100644 --- a/crates/fgumi-bam-io/src/lib.rs +++ b/crates/fgumi-bam-io/src/lib.rs @@ -40,6 +40,6 @@ pub use reader::{ pub use reorder::{DrainReady, ReorderBuffer}; pub use writer::{ BamWriter, BgzfWriterEnum, RawBamWriter, bai_sidecar_path, create_bam_writer, - create_optional_bam_writer, create_raw_bam_writer, write_bai_index, write_bai_sidecar, - write_bam_header, + create_bgzf_writer, create_optional_bam_writer, create_raw_bam_writer, write_bai_index, + write_bai_sidecar, write_bam_header, }; diff --git a/crates/fgumi-bam-io/src/writer.rs b/crates/fgumi-bam-io/src/writer.rs index 089859ac8..3cfbbc8a9 100644 --- a/crates/fgumi-bam-io/src/writer.rs +++ b/crates/fgumi-bam-io/src/writer.rs @@ -112,6 +112,25 @@ pub(crate) fn make_bgzf_writer( } } +/// Create a BGZF-compressed writer over an arbitrary `output`, selecting a +/// multi-threaded encoder when `threads > 1`. +/// +/// This exposes the same BGZF backend fgumi uses for BAM output so other output +/// formats (e.g. gzipped FASTQ) get identical, block-gzip (BGZF) framing and +/// multi-threaded compression rather than a separate implementation. +/// +/// `compression_level` is the BGZF/DEFLATE level (0–12, where 0 means +/// uncompressed); values outside the valid range fall back to the backend +/// default. +#[must_use] +pub fn create_bgzf_writer( + output: Box, + threads: usize, + compression_level: u32, +) -> BgzfWriterEnum { + make_bgzf_writer(output, threads, compression_level) +} + /// Write a BAM header (magic, SAM text, reference sequences) to any writer. /// /// Shared implementation used by the BAM writer factories (e.g. [`RawBamWriter`]). diff --git a/crates/fgumi-pipeline-core/src/handles.rs b/crates/fgumi-pipeline-core/src/handles.rs index 053d42b92..7804fdca7 100644 --- a/crates/fgumi-pipeline-core/src/handles.rs +++ b/crates/fgumi-pipeline-core/src/handles.rs @@ -1031,6 +1031,32 @@ where } } +impl OutputHandles> +where + A: Send + HeapSize + Ordered + 'static, + B: Send + HeapSize + Ordered + 'static, + C: Send + HeapSize + Ordered + 'static, +{ + /// Borrow the typed per-branch view for an ordered + byte-bounded 3-tuple + /// output. Same view type as plain `(A, B, C)` — the view carries only + /// `BranchOutputHandle`s; the ordering is encoded in the queue transport. + /// + /// # Panics + /// + /// Panics if the type-erased outputs view doesn't downcast to + /// `Tuple3OutputsView` (a framework invariant violation). + #[must_use] + #[inline] + pub fn view(&self) -> Tuple3View<'_, A, B, C> { + let v = self + .inner + .inner + .downcast_ref::>() + .expect("OrderedBytesTuple3 outputs view downcast failed"); + Tuple3View { a: &v.a, b: &v.b, c: &v.c } + } +} + pub struct Tuple2View<'a, A: Send + HeapSize + 'static, B: Send + HeapSize + 'static> { pub a: &'a BranchOutputHandle, pub b: &'a BranchOutputHandle, @@ -1196,6 +1222,29 @@ where } } +impl OutputHandles> +where + A: Send + HeapSize + Ordered + 'static, + B: Send + HeapSize + Ordered + 'static, + C: Send + HeapSize + Ordered + 'static, +{ + /// Mark all output branches drained. + /// + /// # Panics + /// + /// Panics if the type-erased outputs view doesn't downcast to + /// `Tuple3OutputsView` (the view shape is shared with plain + /// `(A, B, C)`). + pub fn mark_all_drained(&self) { + let view = self + .inner + .inner + .downcast_ref::>() + .expect("OrderedBytesTuple3 outputs view downcast failed"); + view.mark_all_drained(); + } +} + impl OutputHandles<(A, B, C)> where A: Send + HeapSize + 'static, @@ -1402,6 +1451,49 @@ where (queue_set, outputs_view) } +/// Build queues for `OrderedBytesTuple3` outputs. All three branches +/// support the full `QueueSpec × BranchOrdering` cross-product because each is +/// `Ordered + HeapSize`. The `Tuple3OutputsView` carries only +/// `BranchOutputHandle`s (no ordering metadata), so the view type is the same +/// as for plain `(A, B, C)`. +pub(crate) fn build_tuple3_queues_ordered_bytes( + specs: &[QueueSpec], + ordering: &[BranchOrdering], + level: crate::builder::InstrumentationLevel, +) -> (OutputQueueSet, OutputsViewAny) +where + A: Send + HeapSize + Ordered + 'static, + B: Send + HeapSize + Ordered + 'static, + C: Send + HeapSize + Ordered + 'static, +{ + assert_eq!(specs.len(), 3, "OrderedBytesTuple3::build_queues requires 3 specs"); + assert_eq!(ordering.len(), 3, "OrderedBytesTuple3::build_queues requires 3 orderings"); + + let ba = build_branch_ordered_bytes::(specs[0], ordering[0], level); + let bb = build_branch_ordered_bytes::(specs[1], ordering[1], level); + let bc = build_branch_ordered_bytes::(specs[2], ordering[2], level); + let view = Tuple3OutputsView { a: ba.output, b: bb.output, c: bc.output }; + let outputs_view = OutputsViewAny { inner: Box::new(view) }; + let queue_set = OutputQueueSet::new(vec![ + BranchEntry { + input_handle: Box::new(ba.input), + bounded_queue_handle: ba.bounded_queue_handle, + metrics: ba.metrics, + }, + BranchEntry { + input_handle: Box::new(bb.input), + bounded_queue_handle: bb.bounded_queue_handle, + metrics: bb.metrics, + }, + BranchEntry { + input_handle: Box::new(bc.input), + bounded_queue_handle: bc.bounded_queue_handle, + metrics: bc.metrics, + }, + ]); + (queue_set, outputs_view) +} + pub(crate) fn build_tuple4_queues( specs: &[QueueSpec], ordering: &[BranchOrdering], diff --git a/crates/fgumi-pipeline-core/src/outputs.rs b/crates/fgumi-pipeline-core/src/outputs.rs index bce16e4b0..4e86adca2 100644 --- a/crates/fgumi-pipeline-core/src/outputs.rs +++ b/crates/fgumi-pipeline-core/src/outputs.rs @@ -171,6 +171,46 @@ where } } +/// Ordered + byte-bounded tuple-3 outputs. The 3-branch analog of +/// [`OrderedBytesTuple2`]: use when a step fans out to three branches that all +/// need `BranchOrdering::ByItemOrdinal` + byte-bounded queues — e.g. paired +/// FASTQ output splitting one record batch into R1 / R2 / other byte streams, +/// each feeding an ordered `WriteRawFile` sink. +/// +/// The plain `(A, B, C)` `StepOutputs` impl only bounds each branch `HeapSize`. +/// To use `ByItemOrdinal` on any branch, all three must satisfy +/// `Ordered + HeapSize` — that's what this shape encodes at the type level. +#[allow(clippy::type_complexity)] +pub struct OrderedBytesTuple3(PhantomData (A, B, C)>) +where + A: Send + HeapSize + Ordered + 'static, + B: Send + HeapSize + Ordered + 'static, + C: Send + HeapSize + Ordered + 'static; + +impl StepOutputs for OrderedBytesTuple3 +where + A: Send + HeapSize + Ordered + 'static, + B: Send + HeapSize + Ordered + 'static, + C: Send + HeapSize + Ordered + 'static, +{ + #[inline] + fn arity() -> usize { + 3 + } + + fn build_queues( + specs: &[QueueSpec], + ordering: &[BranchOrdering], + level: crate::builder::InstrumentationLevel, + ) -> (OutputQueueSet, OutputsViewAny) { + super::handles::build_tuple3_queues_ordered_bytes::(specs, ordering, level) + } + + fn mark_all_drained(handles: &super::step::OutputHandles) { + handles.mark_all_drained(); + } +} + impl StepOutputs for (A, B, C) where A: Send + HeapSize + 'static, @@ -258,6 +298,26 @@ mod tests { assert_eq!(<(u32, u64, String) as StepOutputs>::arity(), 3); } + #[derive(Clone, Copy)] + struct OrdU64(u64); + impl crate::item::HeapSize for OrdU64 {} + impl crate::item::Ordered for OrdU64 { + fn ordinal(&self) -> u64 { + self.0 + } + } + + #[test] + fn ordered_bytes_tuple_3_arity_is_three() { + // OrderedBytesTuple3 is the ordered + byte-bounded 3-way fan-out shape. + // Its three branches carry Ordered + HeapSize items (here the u64/u32 + // stand-ins just need to satisfy the bounds at the type level). + fn assert_arity() -> usize { + O::arity() + } + assert_eq!(assert_arity::>(), 3); + } + #[test] fn tuple_4_arity_is_four() { assert_eq!(<(u32, u64, String, Vec) as StepOutputs>::arity(), 4); diff --git a/crates/fgumi-pipeline-io/src/sink/mod.rs b/crates/fgumi-pipeline-io/src/sink/mod.rs index 201bc6e29..56d141d03 100644 --- a/crates/fgumi-pipeline-io/src/sink/mod.rs +++ b/crates/fgumi-pipeline-io/src/sink/mod.rs @@ -1 +1,2 @@ pub mod write_bgzf; +pub mod write_raw; diff --git a/crates/fgumi-pipeline-io/src/sink/write_raw.rs b/crates/fgumi-pipeline-io/src/sink/write_raw.rs new file mode 100644 index 000000000..4f01b08c2 --- /dev/null +++ b/crates/fgumi-pipeline-io/src/sink/write_raw.rs @@ -0,0 +1,187 @@ +//! `WriteRawFile` sink step: writes a byte stream verbatim to a file or stdout, +//! with **no** container header and **no** BGZF EOF marker. +//! +//! Unlike [`super::write_bgzf::WriteBgzfFile`] (which is BAM-specific — it emits +//! a BAM header on open and a BGZF EOF on drain), this sink just concatenates +//! the `bytes` of each block it receives. It backs FASTQ output: the chain's +//! FASTQ-encode step produces `DecompressedBlock`s of FASTQ text, which either +//! go straight here (plain output / stdout) or through `BgzfCompress` first +//! (`.gz`/`.bgz` output, producing `BgzfBlock`s — still just bytes to write). +//! +//! `Serial + Affinity::Writer + sticky`, matching `WriteBgzfFile`: exactly one +//! shared instance drains the (reorder-ordered) block stream to the sink. + +use std::fs::File; +use std::io::{self, BufWriter, Write}; +use std::path::Path; + +use parking_lot::Mutex; + +use crate::types::{BgzfBlock, DecompressedBlock}; +use fgumi_pipeline_core::{ + item::HeapSize, + step::{Affinity, Step, StepCtx, StepKind, StepOutcome, StepProfile}, +}; + +/// A pipeline block whose payload is a run of bytes to write verbatim. +pub trait RawBytesBlock: Send + HeapSize + 'static { + /// The bytes to write for this block. + fn bytes(&self) -> &[u8]; +} + +impl RawBytesBlock for DecompressedBlock { + fn bytes(&self) -> &[u8] { + &self.bytes + } +} + +impl RawBytesBlock for BgzfBlock { + fn bytes(&self) -> &[u8] { + &self.bytes + } +} + +/// `Serial + sticky` sink that writes each block's bytes verbatim to a file or +/// stdout. Generic over the block type so it serves both the plain +/// (`DecompressedBlock`) and BGZF-compressed (`BgzfBlock`) FASTQ tails. +pub struct WriteRawFile { + state: Mutex>>>, + /// Bytes appended once, after the last block, on a clean drain. Empty for + /// plain output; the 28-byte BGZF EOF marker for BGZF output so the `.gz` + /// stream is a complete, non-truncated BGZF file. + trailer: &'static [u8], + _marker: std::marker::PhantomData, +} + +impl WriteRawFile { + /// Open `path` for writing (`-` selects stdout), appending `trailer` once + /// after the final block on clean completion. No header is written. + /// + /// # Errors + /// + /// Returns I/O errors from opening the file. + pub fn new>(path: P, trailer: &'static [u8]) -> io::Result { + let inner: Box = if path.as_ref().as_os_str() == "-" { + Box::new(io::stdout()) + } else { + Box::new(File::create(path.as_ref())?) + }; + Ok(Self { + state: Mutex::new(Some(BufWriter::with_capacity(256 * 1024, inner))), + trailer, + _marker: std::marker::PhantomData, + }) + } +} + +impl Step for WriteRawFile { + type Input = B; + type Outputs = (); + + fn profile(&self) -> StepProfile { + StepProfile { + name: "WriteRawFile", + kind: StepKind::Serial, + sticky: true, + output_queues: vec![], + branch_ordering: vec![], + } + } + + fn affinity(&self) -> Affinity { + Affinity::Writer + } + + fn try_run(&mut self, ctx: &mut StepCtx<'_, Self>) -> io::Result { + let mut guard = self.state.lock(); + let Some(out) = guard.as_mut() else { + return Ok(StepOutcome::Finished); + }; + + if let Some(block) = ctx.input.pop() { + out.write_all(block.bytes())?; + return Ok(StepOutcome::Progress); + } + + if ctx.input.is_drained() { + // Clean end-of-stream: append the trailer (e.g. the BGZF EOF marker) + // exactly once, then flush and retire the sink. + if !self.trailer.is_empty() { + out.write_all(self.trailer)?; + } + out.flush()?; + let _ = guard.take(); + return Ok(StepOutcome::Finished); + } + Ok(StepOutcome::NoProgress) + } +} + +impl Drop for WriteRawFile { + /// Cleanup-only drop: flush buffered bytes but deliberately do **not** write + /// `trailer`. The trailer (e.g. the BGZF EOF marker for `.gz` output) is + /// written exclusively by the drained-finish path in `try_run`, which then + /// takes the state so this drop is a no-op for a cleanly-finished sink. If + /// state is still present here the stream was aborted mid-way, and stamping + /// the BGZF EOF onto a truncated `.gz` would make it look complete and hide + /// the truncation from readers — so we withhold it, exactly as + /// `WriteBgzfFile::drop` does. + fn drop(&mut self) { + if let Some(mut out) = self.state.lock().take() { + let _ = out.flush(); + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn raw_bytes_block_exposes_payload() { + let plain = DecompressedBlock { batch_serial: 0, bytes: b"ACGT".to_vec() }; + assert_eq!(RawBytesBlock::bytes(&plain), b"ACGT"); + let bgzf = BgzfBlock { batch_serial: 0, bytes: b"\x1f\x8b".to_vec(), uncompressed_size: 4 }; + assert_eq!(RawBytesBlock::bytes(&bgzf), b"\x1f\x8b"); + } + + #[test] + fn profile_is_serial_writer_sink() { + let path = tempfile::NamedTempFile::new().unwrap().into_temp_path(); + let step = WriteRawFile::::new(&path, b"").unwrap(); + let profile = step.profile(); + assert_eq!(profile.name, "WriteRawFile"); + assert_eq!(profile.kind, StepKind::Serial); + assert!(profile.sticky); + assert_eq!(step.affinity(), Affinity::Writer); + } + + /// A sink dropped before the drained-finish path (an aborted stream) must + /// NOT append its trailer — stamping the BGZF EOF onto a truncated `.gz` + /// would hide the truncation. Mirrors + /// `WriteBgzfFile::drop_before_finish_does_not_append_eof_marker`. + #[test] + fn drop_before_finish_omits_trailer() { + let tmp = tempfile::NamedTempFile::new().unwrap(); + let path = tmp.path().to_path_buf(); + let trailer = b"\x1f\x8bTRAILER"; + let step = WriteRawFile::::new(&path, trailer).unwrap(); + // Write some payload bytes directly, then drop without draining (abort). + { + let mut guard = step.state.lock(); + guard.as_mut().unwrap().write_all(b"PAYLOAD").unwrap(); + } + drop(step); + + let bytes = std::fs::read(&path).unwrap(); + assert_eq!(bytes, b"PAYLOAD", "aborted stream must contain payload only, no trailer"); + assert!(!bytes.ends_with(trailer), "aborted stream must not end with the trailer"); + } + + #[test] + fn dash_path_selects_stdout_without_error() { + // `-` must construct a stdout-backed sink without touching the filesystem. + let step = WriteRawFile::::new("-", b"").unwrap(); + assert!(step.state.lock().is_some()); + } +} diff --git a/docs/src/guide/best-practices.md b/docs/src/guide/best-practices.md index af6fc6d65..528091f9b 100644 --- a/docs/src/guide/best-practices.md +++ b/docs/src/guide/best-practices.md @@ -525,6 +525,63 @@ fgumi dedup --input sorted.bam --output deduped.bam \ --- +## Alternative: UMI-Extracted FASTQ Output (for third-party tools) + +Some downstream tools expect the UMI in the **FASTQ read name** rather than in a BAM tag — for example Illumina DRAGEN, which reads the UMI from the read header and, for duplex UMIs, requires the two UMIs joined with a `+`. Because `fgumi extract` already parses the UMIs, trims them off the reads, and stores them in the standard `RX` tag, `fgumi fastq --annotate-read-names` can write them straight back into the read name. This replaces the slow `umi-tools` + `samtools fastq` combination with a single fast pipeline that yields both a UMI-tagged BAM and UMI-in-header FASTQs. + +```mermaid +graph TD; +A["fgumi extract"]-->B["extracted.bam (UMI in RX tag, reads trimmed)"]; +B-->C["fgumi fastq --annotate-read-names"]; +C-->D["R1.fastq.gz / R2.fastq.gz (UMI in read name)"]; +``` + +```bash +# Step 1: Extract UMIs from FASTQ (UMIs trimmed off, stored in the RX tag) +fgumi extract \ + --inputs r1.fq.gz r2.fq.gz \ + --read-structures 8M+T 8M+T \ + --sample "sample_name" \ + --library "library_name" \ + --output extracted.bam + +# Step 2: Write paired, gzipped FASTQs with the UMI in the read name +fgumi fastq \ + --input extracted.bam \ + --annotate-read-names \ + --threads 8 \ + --out1 out_R1.fastq.gz \ + --out2 out_R2.fastq.gz \ + --out0 /dev/null +``` + +Key points: +- Set `--read-structures` in Step 1 to match where your UMIs actually live (e.g. `8M+T +T` for a UMI on R1 only). See [Read Structures](read-structures.md). +- **DRAGEN-ready by default**: a duplex UMI stored as `RX:Z:AAAAAAAA-CCCCCCCC` becomes the read name `@readname:AAAAAAAA+CCCCCCCC` — identical to `samtools fastq -U`, and the `+`-joined format DRAGEN expects. +- In split (`--out1`/`--out2`) mode the `/1` `/2` suffixes are omitted so R1 and R2 read names match, as in `samtools fastq -1/-2`. Reads that are neither cleanly R1 nor R2 go to `--out0` (or stdout if `--out0` is omitted). +- Any `.gz`/`.bgz` output path is written as **BGZF** (gzip-compatible for downstream tools, yet block-indexable). Pass `--threads` to parallelize BGZF compression on compressed output. + +### Optional: Interleaved output and custom delimiters + +```bash +# Interleaved FASTQ to stdout (e.g. to pipe elsewhere), UMI in the read name. +# `-U` is a short alias for `--annotate-read-names`, matching `samtools fastq -U`. +fgumi fastq --input extracted.bam -U > umi_in_name.fastq + +# Custom delimiters for non-DRAGEN platforms (e.g. MGI, Element, ONT): +# underscore between the read name and UMI, and no delimiter between duplex UMIs. +fgumi fastq --input extracted.bam --annotate-read-names \ + --umi-name-delim _ --umi-sep '' \ + --out1 out_R1.fastq.gz --out2 out_R2.fastq.gz +``` + +Options: +- `--umi-tag`: tag(s) to read the UMI from, first present wins (default `RX,OX`) +- `--umi-name-delim`: delimiter between the read name and the UMI (default `:`) +- `--umi-sep`: delimiter joining duplex UMIs, replacing the stored `-` (default `+`) + +--- + ## Recommended Parameters by Application Pick a starting point by how you trade sensitivity against specificity for your assay, then tune from there. diff --git a/docs/src/index.md b/docs/src/index.md index 4fed41ad7..0b4a5c4f0 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -60,7 +60,7 @@ The commands below are grouped by pipeline stage. |-------|---------|-------------| | UMI extraction | `extract` | Extract UMIs from FASTQ and create an unmapped BAM | | UMI extraction | `correct` | Correct UMIs in a BAM to a fixed set of UMIs | -| Alignment support | `fastq` | Convert BAM to interleaved FASTQ for an aligner | +| Alignment support | `fastq` | Convert BAM to FASTQ (interleaved or paired), optionally embedding the UMI in the read name | | Alignment support | `zipper` | Merge an aligned BAM with its unmapped BAM, restoring tags | | Alignment support | `sort` | Sort by coordinate, queryname, or template-coordinate | | Alignment support | `merge` | Merge pre-sorted BAM files into one sorted BAM | diff --git a/src/lib/commands/fastq.rs b/src/lib/commands/fastq.rs index 51fe5c92f..698e57fd8 100644 --- a/src/lib/commands/fastq.rs +++ b/src/lib/commands/fastq.rs @@ -1,20 +1,23 @@ //! Convert BAM to FASTQ format. //! -//! This tool reads a BAM file and outputs interleaved FASTQ to stdout for piping to aligners. -//! Input should be queryname-sorted or template-coordinate sorted. - -use crate::commands::common::parse_bool; -use crate::logging::OperationTimer; +//! Reads a BAM and writes FASTQ: interleaved to stdout (for piping to an +//! aligner) or a file, paired split files (`--out1`/`--out2`/`--out0`), and +//! optionally BGZF-compressed (`.gz`/`.bgz`) with the UMI embedded in the read +//! name (`--annotate-read-names`). Both interleaved and paired-split output run +//! on the typed-step pipeline. Input should be queryname-sorted or +//! template-coordinate sorted. + +use crate::commands::common::{QueueMemoryOptions, SchedulerOptions, parse_bool}; +use crate::sam::SamTag; use crate::validation::validate_file_exists; use anyhow::Result; use clap::Parser; -use fgumi_bam_io::create_raw_bam_reader; use fgumi_raw_bam::{ RawRecord, extract_sequence_into, quality_scores_slice, read_name as raw_read_name, }; use log::info; -use std::io::{BufWriter, Write, stdout}; -use std::path::PathBuf; +use std::io::Write; +use std::path::{Path, PathBuf}; use crate::commands::command::Command; @@ -46,16 +49,139 @@ static COMPLEMENT: [u8; 256] = { table }; +/// Parse a comma-separated list of two-character SAM tags (e.g. `"RX,OX"`). +/// +/// Whitespace around each tag is trimmed. Returns an error if any token is not +/// a valid SAM aux-tag identifier, which also rejects empty tokens (`""`, +/// `"RX,"`). +fn parse_umi_tags(spec: &str) -> Result> { + spec.split(',').map(|token| token.trim().parse::()).collect() +} + +/// Configuration for embedding a UMI into the FASTQ read header, mirroring +/// `samtools fastq -U`. +/// +/// The UMI is read from the first present tag in `tags` (default `RX` then +/// `OX`) and appended to the read name as ``, with each +/// `-` in the stored UMI (fgumi's duplex join) replaced by `umi_sep` (default +/// `+`, the delimiter Illumina DRAGEN expects between duplex UMIs). +#[derive(Clone)] +pub struct UmiHeader { + /// Tags to search, in priority order; the first present wins. + tags: Vec, + /// Bytes inserted between the read name and the UMI. + name_delim: Vec, + /// Bytes that replace each `-` (duplex join) in the stored UMI. + umi_sep: Vec, +} + +impl UmiHeader { + /// Build from raw CLI strings, validating the tag list. + pub(crate) fn from_args(umi_tag: &str, name_delim: &str, umi_sep: &str) -> Result { + Ok(Self { + tags: parse_umi_tags(umi_tag)?, + name_delim: name_delim.as_bytes().to_vec(), + umi_sep: umi_sep.as_bytes().to_vec(), + }) + } + + /// Look up the UMI value from the record's tags (first present, in priority order). + fn lookup<'a>(&self, record: &'a RawRecord) -> Option<&'a [u8]> { + let tags = record.tags(); + self.tags.iter().find_map(|tag| tags.find_string(*tag)) + } + + /// Write `` to `writer`, replacing each `-` in `umi` with + /// `umi_sep`. Maximal runs between separators are written in one call. + fn write_annotation(&self, writer: &mut W, umi: &[u8]) -> std::io::Result<()> { + writer.write_all(&self.name_delim)?; + let mut start = 0; + for (i, &b) in umi.iter().enumerate() { + if b == b'-' { + writer.write_all(&umi[start..i])?; + writer.write_all(&self.umi_sep)?; + start = i + 1; + } + } + writer.write_all(&umi[start..])?; + Ok(()) + } +} + +/// Which FASTQ stream a record belongs to, based on its segment flags. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub(crate) enum Segment { + /// First segment of a pair (`FIRST_SEGMENT` set, `LAST_SEGMENT` clear) → R1. + Read1, + /// Last segment of a pair (`LAST_SEGMENT` set, `FIRST_SEGMENT` clear) → R2. + Read2, + /// Neither or both segment bits set (single-end or ambiguous) → "other". + Other, +} + +/// Classify a record into R1 / R2 / other from its flags, matching how +/// `samtools fastq` routes reads to `-1` / `-2` / `-0`. +pub(crate) fn classify_segment(flags: u16) -> Segment { + use fgumi_raw_bam::flags as flag_bits; + let is_first = (flags & flag_bits::FIRST_SEGMENT) != 0; + let is_last = (flags & flag_bits::LAST_SEGMENT) != 0; + match (is_first, is_last) { + (true, false) => Segment::Read1, + (false, true) => Segment::Read2, + _ => Segment::Other, + } +} + +/// Returns `true` if `path` should be written compressed (ends in `.gz`/`.bgz`). +pub(crate) fn path_is_gzip(path: &Path) -> bool { + matches!(path.extension().and_then(|e| e.to_str()), Some("gz" | "bgz" | "bgzf")) +} + +/// Per-stage options for [`crate::pipeline::chains::Stage::Fastq`] — the +/// BAM→FASTQ encode. Carries the flag filters, read-name suffix behavior, and +/// the optional UMI-in-read-name config the encode step needs. The output +/// destination(s) and compression are carried by the chain's `SinkSpec`, not +/// here. +#[derive(Clone)] +pub struct FastqOptions { + /// Exclude reads with any of these flag bits set. + pub exclude_flags: u16, + /// Only include reads with all of these flag bits set. + pub require_flags: u16, + /// When `true`, omit the `/1` `/2` read-name suffix. + pub no_suffix: bool, + /// When `Some`, append the UMI (from the configured tag) to the read name. + pub umi_header: Option, +} + +impl FastqOptions { + /// Returns `true` if a record with `flags` passes the include/exclude filters. + #[must_use] + pub(crate) fn passes_filters(&self, flags: u16) -> bool { + (flags & self.exclude_flags) == 0 && (flags & self.require_flags) == self.require_flags + } +} + +/// Default BGZF/DEFLATE compression level for `.gz` FASTQ output (matches the +/// samtools default). +const FASTQ_GZIP_LEVEL: u32 = 6; + +/// Default for the deprecated, now-ignored `--bwa-chunk-size` flag. +const DEFAULT_BWA_CHUNK_SIZE: u64 = 150_000_000; + /// Convert BAM to FASTQ format. #[derive(Debug, Parser)] #[command( name = "fastq", about = "\x1b[38;5;72m[ALIGNMENT]\x1b[0m \x1b[36mConvert BAM to FASTQ format\x1b[0m", long_about = r#" -Convert a BAM file to interleaved FASTQ format. +Convert a BAM file to FASTQ. -Reads BAM records and outputs FASTQ to stdout for piping to aligners. -Input should be queryname-sorted or template-coordinate sorted. +By default, reads BAM records and writes interleaved FASTQ to stdout for piping +to aligners. Can also write paired split files (--out1/--out2, plus --out0), +gzip (BGZF) output for any .gz/.bgz path, and embed the UMI in the read name +(--annotate-read-names). Input should be queryname-sorted or template-coordinate +sorted. EXAMPLES: @@ -67,6 +193,13 @@ EXAMPLES: # Exclude secondary and supplementary alignments (default) fgumi fastq -i aligned.bam -F 0x900 | bwa mem ... + + # Embed the UMI (RX tag) in the read name for DRAGEN, like `samtools fastq -U`. + # A duplex UMI stored as `AAA-CCC` is emitted as `readname:AAA+CCC`. + fgumi fastq -i extracted.bam --annotate-read-names -n > umi_in_name.fq + + # Custom delimiters for other platforms (underscore separator, no duplex join) + fgumi fastq -i extracted.bam --annotate-read-names --umi-name-delim _ --umi-sep '' ... "# )] pub struct Fastq { @@ -91,16 +224,73 @@ pub struct Fastq { #[arg(short = 'f', long = "require-flags", default_value_t = 0, value_parser = parse_flags)] pub require_flags: u16, - /// Number of threads for BAM decompression. + /// Number of threads for BAM decompression and, for `.gz`/`.bgz` outputs, + /// BGZF compression (per output). Use several threads for compressed output. #[arg(short = '@', short_alias = 't', long = "threads", default_value = "1")] pub threads: usize, - /// BWA -K parameter value (bases per batch). Sizes output buffer to match bwa's - /// batch size for optimal pipe throughput. Default matches common bwa mem usage. - #[arg(short = 'K', long = "bwa-chunk-size", default_value = "150000000")] + /// Deprecated and ignored. Previously sized the stdout flush batch to bwa's + /// `-K`; the typed-step pipeline now manages output batching via its queue + /// budget (`--max-memory`), so this flag has no effect. Kept for backward + /// compatibility; overriding it logs a warning. + #[arg(short = 'K', long = "bwa-chunk-size", default_value_t = DEFAULT_BWA_CHUNK_SIZE)] pub bwa_chunk_size: u64, + + /// Append the UMI (from a BAM tag) to the read name, matching + /// `fgumi extract --annotate-read-names` and `samtools fastq -U`. The UMI is + /// inserted before any /1 or /2 suffix. + #[arg(short = 'a', visible_short_alias = 'U', long = "annotate-read-names", default_value = "false", num_args = 0..=1, default_missing_value = "true", action = clap::ArgAction::Set, value_parser = parse_bool)] + pub annotate_read_names: bool, + + /// Comma-separated tag(s) to source the UMI from; the first present wins. + #[arg(long = "umi-tag", default_value = "RX,OX")] + pub umi_tag: String, + + /// Delimiter inserted between the read name and the UMI. + #[arg(long = "umi-name-delim", default_value = ":")] + pub umi_name_delim: String, + + /// Delimiter joining multiple (duplex) UMIs, replacing the stored '-' + /// [default '+' matches Illumina DRAGEN]. + #[arg(long = "umi-sep", default_value = "+")] + pub umi_sep: String, + + /// Write read 1 (R1) to this file instead of the interleaved stream; requires + /// --out2. A `.gz`/`.bgz` path is written as BGZF. + #[arg(short = '1', long = "out1", requires = "out2", conflicts_with = "output")] + pub out1: Option, + + /// Write read 2 (R2) to this file; requires --out1. + /// A `.gz`/`.bgz` path is written as BGZF. + #[arg(short = '2', long = "out2", requires = "out1", conflicts_with = "output")] + pub out2: Option, + + /// Write reads that are neither cleanly R1 nor R2 (single-end / ambiguous) + /// here; requires --out1. If omitted, such reads go to stdout, matching + /// `samtools fastq` without -0. A `.gz`/`.bgz` path is written as BGZF. + #[arg(short = '0', long = "out0", requires = "out1")] + pub out0: Option, + + /// Pipeline memory limits (`--max-memory`, …). The conversion runs on the + /// typed-step pipeline; these cap the in-flight queue memory. Defaults to a + /// lean fastq budget (see `FASTQ_DEFAULT_QUEUE_MEMORY_MB`); pass + /// `--max-memory` to override. + #[command(flatten)] + pub queue_memory: QueueMemoryOptions, + + /// Pipeline scheduler diagnostics (`--pipeline-stats`, `--pipeline-trace`). + #[command(flatten)] + pub scheduler: SchedulerOptions, } +/// Default total in-flight queue-memory budget for `fgumi fastq` when the user +/// passes no `--max-memory`/`--queue-memory` flag. FASTQ blocks are small and +/// the chain is a short linear source→encode→write pipeline, so the +/// sort/consensus-sized default (768 MiB/thread) is wildly oversized; a modest +/// fixed total keeps RSS flat as `--threads` scales while staying well above +/// the pipeline's steady-state working set. +const FASTQ_DEFAULT_QUEUE_MEMORY_MB: u64 = 256; + /// Parse flag values supporting both decimal and hex (0x) notation. fn parse_flags(s: &str) -> Result { if s.starts_with("0x") || s.starts_with("0X") { @@ -111,118 +301,199 @@ fn parse_flags(s: &str) -> Result { } impl Fastq { - /// Run the BAM-to-FASTQ conversion loop against the given writer. - /// - /// Extracted so `execute` can dispatch between stdout (the default piping - /// path) and a file (`--output`) without duplicating the conversion loop. - fn run_with_writer(&self, mut writer: W) -> Result<()> { - let timer = OperationTimer::new("Converting BAM to FASTQ"); + /// Build the optional UMI-in-read-name config, validating the tag list up front. + fn build_umi_header(&self) -> Result> { + if self.annotate_read_names { + info!("UMI in read name: from tag(s) {}", self.umi_tag); + Ok(Some(UmiHeader::from_args(&self.umi_tag, &self.umi_name_delim, &self.umi_sep)?)) + } else { + Ok(None) + } + } + /// Log the shared conversion configuration. + fn log_config(&self) { info!("Input: {}", self.input.display()); info!("Threads: {}", self.threads); info!("Exclude flags: 0x{:X}", self.exclude_flags); info!("Require flags: 0x{:X}", self.require_flags); info!("Read name suffix: {}", if self.no_suffix { "disabled" } else { "enabled" }); - info!("BWA chunk size: {} bases", self.bwa_chunk_size); - - let (mut reader, _header) = create_raw_bam_reader(&self.input, self.threads)?; - - let mut total_records: u64 = 0; - let mut written_records: u64 = 0; - - // Batch tracking (mirrors bwa's logic) - let mut bases_this_batch: u64 = 0; - let mut records_this_batch: usize = 0; - - // Reusable buffers to avoid per-record allocations - let mut seq_buf: Vec = Vec::with_capacity(512); - let mut qual_buf: Vec = Vec::with_capacity(512); - let mut record = RawRecord::new(); + } - loop { - let n = reader.read_record(&mut record)?; - if n == 0 { - break; // EOF - } - total_records += 1; + /// Build the per-stage [`FastqOptions`] — the single source of truth for the + /// flag filters and UMI-header config, shared by the interleaved and + /// paired-split chain specs. + fn to_fastq_options(&self) -> Result { + Ok(FastqOptions { + exclude_flags: self.exclude_flags, + require_flags: self.require_flags, + no_suffix: self.no_suffix, + umi_header: self.build_umi_header()?, + }) + } - let flags = record.flags(); + /// Assemble the typed-step [`ChainSpec`] for a BAM → FASTQ conversion: + /// `SourceSpec::Bam` → `Stage::Fastq` → the given `sink`. + fn build_chain_spec( + &self, + sink: crate::pipeline::chains::SinkSpec, + command_line: &str, + ) -> Result { + use crate::commands::common::ThreadingOptions; + use crate::pipeline::chains::{ChainSpec, SourceSpec, Stage, StageOptionsBag}; + use fgumi_cli_common::CompressionOptions; + + Ok(ChainSpec { + stages: vec![Stage::Fastq], + source: SourceSpec::Bam(self.input.clone()), + sink, + stage_opts: StageOptionsBag { + fastq: Some(self.to_fastq_options()?), + ..Default::default() + }, + threading: ThreadingOptions::new(self.threads), + compression: CompressionOptions { compression_level: FASTQ_GZIP_LEVEL }, + scheduler: self.scheduler.clone(), + queue_memory: self.resolve_queue_memory(), + async_reader: false, + command_line: command_line.to_string(), + }) + } - // Filter by flags - if (flags & self.exclude_flags) != 0 { - continue; - } - if (flags & self.require_flags) != self.require_flags { - continue; - } + /// Assemble the [`ChainSpec`] for paired-split BAM → FASTQ: + /// `SourceSpec::Bam` → `Stage::Fastq` → `SinkSpec::FastqPaired`. The only + /// difference from [`Self::build_chain_spec`] is the sink; everything else + /// (stage options, threading, compression, scheduler, queue memory) is + /// shared. + fn build_paired_chain_spec( + &self, + command_line: &str, + ) -> Result { + // clap `requires` guarantees out2 is present whenever out1 is (and + // execute() only calls this when out1.is_some()). + let out1 = self.out1.clone().expect("out1 present in paired mode"); + let out2 = self.out2.clone().expect("out2 present in paired mode"); + let sink = + crate::pipeline::chains::SinkSpec::FastqPaired { out1, out2, out0: self.out0.clone() }; + self.build_chain_spec(sink, command_line) + } - // Get sequence length for batch tracking - let seq_len = record.l_seq() as usize; - - // Write FASTQ record - write_fastq_record( - &mut writer, - &record, - flags, - self.no_suffix, - &mut seq_buf, - &mut qual_buf, - )?; - written_records += 1; - - // Track batch progress - bases_this_batch += seq_len as u64; - records_this_batch += 1; - - // Flush at batch boundary (like bwa: bases >= K AND record count is even) - if bases_this_batch >= self.bwa_chunk_size && records_this_batch.is_multiple_of(2) { - writer.flush()?; - bases_this_batch = 0; - records_this_batch = 0; - } + /// The queue-memory options for the chain: the user's `--max-memory` / + /// `--queue-memory` flags if given, otherwise the lean fastq default + /// ([`FASTQ_DEFAULT_QUEUE_MEMORY_MB`]) so a large `--threads` doesn't + /// inflate RSS to the sort/consensus-sized 768 MiB/thread default. + /// + /// The default is applied through the modern `max_memory` field (a fixed + /// total, not per-thread) so it does not trip the deprecation warning the + /// legacy `queue_memory_limit_mb` field would. + fn resolve_queue_memory(&self) -> QueueMemoryOptions { + use crate::commands::common::MemoryLimit; + + let mut qm = self.queue_memory.clone(); + let user_set_limit = qm.max_memory.is_some() + || qm.queue_memory.is_some() + || qm.queue_memory_limit_mb.is_some(); + if !user_set_limit { + let bytes = FASTQ_DEFAULT_QUEUE_MEMORY_MB as usize * 1024 * 1024; + qm.max_memory = Some(MemoryLimit::Fixed(bytes)); + qm.memory_per_thread = false; } - - // Final flush - writer.flush()?; - - info!("Read {total_records} records, wrote {written_records} FASTQ records"); - timer.log_completion(written_records); - Ok(()) + qm } } impl Command for Fastq { - fn execute(&self, _command_line: &str) -> Result<()> { + fn execute(&self, command_line: &str) -> Result<()> { validate_file_exists(&self.input, "Input BAM")?; - // Refuse to clobber the input BAM. `File::create(path)` truncates - // before the input is opened in `run_with_writer`, so an `--output` - // pointing at the same file silently destroys the input data. - // Lexical equality catches the obvious case; canonicalising both - // sides also catches `./in.bam` vs `in.bam`, symlinks, and absolute - // vs relative paths pointing at the same file. - if let Some(output) = &self.output { - let same_path = output == &self.input - || (output.exists() - && std::fs::canonicalize(output)? == std::fs::canonicalize(&self.input)?); - if same_path { + if self.bwa_chunk_size != DEFAULT_BWA_CHUNK_SIZE { + info!( + "--bwa-chunk-size ({}) is deprecated and ignored; output batching is managed \ + by the pipeline (tune with --max-memory)", + self.bwa_chunk_size + ); + } + + // Refuse to clobber the input BAM or write two streams to the same file. + // `File::create(path)` truncates before the input is opened, so an output + // path pointing at the input silently destroys it, and two outputs sharing + // a path silently corrupt each other (two writer threads appending to one + // truncated file). Compare by canonical path where the file exists (which + // also catches `./x` vs `x`, symlinks, and absolute vs relative). For + // outputs that don't exist yet, canonicalize the (existing) parent + // directory and re-join the file name so aliases that only differ through + // a symlinked intermediate directory (or a leading `./`) still collide in + // the `seen` dedup — otherwise two sink threads could truncate the same + // not-yet-existing file. + let output_paths: Vec<&PathBuf> = + [&self.output, &self.out1, &self.out2, &self.out0].into_iter().flatten().collect(); + // Nested `if let` (not a let-chain) to stay compatible with the 1.87.0 MSRV. + let canonical = |p: &Path| -> Result { + if p.exists() { + return Ok(std::fs::canonicalize(p)?); + } + if let Some(parent) = p.parent().filter(|d| !d.as_os_str().is_empty()) { + if let Ok(canon_parent) = std::fs::canonicalize(parent) { + if let Some(file_name) = p.file_name() { + return Ok(canon_parent.join(file_name)); + } + } + } + Ok(p.to_path_buf()) + }; + let input_key = canonical(&self.input)?; + let mut seen: std::collections::HashSet = std::collections::HashSet::new(); + for output in &output_paths { + // The literal `/dev/null` is exempt — routing several streams to it + // is intentional, as with `samtools fastq` (e.g. `--out0 /dev/null`). + // Only the exact path is exempt, not symlinks to it or `/dev/stdout`. + if output.as_os_str() == "/dev/null" { + continue; + } + let key = canonical(output)?; + if key == input_key { anyhow::bail!( - "--output {} must differ from --input {} (would truncate the input BAM)", + "output path {} must differ from --input {} (would truncate the input BAM)", output.display(), self.input.display() ); } + if !seen.insert(key) { + anyhow::bail!( + "output path {} is used more than once; each FASTQ output must be a distinct file", + output.display() + ); + } } - // Use 64MB buffer for efficient pipe throughput. - const BUF_CAPACITY: usize = 64 * 1024 * 1024; - match &self.output { - Some(path) => { - let file = std::fs::File::create(path)?; - self.run_with_writer(BufWriter::with_capacity(BUF_CAPACITY, file)) + // Paired split output (`-1`/`-2`, optionally `-0`) runs through the + // typed-step chain: BAM source → Stage::Fastq 3-way encode → three + // per-file raw writers (BGZF for `.gz`), all sharing one work-stealing + // pool. clap guarantees `--out1`/`--out2` come together and conflict + // with `--output`. + if let Some(out1) = &self.out1 { + self.log_config(); + // clap `requires` guarantees out2 is present whenever out1 is. + let out2 = self.out2.as_ref().expect("out2 present in paired mode"); + info!("R1 -> {}", out1.display()); + info!("R2 -> {}", out2.display()); + if let Some(out0) = &self.out0 { + info!("other -> {}", out0.display()); } - None => self.run_with_writer(BufWriter::with_capacity(BUF_CAPACITY, stdout().lock())), + return crate::pipeline::chains::build_for( + self.build_paired_chain_spec(command_line)?, + )? + .run(); } + + // Interleaved output runs through the typed-step chain: BAM source → + // Stage::Fastq encode (pool-spread) → raw writer (BGZF for `.gz`). `-` + // (or omitted `--output`) writes to stdout. + self.log_config(); + let sink = crate::pipeline::chains::SinkSpec::Fastq( + self.output.clone().unwrap_or_else(|| PathBuf::from("-")), + ); + crate::pipeline::chains::build_for(self.build_chain_spec(sink, command_line)?)?.run() } } @@ -231,12 +502,16 @@ impl Command for Fastq { /// `pub(crate)` so the runall AAM chain (`crate::align_and_merge`) /// can reuse the same BAM→FASTQ logic the standalone `fgumi fastq` /// command uses, instead of duplicating the encoding. +/// +/// When `umi_header` is `Some`, the UMI from the configured tag is appended to +/// the read name (before any `/1`/`/2` suffix), mirroring `samtools fastq -U`. #[inline] pub(crate) fn write_fastq_record( writer: &mut W, record: &RawRecord, flags: u16, no_suffix: bool, + umi_header: Option<&UmiHeader>, seq_buf: &mut Vec, qual_buf: &mut Vec, ) -> Result<()> { @@ -278,6 +553,13 @@ pub(crate) fn write_fastq_record( // Write all parts writer.write_all(b"@")?; writer.write_all(name)?; + // Append the UMI to the read name (before the /1 /2 suffix), if requested. + // Nested `if let` (not a let-chain) to stay compatible with the 1.87.0 MSRV. + if let Some(uh) = umi_header { + if let Some(umi) = uh.lookup(record) { + uh.write_annotation(writer, umi)?; + } + } writer.write_all(suffix)?; writer.write_all(b"\n")?; writer.write_all(seq_buf)?; @@ -441,4 +723,73 @@ mod tests { .expect("write_quality_bytes should succeed"); assert_eq!(output, vec![126, 126, 126]); } + + // ── UMI-in-header (mirrors `samtools fastq -U`) ───────────────────────── + + use fgumi_raw_bam::flags as fb; + use rstest::rstest; + use std::path::Path; + + #[rstest] + #[case::default_list("RX,OX", vec![SamTag::RX, SamTag::OX])] + #[case::single("RX", vec![SamTag::RX])] + #[case::trims_whitespace(" RX , OX ", vec![SamTag::RX, SamTag::OX])] + fn parse_umi_tags_accepts_valid(#[case] spec: &str, #[case] expected: Vec) { + assert_eq!(parse_umi_tags(spec).expect("valid"), expected); + } + + #[rstest] + #[case::one_char("R")] + #[case::empty("")] + #[case::trailing_comma("RX,")] + fn parse_umi_tags_rejects_malformed(#[case] spec: &str) { + assert!(parse_umi_tags(spec).is_err(), "spec {spec:?} must be rejected"); + } + + #[rstest] + // name_delim, umi_sep, stored UMI -> annotation + #[case::single_umi(":", "+", b"AAAA".as_slice(), b":AAAA".as_slice())] + #[case::duplex_hyphen_to_plus(":", "+", b"AAAA-GGGG", b":AAAA+GGGG")] + #[case::multi_segment(":", "+", b"A-B-C", b":A+B+C")] + // Underscore name separator (umi-tools style) with an empty UMI separator + // (platforms that want the duplex UMIs concatenated with no delimiter). + #[case::custom_delims_empty_sep("_", "", b"AAAA-GGGG", b"_AAAAGGGG")] + fn umi_annotation_writes_expected( + #[case] name_delim: &str, + #[case] umi_sep: &str, + #[case] umi: &[u8], + #[case] expected: &[u8], + ) { + let uh = UmiHeader::from_args("RX", name_delim, umi_sep).expect("valid config"); + let mut out = Vec::new(); + uh.write_annotation(&mut out, umi).expect("write"); + assert_eq!(out, expected); + } + + #[test] + fn umi_header_rejects_bad_tag() { + assert!(UmiHeader::from_args("nope", ":", "+").is_err()); + } + + // ── Paired split output (mirrors `samtools fastq -1/-2/-0`) ───────────── + + #[rstest] + #[case::read1(fb::PAIRED | fb::FIRST_SEGMENT, Segment::Read1)] + #[case::read2(fb::PAIRED | fb::LAST_SEGMENT, Segment::Read2)] + // Both segment bits set, or neither, is "other" (single-end / ambiguous). + #[case::both_bits(fb::PAIRED | fb::FIRST_SEGMENT | fb::LAST_SEGMENT, Segment::Other)] + #[case::no_bits(0, Segment::Other)] + fn classify_segment_by_flags(#[case] flags: u16, #[case] expected: Segment) { + assert_eq!(classify_segment(flags), expected); + } + + #[rstest] + #[case::fq_gz("out_R1.fq.gz", true)] + #[case::gz("out.gz", true)] + #[case::bgz("out.bgz", true)] + #[case::plain_fq("out_R1.fq", false)] + #[case::fastq("out.fastq", false)] + fn gz_extension_is_detected(#[case] path: &str, #[case] expected: bool) { + assert_eq!(path_is_gzip(Path::new(path)), expected); + } } diff --git a/src/lib/commands/runall.rs b/src/lib/commands/runall.rs index e2103f61f..7d3b7927a 100644 --- a/src/lib/commands/runall.rs +++ b/src/lib/commands/runall.rs @@ -1694,7 +1694,7 @@ impl RunAll { // Stages that runall does not use. If any of these appear in // the stages list it is a programming error (derive_stages // never emits them for a runall chain). - Stage::Clip | Stage::Dedup | Stage::Downsample => { + Stage::Clip | Stage::Dedup | Stage::Downsample | Stage::Fastq => { bail!( "internal error: build_stage_options_bag encountered unexpected \ stage {stage:?} in a runall chain; this is a bug in derive_stages" diff --git a/src/lib/pipeline/chains/builder.rs b/src/lib/pipeline/chains/builder.rs index 60ef5db48..342a75277 100644 --- a/src/lib/pipeline/chains/builder.rs +++ b/src/lib/pipeline/chains/builder.rs @@ -1230,12 +1230,76 @@ impl<'a> ChainBuilder<'a> { Stage::Zipper => self.add_zipper(position), Stage::Align => self.add_align(position), Stage::Extract => self.add_extract(position), + Stage::Fastq => self.add_fastq(position), Stage::Downsample => { Err(anyhow!("Stage::Downsample is out of scope — not on the typed-step pipeline")) } } } + /// Add the terminal BAM → FASTQ encode stage. + /// + /// Consumes the `DecodedRecordBatch` tail from the BAM source preamble. + /// For a single-output sink, appends the `Parallel` FASTQ-encode step + /// (pool-spread), leaving a `DecompressedBlock` byte tail for `add_sink` + /// to route to a raw writer (optionally through `BgzfCompress` for + /// `.gz`/`.bgz` output). For a `SinkSpec::FastqPaired` sink, appends the + /// 3-output `Process3` encode step and wires branch b (R2) to `out2` and + /// branch c (other) to `out0`/stdout in-method (mirroring how + /// `add_filter`/`wire_consensus_rejects_branch` wire secondary branches), + /// leaving branch a (R1) as `current_tail` for `add_sink` to route to + /// `out1`. + fn add_fastq(&mut self, position: StagePosition) -> Result<()> { + use std::sync::atomic::AtomicU64; + + use crate::logging::OperationTimer; + use crate::pipeline::chains::commands::fastq::{ + FastqFinalizeHook, build_fastq_encode_step, build_fastq_paired_encode_step, + }; + use crate::pipeline::core::topology::BranchIdx; + + debug_assert_eq!(position, StagePosition::Terminal, "Stage::Fastq is always terminal"); + + let opts = self + .spec + .stage_opts + .fastq + .as_ref() + .ok_or_else(|| anyhow!("Stage::Fastq options missing from StageOptionsBag"))?; + + let tail = self.current_tail.expect("add_fastq called before add_source"); + let timer = OperationTimer::new("Converting BAM to FASTQ"); + let records_written = Arc::new(AtomicU64::new(0)); + + if let SinkSpec::FastqPaired { out2, out0, .. } = &self.spec.sink { + let out2 = out2.clone(); + let out0 = out0.clone(); + let step = build_fastq_paired_encode_step( + Arc::new(opts.clone()), + self.tuning.per_step_byte_limit, + Arc::clone(&records_written), + ); + let process_tail = self.pipeline.append_step(step, tail); + // Branch 1 = R2 → out2; branch 2 = "other" → out0 (or stdout via `-`). + self.wire_fastq_output((process_tail.0, BranchIdx(1)), &out2)?; + let out0_path = out0.unwrap_or_else(|| std::path::PathBuf::from("-")); + self.wire_fastq_output((process_tail.0, BranchIdx(2)), &out0_path)?; + // Branch 0 = R1 → current_tail for add_sink → out1. + self.current_tail = Some(process_tail); + } else { + let step = build_fastq_encode_step( + Arc::new(opts.clone()), + self.tuning.per_step_byte_limit, + Arc::clone(&records_written), + ); + self.current_tail = Some(self.pipeline.append_step(step, tail)); + } + self.chain_tail_kind = ChainTailKind::SerializedBytes; + + self.finalize.push(Box::new(FastqFinalizeHook { records_written, timer })); + Ok(()) + } + /// Add the output sink step(s) to the pipeline. /// /// Appends `BgzfCompress` + `WriteBgzfFile`. Reads `spec.sink` for the @@ -1252,6 +1316,20 @@ impl<'a> ChainBuilder<'a> { use crate::pipeline::steps::bgzf::compress::BgzfCompress; use crate::pipeline::steps::sink::write_bgzf::WriteBgzfFile; + // FASTQ sinks write raw bytes (no BAM header, no BGZF EOF) via a + // dedicated writer, optionally through `BgzfCompress` for `.gz`/`.bgz`. + if let SinkSpec::Fastq(path) = &self.spec.sink { + let path = path.clone(); + return self.add_fastq_sink(&path); + } + if let SinkSpec::FastqPaired { out1, .. } = &self.spec.sink { + let out1 = out1.clone(); + // Branches b (R2) and c (other) were wired in add_fastq; branch a + // (R1) is current_tail here. + let tail = self.current_tail.expect("add_sink called before add_source"); + return self.wire_fastq_output(tail, &out1); + } + let tail = self.current_tail.expect("add_sink called before add_source"); let output_path = self.spec.sink.path(); @@ -1288,6 +1366,43 @@ impl<'a> ChainBuilder<'a> { Ok(()) } + /// Route the FASTQ byte tail to a single (interleaved) raw writer. + fn add_fastq_sink(&mut self, path: &std::path::Path) -> Result<()> { + let tail = self.current_tail.expect("add_sink called before add_source"); + self.wire_fastq_output(tail, path)?; + // add_fastq_sink terminates the chain; no downstream tail. + Ok(()) + } + + /// Wire one FASTQ byte branch (`tail`) to its own writer: a `.gz`/`.bgz` + /// `path` is compressed via the shared-pool `BgzfCompress` step first + /// (appending the BGZF EOF marker so the stream is a complete BGZF file); + /// otherwise the bytes are written verbatim (file, or stdout via `-`). + fn wire_fastq_output( + &mut self, + tail: (StepIdx, BranchIdx), + path: &std::path::Path, + ) -> Result<()> { + use crate::commands::fastq::path_is_gzip; + use crate::pipeline::steps::bgzf::compress::BgzfCompress; + use crate::pipeline::steps::sink::write_raw::WriteRawFile; + use crate::pipeline::steps::types::{BgzfBlock, DecompressedBlock}; + + if path_is_gzip(path) { + let compress = + BgzfCompress::new(self.tuning.compression_level, self.tuning.per_step_byte_limit); + let compress_tail = self.pipeline.append_step(compress, tail); + let writer = WriteRawFile::::new(path, &fgumi_bgzf::BGZF_EOF) + .map_err(|e| anyhow!("WriteRawFile: {e}"))?; + self.pipeline.append_step(writer, compress_tail); + } else { + let writer = WriteRawFile::::new(path, b"") + .map_err(|e| anyhow!("WriteRawFile: {e}"))?; + self.pipeline.append_step(writer, tail); + } + Ok(()) + } + /// Replace `self.header` with `header`, invalidating any pending /// [`HeaderHandle`] stashed by [`Self::add_align`]. /// @@ -1440,9 +1555,7 @@ impl<'a> ChainBuilder<'a> { let tail = self.current_tail.expect("add_extract called before add_source"); - let output_path = match &self.spec.sink { - SinkSpec::Bam(p) | SinkSpec::BamWithIndex(p) => p.clone(), - }; + let output_path = self.spec.sink.path().clone(); let timer = OperationTimer::new("Extracting UMIs"); @@ -1569,9 +1682,7 @@ impl<'a> ChainBuilder<'a> { // Resolve source path for log messages only. let input_path = self.resolve_log_input_path(); - let output_path = match &self.spec.sink { - SinkSpec::Bam(p) | SinkSpec::BamWithIndex(p) => p.clone(), - }; + let output_path = self.spec.sink.path().clone(); let timer = OperationTimer::new("Correcting UMIs (new pipeline)"); @@ -2485,9 +2596,7 @@ impl<'a> ChainBuilder<'a> { // summary still lands before any "Indexing BAM:" line, mirroring // the legacy SortFinalizeHook ordering. if let Some(slot) = sort_stats_slot { - let output_path = match &self.spec.sink { - SinkSpec::Bam(p) | SinkSpec::BamWithIndex(p) => p.clone(), - }; + let output_path = self.spec.sink.path().clone(); self.finalize_on_success.push(Box::new( crate::pipeline::chains::commands::sort::SortSummaryFinalizeHook { stats_slot: slot, @@ -2624,9 +2733,7 @@ impl<'a> ChainBuilder<'a> { // Resolve source path for log messages only. let input_path = self.resolve_log_input_path(); - let output_path = match &self.spec.sink { - SinkSpec::Bam(p) | SinkSpec::BamWithIndex(p) => p.clone(), - }; + let output_path = self.spec.sink.path().clone(); let timer = OperationTimer::new("Grouping reads by UMI"); @@ -2891,9 +2998,7 @@ impl<'a> ChainBuilder<'a> { // Resolve source path for log messages only. let input_path = self.resolve_log_input_path(); - let output_path = match &self.spec.sink { - SinkSpec::Bam(p) | SinkSpec::BamWithIndex(p) => p.clone(), - }; + let output_path = self.spec.sink.path().clone(); let timer = OperationTimer::new("Calling simplex consensus"); @@ -3142,9 +3247,7 @@ impl<'a> ChainBuilder<'a> { // Resolve source path for log messages only. let input_path = self.resolve_log_input_path(); - let output_path = match &self.spec.sink { - SinkSpec::Bam(p) | SinkSpec::BamWithIndex(p) => p.clone(), - }; + let output_path = self.spec.sink.path().clone(); let timer = OperationTimer::new("Calling duplex consensus"); @@ -3385,9 +3488,7 @@ impl<'a> ChainBuilder<'a> { // Resolve source path for log messages only. let input_path = self.resolve_log_input_path(); - let output_path = match &self.spec.sink { - SinkSpec::Bam(p) | SinkSpec::BamWithIndex(p) => p.clone(), - }; + let output_path = self.spec.sink.path().clone(); let timer = OperationTimer::new("Calling CODEC consensus"); @@ -3608,9 +3709,7 @@ impl<'a> ChainBuilder<'a> { // Resolve source path for log messages only. let input_path = self.resolve_log_input_path(); - let output_path = match &self.spec.sink { - SinkSpec::Bam(p) | SinkSpec::BamWithIndex(p) => p.clone(), - }; + let output_path = self.spec.sink.path().clone(); // Reference is always required for clip. crate::validation::validate_file_exists(&clip.reference, "Reference FASTA")?; @@ -3773,9 +3872,7 @@ impl<'a> ChainBuilder<'a> { let timer = OperationTimer::new("Filtering consensus reads"); // Resolve output path for log. - let output_path = match &self.spec.sink { - SinkSpec::Bam(p) | SinkSpec::BamWithIndex(p) => p.clone(), - }; + let output_path = self.spec.sink.path().clone(); info!("Starting Filter"); info!("Input: {}", input_path.display()); @@ -3976,9 +4073,7 @@ impl<'a> ChainBuilder<'a> { let timer = OperationTimer::new("Marking duplicates"); - let output_path = match &self.spec.sink { - SinkSpec::Bam(p) | SinkSpec::BamWithIndex(p) => p.clone(), - }; + let output_path = self.spec.sink.path().clone(); info!("Starting dedup"); info!("Output: {}", output_path.display()); diff --git a/src/lib/pipeline/chains/commands/fastq.rs b/src/lib/pipeline/chains/commands/fastq.rs new file mode 100644 index 000000000..425b03057 --- /dev/null +++ b/src/lib/pipeline/chains/commands/fastq.rs @@ -0,0 +1,346 @@ +//! Chain-builder support for [`crate::pipeline::chains::Stage::Fastq`] — the +//! terminal BAM → FASTQ encode. +//! +//! The encode step consumes the [`DecodedRecordBatch`] tail produced by the +//! BAM source preamble and emits a [`DecompressedBlock`] of FASTQ text per +//! batch, reusing [`write_fastq_record`] (the same encoder the standalone +//! command used). It is a `Parallel` step, so the work-stealing pool spreads +//! the encode across N workers; `add_sink` then routes the byte tail to a +//! (optionally BGZF-compressing) raw-file/stdout writer. + +use std::io; +use std::sync::Arc; +use std::sync::atomic::{AtomicU64, Ordering}; + +use anyhow::Result; +use log::info; + +use crate::commands::fastq::{FastqOptions, Segment, classify_segment, write_fastq_record}; +use crate::logging::OperationTimer; +use crate::pipeline::chains::FinalizeHook; +use crate::pipeline::steps::process::{ + Process3Output, Process3WithWorkerState, ProcessWithWorkerState, process_with_worker_state, + process3_with_worker_state, +}; +use crate::pipeline::steps::types::{DecodedRecordBatch, DecompressedBlock}; + +/// Post-pipeline finalize hook for the fastq stage: logs the record count and +/// the per-stage wall-time line. +pub(crate) struct FastqFinalizeHook { + pub(crate) records_written: Arc, + pub(crate) timer: OperationTimer, +} + +impl FinalizeHook for FastqFinalizeHook { + fn finalize(self: Box) -> Result<()> { + let FastqFinalizeHook { records_written, timer } = *self; + let n = records_written.load(Ordering::Relaxed); + info!("Fastq: completed ({n} FASTQ records written)"); + timer.log_completion(n); + Ok(()) + } +} + +/// Per-worker scratch for the FASTQ encoder: the reusable sequence and quality +/// decode buffers `write_fastq_record` fills per record. +pub(crate) struct FastqEncodeState { + seq_buf: Vec, + qual_buf: Vec, +} + +impl crate::pipeline::core::item::HeapSize for FastqEncodeState {} + +fn new_encode_state() -> FastqEncodeState { + FastqEncodeState { seq_buf: Vec::with_capacity(512), qual_buf: Vec::with_capacity(512) } +} + +/// Bump the shared record counter and log a progress line every million records. +fn bump_progress(counter: &AtomicU64, written: u64) { + if written == 0 { + return; + } + let prev = counter.fetch_add(written, Ordering::Relaxed); + if (prev + written) / 1_000_000 > prev / 1_000_000 { + info!("Wrote {} FASTQ records", prev + written); + } +} + +/// Interleaved FASTQ encode step: [`DecodedRecordBatch`] → one +/// [`DecompressedBlock`] of FASTQ text (all passing records, in input order). +pub(crate) fn build_fastq_encode_step( + opts: Arc, + per_step_byte_limit: u64, + records_written: Arc, +) -> ProcessWithWorkerState< + DecodedRecordBatch, + DecompressedBlock, + impl Fn(&mut FastqEncodeState, DecodedRecordBatch) -> io::Result + + Send + + Sync + + 'static, + FastqEncodeState, + impl Fn() -> FastqEncodeState + Send + Sync + 'static, +> { + process_with_worker_state::( + "FastqEncode", + per_step_byte_limit, + new_encode_state, + move |state: &mut FastqEncodeState, + batch: DecodedRecordBatch| + -> io::Result { + let DecodedRecordBatch { batch_serial, records, .. } = batch; + let mut out = Vec::with_capacity(records.len().saturating_mul(512)); + let mut written = 0u64; + for record in &records { + let flags = record.record().flags(); + if !opts.passes_filters(flags) { + continue; + } + write_fastq_record( + &mut out, + record.record(), + flags, + opts.no_suffix, + opts.umi_header.as_ref(), + &mut state.seq_buf, + &mut state.qual_buf, + ) + .map_err(|e| io::Error::other(format!("fastq encode: {e:#}")))?; + written += 1; + } + bump_progress(&records_written, written); + Ok(DecompressedBlock { batch_serial, bytes: out }) + }, + ) +} + +/// Routes and encodes a single [`DecodedRecordBatch`] into R1/R2/other FASTQ +/// byte blocks. +/// +/// Routing matches `samtools fastq -1/-2/-0`: R1 = `FIRST_SEGMENT` set / +/// `LAST_SEGMENT` clear, R2 = the reverse, "other" = single-end or ambiguous +/// (see [`classify_segment`]). Split mode always omits the `/1`/`/2` suffix. +/// +/// **Emits a block on all three branches for every batch** — an empty +/// `Vec` where a branch has no records this batch — to keep each ordered +/// branch's `batch_serial` sequence gap-free (a skipped serial would wedge that +/// branch's reorder stage). Never returns `None` on a branch. +/// +/// Extracted from [`build_fastq_paired_encode_step`]'s closure so the routing +/// logic is directly unit-testable without going through the pipeline +/// machinery (see `mod tests`). +fn encode_paired_batch( + state: &mut FastqEncodeState, + batch: DecodedRecordBatch, + opts: &FastqOptions, + records_written: &AtomicU64, +) -> io::Result> { + let DecodedRecordBatch { batch_serial, records, .. } = batch; + let mut r1: Vec = Vec::new(); + let mut r2: Vec = Vec::new(); + let mut other: Vec = Vec::new(); + let mut written = 0u64; + for record in &records { + let flags = record.record().flags(); + if !opts.passes_filters(flags) { + continue; + } + // Split mode: the R1/R2 file already identifies the mate, so the + // `/1`/`/2` suffix is always omitted (matches `samtools fastq -1/-2`). + let out = match classify_segment(flags) { + Segment::Read1 => &mut r1, + Segment::Read2 => &mut r2, + Segment::Other => &mut other, + }; + write_fastq_record( + out, + record.record(), + flags, + /* no_suffix */ true, + opts.umi_header.as_ref(), + &mut state.seq_buf, + &mut state.qual_buf, + ) + .map_err(|e| io::Error::other(format!("fastq encode: {e:#}")))?; + written += 1; + } + bump_progress(records_written, written); + // Always emit all three branches (empty Vec where none) — no gaps. + Ok(Process3Output::both3( + DecompressedBlock { batch_serial, bytes: r1 }, + DecompressedBlock { batch_serial, bytes: r2 }, + DecompressedBlock { batch_serial, bytes: other }, + )) +} + +/// Paired-split FASTQ encode step: [`DecodedRecordBatch`] → three +/// [`DecompressedBlock`]s of FASTQ text (R1, R2, other), one per output file. +/// +/// See [`encode_paired_batch`] for the per-batch routing/encode logic; this +/// function only wires that logic into the `Process3WithWorkerState` step +/// shape. +#[allow(clippy::type_complexity)] +pub(crate) fn build_fastq_paired_encode_step( + opts: Arc, + per_step_byte_limit: u64, + records_written: Arc, +) -> Process3WithWorkerState< + DecodedRecordBatch, + DecompressedBlock, + DecompressedBlock, + DecompressedBlock, + FastqEncodeState, + impl Fn() -> FastqEncodeState + Send + Sync + 'static, + impl Fn( + &mut FastqEncodeState, + DecodedRecordBatch, + ) -> io::Result> + + Send + + Sync + + 'static, +> { + process3_with_worker_state::< + DecodedRecordBatch, + DecompressedBlock, + DecompressedBlock, + DecompressedBlock, + FastqEncodeState, + _, + _, + >( + "FastqEncodePaired", + per_step_byte_limit, + per_step_byte_limit, + per_step_byte_limit, + new_encode_state, + move |state: &mut FastqEncodeState, batch: DecodedRecordBatch| { + encode_paired_batch(state, batch, &opts, &records_written) + }, + ) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn paired_encode_step_profile_is_parallel_three_branches() { + use crate::commands::fastq::FastqOptions; + use crate::pipeline::core::step::{Step, StepKind}; + use std::sync::Arc; + use std::sync::atomic::AtomicU64; + + let opts = Arc::new(FastqOptions { + exclude_flags: 0, + require_flags: 0, + no_suffix: true, + umi_header: None, + }); + let step = build_fastq_paired_encode_step(opts, 1 << 20, Arc::new(AtomicU64::new(0))); + let profile = step.profile(); + assert_eq!(profile.kind, StepKind::Parallel); + assert_eq!(profile.output_queues.len(), 3); + assert_eq!(profile.branch_ordering.len(), 3); + } + + /// Build a `DecodedRecord` for a single-segment (unpaired-flags-wise + /// simple) paired read: `PAIRED` plus exactly one of `FIRST_SEGMENT` / + /// `LAST_SEGMENT`, matching how `classify_segment` routes R1 vs R2. + fn decoded_record(name: &str, first_segment: bool) -> fgumi_bam_io::DecodedRecord { + use fgumi_bam_io::{DecodedRecord, GroupKey}; + use fgumi_raw_bam::{SamBuilder as RawSamBuilder, flags as raw_flags}; + + let mut flags = raw_flags::PAIRED; + flags |= if first_segment { raw_flags::FIRST_SEGMENT } else { raw_flags::LAST_SEGMENT }; + let mut b = RawSamBuilder::new(); + b.read_name(name.as_bytes()) + .sequence(b"ACGTACGT") + .qualities(&[30; 8]) + .flags(flags) + .ref_id(0) + .pos(99) + .mapq(60); + DecodedRecord::from_raw_bytes(b.build(), GroupKey::default()) + } + + fn no_filter_opts() -> FastqOptions { + FastqOptions { exclude_flags: 0, require_flags: 0, no_suffix: true, umi_header: None } + } + + /// Direct, deterministic regression test for the no-ordinal-gaps rule + /// (see the doc comment on `encode_paired_batch`): a batch containing + /// ONLY R1 records must still return `Some(empty block)` on the R2 and + /// "other" branches, never `None`. Unlike the integration-level + /// `test_fastq_paired_empty_r2_batch_no_hang` (which depends on the + /// fixture fitting in a single BGZF block to avoid a real hang), this + /// calls the routing/encode logic directly and asserts on `Process3Output` + /// itself, so it fails immediately — no subprocess, no timing, no + /// block-size heuristics — if a regression ever changes `b`/`c` to `None`. + #[test] + fn encode_paired_batch_emits_empty_block_not_none_when_branch_has_no_records() { + use std::sync::atomic::AtomicU64; + + let opts = no_filter_opts(); + let mut state = new_encode_state(); + let records_written = AtomicU64::new(0); + + let batch = DecodedRecordBatch::new(0, vec![decoded_record("read1", true)]); + let out = encode_paired_batch(&mut state, batch, &opts, &records_written) + .expect("encode_paired_batch should succeed"); + + let r1 = out.a.expect("R1 branch must be Some: this batch has an R1 record"); + assert!(!r1.bytes.is_empty(), "R1 branch must contain the encoded R1 record"); + + let r2 = out.b.expect( + "R2 branch must be Some (empty block) even though this batch has no R2 \ + records -- returning None here is exactly the regression that wedges \ + the R2 reorder stage", + ); + assert!(r2.bytes.is_empty(), "R2 branch must be an empty block, not populated"); + + let other = out.c.expect( + "other branch must be Some (empty block) even though this batch has no \ + single-end/ambiguous records", + ); + assert!(other.bytes.is_empty(), "other branch must be an empty block, not populated"); + } + + /// Sanity check for the routing side of `encode_paired_batch`: a batch + /// with one record on each branch produces non-empty bytes on all three, + /// so the "empty means absent" assertions above aren't trivially true for + /// every branch regardless of content. + #[test] + fn encode_paired_batch_routes_mixed_batch_to_all_three_branches() { + use std::sync::atomic::AtomicU64; + + let opts = no_filter_opts(); + let mut state = new_encode_state(); + let records_written = AtomicU64::new(0); + + // "Other" = neither/both segment flags set; use an unpaired single-end read. + let mut other_builder = fgumi_raw_bam::SamBuilder::new(); + other_builder + .read_name(b"single1") + .sequence(b"GGGGCCCC") + .qualities(&[30; 8]) + .flags(0) + .ref_id(0) + .pos(199) + .mapq(60); + let other_record = fgumi_bam_io::DecodedRecord::from_raw_bytes( + other_builder.build(), + fgumi_bam_io::GroupKey::default(), + ); + + let batch = DecodedRecordBatch::new( + 0, + vec![decoded_record("read1", true), decoded_record("read1", false), other_record], + ); + let out = encode_paired_batch(&mut state, batch, &opts, &records_written) + .expect("encode_paired_batch should succeed"); + + assert!(out.a.is_some_and(|b| !b.bytes.is_empty()), "R1 branch should be non-empty"); + assert!(out.b.is_some_and(|b| !b.bytes.is_empty()), "R2 branch should be non-empty"); + assert!(out.c.is_some_and(|b| !b.bytes.is_empty()), "other branch should be non-empty"); + } +} diff --git a/src/lib/pipeline/chains/commands/mod.rs b/src/lib/pipeline/chains/commands/mod.rs index b6f7f5ce4..b0927dbbb 100644 --- a/src/lib/pipeline/chains/commands/mod.rs +++ b/src/lib/pipeline/chains/commands/mod.rs @@ -16,6 +16,7 @@ pub mod dedup; #[cfg(feature = "consensus")] pub mod duplex; pub mod extract; +pub mod fastq; pub mod filter; pub mod group; #[cfg(feature = "consensus")] diff --git a/src/lib/pipeline/chains/options_bag.rs b/src/lib/pipeline/chains/options_bag.rs index d95833762..7b55e7b15 100644 --- a/src/lib/pipeline/chains/options_bag.rs +++ b/src/lib/pipeline/chains/options_bag.rs @@ -14,6 +14,7 @@ pub use crate::commands::dedup::MarkDuplicates; #[cfg(feature = "consensus")] pub use crate::commands::duplex::DuplexOptions; pub use crate::commands::extract::ExtractOptions; +pub use crate::commands::fastq::FastqOptions; pub use crate::commands::filter::FilterOptions; pub use crate::commands::group::GroupOptions; #[cfg(feature = "consensus")] @@ -92,6 +93,9 @@ pub struct StageOptionsBag { /// Extract options. Carries header-synthesis fields and behavior knobs /// for the FASTQ→unmapped-BAM extract stage (Phase 5 T5.3). pub extract: Option, + /// Fastq options. Carries flag filters, read-name suffix, and UMI-in-name + /// config for the terminal BAM→FASTQ encode stage. + pub fastq: Option, // Additional fields added during command migrations: // - downsample (T2.13). // Cross-stage fields (read_group, methylation_mode, reference, …) diff --git a/src/lib/pipeline/chains/sink_spec.rs b/src/lib/pipeline/chains/sink_spec.rs index baaac985c..24c7e1c02 100644 --- a/src/lib/pipeline/chains/sink_spec.rs +++ b/src/lib/pipeline/chains/sink_spec.rs @@ -31,14 +31,27 @@ pub enum SinkSpec { /// BAM file write + `.bai` index write after the BAM finishes. /// Only valid when the final stage emits coordinate-sorted output. BamWithIndex(PathBuf), + /// Interleaved FASTQ file write (or `-` for stdout). A `.gz`/`.bgz` + /// path is written as BGZF. Only valid on a `Stage::Fastq`-terminal chain. + /// + /// Paired split output uses the [`SinkSpec::FastqPaired`] variant. + Fastq(PathBuf), + /// Paired split FASTQ output: R1 → `out1`, R2 → `out2`, and "other" + /// (single-end / ambiguous) reads → `out0` if given, else stdout. Each + /// `.gz`/`.bgz` path is written as BGZF. Only valid on a `Stage::Fastq`- + /// terminal chain. The three streams are fanned out by the paired FASTQ + /// encode step (a 3-output `Process3WithWorkerState`). + FastqPaired { out1: PathBuf, out2: PathBuf, out0: Option }, } impl SinkSpec { - /// The output path, regardless of variant. + /// A representative output path for the sink, used for logging and the + /// input-clobber guard. #[must_use] pub fn path(&self) -> &PathBuf { match self { - SinkSpec::Bam(p) | SinkSpec::BamWithIndex(p) => p, + SinkSpec::Bam(p) | SinkSpec::BamWithIndex(p) | SinkSpec::Fastq(p) => p, + SinkSpec::FastqPaired { out1, .. } => out1, } } } diff --git a/src/lib/pipeline/chains/stage.rs b/src/lib/pipeline/chains/stage.rs index 76f676008..c93b72eb0 100644 --- a/src/lib/pipeline/chains/stage.rs +++ b/src/lib/pipeline/chains/stage.rs @@ -27,6 +27,8 @@ pub enum Stage { Filter, Dedup, Downsample, + /// Terminal BAM → FASTQ encode (interleaved or paired split output). + Fastq, } impl Stage { diff --git a/src/lib/pipeline/chains/validate.rs b/src/lib/pipeline/chains/validate.rs index 4b075152d..f81d069a5 100644 --- a/src/lib/pipeline/chains/validate.rs +++ b/src/lib/pipeline/chains/validate.rs @@ -40,6 +40,9 @@ fn stage_ord(stage: Stage) -> usize { Stage::Clip | Stage::Dedup | Stage::Downsample => 6, Stage::Simplex | Stage::Duplex | Stage::Codec => 7, Stage::Filter => 8, + // Terminal BAM→FASTQ export. Only ever appears as the sole stage of a + // `fgumi fastq` chain, so its rank is nominal; place it last. + Stage::Fastq => 9, } } @@ -187,6 +190,7 @@ pub fn validate_stage_opts_present(spec: &ChainSpec) -> Result<()> { } Stage::Align => bag.aligner.is_some(), Stage::Extract => bag.extract.is_some(), + Stage::Fastq => bag.fastq.is_some(), // Downsample is rejected by the early `bail!` above (no bag slot // yet). Report "not provided" rather than `unreachable!()` so that // if the early guard is ever removed without adding a slot here, the @@ -308,6 +312,23 @@ pub fn validate_cross_stage_constraints(spec: &ChainSpec) -> Result<()> { bail!("Stage::Extract requires SourceSpec::Fastqs; got {:?}", spec.source); } + // Rule 5: the FASTQ sink and `Stage::Fastq` must go together. A `Fastq` + // sink pushes raw bytes through `add_fastq_sink` with no BAM header, and + // `Stage::Fastq` emits FASTQ text, not serialized BAM records — pairing + // either with the other format's counterpart would silently write a + // corrupt file (FASTQ text wrapped in a BAM header, or BAM bytes with no + // header). Require the biconditional. + let sink_is_fastq = matches!(spec.sink, SinkSpec::Fastq(_) | SinkSpec::FastqPaired { .. }); + let terminal_is_fastq = spec.stages.last() == Some(&Stage::Fastq); + if sink_is_fastq != terminal_is_fastq { + bail!( + "SinkSpec::Fastq requires Stage::Fastq as the terminal chain stage and vice versa; \ + got sink={:?}, terminal stage={:?}", + spec.sink, + spec.stages.last() + ); + } + Ok(()) } @@ -668,4 +689,31 @@ mod tests { validate_stage_opts_present(&spec) .expect("Stage::Extract with options populated should pass"); } + + // ── Rule 5: SinkSpec::{Fastq,FastqPaired} <-> Stage::Fastq biconditional ─ + + #[test] + fn cross_stage_fastq_paired_accepted_on_terminal_fastq() { + let mut spec = empty_spec(vec![Stage::Fastq]); + spec.sink = SinkSpec::FastqPaired { + out1: std::path::PathBuf::from("r1.fq.gz"), + out2: std::path::PathBuf::from("r2.fq.gz"), + out0: Some(std::path::PathBuf::from("other.fq.gz")), + }; + validate_cross_stage_constraints(&spec) + .expect("FastqPaired + terminal Stage::Fastq must be valid"); + } + + #[test] + fn cross_stage_fastq_paired_rejected_on_non_fastq_terminal() { + // Sort -> Group (Group is terminal, not Fastq). + let mut spec = empty_spec(vec![Stage::Sort, Stage::Group]); + spec.sink = SinkSpec::FastqPaired { + out1: std::path::PathBuf::from("r1.fq"), + out2: std::path::PathBuf::from("r2.fq"), + out0: None, + }; + let err = validate_cross_stage_constraints(&spec).unwrap_err(); + assert!(err.to_string().contains("Stage::Fastq"), "got: {err}"); + } } diff --git a/src/lib/pipeline/steps/align_and_merge.rs b/src/lib/pipeline/steps/align_and_merge.rs index 695f89c30..58a908403 100644 --- a/src/lib/pipeline/steps/align_and_merge.rs +++ b/src/lib/pipeline/steps/align_and_merge.rs @@ -666,6 +666,7 @@ fn writer_loop( record, flags, /* no_suffix */ true, + /* umi_header */ None, &mut seq_buf, &mut qual_buf, ) { diff --git a/src/lib/pipeline/steps/process.rs b/src/lib/pipeline/steps/process.rs index 0abbf4e0f..0df59da48 100644 --- a/src/lib/pipeline/steps/process.rs +++ b/src/lib/pipeline/steps/process.rs @@ -599,6 +599,238 @@ where Process2WithWorkerState::new(name, limit_bytes_a, limit_bytes_b, init, f) } +// ───────────────────────────────────────────────────────────────────────────── +// Process3WithWorkerState — Parallel + lazy per-worker +// state + 3-output fan-out. Mirror of `Process2WithWorkerState` widened to +// three ordered + byte-bounded output branches. Used by the paired FASTQ +// encode (R1 / R2 / other), where each worker reuses one scratch buffer set +// across all batches. Each branch carries its own `HeldSlot` for independent +// backpressure, exactly as the 2-output sibling. +// ───────────────────────────────────────────────────────────────────────────── + +/// Output of a [`Process3WithWorkerState`] closure: one optional value per +/// output branch. For an *ordered* consumer every branch that participates in +/// the run must receive an item for every input ordinal (an empty payload is +/// fine); `None` leaves a gap that wedges that branch's reorder stage. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct Process3Output { + pub a: Option, + pub b: Option, + pub c: Option, +} + +impl Process3Output { + #[must_use] + pub fn both3(a: A, b: B, c: C) -> Self { + Self { a: Some(a), b: Some(b), c: Some(c) } + } +} + +/// `Parallel` + lazy per-worker state + 3-output fan-out. `Clone` resets +/// `state` to `None`; first `try_run` initializes via `init()`; subsequent +/// calls reuse the stored `S` via `&mut`. Each worker owns its clone +/// exclusively (Parallel kind), so no synchronization is needed. Mirror of +/// [`Process2WithWorkerState`] with a third branch. +pub struct Process3WithWorkerState +where + I: Send + HeapSize + 'static, + A: Send + HeapSize + Ordered + 'static, + B: Send + HeapSize + Ordered + 'static, + C: Send + HeapSize + Ordered + 'static, + S: Send + HeapSize + 'static, + Init: Fn() -> S + Send + Sync + 'static, + F: Fn(&mut S, I) -> io::Result> + Send + Sync + 'static, +{ + f: Arc, + init: Arc, + state: Option, + held_a: HeldSlot>, + held_b: HeldSlot>, + held_c: HeldSlot>, + name: &'static str, + limit_bytes_a: u64, + limit_bytes_b: u64, + limit_bytes_c: u64, + _phantom: PhantomData I>, +} + +impl Process3WithWorkerState +where + I: Send + HeapSize + 'static, + A: Send + HeapSize + Ordered + 'static, + B: Send + HeapSize + Ordered + 'static, + C: Send + HeapSize + Ordered + 'static, + S: Send + HeapSize + 'static, + Init: Fn() -> S + Send + Sync + 'static, + F: Fn(&mut S, I) -> io::Result> + Send + Sync + 'static, +{ + pub fn new( + name: &'static str, + limit_bytes_a: u64, + limit_bytes_b: u64, + limit_bytes_c: u64, + init: Init, + f: F, + ) -> Self { + Self { + f: Arc::new(f), + init: Arc::new(init), + state: None, + held_a: HeldSlot::new(), + held_b: HeldSlot::new(), + held_c: HeldSlot::new(), + name, + limit_bytes_a, + limit_bytes_b, + limit_bytes_c, + _phantom: PhantomData, + } + } +} + +impl Clone for Process3WithWorkerState +where + I: Send + HeapSize + 'static, + A: Send + HeapSize + Ordered + 'static, + B: Send + HeapSize + Ordered + 'static, + C: Send + HeapSize + Ordered + 'static, + S: Send + HeapSize + 'static, + Init: Fn() -> S + Send + Sync + 'static, + F: Fn(&mut S, I) -> io::Result> + Send + Sync + 'static, +{ + fn clone(&self) -> Self { + Self { + f: Arc::clone(&self.f), + init: Arc::clone(&self.init), + state: None, // per-worker fresh; lazy on first try_run + held_a: HeldSlot::new(), + held_b: HeldSlot::new(), + held_c: HeldSlot::new(), + name: self.name, + limit_bytes_a: self.limit_bytes_a, + limit_bytes_b: self.limit_bytes_b, + limit_bytes_c: self.limit_bytes_c, + _phantom: PhantomData, + } + } +} + +impl Step for Process3WithWorkerState +where + I: Send + HeapSize + 'static, + A: Send + HeapSize + Ordered + 'static, + B: Send + HeapSize + Ordered + 'static, + C: Send + HeapSize + Ordered + 'static, + S: Send + HeapSize + 'static, + Init: Fn() -> S + Send + Sync + 'static, + F: Fn(&mut S, I) -> io::Result> + Send + Sync + 'static, +{ + type Input = I; + type Outputs = crate::pipeline::core::outputs::OrderedBytesTuple3; + + fn profile(&self) -> StepProfile { + StepProfile { + name: self.name, + kind: StepKind::Parallel, + sticky: false, + output_queues: vec![ + QueueSpec::ByteBounded { limit_bytes: self.limit_bytes_a }, + QueueSpec::ByteBounded { limit_bytes: self.limit_bytes_b }, + QueueSpec::ByteBounded { limit_bytes: self.limit_bytes_c }, + ], + branch_ordering: vec![ + BranchOrdering::ByItemOrdinal, + BranchOrdering::ByItemOrdinal, + BranchOrdering::ByItemOrdinal, + ], + } + } + + fn try_run(&mut self, ctx: &mut StepCtx<'_, Self>) -> io::Result { + let view = ctx.outputs.view(); + + if let Some(unpushed) = self.held_a.take() { + match view.a.retry(unpushed) { + Ok(()) => {} + Err(again) => { + self.held_a.put(again); + return Ok(StepOutcome::Contention); + } + } + } + if let Some(unpushed) = self.held_b.take() { + match view.b.retry(unpushed) { + Ok(()) => {} + Err(again) => { + self.held_b.put(again); + return Ok(StepOutcome::Contention); + } + } + } + if let Some(unpushed) = self.held_c.take() { + match view.c.retry(unpushed) { + Ok(()) => {} + Err(again) => { + self.held_c.put(again); + return Ok(StepOutcome::Contention); + } + } + } + + let Some(item) = ctx.input.pop() else { + if ctx.input.is_drained() { + return Ok(StepOutcome::Finished); + } + return Ok(StepOutcome::NoProgress); + }; + let state = self.state.get_or_insert_with(|| (self.init)()); + let Process3Output { a, b, c } = (self.f)(state, item)?; + + if let Some(av) = a { + if let Err(unpushed) = view.a.push(av) { + self.held_a.put(unpushed); + } + } + if let Some(bv) = b { + if let Err(unpushed) = view.b.push(bv) { + self.held_b.put(unpushed); + } + } + if let Some(cv) = c { + if let Err(unpushed) = view.c.push(cv) { + self.held_c.put(unpushed); + } + } + Ok(StepOutcome::Progress) + } + + fn new_worker_copy(&self) -> Self { + self.clone() + } +} + +/// Convenience constructor for [`Process3WithWorkerState`]. +#[allow(clippy::too_many_arguments)] +pub fn process3_with_worker_state( + name: &'static str, + limit_bytes_a: u64, + limit_bytes_b: u64, + limit_bytes_c: u64, + init: Init, + f: F, +) -> Process3WithWorkerState +where + I: Send + HeapSize + 'static, + A: Send + HeapSize + Ordered + 'static, + B: Send + HeapSize + Ordered + 'static, + C: Send + HeapSize + Ordered + 'static, + S: Send + HeapSize + 'static, + Init: Fn() -> S + Send + Sync + 'static, + F: Fn(&mut S, I) -> io::Result> + Send + Sync + 'static, +{ + Process3WithWorkerState::new(name, limit_bytes_a, limit_bytes_b, limit_bytes_c, init, f) +} + // ───────────────────────────────────────────────────────────────────────────── // MiAssign — Serial + monotonic MI counter. Used for deterministic // MoleculeId numbering across input batches. @@ -1398,6 +1630,79 @@ mod tests { ); } + #[test] + fn process3_with_worker_state_clone_resets_state() { + #[derive(Default)] + struct WS { + counter: u32, + } + impl HeapSize for WS { + fn heap_size(&self) -> usize { + 0 + } + } + + let step = process3_with_worker_state::( + "stateful3", + 1024, + 1024, + 1024, + WS::default, + |state, item| { + state.counter += 1; + Ok(Process3Output::both3( + Item { ord: item.ord, v: state.counter }, + Item { ord: item.ord, v: state.counter }, + Item { ord: item.ord, v: state.counter }, + )) + }, + ); + let cloned = step.clone(); + assert!(step.state.is_none()); + assert!(cloned.state.is_none()); + } + + #[test] + fn process3_with_worker_state_profile_is_parallel_byitem_three_branches() { + #[derive(Default)] + struct WS; + impl HeapSize for WS { + fn heap_size(&self) -> usize { + 0 + } + } + + let step = process3_with_worker_state::( + "stateful3-profile", + 256, + 512, + 1024, + WS::default, + |_state, item| { + Ok(Process3Output::both3(item, Item { ord: 0, v: 0 }, Item { ord: 0, v: 0 })) + }, + ); + let profile = step.profile(); + assert_eq!(profile.kind, StepKind::Parallel); + assert_eq!(profile.output_queues.len(), 3); + assert_eq!( + profile.output_queues, + vec![ + QueueSpec::ByteBounded { limit_bytes: 256 }, + QueueSpec::ByteBounded { limit_bytes: 512 }, + QueueSpec::ByteBounded { limit_bytes: 1024 }, + ] + ); + assert_eq!( + profile.branch_ordering, + vec![ + BranchOrdering::ByItemOrdinal, + BranchOrdering::ByItemOrdinal, + BranchOrdering::ByItemOrdinal, + ] + ); + } + #[test] fn mi_assign_profile_is_serial() { let s: MiAssign = mi_assign("mi", 1024, |next_mi, item: Item| { diff --git a/src/lib/pipeline/steps/sink/mod.rs b/src/lib/pipeline/steps/sink/mod.rs index 7a40a0704..6cda04b6b 100644 --- a/src/lib/pipeline/steps/sink/mod.rs +++ b/src/lib/pipeline/steps/sink/mod.rs @@ -1,3 +1,4 @@ //! Sink steps (BAM file writers, future FASTQ writers). pub mod write_bgzf; +pub mod write_raw; diff --git a/src/lib/pipeline/steps/sink/write_raw.rs b/src/lib/pipeline/steps/sink/write_raw.rs new file mode 100644 index 000000000..d2e670a12 --- /dev/null +++ b/src/lib/pipeline/steps/sink/write_raw.rs @@ -0,0 +1,3 @@ +//! Re-export shim. The `WriteRawFile` sink step lives in the +//! `fgumi-pipeline-io` crate. +pub use fgumi_pipeline_io::sink::write_raw::{RawBytesBlock, WriteRawFile}; diff --git a/tests/integration/test_fastq_command.rs b/tests/integration/test_fastq_command.rs index 98a06c541..be55687c7 100644 --- a/tests/integration/test_fastq_command.rs +++ b/tests/integration/test_fastq_command.rs @@ -6,15 +6,61 @@ use clap::Parser; use fgumi_lib::commands::command::Command; use fgumi_lib::commands::fastq::Fastq; use fgumi_raw_bam::{SamBuilder, flags}; +use fgumi_tag::SamTag; use noodles::bam; use noodles::sam::alignment::io::Write as AlignmentWrite; use std::fs; -use std::io::{BufRead, BufReader}; -use std::path::PathBuf; +use std::io::{BufRead, BufReader, Read}; +use std::path::{Path, PathBuf}; +use std::process::{Command as ProcessCommand, Stdio}; +use std::thread; +use std::time::{Duration, Instant}; use tempfile::TempDir; use crate::helpers::bam_generator::{create_minimal_header, to_record_buf}; +/// Create a BAM with a single read pair carrying a UMI in the given tag. +/// +/// `umi_tag` is `None` to write no UMI tag at all (passthrough test). +fn create_pair_with_umi_tag(path: &Path, name: &str, umi_tag: Option<(SamTag, &str)>) { + let header = create_minimal_header("chr1", 10000); + let mut writer = + bam::io::Writer::new(fs::File::create(path).expect("Failed to create BAM file")); + writer.write_header(&header).expect("Failed to write header"); + + for (seq, segment_flag) in + [("ACGTACGT", flags::FIRST_SEGMENT), ("TGCATGCA", flags::LAST_SEGMENT)] + { + let mut b = SamBuilder::new(); + b.read_name(name.as_bytes()) + .sequence(seq.as_bytes()) + .qualities(&[40, 40, 40, 40, 40, 40, 40, 40]) + .flags(flags::PAIRED | segment_flag) + .ref_id(0) + .pos(99) + .mapq(60); + if let Some((tag, value)) = umi_tag { + b.add_string_tag(tag, value.as_bytes()); + } + writer + .write_alignment_record(&header, &to_record_buf(&b.build())) + .expect("Failed to write record"); + } + writer.try_finish().expect("Failed to finish BAM"); +} + +/// Run `fgumi fastq` with the given extra args and return parsed FASTQ records. +fn run_fastq_with_args(input_bam: &Path, extra: &[&str]) -> Vec<(String, String, String)> { + let temp = input_bam.parent().expect("parent"); + let output_fq = temp.join("umi_out.fq"); + let mut args = + vec!["fastq", "-i", input_bam.to_str().unwrap(), "-o", output_fq.to_str().unwrap()]; + args.extend_from_slice(extra); + let cmd = Fastq::try_parse_from(args).expect("failed to parse fastq args"); + cmd.execute("fgumi fastq").expect("fastq command failed"); + parse_fastq_records(&output_fq) +} + /// Create a BAM file with paired-end reads for testing. fn create_paired_bam(path: &PathBuf, read_pairs: Vec<(&str, &str, &str, &str, &str, bool)>) { let header = create_minimal_header("chr1", 10000); @@ -370,6 +416,387 @@ fn test_fastq_hex_flags() { assert_eq!(records.len(), 2); } +/// `--umi-in-header` with defaults reproduces `samtools fastq -U`: the duplex +/// UMI stored as `AAAAAAAA-CCCCCCCC` in `RX` becomes `:AAAAAAAA+CCCCCCCC` +/// appended before the /1 /2 suffix. +#[test] +fn test_fastq_umi_in_header_dragen_format() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + create_pair_with_umi_tag(&input_bam, "read1", Some((SamTag::RX, "AAAAAAAA-CCCCCCCC"))); + + let records = run_fastq_with_args(&input_bam, &["--annotate-read-names"]); + assert_eq!(records.len(), 2); + assert_eq!(records[0].0, "read1:AAAAAAAA+CCCCCCCC/1"); + assert_eq!(records[1].0, "read1:AAAAAAAA+CCCCCCCC/2"); +} + +/// `-n` drops the /1 /2 suffix but keeps the UMI annotation. +#[test] +fn test_fastq_umi_in_header_no_suffix() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + create_pair_with_umi_tag(&input_bam, "read1", Some((SamTag::RX, "AAAAAAAA-CCCCCCCC"))); + + let records = run_fastq_with_args(&input_bam, &["--annotate-read-names", "-n"]); + assert_eq!(records[0].0, "read1:AAAAAAAA+CCCCCCCC"); + assert_eq!(records[1].0, "read1:AAAAAAAA+CCCCCCCC"); +} + +/// Custom delimiters: underscore name separator, empty UMI separator. +#[test] +fn test_fastq_umi_in_header_custom_delims() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + create_pair_with_umi_tag(&input_bam, "read1", Some((SamTag::RX, "AAAA-GGGG"))); + + let records = run_fastq_with_args( + &input_bam, + &["--annotate-read-names", "-n", "--umi-name-delim", "_", "--umi-sep", ""], + ); + assert_eq!(records[0].0, "read1_AAAAGGGG"); +} + +/// When `RX` is absent the fallback `OX` tag is used (default tag list `RX,OX`). +#[test] +fn test_fastq_umi_in_header_ox_fallback() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + create_pair_with_umi_tag(&input_bam, "read1", Some((SamTag::OX, "TTTT-GGGG"))); + + let records = run_fastq_with_args(&input_bam, &["--annotate-read-names", "-n"]); + assert_eq!(records[0].0, "read1:TTTT+GGGG"); +} + +/// With the default tag list `RX,OX`, `RX` takes priority when both are present. +#[test] +fn test_fastq_umi_in_header_rx_preferred_over_ox() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + // Write both RX and OX; RX must win. + let header = create_minimal_header("chr1", 10000); + let mut writer = bam::io::Writer::new(fs::File::create(&input_bam).expect("create BAM")); + writer.write_header(&header).expect("write header"); + let mut b = SamBuilder::new(); + b.read_name(b"read1") + .sequence(b"ACGT") + .qualities(&[40, 40, 40, 40]) + .flags(flags::PAIRED | flags::FIRST_SEGMENT) + .ref_id(0) + .pos(99) + .mapq(60); + b.add_string_tag(SamTag::RX, b"RRRR"); + b.add_string_tag(SamTag::OX, b"OOOO"); + writer.write_alignment_record(&header, &to_record_buf(&b.build())).expect("write"); + writer.try_finish().expect("finish"); + + let records = run_fastq_with_args(&input_bam, &["--annotate-read-names", "-n"]); + assert_eq!(records[0].0, "read1:RRRR"); +} + +/// A record with no UMI tag passes through with the name unchanged. +#[test] +fn test_fastq_umi_in_header_missing_tag_passthrough() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + create_pair_with_umi_tag(&input_bam, "read1", None); + + let records = run_fastq_with_args(&input_bam, &["--annotate-read-names", "-n"]); + assert_eq!(records[0].0, "read1"); +} + +/// Without `--umi-in-header`, the tag is ignored and the name is unchanged. +#[test] +fn test_fastq_umi_not_in_header_by_default() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + create_pair_with_umi_tag(&input_bam, "read1", Some((SamTag::RX, "AAAA-GGGG"))); + + let records = run_fastq_with_args(&input_bam, &["-n"]); + assert_eq!(records[0].0, "read1"); +} + +/// Read a FASTQ file, transparently decompressing BGZF/gzip (`.gz`/`.bgz`) via +/// multi-member gzip decoding (BGZF is a series of gzip blocks). Both extensions +/// trigger BGZF writing on the output side, so the reader honors both. +fn read_fastq_maybe_gz(path: &Path) -> Vec<(String, String, String)> { + let file = fs::File::open(path).expect("open fastq"); + let is_gz = matches!(path.extension().and_then(|e| e.to_str()), Some("gz" | "bgz")); + let reader: Box = if is_gz { + Box::new(BufReader::new(flate2::read::MultiGzDecoder::new(file))) + } else { + Box::new(BufReader::new(file)) + }; + let lines: Vec = reader.lines().map(|l| l.unwrap()).collect(); + let mut records = Vec::new(); + for chunk in lines.chunks(4) { + if chunk.len() == 4 { + records.push(( + chunk[0].trim_start_matches('@').to_string(), + chunk[1].clone(), + chunk[3].clone(), + )); + } + } + records +} + +/// Paired `-1`/`-2` output routes R1 to out1 and R2 to out2, with the UMI +/// annotation applied to both. +#[test] +fn test_fastq_paired_split_output() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + create_pair_with_umi_tag(&input_bam, "read1", Some((SamTag::RX, "AAAA-GGGG"))); + let r1 = temp_dir.path().join("R1.fq"); + let r2 = temp_dir.path().join("R2.fq"); + + let cmd = Fastq::try_parse_from([ + "fastq", + "-i", + input_bam.to_str().unwrap(), + "--annotate-read-names", + "-1", + r1.to_str().unwrap(), + "-2", + r2.to_str().unwrap(), + ]) + .expect("parse"); + cmd.execute("fgumi fastq").expect("run"); + + let recs1 = read_fastq_maybe_gz(&r1); + let recs2 = read_fastq_maybe_gz(&r2); + assert_eq!(recs1.len(), 1, "one R1 read"); + assert_eq!(recs2.len(), 1, "one R2 read"); + // Split mode omits the /1 /2 suffix (like samtools -1/-2); R1 and R2 names match. + assert_eq!(recs1[0].0, "read1:AAAA+GGGG"); + assert_eq!(recs2[0].0, "read1:AAAA+GGGG"); + // R1 sequence is the FIRST_SEGMENT read; R2 the LAST_SEGMENT read. + assert_eq!(recs1[0].1, "ACGTACGT"); + assert_eq!(recs2[0].1, "TGCATGCA"); +} + +/// A `.gz` output path is written as BGZF (gzip magic bytes + BGZF `BC` subfield) +/// and round-trips to the expected records. +#[test] +fn test_fastq_paired_output_is_bgzf() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + create_pair_with_umi_tag(&input_bam, "read1", None); + let r1 = temp_dir.path().join("R1.fq.gz"); + let r2 = temp_dir.path().join("R2.fq.gz"); + + let cmd = Fastq::try_parse_from([ + "fastq", + "-i", + input_bam.to_str().unwrap(), + "-1", + r1.to_str().unwrap(), + "-2", + r2.to_str().unwrap(), + ]) + .expect("parse"); + cmd.execute("fgumi fastq").expect("run"); + + // Verify BGZF framing: gzip magic 1f 8b, and the BGZF `BC` extra subfield. + let raw = fs::read(&r1).expect("read gz"); + assert!(raw.len() > 18, "non-empty bgzf"); + assert_eq!(&raw[0..2], &[0x1f, 0x8b], "gzip magic"); + assert_eq!(raw[3] & 0x04, 0x04, "FEXTRA flag set (BGZF has an extra field)"); + assert_eq!([raw[12], raw[13]], [b'B', b'C'], "BGZF BC subfield identifier"); + + // And it decodes to the expected read (split mode omits the /1 suffix). + let recs1 = read_fastq_maybe_gz(&r1); + assert_eq!(recs1.len(), 1); + assert_eq!(recs1[0].0, "read1"); +} + +/// The `.bgz` extension triggers BGZF writing on the same footing as `.gz`, +/// so a `.bgz`-suffixed paired output is BGZF-framed and round-trips. Guards +/// against a writer bug specific to `.bgz` extension detection. +#[test] +fn test_fastq_paired_output_bgz_extension_is_bgzf() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + create_pair_with_umi_tag(&input_bam, "read1", None); + let r1 = temp_dir.path().join("R1.fq.bgz"); + let r2 = temp_dir.path().join("R2.fq.bgz"); + + let cmd = Fastq::try_parse_from([ + "fastq", + "-i", + input_bam.to_str().unwrap(), + "-1", + r1.to_str().unwrap(), + "-2", + r2.to_str().unwrap(), + ]) + .expect("parse"); + cmd.execute("fgumi fastq").expect("run"); + + // Verify BGZF framing on the `.bgz` output: gzip magic + BGZF `BC` subfield. + let raw = fs::read(&r1).expect("read bgz"); + assert!(raw.len() > 18, "non-empty bgzf"); + assert_eq!(&raw[0..2], &[0x1f, 0x8b], "gzip magic"); + assert_eq!(raw[3] & 0x04, 0x04, "FEXTRA flag set (BGZF has an extra field)"); + assert_eq!([raw[12], raw[13]], [b'B', b'C'], "BGZF BC subfield identifier"); + + // And both `.bgz` mates decode to the expected read (split mode omits /1 /2). + let recs1 = read_fastq_maybe_gz(&r1); + let recs2 = read_fastq_maybe_gz(&r2); + assert_eq!(recs1.len(), 1); + assert_eq!(recs2.len(), 1); + assert_eq!(recs1[0].0, "read1"); + assert_eq!(recs2[0].0, "read1"); +} + +/// An interleaved `--output` ending in `.gz` is also written as BGZF. +#[test] +fn test_fastq_interleaved_gz_is_bgzf() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + create_pair_with_umi_tag(&input_bam, "read1", None); + let out = temp_dir.path().join("out.fq.gz"); + + let cmd = Fastq::try_parse_from([ + "fastq", + "-i", + input_bam.to_str().unwrap(), + "-o", + out.to_str().unwrap(), + ]) + .expect("parse"); + cmd.execute("fgumi fastq").expect("run"); + + let raw = fs::read(&out).expect("read gz"); + assert_eq!(&raw[0..2], &[0x1f, 0x8b], "gzip magic"); + assert_eq!([raw[12], raw[13]], [b'B', b'C'], "BGZF BC subfield"); + let recs = read_fastq_maybe_gz(&out); + assert_eq!(recs.len(), 2, "both mates interleaved"); + // Assert record identity, not just count: R1 then R2, with the /1 /2 suffix. + assert_eq!(recs[0].0, "read1/1"); + assert_eq!(recs[1].0, "read1/2"); + assert_eq!(recs[0].1, "ACGTACGT"); + assert_eq!(recs[1].1, "TGCATGCA"); +} + +/// `--out1` without `--out2` is rejected by clap (`requires`). +#[test] +fn test_fastq_out1_requires_out2() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + let r1 = temp_dir.path().join("R1.fq"); + let parsed = Fastq::try_parse_from([ + "fastq", + "-i", + input_bam.to_str().unwrap(), + "-1", + r1.to_str().unwrap(), + ]); + assert!(parsed.is_err(), "--out1 without --out2 must be rejected"); +} + +/// `-U` remains a working alias for `--annotate-read-names`. +#[test] +fn test_fastq_annotate_u_alias() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + create_pair_with_umi_tag(&input_bam, "read1", Some((SamTag::RX, "AAAA-GGGG"))); + let records = run_fastq_with_args(&input_bam, &["-U", "-n"]); + assert_eq!(records[0].0, "read1:AAAA+GGGG"); +} + +/// Two outputs pointing at the same file must be rejected before any writing, +/// so two writer streams cannot silently truncate/corrupt one file. +#[test] +fn test_fastq_duplicate_output_paths_rejected() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + create_pair_with_umi_tag(&input_bam, "read1", None); + let same = temp_dir.path().join("both.fq"); + + let cmd = Fastq::try_parse_from([ + "fastq", + "-i", + input_bam.to_str().unwrap(), + "-1", + same.to_str().unwrap(), + "-2", + same.to_str().unwrap(), + ]) + .expect("parse"); + let err = cmd.execute("fgumi fastq").expect_err("must reject --out1 == --out2"); + assert!(err.to_string().contains("used more than once"), "unexpected error: {err}"); + assert!(!same.exists(), "no output file should be created when paths collide"); +} + +/// Two lexically-different but equivalent output paths that don't exist yet — +/// here the same file reached through a real directory and a symlink to it — +/// must still be rejected. The dedup canonicalizes the (existing) parent +/// directory, so the aliases collide before two sink threads can truncate the +/// same not-yet-existing target. Without parent-canonicalization the raw +/// lexical paths differ and the collision slips through. +#[cfg(unix)] +#[test] +fn test_fastq_aliased_nonexistent_output_paths_rejected() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + create_pair_with_umi_tag(&input_bam, "read1", None); + + // `real/` exists; `link/` is a symlink to it. Same not-yet-existing file + // `r.fq` reached two lexically-distinct ways that resolve to one inode. + let real_dir = temp_dir.path().join("real"); + fs::create_dir(&real_dir).expect("mkdir real"); + let link_dir = temp_dir.path().join("link"); + std::os::unix::fs::symlink(&real_dir, &link_dir).expect("symlink"); + let out1 = real_dir.join("r.fq"); + let out2 = link_dir.join("r.fq"); + assert!(!out1.exists(), "target must not exist yet for this to test the non-existent path"); + assert_ne!(out1, out2, "the two paths must be lexically distinct"); + + let cmd = Fastq::try_parse_from([ + "fastq", + "-i", + input_bam.to_str().unwrap(), + "-1", + out1.to_str().unwrap(), + "-2", + out2.to_str().unwrap(), + ]) + .expect("parse"); + let err = cmd.execute("fgumi fastq").expect_err("must reject aliased non-existent outputs"); + assert!(err.to_string().contains("used more than once"), "unexpected error: {err}"); + assert!(!out1.exists(), "no output file should be created when paths collide"); +} + +/// `/dev/null` is exempt from the duplicate-path check: routing several streams +/// to it (e.g. `--out0 /dev/null`) is intentional, as with `samtools fastq`. +#[cfg(unix)] +#[test] +fn test_fastq_dev_null_out0_allowed() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + create_pair_with_umi_tag(&input_bam, "read1", None); + let r1 = temp_dir.path().join("R1.fq"); + let r2 = temp_dir.path().join("R2.fq"); + + let cmd = Fastq::try_parse_from([ + "fastq", + "-i", + input_bam.to_str().unwrap(), + "-1", + r1.to_str().unwrap(), + "-2", + r2.to_str().unwrap(), + "-0", + "/dev/null", + ]) + .expect("parse"); + cmd.execute("fgumi fastq").expect("run with --out0 /dev/null"); + assert_eq!(read_fastq_maybe_gz(&r1).len(), 1); + assert_eq!(read_fastq_maybe_gz(&r2).len(), 1); +} + /// `--output` pointing at the same file as `--input` must fail before any /// truncation, since `File::create` would otherwise clobber the BAM data. #[test] @@ -433,3 +860,334 @@ fn test_fastq_output_symlink_to_input_rejected() { "input BAM was truncated/clobbered through symlinked --output" ); } + +// ── Paired output through the typed-step chain (Task 6) ───────────────── +// +// `test_fastq_paired_split_output` and `test_fastq_paired_output_is_bgzf` +// above already exercise the chain path (paired output no longer runs on the +// standalone `run_paired` loop); the tests below extend paired coverage to +// `--out0` routing, both `.gz` outputs, multi-pair UMI ordering, stdout +// fallback, and the empty-R2-branch no-ordinal-gaps regression. + +/// Create a BAM with one R1/R2 pair plus one single-end ("other") read — +/// exercises `--out0` routing (`samtools fastq -0` semantics). +fn create_paired_bam_with_other(path: &Path) { + let header = create_minimal_header("chr1", 10000); + let mut writer = + bam::io::Writer::new(fs::File::create(path).expect("Failed to create BAM file")); + writer.write_header(&header).expect("Failed to write header"); + + for (name, seq, flag_bits) in [ + ("pair1", "ACGTACGT", flags::PAIRED | flags::FIRST_SEGMENT), + ("pair1", "TGCATGCA", flags::PAIRED | flags::LAST_SEGMENT), + // No segment bits set: classify_segment routes this to "other". + ("single1", "GGGGCCCC", 0), + ] { + let mut b = SamBuilder::new(); + b.read_name(name.as_bytes()) + .sequence(seq.as_bytes()) + .qualities(&[40; 8]) + .flags(flag_bits) + .ref_id(0) + .pos(99) + .mapq(60); + writer + .write_alignment_record(&header, &to_record_buf(&b.build())) + .expect("Failed to write record"); + } + writer.try_finish().expect("Failed to finish BAM"); +} + +/// Parse FASTQ records from an in-memory byte buffer (e.g. captured stdout). +fn parse_fastq_bytes(bytes: &[u8]) -> Vec<(String, String, String)> { + let lines: Vec = + BufReader::new(bytes).lines().map(|l| l.expect("read stdout line")).collect(); + let mut records = Vec::new(); + for chunk in lines.chunks(4) { + if chunk.len() == 4 { + records.push(( + chunk[0].trim_start_matches('@').to_string(), + chunk[1].clone(), + chunk[3].clone(), + )); + } + } + records +} + +/// Run `fgumi` as a real subprocess with the given args and a watchdog +/// timeout, returning `(exited_successfully, stdout, stderr)`. +/// +/// Needed for cases that a plain in-process `Command::execute()` call can't +/// observe or guard: capturing real process stdout (the paired chain's +/// `--out0`-omitted branch writes straight to `io::stdout()` for path `-`), +/// and detecting a pipeline deadlock without hanging the test suite. +/// +/// stdout/stderr are drained on worker threads while the main thread polls +/// `try_wait`, mirroring `test_simulate_aligner.rs`'s pattern: a deadlocked +/// child never closes its pipes, so a synchronous `read_to_end` on the main +/// thread would block forever and the watchdog would never get to run. +fn run_fastq_subprocess_with_timeout(args: &[&str], timeout: Duration) -> (bool, Vec, Vec) { + let mut cmd = ProcessCommand::new(env!("CARGO_BIN_EXE_fgumi")); + cmd.args(args).stdin(Stdio::null()).stdout(Stdio::piped()).stderr(Stdio::piped()); + let mut child = cmd.spawn().expect("spawn fgumi fastq"); + + let mut child_stdout = child.stdout.take().expect("child stdout"); + let stdout_reader = thread::spawn(move || { + let mut out = Vec::new(); + child_stdout.read_to_end(&mut out).expect("read child stdout"); + out + }); + let mut child_stderr = child.stderr.take().expect("child stderr"); + let stderr_reader = thread::spawn(move || { + let mut out = Vec::new(); + child_stderr.read_to_end(&mut out).expect("read child stderr"); + out + }); + + let deadline = Instant::now() + timeout; + let success = loop { + if let Some(status) = child.try_wait().expect("try_wait") { + break status.success(); + } + if Instant::now() >= deadline { + let _ = child.kill(); + panic!("fgumi fastq did not exit within {timeout:?} (deadlock?)"); + } + thread::sleep(Duration::from_millis(20)); + }; + let stdout = stdout_reader.join().expect("stdout reader thread"); + let stderr = stderr_reader.join().expect("stderr reader thread"); + (success, stdout, stderr) +} + +/// `--out0` routes single-end / "other" reads to their own file, matching +/// `samtools fastq -0`; R1/R2 files contain only the paired reads. +#[test] +fn test_fastq_paired_other_reads_route_to_out0() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + create_paired_bam_with_other(&input_bam); + let r1 = temp_dir.path().join("R1.fq"); + let r2 = temp_dir.path().join("R2.fq"); + let r0 = temp_dir.path().join("R0.fq"); + + let cmd = Fastq::try_parse_from([ + "fastq", + "-i", + input_bam.to_str().unwrap(), + "-1", + r1.to_str().unwrap(), + "-2", + r2.to_str().unwrap(), + "-0", + r0.to_str().unwrap(), + ]) + .expect("parse"); + cmd.execute("fgumi fastq").expect("run"); + + let recs1 = read_fastq_maybe_gz(&r1); + let recs2 = read_fastq_maybe_gz(&r2); + let recs0 = read_fastq_maybe_gz(&r0); + assert_eq!(recs1.len(), 1, "R1 file has only the paired R1 read"); + assert_eq!(recs2.len(), 1, "R2 file has only the paired R2 read"); + assert_eq!(recs0.len(), 1, "out0 has the single-end read"); + assert_eq!(recs1[0].0, "pair1"); + assert_eq!(recs1[0].1, "ACGTACGT"); + assert_eq!(recs2[0].0, "pair1"); + assert_eq!(recs2[0].1, "TGCATGCA"); + assert_eq!(recs0[0].0, "single1"); + assert_eq!(recs0[0].1, "GGGGCCCC"); +} + +/// Both `-1`/`-2` outputs as `.fastq.gz` decompress to the expected R1/R2 +/// records and are each valid BGZF. Extends `test_fastq_paired_output_is_bgzf` +/// (which only checks R1) to both files. +#[test] +fn test_fastq_paired_gz_output_round_trips_both_files() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + create_pair_with_umi_tag(&input_bam, "read1", None); + let r1 = temp_dir.path().join("R1.fastq.gz"); + let r2 = temp_dir.path().join("R2.fastq.gz"); + + let cmd = Fastq::try_parse_from([ + "fastq", + "-i", + input_bam.to_str().unwrap(), + "-1", + r1.to_str().unwrap(), + "-2", + r2.to_str().unwrap(), + ]) + .expect("parse"); + cmd.execute("fgumi fastq").expect("run"); + + for path in [&r1, &r2] { + let raw = fs::read(path).expect("read gz"); + assert!(raw.len() > 18, "non-empty bgzf: {}", path.display()); + assert_eq!(&raw[0..2], &[0x1f, 0x8b], "gzip magic: {}", path.display()); + assert_eq!(raw[3] & 0x04, 0x04, "FEXTRA flag set: {}", path.display()); + assert_eq!([raw[12], raw[13]], [b'B', b'C'], "BGZF BC subfield: {}", path.display()); + } + + let recs1 = read_fastq_maybe_gz(&r1); + let recs2 = read_fastq_maybe_gz(&r2); + assert_eq!(recs1.len(), 1); + assert_eq!(recs2.len(), 1); + assert_eq!(recs1[0].0, "read1"); + assert_eq!(recs2[0].0, "read1"); + assert_eq!(recs1[0].1, "ACGTACGT"); + assert_eq!(recs2[0].1, "TGCATGCA"); +} + +/// Paired split output with `--annotate-read-names` across multiple pairs: +/// each pair's DRAGEN-format UMI annotation lands on both R1 and R2, and the +/// per-pair records line up positionally between the two files. +#[test] +fn test_fastq_paired_umi_header_multiple_pairs() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + let header = create_minimal_header("chr1", 10000); + let mut writer = bam::io::Writer::new(fs::File::create(&input_bam).expect("create BAM")); + writer.write_header(&header).expect("write header"); + + for (name, umi) in [("readA", "AAAA-CCCC"), ("readB", "GGGG-TTTT")] { + for (seq, segment_flag) in + [("ACGTACGT", flags::FIRST_SEGMENT), ("TGCATGCA", flags::LAST_SEGMENT)] + { + let mut b = SamBuilder::new(); + b.read_name(name.as_bytes()) + .sequence(seq.as_bytes()) + .qualities(&[40; 8]) + .flags(flags::PAIRED | segment_flag) + .ref_id(0) + .pos(99) + .mapq(60); + b.add_string_tag(SamTag::RX, umi.as_bytes()); + writer + .write_alignment_record(&header, &to_record_buf(&b.build())) + .expect("Failed to write record"); + } + } + writer.try_finish().expect("Failed to finish BAM"); + + let r1 = temp_dir.path().join("R1.fq"); + let r2 = temp_dir.path().join("R2.fq"); + let cmd = Fastq::try_parse_from([ + "fastq", + "-i", + input_bam.to_str().unwrap(), + "--annotate-read-names", + "-1", + r1.to_str().unwrap(), + "-2", + r2.to_str().unwrap(), + ]) + .expect("parse"); + cmd.execute("fgumi fastq").expect("run"); + + let recs1 = read_fastq_maybe_gz(&r1); + let recs2 = read_fastq_maybe_gz(&r2); + assert_eq!(recs1.len(), 2); + assert_eq!(recs2.len(), 2); + assert_eq!(recs1[0].0, "readA:AAAA+CCCC"); + assert_eq!(recs2[0].0, "readA:AAAA+CCCC"); + assert_eq!(recs1[1].0, "readB:GGGG+TTTT"); + assert_eq!(recs2[1].0, "readB:GGGG+TTTT"); +} + +/// Without `--out0`, single-end/"other" reads go to stdout, matching +/// `samtools fastq` without `-0`; R1/R2 files contain only the paired reads. +/// Runs the real binary as a subprocess: the paired chain writes the "other" +/// branch straight to `io::stdout()` for path `-`, which an in-process +/// `Command::execute()` call in this test binary cannot observe. +#[test] +fn test_fastq_paired_no_out0_sends_other_to_stdout() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + create_paired_bam_with_other(&input_bam); + let r1 = temp_dir.path().join("R1.fq"); + let r2 = temp_dir.path().join("R2.fq"); + + let (success, stdout, stderr) = run_fastq_subprocess_with_timeout( + &[ + "fastq", + "-i", + input_bam.to_str().unwrap(), + "-1", + r1.to_str().unwrap(), + "-2", + r2.to_str().unwrap(), + ], + Duration::from_secs(30), + ); + assert!(success, "fgumi fastq failed: {}", String::from_utf8_lossy(&stderr)); + + let stdout_recs = parse_fastq_bytes(&stdout); + assert_eq!(stdout_recs.len(), 1, "single-end read on stdout"); + assert_eq!(stdout_recs[0].0, "single1"); + assert_eq!(stdout_recs[0].1, "GGGGCCCC"); + + let recs1 = read_fastq_maybe_gz(&r1); + let recs2 = read_fastq_maybe_gz(&r2); + assert_eq!(recs1.len(), 1, "R1 file has only the paired read"); + assert_eq!(recs2.len(), 1, "R2 file has only the paired read"); +} + +/// A BAM whose reads are ALL R1 (no R2 at all) must not deadlock. The paired +/// encode step always emits a block on every branch per batch — an empty +/// `Vec` where a branch has no records — precisely so the R2 branch's +/// `batch_serial` sequence has no gaps; a gap would wedge the R2 writer's +/// reorder stage waiting on an ordinal that never arrives. This is the direct +/// regression guard for that rule: it asserts the command actually completes +/// (via the subprocess watchdog) rather than hanging, and that the R2 output +/// is a valid empty file rather than missing or truncated. +#[test] +fn test_fastq_paired_empty_r2_batch_no_hang() { + let temp_dir = TempDir::new().expect("temp dir"); + let input_bam = temp_dir.path().join("input.bam"); + let header = create_minimal_header("chr1", 10000); + let mut writer = bam::io::Writer::new(fs::File::create(&input_bam).expect("create BAM")); + writer.write_header(&header).expect("write header"); + for i in 0..8 { + let name = format!("read{i}"); + let mut b = SamBuilder::new(); + b.read_name(name.as_bytes()) + .sequence(b"ACGTACGT") + .qualities(&[40; 8]) + .flags(flags::PAIRED | flags::FIRST_SEGMENT) + .ref_id(0) + .pos(99) + .mapq(60); + writer + .write_alignment_record(&header, &to_record_buf(&b.build())) + .expect("Failed to write record"); + } + writer.try_finish().expect("Failed to finish BAM"); + + let r1 = temp_dir.path().join("R1.fq"); + let r2 = temp_dir.path().join("R2.fq"); + let r0 = temp_dir.path().join("R0.fq"); + + let (success, _stdout, stderr) = run_fastq_subprocess_with_timeout( + &[ + "fastq", + "-i", + input_bam.to_str().unwrap(), + "-1", + r1.to_str().unwrap(), + "-2", + r2.to_str().unwrap(), + "-0", + r0.to_str().unwrap(), + ], + Duration::from_secs(30), + ); + assert!(success, "fgumi fastq did not complete: {}", String::from_utf8_lossy(&stderr)); + + assert_eq!(read_fastq_maybe_gz(&r1).len(), 8, "all reads are R1"); + let r2_len = fs::metadata(&r2).expect("stat R2").len(); + assert_eq!(r2_len, 0, "R2 file must be empty (zero-length), not missing/corrupt"); + assert_eq!(read_fastq_maybe_gz(&r0).len(), 0, "no other reads"); +}