From 84bee5897123156dd569aa689ebb5b39e7be8b30 Mon Sep 17 00:00:00 2001 From: alejandrogzi Date: Thu, 6 Aug 2026 15:16:26 +0200 Subject: [PATCH 1/2] BREAKING CHANGE: v0.0.16 -> CBQ input/output implementation --- CHANGELOG.md | 3 + Cargo.lock | 88 ++++- Cargo.toml | 5 +- README.md | 6 +- deacon-py/Cargo.toml | 2 +- deacon-py/src/lib.rs | 2 +- src/filter.rs | 781 +++++++++++++++++++++++++++++++----------- src/index.rs | 8 +- src/main.rs | 8 +- tests/filter_tests.rs | 273 +++++++++++++++ 10 files changed, 957 insertions(+), 219 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3ffbda4..a774c0a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,9 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.16.0] - 2026-08-05 + ### Added - Support for filtering using smaller static fuse filter indexes (`.pidx`). These can be generated from standard Deacon indexes using `deacon index freeze`. Fuse filters reduce index size from ~64 bits/k-mer to ~36 bits/k-mer using default 32 bit fingerprints, or to ~18 bits/k-mer using 16 bit fingerprints. Default 32 bit fingerprints deliver a negligible false positive rate of ~2-32. +- `deacon filter` accepts CBQ (BINSEQ) input, detected automatically, and writes CBQ when the output path ends in `.cbq`; paired reads are stored as native paired records in a single file. - Browser-based streaming FASTA/FASTQ filtering with `deacon-wasm`. - `deacon index filter` keeps or discards index minimizers based on a complexity threshold (`-c`), with a choice of `--algorithm` (`-a`). - `deacon filter --complexity-threshold` (`-c`) ignores minimizers below a kdust complexity threshold at filter time. diff --git a/Cargo.lock b/Cargo.lock index 32eeadf..86810f3 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -94,6 +94,17 @@ dependencies = [ "wait-timeout", ] +[[package]] +name = "auto_impl" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ffdcb70bdbc4d478427380519163274ac86e52916e10f0a8889adf0f96d3fee7" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + [[package]] name = "autocfg" version = "1.5.1" @@ -155,6 +166,27 @@ dependencies = [ "virtue", ] +[[package]] +name = "binseq" +version = "0.9.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c03ce47e3c043402933b71d00e26586b91f70d1cc237251097e2eb1ee9621d83" +dependencies = [ + "anyhow", + "auto_impl", + "bitnuc", + "bytemuck", + "byteorder", + "itoa", + "memchr", + "memmap2", + "num_cpus", + "rand 0.9.5", + "sucds", + "thiserror 2.0.18", + "zstd", +] + [[package]] name = "bitflags" version = "1.3.2" @@ -167,6 +199,12 @@ version = "2.13.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b4388bee8683e3d04af747c73422af53102d2bd24d9eadb6cbc100baef4b43f8" +[[package]] +name = "bitnuc" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d7871516cbb4e097e220623917e1491c995ee58c8018d929d4b1e76c378002bf" + [[package]] name = "bstr" version = "1.12.3" @@ -186,9 +224,23 @@ checksum = "72f5acc6cb2ba439de613abc23857ec3d78374d8ed5ac84e9d11336e87da8649" [[package]] name = "bytemuck" -version = "1.25.0" +version = "1.25.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c8efb64bd706a16a1bdde310ae86b351e4d21550d98d056f22f8a7f7a2183fec" +checksum = "95832e849adfb21180ccb6826a99da14e5d266ae5c2e668e1602cf234f153797" +dependencies = [ + "bytemuck_derive", +] + +[[package]] +name = "bytemuck_derive" +version = "1.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f65693059b6b9c588b9f62fed1cedbf0a8b805631457ea162d68f0de186f3de5" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] [[package]] name = "byteorder" @@ -425,6 +477,7 @@ dependencies = [ "anyhow", "assert_cmd", "bincode", + "binseq", "clap", "ensure_simd", "flate2", @@ -438,7 +491,7 @@ dependencies = [ "paraseq", "parking_lot", "predicates", - "rand 0.9.4", + "rand 0.9.5", "rapidhash", "rayon", "rstest", @@ -453,7 +506,7 @@ dependencies = [ [[package]] name = "deacon-py" -version = "0.16.0-rc4" +version = "0.16.0" dependencies = [ "anyhow", "deacon", @@ -1009,9 +1062,18 @@ dependencies = [ [[package]] name = "memchr" -version = "2.8.2" +version = "2.8.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "88904434abc2901f197fe8cc55f0445e7ded921dba5911dad2e2b39b48e663c4" +checksum = "cf8baf1c55e62ffcace7a9f06f4bd9cd3f0c4beb022d3b367256b91b87513d98" + +[[package]] +name = "memmap2" +version = "0.9.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d1219ed1b7f229ee7104d281dd01d6802fe28bb6e95d292942c4daacdeb798c0" +dependencies = [ + "libc", +] [[package]] name = "miniz_oxide" @@ -1382,9 +1444,9 @@ checksum = "f8dcc9c7d52a811697d2151c701e0d08956f92b0e24136cf4cf27b57a6a0d9bf" [[package]] name = "rand" -version = "0.9.4" +version = "0.9.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "44c5af06bb1b7d3216d91932aed5265164bf384dc89cd6ba05cf59a35f5f76ea" +checksum = "b9ef1d0d795eb7d84685bca4f72f3649f064e6641543d3a8c415898726a57b41" dependencies = [ "rand_chacha", "rand_core 0.9.5", @@ -1835,6 +1897,16 @@ version = "2.6.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "13c2bddecc57b384dee18652358fb23172facb8a2c51ccc10d74c157bdea3292" +[[package]] +name = "sucds" +version = "0.8.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cd324eaa05be64f105ea5269bb8aabd70e5dd57fa5c673b167f451b07d6c0dcd" +dependencies = [ + "anyhow", + "num-traits", +] + [[package]] name = "syn" version = "2.0.118" diff --git a/Cargo.toml b/Cargo.toml index b7e3125..6edb685 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -42,6 +42,8 @@ rand = { version = "0.9", optional = true } serde_json = { version = "1.0", optional = true } niffler = { version = "3.0.1", default-features = false, optional = true } paraseq = { version = "0.4.3", default-features = false, features = ["anyhow", "niffler"], optional = true } +# Binary sequence format (CBQ) support. Defaults disabled: binseq's own paraseq integration is not needed. +binseq = { version = "0.9.4", default-features = false, features = ["anyhow"], optional = true } ensure_simd = { version = "0.1.0", optional = true } # Compression deps (optional) @@ -60,10 +62,11 @@ predicates = "3.0" tempfile = "3.20" rstest = "0.26" nix = { version = "0.31", features = ["fs"] } +binseq = { version = "0.9.4", default-features = false, features = ["anyhow"] } [features] scalar = ["simd-minimizers/scalar", "packed-seq/scalar"] -cli = ["rayon", "parking_lot", "indicatif", "paraseq", "ensure_simd", "clap", "rand", "niffler", "serde_json"] +cli = ["rayon", "parking_lot", "indicatif", "paraseq", "binseq", "ensure_simd", "clap", "rand", "niffler", "serde_json"] # Disable for faster builds. compression = ["zstd", "liblzma", "flate2", "gzp", "paraseq/default"] # Use to still handle .gz files when "compression" is not enabled. diff --git a/README.md b/README.md index 0f618b0..51bdf42 100644 --- a/README.md +++ b/README.md @@ -89,7 +89,7 @@ Prebuilt pangenome indexes are provided. These can be downloaded using the links ### Filtering -The main command `deacon filter` accepts an index path followed by up to two FASTA/FASTQ file paths, depending on whether input sequences originate from stdin, a single file, or paired input files. Indexes are built with `deacon index build`. Paired inputs are supported as either two separate or one interleaved file/stream when using `--interleaved`, and may be written either to separate paired output files or one interleaved file. For paired sequences, distinct minimizer hits originating from either mate are counted. By default, input sequences must meet both an absolute threshold of 2 minimizer hits (`-a 2`) and a relative threshold of 1% of minimizers (`-r 0.01`) to pass the filter. Filtering can be inverted for e.g. host depletion using the `--deplete` (`-d`) flag. Gzip, Zstandard, and xz compression formats are detected automatically by file extension. Paired read headers can be validated using `--check-pairs`. +The main command `deacon filter` accepts an index path followed by up to two FASTA/FASTQ file paths, depending on whether input sequences originate from stdin, a single file, or paired input files. Indexes are built with `deacon index build`. Paired inputs are supported as either two separate or one interleaved file/stream when using `--interleaved`, and may be written either to separate paired output files or one interleaved file. For paired sequences, distinct minimizer hits originating from either mate are counted. By default, input sequences must meet both an absolute threshold of 2 minimizer hits (`-a 2`) and a relative threshold of 1% of minimizers (`-r 0.01`) to pass the filter. Filtering can be inverted for e.g. host depletion using the `--deplete` (`-d`) flag. Gzip, Zstandard, and xz compression formats are detected automatically by file extension. Paired read headers can be validated using `--check-pairs`. CBQ files (the columnar binary format of the [BINSEQ](https://www.biorxiv.org/content/10.1101/2025.04.08.647863v2) family) are detected on input by magic bytes, and written when the output path ends in `.cbq`; CBQ pairing is native, so paired reads stay in one file (do not pass `--output2`). #### Examples @@ -131,6 +131,10 @@ zcat r12.fq.gz | deacon filter -d panhuman-1.k31w15.idx - - > filt12.fq # Save summary JSON deacon filter -d panhuman-1.k31w15.idx reads.fq.gz -o filt.fq.gz -s summary.json +# CBQ (BINSEQ) input/output; paired reads stay in one CBQ file +deacon filter panhuman-1.k31w15.idx reads.fq.gz -o filt.cbq +deacon filter panhuman-1.k31w15.idx filt.cbq -o filt.fq.gz + # Replace read headers with incrementing integers deacon filter -d -R panhuman-1.k31w15.idx reads.fq.gz > filt.fq diff --git a/deacon-py/Cargo.toml b/deacon-py/Cargo.toml index 6b94e07..f2f0193 100644 --- a/deacon-py/Cargo.toml +++ b/deacon-py/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "deacon-py" -version = "0.16.0-rc4" +version = "0.16.0" edition = "2024" license = "MIT" publish = false diff --git a/deacon-py/src/lib.rs b/deacon-py/src/lib.rs index 58b21d0..85063ef 100644 --- a/deacon-py/src/lib.rs +++ b/deacon-py/src/lib.rs @@ -157,7 +157,7 @@ impl Index { let mins = Arc::clone(&self.minimizers); let (k, w) = (self.k, self.w); let summary = py - .detach(|| run_with_index(&mins, &IndexHeader::new(k, w), &cfg)) + .detach(|| run_with_index(mins, &IndexHeader::new(k, w), &cfg)) .map_err(to_pyerr)?; Ok(pythonize::pythonize(py, &summary)?.unbind()) } diff --git a/src/filter.rs b/src/filter.rs index eb60195..fee5099 100644 --- a/src/filter.rs +++ b/src/filter.rs @@ -4,6 +4,9 @@ use crate::{ MinimizerSet, }; use anyhow::{Context, Result}; +use binseq::cbq; +use binseq::write::{BinseqWriterBuilder, Format as BinseqFormat}; +use binseq::{BinseqRecord, ParallelReader as BinseqParallelReader, SequencingRecordBuilder}; use indicatif::{ProgressBar, ProgressDrawTarget, ProgressStyle}; use paraseq::Record; use paraseq::fastx::Reader; @@ -11,8 +14,9 @@ use paraseq::parallel::{PairedParallelProcessor, ParallelProcessor, ParallelRead use parking_lot::Mutex; use rand::Rng; use serde::{Deserialize, Serialize}; +use std::borrow::Cow; use std::fs::{File, OpenOptions}; -use std::io::{self, BufWriter, Write}; +use std::io::{self, BufWriter, Read, Write}; use std::path::PathBuf; use std::sync::Arc; use std::sync::atomic::{AtomicU64, Ordering}; @@ -23,6 +27,109 @@ const DEFAULT_BUFFER_SIZE: usize = 64 * 1024; type BoxedWriter = Box; +/// Input file format +#[derive(Clone, Copy, PartialEq, Eq)] +enum InputFormat { + Fastx, + Cbq, +} + +/// Output file format +#[derive(Clone, Copy, PartialEq, Eq)] +enum OutputFormat { + Fastx, + Cbq, +} + +/// Normalized input metadata, populated before any output is opened +struct InputLayout { + format: InputFormat, + paired: bool, + qualities: bool, + headers: bool, +} + +/// Shared output state; every processor clone merges into this under a lock +#[derive(Clone)] +enum SharedOutput { + Fastx(Arc>), + Cbq(Arc>>), +} + +/// Thread-local output buffer merged into the shared output on batch completion +#[allow(clippy::large_enum_variant)] // mirrors binseq's BinseqWriter, which is not boxed upstream +#[derive(Clone)] +enum LocalOutput { + Fastx(Vec), + Cbq(binseq::BinseqWriter>), +} + +impl SharedOutput { + fn is_cbq(&self) -> bool { + matches!(self, Self::Cbq(_)) + } +} + +impl LocalOutput { + fn fastx_mut(&mut self) -> &mut Vec { + match self { + Self::Fastx(v) => v, + Self::Cbq(_) => unreachable!("CBQ local buffer used in FASTX path"), + } + } + + fn fastx(&self) -> &[u8] { + match self { + Self::Fastx(v) => v, + Self::Cbq(_) => unreachable!("CBQ local buffer used in FASTX path"), + } + } + + fn cbq_mut(&mut self) -> &mut binseq::BinseqWriter> { + match self { + Self::Cbq(w) => w, + Self::Fastx(_) => unreachable!("FASTX local buffer used in CBQ path"), + } + } +} + +/// Input processing plan: every reader opens during layout resolution, before +/// any output file is created. +#[allow(clippy::large_enum_variant)] // MmapReader holds its mmap; boxed variants would add noise for no measurable gain +enum InputPrep { + Cbq(cbq::MmapReader), + FastxSingle(Reader>), + FastxInterleaved(Reader>), + FastxPaired(Reader>, Reader>), + /// Empty input: create empty output without processing + Empty, +} + +/// Adapter exposing the primary/secondary halves of a CBQ record as paraseq records +struct CbqRecord<'a> { + id: &'a [u8], + seq: &'a [u8], + qual: Option<&'a [u8]>, +} + +impl Record for CbqRecord<'_> { + fn id(&self) -> &[u8] { + self.id + } + + fn seq(&self) -> Cow<'_, [u8]> { + Cow::Borrowed(self.seq) + } + + fn seq_raw(&self) -> &[u8] { + self.seq + } + + fn qual(&self) -> Option<&[u8]> { + self.qual + } +} + /// Filtering config for an already-loaded index (no index path; see [`FilterConfig`]). pub struct FilterRunConfig { /// Path to input fastx file (or - for stdin) @@ -197,6 +304,229 @@ fn create_paraseq_reader(path: Option<&str>) -> Result Result { + let Some(path) = config.output_path.as_deref() else { + return Ok(OutputFormat::Fastx); + }; + let name = path.to_string_lossy(); + if name.ends_with(".cbq") { + Ok(OutputFormat::Cbq) + } else if name.ends_with(".cbq.gz") + || name.ends_with(".cbq.zst") + || name.ends_with(".cbq.xz") + { + anyhow::bail!( + "Compressed CBQ output is not supported (CBQ performs its own block compression): {}", + name + ); + } else { + Ok(OutputFormat::Fastx) + } +} + +/// Resolve the input format and layout, opening all readers up front so input +/// errors surface before any output file is created. +fn prepare_input( + config: &FilterRunConfig, + interleaved_input: bool, + output_format: OutputFormat, +) -> Result<(InputLayout, InputPrep)> { + // ponytail: CBQ input is file-only so the mmap reader stays the only CBQ + // path; add binseq's streaming reader if stdin support is ever needed. + if config.input_path == "-" || is_special_input_path(&config.input_path) { + return prepare_fastx(config, interleaved_input, output_format); + } + + // Regular files: sniff the CBQ magic; everything else is FASTX + let mut file = File::open(&config.input_path) + .map_err(|e| anyhow::anyhow!("Failed to open file {}: {}", config.input_path, e))?; + let mut magic = [0u8; 64]; + let n = file.read(&mut magic)?; + match BinseqFormat::sniff(&magic[..n]) { + Some(BinseqFormat::Cbq) => { + match cbq::MmapReader::new(&config.input_path) { + Ok(reader) => { + let header = reader.header(); + let layout = InputLayout { + format: InputFormat::Cbq, + paired: header.is_paired(), + qualities: header.has_qualities(), + headers: header.has_headers(), + }; + Ok((layout, InputPrep::Cbq(reader))) + } + // binseq 0.9.4 cannot reopen zero-record CBQ files (its index + // cast fails on zero entries); treat them as valid empty inputs. + Err(binseq::Error::CbqError(binseq::error::CbqError::IndexCastingError)) + if n >= std::mem::size_of::() => + { + let header = cbq::FileHeader::from_bytes( + &magic[..std::mem::size_of::()], + )?; + Ok(( + InputLayout { + format: InputFormat::Cbq, + paired: header.is_paired(), + qualities: header.has_qualities(), + headers: header.has_headers(), + }, + InputPrep::Empty, + )) + } + Err(e) => Err(e).context("Failed to open CBQ input"), + } + } + Some(f) => anyhow::bail!("{f:?} input is not supported; convert it to FASTQ or CBQ"), + None => prepare_fastx(config, interleaved_input, output_format), + } +} + +/// Prepare a FASTX input: resolve the layout and open the reader(s) up front. +fn prepare_fastx( + config: &FilterRunConfig, + interleaved_input: bool, + output_format: OutputFormat, +) -> Result<(InputLayout, InputPrep)> { + let layout = InputLayout { + format: InputFormat::Fastx, + paired: interleaved_input || config.input2_path.is_some(), + qualities: false, + headers: true, + }; + + let input1_empty = is_empty_file(&config.input_path)?; + let input2_empty = config + .input2_path + .as_deref() + .map(is_empty_file) + .transpose()? + .unwrap_or(false); + + if interleaved_input { + if input1_empty { + return Ok((layout, InputPrep::Empty)); + } + return match create_paraseq_reader(Some(config.input_path.as_str())) { + Ok(reader) => { + let qualities = reader.format() == paraseq::fastx::Format::Fastq; + Ok(( + InputLayout { qualities, ..layout }, + InputPrep::FastxInterleaved(reader), + )) + } + Err(e) if is_empty_input_error(&e) => Ok((layout, InputPrep::Empty)), + Err(e) => Err(e), + }; + } + + if let Some(input2_path) = config.input2_path.as_deref() { + if input1_empty && input2_empty { + return Ok((layout, InputPrep::Empty)); + } + if input1_empty || input2_empty { + return Err(anyhow::anyhow!( + "One paired file is empty but the other is not" + )); + } + let r1 = create_paraseq_reader(Some(config.input_path.as_str())); + let r2 = create_paraseq_reader(Some(input2_path)); + return match (r1, r2) { + (Ok(reader1), Ok(reader2)) => { + if output_format == OutputFormat::Cbq + && !config.output_fasta + && reader1.format() != reader2.format() + { + anyhow::bail!( + "Mixed FASTA/FASTQ paired input cannot be written to a quality CBQ output" + ); + } + let qualities = reader1.format() == paraseq::fastx::Format::Fastq; + Ok(( + InputLayout { qualities, ..layout }, + InputPrep::FastxPaired(reader1, reader2), + )) + } + (Err(e1), Err(e2)) if is_empty_input_error(&e1) && is_empty_input_error(&e2) => { + Ok((layout, InputPrep::Empty)) + } + (Err(e), _) if is_empty_input_error(&e) => Err(anyhow::anyhow!( + "First paired file appears empty while second is not" + )), + (_, Err(e)) if is_empty_input_error(&e) => Err(anyhow::anyhow!( + "Second paired file appears empty while first is not" + )), + (Err(e), _) => Err(e), + (_, Err(e)) => Err(e), + }; + } + + if input1_empty { + return Ok((layout, InputPrep::Empty)); + } + match create_paraseq_reader(Some(config.input_path.as_str())) { + Ok(reader) => { + let qualities = reader.format() == paraseq::fastx::Format::Fastq; + Ok(( + InputLayout { qualities, ..layout }, + InputPrep::FastxSingle(reader), + )) + } + Err(e) if is_empty_input_error(&e) => Ok((layout, InputPrep::Empty)), + Err(e) => Err(e), + } +} + +/// Validate format combinations before any output file is opened or truncated +fn validate_input_output( + layout: &InputLayout, + output_format: OutputFormat, + config: &FilterRunConfig, +) -> Result<()> { + if layout.format == InputFormat::Cbq && config.input2_path.is_some() { + anyhow::bail!("CBQ input does not support INPUT2"); + } + if layout.format == InputFormat::Cbq && config.interleaved { + anyhow::bail!("CBQ input does not support --interleaved"); + } + if output_format == OutputFormat::Cbq && config.output2_path.is_some() { + anyhow::bail!("CBQ output does not support OUTPUT2; CBQ pairing is native"); + } + if output_format == OutputFormat::Cbq && config.output_path.is_none() { + anyhow::bail!("CBQ output to stdout is not supported"); + } + if config.check_pairs && layout.format == InputFormat::Cbq && !layout.headers { + anyhow::bail!("--check-pairs requires CBQ input with headers"); + } + validate_check_pairs_mode(config.check_pairs, layout.paired)?; + Ok(()) +} + +/// Renamed or original header for a record. Renamed headers are +/// `counter[+ -random][+ suffix]`, matching the existing FASTX output. +fn record_header<'a>( + id: &'a [u8], + counter: u64, + rename: bool, + rename_random: bool, + read_suffix: &[u8], +) -> Cow<'a, [u8]> { + if rename || rename_random { + let mut header = counter.to_string().into_bytes(); + if rename_random { + header.push(b'-'); + header.extend_from_slice(rand::rng().random::().to_string().as_bytes()); + } + if !read_suffix.is_empty() { + header.push(b' '); + header.extend_from_slice(read_suffix); + } + Cow::Owned(header) + } else { + Cow::Borrowed(id) + } +} + /// Format a single record into a buffer (FASTA/FASTQ format) /// `seq` is the newline-stripped sequence corresponding to the record from `record.seq()`. fn format_record_to_buffer( @@ -213,20 +543,8 @@ fn format_record_to_buffer( // Header buffer.write_all(if is_fasta { b">" } else { b"@" })?; - if rename || rename_random { - buffer.extend_from_slice(counter.to_string().as_bytes()); - if rename_random { - buffer.write_all(b"-")?; - let random_suffix = rand::rng().random::(); - buffer.extend_from_slice(random_suffix.to_string().as_bytes()); - } - if !read_suffix.is_empty() { - buffer.write_all(b" ")?; - buffer.write_all(read_suffix)?; - } - } else { - buffer.extend_from_slice(record.id()); - } + let header = record_header(record.id(), counter, rename, rename_random, read_suffix); + buffer.write_all(&header)?; buffer.write_all(b"\n")?; // Sequence @@ -381,9 +699,9 @@ pub struct FilterSummary { } #[derive(Clone)] -struct FilterProcessor<'a> { +struct FilterProcessor { // Minimizer matching parameters - minimizers: &'a MinimizerSet, + minimizers: Arc, rename: bool, rename_random: bool, output_fasta: bool, @@ -392,15 +710,15 @@ struct FilterProcessor<'a> { kernel: FilterKernel, // Local buffers - local_buffer: Vec, - local_buffer2: Vec, // Second buffer for paired output + local_buffer: LocalOutput, + local_buffer2: Vec, // Second buffer for paired FASTX output local_stats: ProcessingStats, // Shared atomic counter for seq renaming rename_counter: Arc, // Global state - global_writer: Arc>, + global_writer: SharedOutput, global_writer2: Option>>, global_stats: Arc>, spinner: Option>>, @@ -417,18 +735,22 @@ pub(crate) struct ProcessingStats { pub last_reported: u64, } -impl<'a> FilterProcessor<'a> { +impl FilterProcessor { fn new( - minimizers: &'a MinimizerSet, + minimizers: Arc, kmer_length: u8, window_size: u8, config: &FilterProcessorConfig, - writer: BoxedWriter, + global_writer: SharedOutput, writer2: Option, spinner: Option>>, filtering_start_time: Instant, - ) -> Self { - Self { + ) -> Result { + let local_buffer = match &global_writer { + SharedOutput::Cbq(writer) => LocalOutput::Cbq(writer.lock().new_headless_buffer()?), + SharedOutput::Fastx(_) => LocalOutput::Fastx(Vec::with_capacity(DEFAULT_BUFFER_SIZE)), + }; + Ok(Self { minimizers, rename: config.rename, rename_random: config.rename_random, @@ -445,25 +767,25 @@ impl<'a> FilterProcessor<'a> { prefix_length: config.prefix_length, }, ), - local_buffer: Vec::with_capacity(DEFAULT_BUFFER_SIZE), + local_buffer, local_buffer2: Vec::with_capacity(DEFAULT_BUFFER_SIZE), - rename_counter: Arc::new(AtomicU64::new(0)), local_stats: ProcessingStats::default(), - global_writer: Arc::new(Mutex::new(writer)), + rename_counter: Arc::new(AtomicU64::new(0)), + global_writer, global_writer2: writer2.map(|w| Arc::new(Mutex::new(w))), global_stats: Arc::new(Mutex::new(ProcessingStats::default())), spinner, filtering_start_time, - } + }) } fn should_keep_sequence(&mut self, seq: &[u8]) -> FilterDecision { - self.kernel.classify_read(self.minimizers, seq, self.debug) + self.kernel .classify_read(&self.minimizers, seq, self.debug) } fn should_keep_pair(&mut self, seq1: &[u8], seq2: &[u8]) -> FilterDecision { self.kernel - .classify_pair(self.minimizers, seq1, seq2, self.debug) + .classify_pair(&self.minimizers, seq1, seq2, self.debug) } fn write_record( @@ -473,6 +795,9 @@ impl<'a> FilterProcessor<'a> { counter: u64, read_suffix: &[u8], ) -> Result<()> { + if self.global_writer.is_cbq() { + return self.push_cbq(record, seq, counter, read_suffix); + } format_record_to_buffer( record, seq, @@ -481,7 +806,7 @@ impl<'a> FilterProcessor<'a> { self.rename_random, read_suffix, self.output_fasta, - &mut self.local_buffer, + self.local_buffer.fastx_mut(), ) } @@ -504,6 +829,47 @@ impl<'a> FilterProcessor<'a> { ) } + /// Push a single record into the thread-local CBQ writer + fn push_cbq( + &mut self, + record: &Rf, + seq: &[u8], + counter: u64, + read_suffix: &[u8], + ) -> Result<()> { + let header = record_header(record.id(), counter, self.rename, self.rename_random, read_suffix); + let seq_record = SequencingRecordBuilder::default() + .s_seq(seq) + .s_header(&header) + .opt_s_qual(record.qual()) + .build()?; + self.local_buffer.cbq_mut().push(seq_record)?; + Ok(()) + } + + /// Push a paired record (both mates) into the thread-local CBQ writer + fn push_cbq_pair( + &mut self, + record1: &Rf, + seq1: &[u8], + record2: &Rf, + seq2: &[u8], + counter: u64, + ) -> Result<()> { + let header1 = record_header(record1.id(), counter, self.rename, self.rename_random, b"/1"); + let header2 = record_header(record2.id(), counter, self.rename, self.rename_random, b"/2"); + let seq_record = SequencingRecordBuilder::default() + .s_seq(seq1) + .s_header(&header1) + .opt_s_qual(record1.qual()) + .x_seq(seq2) + .x_header(&header2) + .opt_x_qual(record2.qual()) + .build()?; + self.local_buffer.cbq_mut().push(seq_record)?; + Ok(()) + } + fn update_spinner(&self) { if let Some(ref spinner) = self.spinner { let stats = self.global_stats.lock(); @@ -538,10 +904,9 @@ impl<'a> FilterProcessor<'a> { )); } } -} -impl<'a, Rf: Record> ParallelProcessor for FilterProcessor<'a> { - fn process_record(&mut self, record: Rf) -> paraseq::parallel::Result<()> { + /// Shared per-record logic for every reader (paraseq FASTX, CBQ mmap, ...) + fn handle_record(&mut self, record: &Rf) -> Result<()> { let seq = record.seq(); self.local_stats.total_seqs += 1; self.local_stats.total_bp += seq.len() as u64; @@ -567,7 +932,7 @@ impl<'a, Rf: Record> ParallelProcessor for FilterProcessor<'a> { } else { 0 }; - self.write_record(&record, &seq, counter, b"")?; + self.write_record(record, &seq, counter, b"")?; } else { self.local_stats.filtered_seqs += 1; self.local_stats.filtered_bp += seq.len() as u64; @@ -576,46 +941,14 @@ impl<'a, Rf: Record> ParallelProcessor for FilterProcessor<'a> { Ok(()) } - fn on_batch_complete(&mut self) -> paraseq::parallel::Result<()> { - // Write buffer to output - if !self.local_buffer.is_empty() { - let mut global_writer = self.global_writer.lock(); - global_writer.write_all(&self.local_buffer)?; - global_writer.flush()?; - } - - // Clear buffer after releasing the lock - self.local_buffer.clear(); - - // Update global stats - { - let mut stats = self.global_stats.lock(); - stats.total_seqs += self.local_stats.total_seqs; - stats.filtered_seqs += self.local_stats.filtered_seqs; - stats.total_bp += self.local_stats.total_bp; - stats.output_bp += self.local_stats.output_bp; - stats.filtered_bp += self.local_stats.filtered_bp; - } - - // Update spinner - self.update_spinner(); - - // Reset local stats - self.local_stats = ProcessingStats::default(); - - Ok(()) - } -} - -impl<'a, Rf: Record> PairedParallelProcessor for FilterProcessor<'a> { - fn process_record_pair(&mut self, record1: Rf, record2: Rf) -> paraseq::parallel::Result<()> { + /// Shared per-pair logic for every reader (paraseq FASTX, CBQ mmap, ...) + fn handle_record_pair(&mut self, record1: &Rf, record2: &Rf) -> Result<()> { if self.check_pairs && !paired_record_names_match(record1.id(), record2.id()) { return Err(anyhow::anyhow!( "Paired record name mismatch: R1='{}', R2='{}'. Expected matching Illumina CASAVA 1: and 2: fields or names suffixed with /1 and /2", String::from_utf8_lossy(record1.id()), String::from_utf8_lossy(record2.id()) - ) - .into()); + )); } let seq1 = record1.seq(); @@ -647,15 +980,17 @@ impl<'a, Rf: Record> PairedParallelProcessor for FilterProcessor<'a> { 0 }; - // Write to appropriate writers - if self.global_writer2.is_some() { + if self.global_writer.is_cbq() { + // Native paired record: both mates in one CBQ record + self.push_cbq_pair(record1, &seq1, record2, &seq2, counter)?; + } else if self.global_writer2.is_some() { // Separate outputs - self.write_record(&record1, &seq1, counter, b"/1")?; - self.write_record_to_buffer2(&record2, &seq2, counter, b"/2")?; + self.write_record(record1, &seq1, counter, b"/1")?; + self.write_record_to_buffer2(record2, &seq2, counter, b"/2")?; } else { // Interleaved output - self.write_record(&record1, &seq1, counter, b"/1")?; - self.write_record(&record2, &seq2, counter, b"/2")?; + self.write_record(record1, &seq1, counter, b"/1")?; + self.write_record(record2, &seq2, counter, b"/2")?; } } else { self.local_stats.filtered_seqs += 2; @@ -665,29 +1000,38 @@ impl<'a, Rf: Record> PairedParallelProcessor for FilterProcessor<'a> { Ok(()) } - fn on_batch_complete(&mut self) -> paraseq::parallel::Result<()> { - if let Some(ref writer2) = self.global_writer2 { - // Atomic paired batch writing - if !self.local_buffer.is_empty() || !self.local_buffer2.is_empty() { - let mut writer1 = self.global_writer.lock(); - let mut writer2 = writer2.lock(); - - writer1.write_all(&self.local_buffer)?; - writer1.flush()?; - writer2.write_all(&self.local_buffer2)?; - writer2.flush()?; + /// Merge thread-local buffers and stats into the shared state (per batch) + fn flush_batch(&mut self) -> Result<()> { + match &self.global_writer { + SharedOutput::Cbq(global) => { + global.lock().ingest_completed(self.local_buffer.cbq_mut())?; } - } else { - // Interleaved output - if !self.local_buffer.is_empty() { - let mut writer = self.global_writer.lock(); - writer.write_all(&self.local_buffer)?; - writer.flush()?; + SharedOutput::Fastx(global) => { + if let Some(writer2) = &self.global_writer2 { + // Atomic paired batch writing + if !self.local_buffer.fastx().is_empty() || !self.local_buffer2.is_empty() { + let mut writer1 = global.lock(); + let mut writer2 = writer2.lock(); + + writer1.write_all(self.local_buffer.fastx())?; + writer1.flush()?; + writer2.write_all(&self.local_buffer2)?; + writer2.flush()?; + } + } else { + // Interleaved output + if !self.local_buffer.fastx().is_empty() { + let mut writer = global.lock(); + writer.write_all(self.local_buffer.fastx())?; + writer.flush()?; + } + } } } - // Clear buffer after releasing the lock for better performance - self.local_buffer.clear(); + if let LocalOutput::Fastx(v) = &mut self.local_buffer { + v.clear(); + } self.local_buffer2.clear(); // Update global stats @@ -708,6 +1052,84 @@ impl<'a, Rf: Record> PairedParallelProcessor for FilterProcessor<'a> { Ok(()) } + + /// Merge any remaining local blocks into the shared state (per thread) + fn flush_thread(&mut self) -> Result<()> { + if let SharedOutput::Cbq(global) = &self.global_writer { + global.lock().ingest(self.local_buffer.cbq_mut())?; + } + Ok(()) + } +} + +impl ParallelProcessor for FilterProcessor { + fn process_record(&mut self, record: Rf) -> paraseq::parallel::Result<()> { + self.handle_record(&record)?; + Ok(()) + } + + fn on_batch_complete(&mut self) -> paraseq::parallel::Result<()> { + self.flush_batch()?; + Ok(()) + } + + fn on_thread_complete(&mut self) -> paraseq::parallel::Result<()> { + self.flush_thread()?; + Ok(()) + } +} + +impl PairedParallelProcessor for FilterProcessor { + fn process_record_pair(&mut self, record1: Rf, record2: Rf) -> paraseq::parallel::Result<()> { + self.handle_record_pair(&record1, &record2)?; + Ok(()) + } + + fn on_batch_complete(&mut self) -> paraseq::parallel::Result<()> { + self.flush_batch()?; + Ok(()) + } + + fn on_thread_complete(&mut self) -> paraseq::parallel::Result<()> { + self.flush_thread()?; + Ok(()) + } +} + +impl binseq::ParallelProcessor for FilterProcessor { + fn process_record(&mut self, record: R) -> binseq::Result<()> { + if record.is_paired() { + let record1 = CbqRecord { + id: record.sheader(), + seq: record.sseq(), + qual: record.has_quality().then(|| record.squal()), + }; + let record2 = CbqRecord { + id: record.xheader(), + seq: record.xseq(), + qual: record.has_quality().then(|| record.xqual()), + }; + self.handle_record_pair(&record1, &record2)?; + } else { + let record1 = CbqRecord { + id: record.sheader(), + seq: record.sseq(), + qual: record.has_quality().then(|| record.squal()), + }; + self.handle_record(&record1)?; + } + Ok(()) + } + + fn on_batch_complete(&mut self) -> binseq::Result<()> { + self.flush_batch()?; + Ok(()) + } + + fn on_thread_complete(&mut self) -> binseq::Result<()> { + self.flush_thread()?; + Ok(()) + } } pub fn run(config: &FilterConfig) -> Result { @@ -776,7 +1198,7 @@ pub fn run(config: &FilterConfig) -> Result { threshold ); } - return run_with_index(&minimizers, &header, &run_config); + return run_with_index(Arc::new(minimizers), &header, &run_config); } let (minimizers, header) = load_minimizers_cached(config.minimizers_path)?; @@ -797,7 +1219,7 @@ pub fn run(config: &FilterConfig) -> Result { /// Reusable entry point behind the Python bindings... load the index once, call repeatedly. /// Does no index-path validation; the index is already in memory. pub fn run_with_index( - minimizers: &MinimizerSet, + minimizers: Arc, header: &IndexHeader, config: &FilterRunConfig, ) -> Result { @@ -851,15 +1273,22 @@ pub fn run_with_index( .context("Failed to initialize thread pool"); } + check_input_paths(config)?; + + // Resolve formats and input metadata before opening or truncating outputs + let interleaved_stdin = config.input_path == "-" && config.input2_path.as_deref() == Some("-"); + let interleaved_input = config.interleaved || interleaved_stdin; + let output_format = resolve_output_format(config)?; + let (layout, prepared) = prepare_input(config, interleaved_input, output_format)?; + validate_input_output(&layout, output_format, config)?; + let mode = if config.deplete { "deplete" } else { "search" }; let mut input_type = String::new(); let mut options = Vec::::new(); - let interleaved_stdin = config.input_path == "-" && config.input2_path.as_deref() == Some("-"); - let interleaved_input = config.interleaved || interleaved_stdin; if interleaved_input { input_type.push_str("interleaved"); - } else if config.input2_path.is_some() { + } else if layout.paired { input_type.push_str("paired"); } else { input_type.push_str("single"); @@ -904,25 +1333,42 @@ pub fn run_with_index( ); } - check_input_paths(config)?; - - let writer = get_writer( - config.output_path.as_deref(), - config.compression_level, - compression_threads_per_output, - )?; - let writer2 = if let Some(output2) = config.output2_path.as_deref() { - if config.input2_path.is_some() || interleaved_input { - Some(get_writer( - Some(std::path::Path::new(output2)), + let (global_writer, writer2) = match output_format { + OutputFormat::Cbq => { + let cbq_writer = BinseqWriterBuilder::new(BinseqFormat::Cbq) + .paired(layout.paired) + .quality(layout.qualities && !config.output_fasta) + .headers(layout.headers || config.rename || config.rename_random) + .flags(false) + .build(get_writer( + config.output_path.as_deref(), + config.compression_level, + compression_threads_per_output, + )?) + .context("Failed to create CBQ writer")?; + (SharedOutput::Cbq(Arc::new(Mutex::new(cbq_writer))), None) + } + OutputFormat::Fastx => { + let writer = get_writer( + config.output_path.as_deref(), config.compression_level, compression_threads_per_output, - )?) - } else { - None + )?; + let writer2 = if let Some(output2) = config.output2_path.as_deref() { + if layout.paired { + Some(get_writer( + Some(std::path::Path::new(output2)), + config.compression_level, + compression_threads_per_output, + )?) + } else { + None + } + } else { + None + }; + (SharedOutput::Fastx(Arc::new(Mutex::new(writer))), writer2) } - } else { - None }; // Progress bar setup if not quiet @@ -958,100 +1404,32 @@ pub fn run_with_index( kmer_length, window_size, &processor_config, - writer, + global_writer, writer2, spinner.clone(), filtering_start_time, - ); - - // Check for empty files via metadata (fast path for uncompressed files <5 bytes) - let input1_empty = is_empty_file(&config.input_path)?; - let input2_empty = config - .input2_path - .as_deref() - .map(is_empty_file) - .transpose()? - .unwrap_or(false); + )?; // Process based on input type - use filtering threads (already calculated above) let num_threads = filtering_threads; - if interleaved_input { - // Interleaved paired input from stdin or a file - if input1_empty { - if !quiet { - eprintln!("Empty input file(s) detected"); - } - } else { - match create_paraseq_reader(Some(config.input_path.as_str())) { - Ok(reader) => { - reader.process_parallel_interleaved(&mut processor, num_threads)?; - } - Err(e) if is_empty_input_error(&e) => { - if !quiet { - eprintln!("Empty input file(s) detected"); - } - } - Err(e) => return Err(e), - } + match prepared { + InputPrep::Cbq(reader) => { + reader.process_parallel(processor.clone(), num_threads)?; } - } else if let Some(input2_path) = config.input2_path.as_deref() { - // Paired files - both must be empty or both non-empty - if input1_empty && input2_empty { - if !quiet { - eprintln!("Empty input file(s) detected"); - } - } else if input1_empty || input2_empty { - return Err(anyhow::anyhow!( - "One paired file is empty but the other is not" - )); - } else { - // Try to create readers, catching empty compressed files - let r1_result = create_paraseq_reader(Some(config.input_path.as_str())); - let r2_result = create_paraseq_reader(Some(input2_path)); - - match (r1_result, r2_result) { - (Ok(r1), Ok(r2)) => { - r1.process_parallel_paired(r2, &mut processor, num_threads)?; - } - (Err(e1), Err(e2)) if is_empty_input_error(&e1) && is_empty_input_error(&e2) => { - if !quiet { - eprintln!("Empty input file(s) detected"); - } - } - (Err(e), _) if is_empty_input_error(&e) => { - return Err(anyhow::anyhow!( - "First paired file appears empty while second is not" - )); - } - (_, Err(e)) if is_empty_input_error(&e) => { - return Err(anyhow::anyhow!( - "Second paired file appears empty while first is not" - )); - } - (Err(e), _) => return Err(e), - (_, Err(e)) => return Err(e), - } + InputPrep::FastxSingle(reader) => { + reader.process_parallel(&mut processor, num_threads)?; } - } else { - // Single file or stdin - if input1_empty { + InputPrep::FastxInterleaved(reader) => { + reader.process_parallel_interleaved(&mut processor, num_threads)?; + } + InputPrep::FastxPaired(reader1, reader2) => { + reader1.process_parallel_paired(reader2, &mut processor, num_threads)?; + } + InputPrep::Empty => { if !quiet { eprintln!("Empty input file(s) detected"); } - } else { - // Try to create reader, catching empty compressed files - match create_paraseq_reader(Some(config.input_path.as_str())) { - Ok(reader) => { - reader.process_parallel(&mut processor, num_threads)?; - } - Err(e) if is_empty_input_error(&e) => { - if !quiet { - eprintln!("Empty input file(s) detected"); - } - } - Err(e) => return Err(e), - } } } @@ -1064,6 +1442,11 @@ pub fn run_with_index( drop(final_stats); // Release lock + // Finish the CBQ stream: flush remaining blocks and write the embedded index + if let SharedOutput::Cbq(global) = &processor.global_writer { + global.lock().finish()?; + } + // Flush writers - they should auto-flush on drop drop(processor.global_writer); if let Some(w2) = processor.global_writer2 { diff --git a/src/index.rs b/src/index.rs index b6f00a1..be0d6aa 100644 --- a/src/index.rs +++ b/src/index.rs @@ -167,7 +167,7 @@ pub fn load_header_and_count>(path: &P) -> Result<(IndexHeader, u } #[cfg(feature = "cli")] -static INDEX: OnceLock<(PathBuf, crate::MinimizerSet, IndexHeader)> = OnceLock::new(); +static INDEX: OnceLock<(PathBuf, Arc, IndexHeader)> = OnceLock::new(); #[cfg(feature = "cli")] pub fn current_index_path() -> Option { @@ -177,18 +177,18 @@ pub fn current_index_path() -> Option { #[cfg(feature = "cli")] pub fn load_minimizers_cached( path: &Path, -) -> Result<(&'static crate::MinimizerSet, &'static IndexHeader)> { +) -> Result<(Arc, &'static IndexHeader)> { let (p, minimizers, header) = INDEX.get_or_init(|| { // Auto-detect exact vs BFF format let (m, h) = load_index_from_path_auto(path).unwrap(); - (path.to_owned(), m, h) + (path.to_owned(), Arc::new(m), h) }); assert_eq!( p, path, "Currently, the server can only have one index loaded." ); - Ok((minimizers, header)) + Ok((Arc::clone(minimizers), header)) } /// Load minimizers from a reader (generic over any Read impl) diff --git a/src/main.rs b/src/main.rs index 0e23fa9..d6f2ba9 100644 --- a/src/main.rs +++ b/src/main.rs @@ -33,11 +33,11 @@ enum Command { /// Path to minimizer index file index: PathBuf, - /// Optional path to fastx file (or - for stdin) + /// Optional path to fastx or CBQ file (or - for stdin) #[arg(default_value = "-")] input: String, - /// Optional path to second paired fastx file + /// Optional path to second paired fastx file (not supported with CBQ input) input2: Option, /// Minimum absolute number of minimizer hits for a match @@ -72,11 +72,11 @@ enum Command { #[arg(short = 'f', long = "fasta", default_value_t = false)] output_fasta: bool, - /// Path to output fastx file (stdout if not specified; detects .gz and .zst) + /// Path to output fastx file (stdout if not specified; detects .gz and .zst; .cbq suffix writes CBQ) #[arg(short = 'o', long = "output")] output: Option, - /// Optional path to second paired output fastx file (detects .gz and .zst) + /// Optional path to second paired output fastx file (detects .gz and .zst; not supported with CBQ output) #[arg(short = 'O', long = "output2")] output2: Option, diff --git a/tests/filter_tests.rs b/tests/filter_tests.rs index 9431b22..361e9b5 100644 --- a/tests/filter_tests.rs +++ b/tests/filter_tests.rs @@ -2576,3 +2576,276 @@ fn test_filter_paired_deplete_with_rename() { ); } } + +/// Write a small paired, headerless, quality-free CBQ file for validation tests +fn write_headerless_paired_cbq(path: &Path) { + use binseq::write::{BinseqWriterBuilder, Format}; + use binseq::SequencingRecordBuilder; + let file = File::create(path).unwrap(); + let mut writer = BinseqWriterBuilder::new(Format::Cbq) + .paired(true) + .build(file) + .unwrap(); + for seq in [b"ACGTACGTACGTACGTACGT".as_slice(), b"TGCAACGTACGTACGTACGT".as_slice()] { + writer + .push( + SequencingRecordBuilder::default() + .s_seq(seq) + .x_seq(seq) + .build() + .unwrap(), + ) + .unwrap(); + } + writer.finish().unwrap(); +} + +/// CBQ is only an I/O concern: FASTX -> CBQ -> FASTX must preserve reads, +/// pairing, headers, qualities, and summary counts exactly. +#[test] +fn cbq_roundtrip_matches_fastx() { + let temp_dir = tempdir().unwrap(); + let fasta_path = temp_dir.path().join("ref.fasta"); + let bin_path = temp_dir.path().join("ref.bin"); + let fastq_path = temp_dir.path().join("reads.fastq"); + let r1_path = temp_dir.path().join("reads_1.fastq"); + let r2_path = temp_dir.path().join("reads_2.fastq"); + + create_test_fasta(&fasta_path); + create_test_fastq(&fastq_path); + create_test_paired_fastq(&r1_path, &r2_path); + build_index(&fasta_path, &bin_path); + + // Single-end: FASTQ -> FASTQ baseline vs FASTQ -> CBQ + let baseline = temp_dir.path().join("baseline.fastq"); + let cbq_path = temp_dir.path().join("reads.cbq"); + let summary1 = temp_dir.path().join("summary1.json"); + let summary2 = temp_dir.path().join("summary2.json"); + + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&fastq_path) + .arg("--output") + .arg(&baseline) + .arg("--summary") + .arg(&summary1) + .assert() + .success(); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&fastq_path) + .arg("--output") + .arg(&cbq_path) + .arg("--summary") + .arg(&summary2) + .assert() + .success(); + + let baseline_content = fs::read_to_string(&baseline).unwrap(); + assert_eq!(count_records(&baseline_content), 2); + + // The CBQ itself: single, with headers and qualities, same record count + let reader = binseq::cbq::MmapReader::new(&cbq_path).unwrap(); + assert!(!reader.is_paired(), "single-end CBQ must be unpaired"); + assert!(reader.header().has_headers(), "CBQ must keep headers"); + assert!(reader.header().has_qualities(), "CBQ must keep qualities"); + assert_eq!(reader.num_records(), 2); + + // CBQ -> FASTQ matches the FASTQ baseline byte-for-byte + let roundtrip = temp_dir.path().join("roundtrip.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&cbq_path) + .arg("--output") + .arg(&roundtrip) + .assert() + .success(); + assert_eq!( + fs::read(&roundtrip).unwrap(), + fs::read(&baseline).unwrap(), + "CBQ round trip must match the FASTQ output" + ); + + // Summary sequence/base counts agree across formats + let s1: serde_json::Value = + serde_json::from_str(&fs::read_to_string(&summary1).unwrap()).unwrap(); + let s2: serde_json::Value = + serde_json::from_str(&fs::read_to_string(&summary2).unwrap()).unwrap(); + for field in ["seqs_in", "seqs_out", "bp_in", "bp_out"] { + assert_eq!(s1[field], s2[field], "summary field {field} differs"); + } + + // Paired: FASTQ -> one paired CBQ (2 native paired records) + let paired_cbq = temp_dir.path().join("paired.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&r1_path) + .arg(&r2_path) + .arg("--output") + .arg(&paired_cbq) + .assert() + .success(); + let reader = binseq::cbq::MmapReader::new(&paired_cbq).unwrap(); + assert!(reader.is_paired(), "paired FASTQ must produce a paired CBQ"); + assert_eq!(reader.num_records(), 2); + + // CBQ -> CBQ preserves pairing, headers, sequences, qualities + let cbq2 = temp_dir.path().join("paired2.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&paired_cbq) + .arg("--output") + .arg(&cbq2) + .assert() + .success(); + let reader = binseq::cbq::MmapReader::new(&cbq2).unwrap(); + assert!(reader.is_paired(), "CBQ -> CBQ must preserve pairing"); + assert_eq!(reader.num_records(), 2); + + // CBQ -> FASTQ matches the paired FASTQ -> FASTQ interleaved baseline + let paired_baseline = temp_dir.path().join("paired_baseline.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&r1_path) + .arg(&r2_path) + .arg("--output") + .arg(&paired_baseline) + .assert() + .success(); + let paired_roundtrip = temp_dir.path().join("paired_roundtrip.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&cbq2) + .arg("--output") + .arg(&paired_roundtrip) + .assert() + .success(); + assert_eq!( + fs::read(&paired_roundtrip).unwrap(), + fs::read(&paired_baseline).unwrap(), + "paired CBQ round trip must match the interleaved FASTQ output" + ); + + // --fasta creates a quality-free CBQ + let fasta_cbq = temp_dir.path().join("fasta.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1", "--fasta"]) + .arg(&bin_path) + .arg(&fastq_path) + .arg("--output") + .arg(&fasta_cbq) + .assert() + .success(); + let reader = binseq::cbq::MmapReader::new(&fasta_cbq).unwrap(); + assert!( + !reader.header().has_qualities(), + "--fasta must produce a quality-free CBQ" + ); + assert_eq!(reader.num_records(), 2); + + // Empty retained output is written as a structurally valid zero-record CBQ + let aaa_path = temp_dir.path().join("aaa.fasta"); + let aaa_bin = temp_dir.path().join("aaa.bin"); + create_test_fasta_aaa(&aaa_path); + build_index(&aaa_path, &aaa_bin); + let empty_cbq = temp_dir.path().join("empty.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-t", "1"]) + .arg(&aaa_bin) + .arg(&fastq_path) + .arg("--output") + .arg(&empty_cbq) + .assert() + .success(); + let empty_bytes = fs::read(&empty_cbq).unwrap(); + assert!( + empty_bytes.starts_with(b"CBQFILE"), + "empty retained output must be a CBQ file" + ); + // binseq 0.9.4 cannot reopen zero-record CBQ files (upstream index-cast + // bug); deacon treats such inputs as valid empty inputs. + let header = binseq::cbq::FileHeader::from_bytes(&empty_bytes[..64]).unwrap(); + assert!(header.has_qualities(), "empty CBQ must keep the writer quality flag"); + + // ... and it reads back as empty FASTX + let empty_out = temp_dir.path().join("empty_out.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&empty_cbq) + .arg("--output") + .arg(&empty_out) + .assert() + .success(); + assert!(fs::read_to_string(&empty_out).unwrap().is_empty()); + + // Invalid CBQ argument combinations fail before output creation + let no_create = temp_dir.path().join("should_not_exist.cbq"); + + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-t", "1"]) + .arg(&bin_path) + .arg(&cbq_path) + .arg(&r2_path) // CBQ input + INPUT2 + .arg("--output") + .arg(&no_create) + .assert() + .failure(); + assert!(!no_create.exists(), "CBQ input + INPUT2 must fail before output creation"); + + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "--interleaved", "-t", "1"]) + .arg(&bin_path) + .arg(&cbq_path) // CBQ input + --interleaved + .arg("--output") + .arg(&no_create) + .assert() + .failure(); + assert!(!no_create.exists(), "CBQ input + --interleaved must fail before output creation"); + + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-t", "1"]) + .arg(&bin_path) + .arg(&fastq_path) + .arg("--output") + .arg(&no_create) + .arg("--output2") + .arg(temp_dir.path().join("x.fastq")) // CBQ output + OUTPUT2 + .assert() + .failure(); + assert!(!no_create.exists(), "CBQ output + OUTPUT2 must fail before output creation"); + + let gz_cbq = temp_dir.path().join("bad.cbq.gz"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-t", "1"]) + .arg(&bin_path) + .arg(&fastq_path) + .arg("--output") + .arg(&gz_cbq) // compressed CBQ output + .assert() + .failure(); + assert!(!gz_cbq.exists(), ".cbq.gz output must fail before output creation"); + + // --check-pairs + headerless CBQ input + let headerless = temp_dir.path().join("headerless.cbq"); + write_headerless_paired_cbq(&headerless); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "--check-pairs", "-t", "1"]) + .arg(&bin_path) + .arg(&headerless) + .arg("--output") + .arg(&no_create) + .assert() + .failure(); + assert!( + !no_create.exists(), + "--check-pairs + headerless CBQ must fail before output creation" + ); +} From 8d657f5e6c89b4cec83f7d5da4b60d21a0982d64 Mon Sep 17 00:00:00 2001 From: alejandrogzi Date: Thu, 6 Aug 2026 15:18:07 +0200 Subject: [PATCH 2/2] chore(filter): drop footprint --- src/filter.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/filter.rs b/src/filter.rs index fee5099..f0e0596 100644 --- a/src/filter.rs +++ b/src/filter.rs @@ -332,7 +332,7 @@ fn prepare_input( interleaved_input: bool, output_format: OutputFormat, ) -> Result<(InputLayout, InputPrep)> { - // ponytail: CBQ input is file-only so the mmap reader stays the only CBQ + // CBQ input is file-only so the mmap reader stays the only CBQ // path; add binseq's streaming reader if stdin support is ever needed. if config.input_path == "-" || is_special_input_path(&config.input_path) { return prepare_fastx(config, interleaved_input, output_format);