From 98b300e8a28b53f665a276d38e4fe1d8599eee43 Mon Sep 17 00:00:00 2001 From: alejandrogzi Date: Thu, 6 Aug 2026 15:16:26 +0200 Subject: [PATCH 01/12] BREAKING CHANGE: v0.0.16 -> CBQ input/output implementation --- CHANGELOG.md | 1 + Cargo.lock | 173 ++++++++- Cargo.toml | 5 +- README.md | 6 +- deacon-py/src/lib.rs | 2 +- src/filter.rs | 817 ++++++++++++++++++++++++++++++------------ src/index.rs | 8 +- src/main.rs | 8 +- tests/filter_tests.rs | 275 ++++++++++++++ 9 files changed, 1059 insertions(+), 236 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ea5ef0e..fc1cef3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added - `deacon filter --ordered` uses paraseq 0.5.0 order preservation for deterministic output. Also exposed as `ordered` in the Python bindings and JSON summary. +- `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. ### Changed diff --git a/Cargo.lock b/Cargo.lock index 2d5dd3f..3e3759c 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 2.0.119", +] + [[package]] name = "autocfg" version = "1.5.1" @@ -155,6 +166,27 @@ dependencies = [ "virtue", ] +[[package]] +name = "binseq" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e468ca3b52411b87624a45d8bbdefd4e510c9121ae7da2a652fa6b90c90181dc" +dependencies = [ + "anyhow", + "auto_impl", + "bitnuc", + "bytemuck", + "byteorder", + "itoa", + "memchr", + "memmap2", + "num_cpus", + "rand 0.9.5", + "sucds", + "thiserror 2.0.20", + "zstd", +] + [[package]] name = "bitflags" version = "1.3.2" @@ -167,6 +199,12 @@ version = "2.13.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b588b76d00fde79687d7646a9b5bdf3cc0f655e0bbd080335a95d7e96f3587da" +[[package]] +name = "bitnuc" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d7871516cbb4e097e220623917e1491c995ee58c8018d929d4b1e76c378002bf" + [[package]] name = "bstr" version = "1.13.1" @@ -189,6 +227,20 @@ name = "bytemuck" version = "1.25.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "95832e849adfb21180ccb6826a99da14e5d266ae5c2e668e1602cf234f153797" +dependencies = [ + "bytemuck_derive", +] + +[[package]] +name = "bytemuck_derive" +version = "1.12.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fc0e56a716f1e132ff6bf4bdac1c944a3fcdc1cae65f70a4a2a1ac3b401d2d1f" +dependencies = [ + "proc-macro2", + "quote", + "syn 3.0.3", +] [[package]] name = "byteorder" @@ -249,7 +301,7 @@ checksum = "d524456ba66e72eb8b115ff89e01e497f8e6d11d78b70b1aa13c0fbd97540a81" dependencies = [ "cfg-if", "cpufeatures", - "rand_core", + "rand_core 0.10.1", ] [[package]] @@ -425,6 +477,7 @@ dependencies = [ "anyhow", "assert_cmd", "bincode", + "binseq", "clap", "ensure_simd", "flate2", @@ -693,6 +746,18 @@ dependencies = [ "wasm-bindgen", ] +[[package]] +name = "getrandom" +version = "0.3.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "899def5c37c4fd7b2664648c28120ecec138e4d395b459e5ca34f9cce2dd77fd" +dependencies = [ + "cfg-if", + "libc", + "r-efi 5.3.0", + "wasip2", +] + [[package]] name = "getrandom" version = "0.4.3" @@ -702,8 +767,8 @@ dependencies = [ "cfg-if", "js-sys", "libc", - "r-efi", - "rand_core", + "r-efi 6.0.0", + "rand_core 0.10.1", "wasm-bindgen", ] @@ -1010,6 +1075,15 @@ version = "2.8.3" source = "registry+https://github.com/rust-lang/crates.io-index" 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" version = "0.8.9" @@ -1118,7 +1192,7 @@ dependencies = [ "cfg-if", "ensure_simd", "mem_dbg", - "rand", + "rand 0.10.2", "wide", ] @@ -1198,6 +1272,15 @@ dependencies = [ "portable-atomic", ] +[[package]] +name = "ppv-lite86" +version = "0.2.21" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "85eae3c4ed2f50dcfe72643da4befc30deadb458a9b590d720cde2f2b1e97da9" +dependencies = [ + "zerocopy", +] + [[package]] name = "predicates" version = "3.1.4" @@ -1334,12 +1417,28 @@ dependencies = [ "proc-macro2", ] +[[package]] +name = "r-efi" +version = "5.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "69cdb34c158ceb288df11e18b4bd39de994f6657d83847bdffdbd7f346754b0f" + [[package]] name = "r-efi" version = "6.0.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "f8dcc9c7d52a811697d2151c701e0d08956f92b0e24136cf4cf27b57a6a0d9bf" +[[package]] +name = "rand" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b9ef1d0d795eb7d84685bca4f72f3649f064e6641543d3a8c415898726a57b41" +dependencies = [ + "rand_chacha", + "rand_core 0.9.5", +] + [[package]] name = "rand" version = "0.10.2" @@ -1348,7 +1447,26 @@ checksum = "c7f5fa3a058cd35567ef9bfa5e75732bee0f9e4c55fa90477bef2dfcdbc4be80" dependencies = [ "chacha20", "getrandom 0.4.3", - "rand_core", + "rand_core 0.10.1", +] + +[[package]] +name = "rand_chacha" +version = "0.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d3022b5f1df60f26e1ffddd6c66e8aa15de382ae63b3a0c1bfc0e4d3e3f325cb" +dependencies = [ + "ppv-lite86", + "rand_core 0.9.5", +] + +[[package]] +name = "rand_core" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "76afc826de14238e6e8c374ddcc1fa19e374fd8dd986b0d2af0d02377261d83c" +dependencies = [ + "getrandom 0.3.4", ] [[package]] @@ -1766,6 +1884,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.119" @@ -1950,6 +2078,15 @@ version = "0.11.1+wasi-snapshot-preview1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "ccf3ec651a847eb01de73ccad15eb7d99f80485de043efb2f370cd654f4ea44b" +[[package]] +name = "wasip2" +version = "1.0.4+wasi-0.2.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b67efb37e106e55ce722a510d6b5f9c17f083e5fc79afc2badeb12cc313d9487" +dependencies = [ + "wit-bindgen", +] + [[package]] name = "wasm-bindgen" version = "0.2.127" @@ -2228,6 +2365,12 @@ dependencies = [ "memchr", ] +[[package]] +name = "wit-bindgen" +version = "0.57.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1ebf944e87a7c253233ad6766e082e3cd714b5d03812acc24c318f549614536e" + [[package]] name = "xorf" version = "0.12.0" @@ -2238,6 +2381,26 @@ dependencies = [ "libm", ] +[[package]] +name = "zerocopy" +version = "0.8.56" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "556764e583adb45a9f8d413c2a147fa7e8d821e48e12b14fd560b607998b75eb" +dependencies = [ + "zerocopy-derive", +] + +[[package]] +name = "zerocopy-derive" +version = "0.8.56" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f2ab42fc20575779bd240faa45f94a74256f755c0fa9e89f0ede20d91d0cdfc1" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.119", +] + [[package]] name = "zeroize" version = "1.9.0" diff --git a/Cargo.toml b/Cargo.toml index 22c6498..8720214 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -41,6 +41,8 @@ indicatif = { version = "0.18", optional = true } serde_json = { version = "1.0", optional = true } niffler = { version = "3.0.1", default-features = false, optional = true } paraseq = { version = "0.5.0", 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.5", default-features = false, features = ["anyhow"], optional = true } ensure_simd = { version = "0.1.0", optional = true } # Compression deps (optional) @@ -59,10 +61,11 @@ predicates = "3.0" tempfile = "3.20" rstest = "0.26" nix = { version = "0.31", features = ["fs"] } +binseq = { version = "0.9.5", default-features = false, features = ["anyhow"] } [features] scalar = ["simd-minimizers/scalar", "packed-seq/scalar"] -cli = ["rayon", "parking_lot", "indicatif", "paraseq", "ensure_simd", "clap", "niffler", "serde_json"] +cli = ["rayon", "parking_lot", "indicatif", "paraseq", "binseq", "ensure_simd", "clap", "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 8d36384..6a37c02 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 @@ -130,6 +130,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/src/lib.rs b/deacon-py/src/lib.rs index a5d5043..c3d87a9 100644 --- a/deacon-py/src/lib.rs +++ b/deacon-py/src/lib.rs @@ -160,7 +160,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 b2c38b3..5041882 100644 --- a/src/filter.rs +++ b/src/filter.rs @@ -4,14 +4,18 @@ use crate::{ MinimizerSet, validate_unit_interval, }; 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; use paraseq::parallel::{PairedParallelProcessor, ParallelProcessor, ParallelReader}; use parking_lot::Mutex; 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}; @@ -22,6 +26,114 @@ 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]>, + index: u64, +} + +impl Record for CbqRecord<'_> { + fn id(&self) -> &[u8] { + self.id + } + + fn index(&self) -> u64 { + self.index + } + + 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) @@ -196,6 +308,190 @@ 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)> { + // 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) => { + let reader = + cbq::MmapReader::new(&config.input_path).context("Failed to open CBQ input")?; + let header = reader.header(); + let layout = InputLayout { + format: InputFormat::Cbq, + paired: header.is_paired(), + qualities: header.has_qualities(), + headers: header.has_headers(), + }; + // binseq's parallel reader rejects an empty record range, so + // zero-record CBQ files skip processing entirely + let prep = if reader.num_records() == 0 { + InputPrep::Empty + } else { + InputPrep::Cbq(reader) + }; + Ok((layout, prep)) + } + 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(()) +} + /// A record buffered without its header, which is written later by [`write_renamed`] /// once its final number is known (numbering depends on earlier batches finishing). #[derive(Clone)] @@ -405,9 +701,9 @@ pub struct FilterSummary { } #[derive(Clone)] -struct FilterProcessor<'a> { +struct FilterProcessor { // Minimizer matching parameters - minimizers: &'a MinimizerSet, + minimizers: Arc, rename: bool, output_fasta: bool, debug: bool, @@ -417,8 +713,8 @@ 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, // Headers awaiting a number, and how many records/pairs this batch kept @@ -429,7 +725,7 @@ struct FilterProcessor<'a> { rename_counter: Arc, // Global state - global_writer: Arc>, + global_writer: SharedOutput, global_writer2: Option>>, global_stats: Arc>, spinner: Option>>, @@ -446,17 +742,21 @@ 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, ) -> 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, @@ -474,14 +774,14 @@ 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), pending_renames: Vec::new(), pending_renames2: Vec::new(), batch_kept: 0, local_stats: ProcessingStats::default(), rename_counter: Arc::new(AtomicU64::new(0)), - global_writer: Arc::new(Mutex::new(writer)), + global_writer, global_writer2: writer2.map(|w| Arc::new(Mutex::new(w))), global_stats: Arc::new(Mutex::new(ProcessingStats::default())), spinner, @@ -490,12 +790,12 @@ impl<'a> FilterProcessor<'a> { } 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( @@ -504,13 +804,16 @@ impl<'a> FilterProcessor<'a> { seq: &[u8], read_suffix: &'static [u8], ) -> Result<()> { - let offset = self.local_buffer.len(); + if self.global_writer.is_cbq() { + return self.push_cbq(record, seq, read_suffix); + } + let offset = self.local_buffer.fastx().len(); let marker = format_record_to_buffer( record, seq, self.rename, self.output_fasta, - &mut self.local_buffer, + self.local_buffer.fastx_mut(), )?; if self.rename { self.pending_renames.push(PendingRename { @@ -558,6 +861,72 @@ impl<'a> FilterProcessor<'a> { + 1 } + /// Header for a CBQ record: the original id, or `number[ suffix]` when + /// renaming. CBQ headers are baked into blocks at push time, so numbers + /// are claimed per record rather than deferred like the FASTX path; + /// they are unique but sequential only with `--ordered` or `-t 1`. + fn cbq_header<'a>(&self, id: &'a [u8], counter: u64, read_suffix: &[u8]) -> Cow<'a, [u8]> { + if self.rename { + let mut header = counter.to_string().into_bytes(); + if !read_suffix.is_empty() { + header.push(b' '); + header.extend_from_slice(read_suffix); + } + Cow::Owned(header) + } else { + Cow::Borrowed(id) + } + } + + /// Push a single record into the thread-local CBQ writer + fn push_cbq( + &mut self, + record: &Rf, + seq: &[u8], + read_suffix: &[u8], + ) -> Result<()> { + let counter = if self.rename { + self.rename_counter.fetch_add(1, Ordering::Relaxed) + 1 + } else { + 0 + }; + let header = self.cbq_header(record.id(), counter, 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], + ) -> Result<()> { + let counter = if self.rename { + self.rename_counter.fetch_add(1, Ordering::Relaxed) + 1 + } else { + 0 + }; + let header1 = self.cbq_header(record1.id(), counter, b"/1"); + let header2 = self.cbq_header(record2.id(), counter, 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(); @@ -592,14 +961,9 @@ impl<'a> FilterProcessor<'a> { )); } } -} -impl<'a, Rf: Record> ParallelProcessor for FilterProcessor<'a> { - fn requires_ordering(&self) -> bool { - self.ordered - } - - 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; @@ -620,7 +984,7 @@ impl<'a, Rf: Record> ParallelProcessor for FilterProcessor<'a> { if decision.keep { self.local_stats.output_bp += seq.len() as u64; - self.write_record(&record, &seq, b"")?; + self.write_record(record, &seq, b"")?; if self.rename { self.batch_kept += 1; } @@ -632,62 +996,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(); - if self.rename { - let base = self.reserve_rename_ids(); - write_renamed( - &mut global_writer, - &self.local_buffer, - &self.pending_renames, - base, - )?; - } else { - global_writer.write_all(&self.local_buffer)?; - } - global_writer.flush()?; - } - - // Clear buffer after releasing the lock - self.local_buffer.clear(); - self.pending_renames.clear(); - self.batch_kept = 0; - - // 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 requires_ordering(&self) -> bool { - self.ordered - } - - 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(); @@ -714,14 +1030,17 @@ impl<'a, Rf: Record> PairedParallelProcessor for FilterProcessor<'a> { if decision.keep { self.local_stats.output_bp += (seq1.len() + seq2.len()) as u64; // Both mates take the pair's number - 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)?; + } else if self.global_writer2.is_some() { // Separate outputs - self.write_record(&record1, &seq1, b"/1")?; - self.write_record_to_buffer2(&record2, &seq2, b"/2")?; + self.write_record(record1, &seq1, b"/1")?; + self.write_record_to_buffer2(record2, &seq2, b"/2")?; } else { // Interleaved output - self.write_record(&record1, &seq1, b"/1")?; - self.write_record(&record2, &seq2, b"/2")?; + self.write_record(record1, &seq1, b"/1")?; + self.write_record(record2, &seq2, b"/2")?; } if self.rename { self.batch_kept += 1; @@ -734,51 +1053,65 @@ 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(); - - if self.rename { - // Both mates share one block of numbers - let base = self.reserve_rename_ids(); - write_renamed( - &mut writer1, - &self.local_buffer, - &self.pending_renames, - base, - )?; - write_renamed( - &mut writer2, - &self.local_buffer2, - &self.pending_renames2, - base, - )?; - } else { - writer1.write_all(&self.local_buffer)?; - writer2.write_all(&self.local_buffer2)?; - } - writer1.flush()?; - 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(); - if self.rename { - let base = self.reserve_rename_ids(); - write_renamed(&mut writer, &self.local_buffer, &self.pending_renames, base)?; + 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(); + + if self.rename { + // Both mates share one block of numbers + let base = self.reserve_rename_ids(); + write_renamed( + &mut writer1, + self.local_buffer.fastx(), + &self.pending_renames, + base, + )?; + write_renamed( + &mut writer2, + &self.local_buffer2, + &self.pending_renames2, + base, + )?; + } else { + writer1.write_all(self.local_buffer.fastx())?; + writer2.write_all(&self.local_buffer2)?; + } + writer1.flush()?; + writer2.flush()?; + } } else { - writer.write_all(&self.local_buffer)?; + // Interleaved output + if !self.local_buffer.fastx().is_empty() { + let mut writer = global.lock(); + if self.rename { + let base = self.reserve_rename_ids(); + write_renamed( + &mut writer, + self.local_buffer.fastx(), + &self.pending_renames, + base, + )?; + } else { + writer.write_all(self.local_buffer.fastx())?; + } + writer.flush()?; + } } - 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(); self.pending_renames.clear(); self.pending_renames2.clear(); @@ -802,6 +1135,89 @@ 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 requires_ordering(&self) -> bool { + self.ordered + } + + 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 requires_ordering(&self) -> bool { + self.ordered + } + + 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<()> { + let record1 = CbqRecord { + id: record.sheader(), + seq: record.sseq(), + qual: record.has_quality().then(|| record.squal()), + index: record.index(), + }; + if record.is_paired() { + let record2 = CbqRecord { + id: record.xheader(), + seq: record.xseq(), + qual: record.has_quality().then(|| record.xqual()), + index: record.index(), + }; + self.handle_record_pair(&record1, &record2)?; + } else { + 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 { @@ -874,7 +1290,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)?; @@ -895,7 +1311,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 { @@ -950,15 +1366,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"); @@ -1003,25 +1426,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) + .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 @@ -1057,100 +1497,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), - } } } @@ -1163,6 +1535,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 7e4426b..5fbf49b 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 449e79e..cdb2c94 100644 --- a/src/main.rs +++ b/src/main.rs @@ -34,11 +34,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 @@ -69,11 +69,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 cb2f2b6..2fe1c2d 100644 --- a/tests/filter_tests.rs +++ b/tests/filter_tests.rs @@ -2783,3 +2783,278 @@ 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" + ); + let reader = binseq::cbq::MmapReader::new(&empty_cbq).unwrap(); + assert_eq!(reader.num_records(), 0, "empty CBQ must hold zero records"); + assert!( + reader.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 777db1389825e39797b0b583e79eed4801556897 Mon Sep 17 00:00:00 2001 From: Bede Constantinides Date: Wed, 12 Aug 2026 15:52:17 +0100 Subject: [PATCH 02/12] CBQ-related changes on top of PR 95 - io abstractions: Output enum, ReadView etc - Fix --ordered, --check-pairs, flags passthrough, compression level validation - customisable default 16 MiB blocks or input block size if cbq, w/e larger - Comments and misc trivial CLI churn for CBQ etc --- CHANGELOG.md | 10 +- README.md | 22 +- deacon-py/src/lib.rs | 5 +- src/filter.rs | 944 ++++++++++++++++++++++-------------------- src/lib.rs | 9 +- src/main.rs | 39 +- tests/filter_tests.rs | 525 ++++++++++++++++++++++- 7 files changed, 1044 insertions(+), 510 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index fc1cef3..658028e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,20 +5,20 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## [0.17.0] - 2026-08-11 +## [Unreleased] ### Added -- `deacon filter --ordered` uses paraseq 0.5.0 order preservation for deterministic output. Also exposed as `ordered` in the Python bindings and JSON summary. -- `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. +- `deacon filter --ordered` preserves output order via paraseq 0.5.0; also exposed as `ordered` in the Python bindings and JSON summary. +- `deacon filter` accepts CBQ (BINSEQ) input, detected automatically, and writes CBQ when given a `.cbq` path; paired reads are stored as native paired records in a single file. Block size is set with `--cbq-block-size` (default 16 MiB). ### Changed -- Python `Index.filter()` options after `fastq2` are now keyword-only. +- Python `Index.filter()` options after `fastq2` are keyword-only. ### Removed -- `deacon filter --rename-random`, and `rename_random` from the JSON summary and Python bindings. +- Removed `--rename-random` CLI arg and `rename_random` from the JSON summary and Python bindings. ## [0.16.0] - 2026-08-09 diff --git a/README.md b/README.md index 6a37c02..58e9046 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`. 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`). +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`). With `--ordered`, CBQ input (and CBQ output combined with `--rename`) falls back to one filtering thread to guarantee deterministic output; without `--ordered`, CBQ `--rename` numbers are unique but not sequential across threads. #### Examples @@ -130,7 +130,7 @@ 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 +# CBQ (BINSEQ) input and output deacon filter panhuman-1.k31w15.idx reads.fq.gz -o filt.cbq deacon filter panhuman-1.k31w15.idx filt.cbq -o filt.fq.gz @@ -185,13 +185,13 @@ A differentiating feature of Deacon is the ease of combining, subtracting and in ```bash $ deacon filter -h -Retain or deplete sequence records with sufficient minimizer hits to an indexed query +Retain or deplete sequence records with sufficient minimizer hits to the index Usage: deacon filter [OPTIONS] [INPUT] [INPUT2] Arguments: Path to minimizer index file - [INPUT] Optional path to fastx file (or - for stdin) [default: -] + [INPUT] Optional path to fastx or binseq cbq file (or - for stdin) [default: -] [INPUT2] Optional path to second paired fastx file Options: @@ -210,9 +210,9 @@ Options: -f, --fasta Output FASTA format regardless of input format -o, --output - Path to output fastx file (stdout if not specified; detects .gz and .zst) + Path to output file (stdout by default; detects fastx with .gz, .zst, .xz, or binseq with .cbq) -O, --output2 - Optional path to second paired output fastx file (detects .gz and .zst) + Optional path to second paired output fastx file (detects .gz, .zst, .xz) -s, --summary Path to JSON summary output file -t, --threads @@ -220,17 +220,19 @@ Options: --compression-threads Number of threads used for output compression (0 = auto) [default: 0] --compression-level - Output compression level (1-9 for gz & xz; 1-22 for zstd) [default: 2] + Output compression level (1-9 for gz & xz; 1-22 for zstd, including cbq) [default: 2] + --cbq-block-size + cbq output block size in MiB (or cbq input block size if higher) [default: 16] --interleaved - Treat INPUT as interleaved paired reads from a file or stdin + Treat INPUT as interleaved paired records from single file or stdin --ordered Preserve input record ordering (deterministic, slightly slower) --check-pairs Validate paired record names (Illumina CASAVA or /1 /2 suffixes) - -q, --quiet - Suppress progress reporting --debug Output sequences with minimizer hits to stderr + -q, --quiet + Suppress progress reporting -h, --help Print help ``` diff --git a/deacon-py/src/lib.rs b/deacon-py/src/lib.rs index c3d87a9..a37fab6 100644 --- a/deacon-py/src/lib.rs +++ b/deacon-py/src/lib.rs @@ -5,8 +5,8 @@ use std::path::{Path, PathBuf}; use std::sync::Arc; use ::deacon::{ - ComplexityAlgorithm, FilterRunConfig, IndexHeader, MinimizerSet, index_fetch, - load_index_from_path_auto, run_with_index, + ComplexityAlgorithm, DEFAULT_CBQ_BLOCK_SIZE_MIB, FilterRunConfig, IndexHeader, MinimizerSet, + index_fetch, load_index_from_path_auto, run_with_index, }; use pyo3::exceptions::{PyRuntimeError, PyValueError}; use pyo3::prelude::*; @@ -151,6 +151,7 @@ impl Index { ordered, threads, compression_level, + cbq_block_size: DEFAULT_CBQ_BLOCK_SIZE_MIB, compression_threads, debug, quiet, diff --git a/src/filter.rs b/src/filter.rs index 5041882..a98bd51 100644 --- a/src/filter.rs +++ b/src/filter.rs @@ -23,114 +23,323 @@ use std::time::Instant; const OUTPUT_BUFFER_SIZE: usize = 8 * 1024 * 1024; // Opt: 8MB output buffer const DEFAULT_BUFFER_SIZE: usize = 64 * 1024; +/// Default CBQ block size in MiB; a block must hold the largest single record +pub const DEFAULT_CBQ_BLOCK_SIZE_MIB: u16 = 16; type BoxedWriter = Box; +/// CBQ output is always a named file, no stdout +type CbqWriter = binseq::BinseqWriter>; -/// Input file format +/// Sequence file format #[derive(Clone, Copy, PartialEq, Eq)] -enum InputFormat { +enum Format { Fastx, Cbq, } -/// Output file format -#[derive(Clone, Copy, PartialEq, Eq)] -enum OutputFormat { - Fastx, - Cbq, -} - -/// Normalized input metadata, populated before any output is opened +/// Input metadata, resolved before any output is opened struct InputLayout { - format: InputFormat, + format: Format, paired: bool, qualities: bool, headers: bool, + /// CBQ per-record flags present (always false for FASTX) + flags: bool, + /// CBQ input block size, a lower bound for CBQ output (None for FASTX) + block_size: Option, } -/// Shared output state; every processor clone merges into this under a lock -#[derive(Clone)] -enum SharedOutput { - Fastx(Arc>), - Cbq(Arc>>), +#[allow(clippy::large_enum_variant)] +enum Input { + Cbq(cbq::MmapReader), + FastxSingle(Reader>), + FastxInterleaved(Reader>), + FastxPaired( + Reader>, + Reader>, + ), + Empty, +} + +/// Borrowed view of one read, shared by every reader and writer +struct ReadView<'a> { + id: &'a [u8], + seq: &'a [u8], + qual: Option<&'a [u8]>, + /// CBQ per-record flag, preserved on CBQ output (None elsewhere) + flag: Option, } -/// Thread-local output buffer merged into the shared output on batch completion +/// Format-specific output state: thread-local buffers plus the shared +/// writers each processor clone merges into #[allow(clippy::large_enum_variant)] // mirrors binseq's BinseqWriter, which is not boxed upstream #[derive(Clone)] -enum LocalOutput { - Fastx(Vec), - Cbq(binseq::BinseqWriter>), +enum Output { + Fastx { + local: Vec, + local2: Vec, + pending: Vec, + pending2: Vec, + batch_kept: u64, + shared: Arc>, + shared2: Option>>, + }, + Cbq { + /// Thread-local block writer; completed blocks merge into the shared writer + local: binseq::BinseqWriter>, + shared: Arc>, + }, } -impl SharedOutput { - fn is_cbq(&self) -> bool { - matches!(self, Self::Cbq(_)) +impl Output { + fn fastx(writer: BoxedWriter, writer2: Option) -> Self { + Output::Fastx { + local: Vec::with_capacity(DEFAULT_BUFFER_SIZE), + local2: Vec::with_capacity(DEFAULT_BUFFER_SIZE), + pending: Vec::new(), + pending2: Vec::new(), + batch_kept: 0, + shared: Arc::new(Mutex::new(writer)), + shared2: writer2.map(|w| Arc::new(Mutex::new(w))), + } } -} -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 cbq(writer: CbqWriter) -> Result { + let shared = Arc::new(Mutex::new(writer)); + let local = shared.lock().new_headless_buffer()?; + Ok(Output::Cbq { local, shared }) + } + + /// Format one read into a FASTX local buffer, deferring the header to + /// [`write_renamed`] when renaming + fn push_fastx( + buffer: &mut Vec, + pending: &mut Vec, + batch_kept: u64, + read: &ReadView, + suffix: &'static [u8], + rename: bool, + output_fasta: bool, + ) -> Result<()> { + let offset = buffer.len(); + let marker = format_record_to_buffer(read, rename, output_fasta, buffer)?; + if rename { + pending.push(PendingRename { + offset, + ordinal: batch_kept, + marker, + suffix, + }); } + Ok(()) } - fn fastx(&self) -> &[u8] { + /// Write one unpaired read + fn push_read( + &mut self, + read: &ReadView, + rename: bool, + output_fasta: bool, + rename_counter: &AtomicU64, + ) -> Result<()> { match self { - Self::Fastx(v) => v, - Self::Cbq(_) => unreachable!("CBQ local buffer used in FASTX path"), + Output::Fastx { + local, + pending, + batch_kept, + .. + } => { + Self::push_fastx(local, pending, *batch_kept, read, b"", rename, output_fasta)?; + if rename { + *batch_kept += 1; + } + } + Output::Cbq { local, .. } => { + let number = rename.then(|| rename_counter.fetch_add(1, Ordering::Relaxed) + 1); + let header = cbq_header(read.id, number, b""); + let seq_record = SequencingRecordBuilder::default() + .s_seq(read.seq) + .s_header(&header) + .opt_s_qual(read.qual) + .opt_flag(read.flag) + .build()?; + local.push(seq_record)?; + } } + Ok(()) } - fn cbq_mut(&mut self) -> &mut binseq::BinseqWriter> { + /// Write a kept pair: one native CBQ record, two FASTX files, or + /// interleaved FASTX. Both mates take the pair's number when renaming. + fn push_pair( + &mut self, + read1: &ReadView, + read2: &ReadView, + rename: bool, + output_fasta: bool, + rename_counter: &AtomicU64, + ) -> Result<()> { match self { - Self::Cbq(w) => w, - Self::Fastx(_) => unreachable!("FASTX local buffer used in CBQ path"), + Output::Fastx { + local, + local2, + pending, + pending2, + batch_kept, + shared2, + .. + } => { + Self::push_fastx( + local, + pending, + *batch_kept, + read1, + b"/1", + rename, + output_fasta, + )?; + if shared2.is_some() { + // Separate outputs + Self::push_fastx( + local2, + pending2, + *batch_kept, + read2, + b"/2", + rename, + output_fasta, + )?; + } else { + // Interleaved output + Self::push_fastx( + local, + pending, + *batch_kept, + read2, + b"/2", + rename, + output_fasta, + )?; + } + if rename { + *batch_kept += 1; + } + } + Output::Cbq { local, .. } => { + let number = rename.then(|| rename_counter.fetch_add(1, Ordering::Relaxed) + 1); + let header1 = cbq_header(read1.id, number, b"/1"); + let header2 = cbq_header(read2.id, number, b"/2"); + let seq_record = SequencingRecordBuilder::default() + .s_seq(read1.seq) + .s_header(&header1) + .opt_s_qual(read1.qual) + .x_seq(read2.seq) + .x_header(&header2) + .opt_x_qual(read2.qual) + .opt_flag(read1.flag) + .build()?; + local.push(seq_record)?; + } } + Ok(()) } -} - -/// 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]>, - index: u64, -} - -impl Record for CbqRecord<'_> { - fn id(&self) -> &[u8] { - self.id - } + /// Merge thread-local buffers into the shared writer(s) (per batch). Rename + /// numbers are claimed under the writer lock: monotone in output order, but + /// deterministic only with `--ordered` or `-t 1`. + fn flush_batch( + &mut self, + rename: bool, + rename_counter: &AtomicU64, + ordered: bool, + ) -> Result<()> { + match self { + Output::Cbq { local, shared } => { + if ordered { + // Drain the incomplete block too, or records linger across batches + shared.lock().ingest(local)?; + } else { + shared.lock().ingest_completed(local)?; + } + } + Output::Fastx { + local, + local2, + pending, + pending2, + batch_kept, + shared, + shared2, + } => { + if let Some(shared2) = shared2 { + // Atomic paired batch writing + if !local.is_empty() || !local2.is_empty() { + let mut writer1 = shared.lock(); + let mut writer2 = shared2.lock(); - fn index(&self) -> u64 { - self.index + if rename { + // Both mates share one block of numbers + let base = rename_counter.fetch_add(*batch_kept, Ordering::Relaxed) + 1; + write_renamed(&mut writer1, local, pending, base)?; + write_renamed(&mut writer2, local2, pending2, base)?; + } else { + writer1.write_all(local)?; + writer2.write_all(local2)?; + } + writer1.flush()?; + writer2.flush()?; + } + } else if !local.is_empty() { + let mut writer = shared.lock(); + if rename { + let base = rename_counter.fetch_add(*batch_kept, Ordering::Relaxed) + 1; + write_renamed(&mut writer, local, pending, base)?; + } else { + writer.write_all(local)?; + } + writer.flush()?; + } + local.clear(); + local2.clear(); + pending.clear(); + pending2.clear(); + *batch_kept = 0; + } + } + Ok(()) } - fn seq(&self) -> Cow<'_, [u8]> { - Cow::Borrowed(self.seq) + /// Merge any remaining local CBQ blocks into the shared writer (per thread) + fn flush_thread(&mut self) -> Result<()> { + if let Output::Cbq { local, shared } = self { + shared.lock().ingest(local)?; + } + Ok(()) } - fn seq_raw(&self) -> &[u8] { - self.seq + /// Finish a CBQ stream: flush remaining blocks and write the embedded index + fn finish(&self) -> Result<()> { + if let Output::Cbq { shared, .. } = self { + shared.lock().finish()?; + } + Ok(()) } +} - fn qual(&self) -> Option<&[u8]> { - self.qual +/// Header for a CBQ record: the original id, or `number[ suffix]` when +/// renaming. Numbers are claimed per record: unique, but sequential only +/// with `--ordered` or `-t 1`. +fn cbq_header<'a>(id: &'a [u8], number: Option, suffix: &[u8]) -> Cow<'a, [u8]> { + match number { + Some(n) => { + let mut header = n.to_string().into_bytes(); + if !suffix.is_empty() { + header.push(b' '); + header.extend_from_slice(suffix); + } + Cow::Owned(header) + } + None => Cow::Borrowed(id), } } @@ -168,6 +377,8 @@ pub struct FilterRunConfig { pub threads: u16, /// Compression level for output files (1-22 for zst, 1-9 for gz) pub compression_level: u8, + /// cbq output block size in MiB (raised to the cbq input's block size if larger) + pub cbq_block_size: u16, /// Number of threads for compression (0 = auto) pub compression_threads: u16, /// Debug mode: output sequences with minimizer hits to stderr @@ -309,80 +520,68 @@ 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) +fn resolve_output_format(config: &FilterRunConfig) -> Format { + match config.output_path.as_deref() { + Some(path) if path.to_string_lossy().ends_with(".cbq") => Format::Cbq, + _ => Format::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)> { - // 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. +fn open_input(config: &FilterRunConfig, interleaved_input: bool) -> Result<(InputLayout, Input)> { + // CBQ input is file-only (mmap reader); stdin would need binseq's streaming reader if config.input_path == "-" || is_special_input_path(&config.input_path) { - return prepare_fastx(config, interleaved_input, output_format); + return open_fastx(config, interleaved_input); } // Regular files: sniff the CBQ magic; everything else is FASTX - let mut file = File::open(&config.input_path) + let 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]) { + let file_len = file.metadata()?.len(); + let mut magic = Vec::with_capacity(64); + file.take(64).read_to_end(&mut magic)?; + match BinseqFormat::sniff(&magic) { Some(BinseqFormat::Cbq) => { + // The mmap reader slices header and footer unchecked, so guard the length + let min_len = (std::mem::size_of::() + + std::mem::size_of::()) as u64; + if file_len < min_len { + anyhow::bail!("Truncated or corrupt CBQ input: {}", config.input_path); + } let reader = cbq::MmapReader::new(&config.input_path).context("Failed to open CBQ input")?; let header = reader.header(); let layout = InputLayout { - format: InputFormat::Cbq, + format: Format::Cbq, paired: header.is_paired(), qualities: header.has_qualities(), headers: header.has_headers(), + flags: header.has_flags(), + block_size: Some(header.block_size as usize), }; - // binseq's parallel reader rejects an empty record range, so - // zero-record CBQ files skip processing entirely - let prep = if reader.num_records() == 0 { - InputPrep::Empty + // binseq's parallel reader rejects an empty record range + let input = if reader.num_records() == 0 { + Input::Empty } else { - InputPrep::Cbq(reader) + Input::Cbq(reader) }; - Ok((layout, prep)) + Ok((layout, input)) } - Some(f) => anyhow::bail!("{f:?} input is not supported; convert it to FASTQ or CBQ"), - None => prepare_fastx(config, interleaved_input, output_format), + Some(f) => anyhow::bail!("{f:?} input is not supported"), + None => open_fastx(config, interleaved_input), } } -/// 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)> { +/// Open a FASTX input and resolve its layout +fn open_fastx(config: &FilterRunConfig, interleaved_input: bool) -> Result<(InputLayout, Input)> { let layout = InputLayout { - format: InputFormat::Fastx, + format: Format::Fastx, paired: interleaved_input || config.input2_path.is_some(), qualities: false, headers: true, + flags: false, + block_size: None, }; let input1_empty = is_empty_file(&config.input_path)?; @@ -395,24 +594,27 @@ fn prepare_fastx( if interleaved_input { if input1_empty { - return Ok((layout, InputPrep::Empty)); + return Ok((layout, Input::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), + InputLayout { + qualities, + ..layout + }, + Input::FastxInterleaved(reader), )) } - Err(e) if is_empty_input_error(&e) => Ok((layout, InputPrep::Empty)), + Err(e) if is_empty_input_error(&e) => Ok((layout, Input::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)); + return Ok((layout, Input::Empty)); } if input1_empty || input2_empty { return Err(anyhow::anyhow!( @@ -423,22 +625,17 @@ fn prepare_fastx( 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), + InputLayout { + qualities, + ..layout + }, + Input::FastxPaired(reader1, reader2), )) } (Err(e1), Err(e2)) if is_empty_input_error(&e1) && is_empty_input_error(&e2) => { - Ok((layout, InputPrep::Empty)) + Ok((layout, Input::Empty)) } (Err(e), _) if is_empty_input_error(&e) => Err(anyhow::anyhow!( "First paired file appears empty while second is not" @@ -452,17 +649,20 @@ fn prepare_fastx( } if input1_empty { - return Ok((layout, InputPrep::Empty)); + return Ok((layout, Input::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), + InputLayout { + qualities, + ..layout + }, + Input::FastxSingle(reader), )) } - Err(e) if is_empty_input_error(&e) => Ok((layout, InputPrep::Empty)), + Err(e) if is_empty_input_error(&e) => Ok((layout, Input::Empty)), Err(e) => Err(e), } } @@ -470,22 +670,32 @@ fn prepare_fastx( /// Validate format combinations before any output file is opened or truncated fn validate_input_output( layout: &InputLayout, - output_format: OutputFormat, + output_format: Format, config: &FilterRunConfig, ) -> Result<()> { - if layout.format == InputFormat::Cbq && config.input2_path.is_some() { + if layout.format == Format::Cbq && config.input2_path.is_some() { anyhow::bail!("CBQ input does not support INPUT2"); } - if layout.format == InputFormat::Cbq && config.interleaved { + if layout.format == Format::Cbq && config.interleaved { anyhow::bail!("CBQ input does not support --interleaved"); } - if output_format == OutputFormat::Cbq && config.output2_path.is_some() { + if output_format == Format::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 + .output2_path + .as_deref() + .is_some_and(|output2| output2.ends_with(".cbq")) + { + anyhow::bail!("OUTPUT2 cannot be CBQ; CBQ pairing is native (use a single --output)"); + } + if output_format == Format::Cbq && !(1..=22).contains(&config.compression_level) { + anyhow::bail!( + "Invalid CBQ compression level {} (must be 1-22)", + config.compression_level + ); } - if config.check_pairs && layout.format == InputFormat::Cbq && !layout.headers { + if config.check_pairs && layout.format == Format::Cbq && !layout.headers { anyhow::bail!("--check-pairs requires CBQ input with headers"); } validate_check_pairs_mode(config.check_pairs, layout.paired)?; @@ -509,32 +719,31 @@ struct PendingRename { /// Format a record into a buffer (FASTA/FASTQ), returning its line prefix /// `seq` is the newline-stripped sequence from `record.seq()`. /// When renaming, the header is left to [`write_renamed`]. -fn format_record_to_buffer( - record: &R, - seq: &[u8], +fn format_record_to_buffer( + read: &ReadView, rename: bool, output_fasta: bool, buffer: &mut Vec, ) -> Result { - let is_fasta = output_fasta || record.qual().is_none(); + let is_fasta = output_fasta || read.qual.is_none(); let marker = if is_fasta { b'>' } else { b'@' }; // Header (omitted when renaming) if !rename { buffer.push(marker); - buffer.extend_from_slice(record.id()); + buffer.extend_from_slice(read.id); buffer.write_all(b"\n")?; } // Sequence - buffer.extend_from_slice(seq); + buffer.extend_from_slice(read.seq); if is_fasta { buffer.write_all(b"\n")?; } else { // FASTQ: plus and qual lines buffer.write_all(b"\n+\n")?; - if let Some(qual) = record.qual() { + if let Some(qual) = read.qual { buffer.extend_from_slice(qual); } buffer.write_all(b"\n")?; @@ -601,7 +810,14 @@ fn count_compressed_outputs(config: &FilterRunConfig) -> u8 { count } -// Return a suitable writer for the output path extension +/// Open a CBQ output file directly, bypassing the compressed writer stack +fn open_cbq_output(path: &std::path::Path) -> Result> { + let file = File::create(path) + .with_context(|| format!("Failed to create output file: {}", path.display()))?; + Ok(BufWriter::with_capacity(OUTPUT_BUFFER_SIZE, file)) +} + +/// Return a suitable writer for the output path extension #[cfg_attr(not(feature = "compression"), allow(unused_variables))] fn get_writer( output_path: Option<&std::path::Path>, @@ -712,21 +928,13 @@ struct FilterProcessor { ordered: bool, kernel: FilterKernel, - // Local buffers - local_buffer: LocalOutput, - local_buffer2: Vec, // Second buffer for paired FASTX output + output: Output, local_stats: ProcessingStats, - // Headers awaiting a number, and how many records/pairs this batch kept - pending_renames: Vec, - pending_renames2: Vec, // Parallel to local_buffer2 - batch_kept: u64, - /// Shared across workers, handing each batch the first number of its block + /// Shared across workers, handing out rename numbers rename_counter: Arc, // Global state - global_writer: SharedOutput, - global_writer2: Option>>, global_stats: Arc>, spinner: Option>>, filtering_start_time: Instant, @@ -748,15 +956,10 @@ impl FilterProcessor { kmer_length: u8, window_size: u8, config: &FilterProcessorConfig, - global_writer: SharedOutput, - writer2: Option, + output: Output, spinner: Option>>, filtering_start_time: Instant, ) -> 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, @@ -774,15 +977,9 @@ impl FilterProcessor { prefix_length: config.prefix_length, }, )?, - local_buffer, - local_buffer2: Vec::with_capacity(DEFAULT_BUFFER_SIZE), - pending_renames: Vec::new(), - pending_renames2: Vec::new(), - batch_kept: 0, + output, local_stats: ProcessingStats::default(), 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, @@ -790,7 +987,7 @@ impl FilterProcessor { } 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 { @@ -798,135 +995,6 @@ impl FilterProcessor { .classify_pair(&self.minimizers, seq1, seq2, self.debug) } - fn write_record( - &mut self, - record: &Rf, - seq: &[u8], - read_suffix: &'static [u8], - ) -> Result<()> { - if self.global_writer.is_cbq() { - return self.push_cbq(record, seq, read_suffix); - } - let offset = self.local_buffer.fastx().len(); - let marker = format_record_to_buffer( - record, - seq, - self.rename, - self.output_fasta, - self.local_buffer.fastx_mut(), - )?; - if self.rename { - self.pending_renames.push(PendingRename { - offset, - ordinal: self.batch_kept, - marker, - suffix: read_suffix, - }); - } - Ok(()) - } - - fn write_record_to_buffer2( - &mut self, - record: &Rf, - seq: &[u8], - read_suffix: &'static [u8], - ) -> Result<()> { - let offset = self.local_buffer2.len(); - let marker = format_record_to_buffer( - record, - seq, - self.rename, - self.output_fasta, - &mut self.local_buffer2, - )?; - if self.rename { - self.pending_renames2.push(PendingRename { - offset, - ordinal: self.batch_kept, - marker, - suffix: read_suffix, - }); - } - Ok(()) - } - - /// Claim this batch's block of numbers, returning the first - /// - /// Called under the writer lock: numbers are monotone in output order, - /// but deterministic only with `--ordered` or `-t 1`. - fn reserve_rename_ids(&self) -> u64 { - self.rename_counter - .fetch_add(self.batch_kept, Ordering::Relaxed) - + 1 - } - - /// Header for a CBQ record: the original id, or `number[ suffix]` when - /// renaming. CBQ headers are baked into blocks at push time, so numbers - /// are claimed per record rather than deferred like the FASTX path; - /// they are unique but sequential only with `--ordered` or `-t 1`. - fn cbq_header<'a>(&self, id: &'a [u8], counter: u64, read_suffix: &[u8]) -> Cow<'a, [u8]> { - if self.rename { - let mut header = counter.to_string().into_bytes(); - if !read_suffix.is_empty() { - header.push(b' '); - header.extend_from_slice(read_suffix); - } - Cow::Owned(header) - } else { - Cow::Borrowed(id) - } - } - - /// Push a single record into the thread-local CBQ writer - fn push_cbq( - &mut self, - record: &Rf, - seq: &[u8], - read_suffix: &[u8], - ) -> Result<()> { - let counter = if self.rename { - self.rename_counter.fetch_add(1, Ordering::Relaxed) + 1 - } else { - 0 - }; - let header = self.cbq_header(record.id(), counter, 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], - ) -> Result<()> { - let counter = if self.rename { - self.rename_counter.fetch_add(1, Ordering::Relaxed) + 1 - } else { - 0 - }; - let header1 = self.cbq_header(record1.id(), counter, b"/1"); - let header2 = self.cbq_header(record2.id(), counter, 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(); @@ -962,19 +1030,18 @@ impl FilterProcessor { } } - /// Shared per-record logic for every reader (paraseq FASTX, CBQ mmap, ...) - fn handle_record(&mut self, record: &Rf) -> Result<()> { - let seq = record.seq(); + /// Shared per-read logic for every reader (paraseq FASTX, CBQ mmap) + fn handle_read(&mut self, read: &ReadView) -> Result<()> { self.local_stats.total_seqs += 1; - self.local_stats.total_bp += seq.len() as u64; + self.local_stats.total_bp += read.seq.len() as u64; - let decision = self.should_keep_sequence(&seq); + let decision = self.should_keep_sequence(read.seq); // Show debug info for sequences with hits if self.debug { eprintln!( "DEBUG: {} hits={}/{} keep={} kmers=[{}]", - String::from_utf8_lossy(record.id()), + String::from_utf8_lossy(read.id), decision.hit_count, decision.total_minimizers, decision.keep, @@ -983,43 +1050,38 @@ impl FilterProcessor { } if decision.keep { - self.local_stats.output_bp += seq.len() as u64; - self.write_record(record, &seq, b"")?; - if self.rename { - self.batch_kept += 1; - } + self.local_stats.output_bp += read.seq.len() as u64; + self.output + .push_read(read, self.rename, self.output_fasta, &self.rename_counter)?; } else { self.local_stats.filtered_seqs += 1; - self.local_stats.filtered_bp += seq.len() as u64; + self.local_stats.filtered_bp += read.seq.len() as u64; } Ok(()) } - /// 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()) { + /// Shared per-pair logic for every reader (paraseq FASTX, CBQ mmap) + fn handle_pair(&mut self, read1: &ReadView, read2: &ReadView) -> Result<()> { + if self.check_pairs && !paired_record_names_match(read1.id, read2.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()) + String::from_utf8_lossy(read1.id), + String::from_utf8_lossy(read2.id) )); } - let seq1 = record1.seq(); - let seq2 = record2.seq(); - self.local_stats.total_seqs += 2; - self.local_stats.total_bp += (seq1.len() + seq2.len()) as u64; + self.local_stats.total_bp += (read1.seq.len() + read2.seq.len()) as u64; - let decision = self.should_keep_pair(&seq1, &seq2); + let decision = self.should_keep_pair(read1.seq, read2.seq); // Debug info for interleaved pairs if self.debug && decision.hit_count > 0 { eprintln!( "DEBUG: {}/{} hits={}/{} keep={} kmers=[{}]", - String::from_utf8_lossy(record1.id()), - String::from_utf8_lossy(record2.id()), + String::from_utf8_lossy(read1.id), + String::from_utf8_lossy(read2.id), decision.hit_count, decision.total_minimizers, decision.keep, @@ -1028,94 +1090,25 @@ impl FilterProcessor { } if decision.keep { - self.local_stats.output_bp += (seq1.len() + seq2.len()) as u64; - // Both mates take the pair's number - if self.global_writer.is_cbq() { - // Native paired record: both mates in one CBQ record - self.push_cbq_pair(record1, &seq1, record2, &seq2)?; - } else if self.global_writer2.is_some() { - // Separate outputs - self.write_record(record1, &seq1, b"/1")?; - self.write_record_to_buffer2(record2, &seq2, b"/2")?; - } else { - // Interleaved output - self.write_record(record1, &seq1, b"/1")?; - self.write_record(record2, &seq2, b"/2")?; - } - if self.rename { - self.batch_kept += 1; - } + self.local_stats.output_bp += (read1.seq.len() + read2.seq.len()) as u64; + self.output.push_pair( + read1, + read2, + self.rename, + self.output_fasta, + &self.rename_counter, + )?; } else { self.local_stats.filtered_seqs += 2; - self.local_stats.filtered_bp += (seq1.len() + seq2.len()) as u64; + self.local_stats.filtered_bp += (read1.seq.len() + read2.seq.len()) as u64; } Ok(()) } - /// 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())?; - } - 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(); - - if self.rename { - // Both mates share one block of numbers - let base = self.reserve_rename_ids(); - write_renamed( - &mut writer1, - self.local_buffer.fastx(), - &self.pending_renames, - base, - )?; - write_renamed( - &mut writer2, - &self.local_buffer2, - &self.pending_renames2, - base, - )?; - } else { - writer1.write_all(self.local_buffer.fastx())?; - writer2.write_all(&self.local_buffer2)?; - } - writer1.flush()?; - writer2.flush()?; - } - } else { - // Interleaved output - if !self.local_buffer.fastx().is_empty() { - let mut writer = global.lock(); - if self.rename { - let base = self.reserve_rename_ids(); - write_renamed( - &mut writer, - self.local_buffer.fastx(), - &self.pending_renames, - base, - )?; - } else { - writer.write_all(self.local_buffer.fastx())?; - } - writer.flush()?; - } - } - } - } - - if let LocalOutput::Fastx(v) = &mut self.local_buffer { - v.clear(); - } - self.local_buffer2.clear(); - self.pending_renames.clear(); - self.pending_renames2.clear(); - self.batch_kept = 0; + self.output + .flush_batch(self.rename, &self.rename_counter, self.ordered)?; // Update global stats { @@ -1136,12 +1129,8 @@ impl FilterProcessor { 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(()) + self.output.flush_thread() } } @@ -1151,7 +1140,13 @@ impl ParallelProcessor for FilterProcessor { } fn process_record(&mut self, record: Rf) -> paraseq::parallel::Result<()> { - self.handle_record(&record)?; + let seq = record.seq(); + self.handle_read(&ReadView { + id: record.id(), + seq: &seq, + qual: record.qual(), + flag: None, + })?; Ok(()) } @@ -1172,7 +1167,22 @@ impl PairedParallelProcessor for FilterProcessor { } fn process_record_pair(&mut self, record1: Rf, record2: Rf) -> paraseq::parallel::Result<()> { - self.handle_record_pair(&record1, &record2)?; + let seq1 = record1.seq(); + let seq2 = record2.seq(); + self.handle_pair( + &ReadView { + id: record1.id(), + seq: &seq1, + qual: record1.qual(), + flag: None, + }, + &ReadView { + id: record2.id(), + seq: &seq2, + qual: record2.qual(), + flag: None, + }, + )?; Ok(()) } @@ -1189,22 +1199,22 @@ impl PairedParallelProcessor for FilterProcessor { impl binseq::ParallelProcessor for FilterProcessor { fn process_record(&mut self, record: R) -> binseq::Result<()> { - let record1 = CbqRecord { + let read1 = ReadView { id: record.sheader(), seq: record.sseq(), qual: record.has_quality().then(|| record.squal()), - index: record.index(), + flag: record.flag(), }; if record.is_paired() { - let record2 = CbqRecord { + let read2 = ReadView { id: record.xheader(), seq: record.xseq(), qual: record.has_quality().then(|| record.xqual()), - index: record.index(), + flag: record.flag(), }; - self.handle_record_pair(&record1, &record2)?; + self.handle_pair(&read1, &read2)?; } else { - self.handle_record(&record1)?; + self.handle_read(&read1)?; } Ok(()) } @@ -1225,10 +1235,16 @@ pub fn run(config: &FilterConfig) -> Result { if let Some(threshold) = config.complexity_threshold { validate_unit_interval("complexity threshold", threshold)?; } - validate_check_pairs_mode( - config.check_pairs, - config.interleaved || config.input2_path.is_some(), - )?; + + // Fail --check-pairs before the index loads; CBQ may be natively paired, + // so skip it there and let validate_input_output recheck once the header is read + if config.check_pairs + && !config.interleaved + && config.input2_path.is_none() + && !config.input_path.ends_with(".cbq") + { + validate_check_pairs_mode(true, false)?; + } // Validate the index path once here; run_with_index never touches it again. if !config.minimizers_path.exists() { @@ -1258,6 +1274,7 @@ pub fn run(config: &FilterConfig) -> Result { ordered: config.ordered, threads: config.threads, compression_level: config.compression_level, + cbq_block_size: config.cbq_block_size, compression_threads: config.compression_threads, debug: config.debug, quiet: config.quiet, @@ -1316,10 +1333,6 @@ pub fn run_with_index( config: &FilterRunConfig, ) -> Result { validate_unit_interval("relative threshold", config.rel_threshold)?; - validate_check_pairs_mode( - config.check_pairs, - config.interleaved || config.input2_path.is_some(), - )?; let start_time = Instant::now(); let version: String = env!("CARGO_PKG_VERSION").to_string(); @@ -1341,7 +1354,7 @@ pub fn run_with_index( let compressed_output_count = count_compressed_outputs(config); // Allocate threads between filtering (rayon) and compression (gzp). - // Rayon pool can only be initialized once, so calculate before build_global(). + // Rayon pool can only be initialised once, so calculate before build_global(). let (filtering_threads, compression_threads_per_output) = if compressed_output_count > 0 { let compression_threads_total = if config.compression_threads > 0 { config.compression_threads as usize @@ -1359,11 +1372,11 @@ pub fn run_with_index( }; if filtering_threads > 0 { - // error is OK here when we initialize a 2nd time in server mode. + // error is OK here when we initialise a 2nd time in server mode. let _ = rayon::ThreadPoolBuilder::new() .num_threads(filtering_threads) .build_global() - .context("Failed to initialize thread pool"); + .context("Failed to initialise thread pool"); } check_input_paths(config)?; @@ -1371,10 +1384,23 @@ pub fn run_with_index( // 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)?; + let output_format = resolve_output_format(config); + let (layout, input) = open_input(config, interleaved_input)?; validate_input_output(&layout, output_format, config)?; + // The binseq reader and CBQ rename numbering cannot honour --ordered + // multithreaded + let ordered_cbq = config.ordered + && (layout.format == Format::Cbq || (output_format == Format::Cbq && config.rename)); + let filtering_threads = if ordered_cbq && filtering_threads > 1 { + if !quiet { + eprintln!("Using 1 filtering thread: --ordered with CBQ input or renamed CBQ output"); + } + 1 + } else { + filtering_threads + }; + let mode = if config.deplete { "deplete" } else { "search" }; let mut input_type = String::new(); @@ -1426,22 +1452,27 @@ pub fn run_with_index( ); } - let (global_writer, writer2) = match output_format { - OutputFormat::Cbq => { + let output = match output_format { + Format::Cbq => { + let path = config + .output_path + .as_deref() + .expect("CBQ output implies a named path"); let cbq_writer = BinseqWriterBuilder::new(BinseqFormat::Cbq) .paired(layout.paired) .quality(layout.qualities && !config.output_fasta) .headers(layout.headers || config.rename) - .flags(false) - .build(get_writer( - config.output_path.as_deref(), - config.compression_level, - compression_threads_per_output, - )?) + .flags(layout.flags) + .block_size( + (config.cbq_block_size as usize * 1024 * 1024) + .max(layout.block_size.unwrap_or(0)), + ) + .compression_level(i32::from(config.compression_level)) + .build(open_cbq_output(path)?) .context("Failed to create CBQ writer")?; - (SharedOutput::Cbq(Arc::new(Mutex::new(cbq_writer))), None) + Output::cbq(cbq_writer)? } - OutputFormat::Fastx => { + Format::Fastx => { let writer = get_writer( config.output_path.as_deref(), config.compression_level, @@ -1460,7 +1491,7 @@ pub fn run_with_index( } else { None }; - (SharedOutput::Fastx(Arc::new(Mutex::new(writer))), writer2) + Output::fastx(writer, writer2) } }; @@ -1497,8 +1528,7 @@ pub fn run_with_index( kmer_length, window_size, &processor_config, - global_writer, - writer2, + output, spinner.clone(), filtering_start_time, )?; @@ -1506,20 +1536,20 @@ pub fn run_with_index( // Process based on input type - use filtering threads (already calculated above) let num_threads = filtering_threads; - match prepared { - InputPrep::Cbq(reader) => { + match input { + Input::Cbq(reader) => { reader.process_parallel(processor.clone(), num_threads)?; } - InputPrep::FastxSingle(reader) => { + Input::FastxSingle(reader) => { reader.process_parallel(&mut processor, num_threads)?; } - InputPrep::FastxInterleaved(reader) => { + Input::FastxInterleaved(reader) => { reader.process_parallel_interleaved(&mut processor, num_threads)?; } - InputPrep::FastxPaired(reader1, reader2) => { + Input::FastxPaired(reader1, reader2) => { reader1.process_parallel_paired(reader2, &mut processor, num_threads)?; } - InputPrep::Empty => { + Input::Empty => { if !quiet { eprintln!("Empty input file(s) detected"); } @@ -1535,16 +1565,10 @@ 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 { - drop(w2); - } + // Finish any CBQ stream (writes the embedded index), then drop the + // processor so the writers flush + processor.output.finish()?; + drop(processor); let total_time = start_time.elapsed(); let filtering_time = filtering_start_time.elapsed(); diff --git a/src/lib.rs b/src/lib.rs index dc12d4f..e17a13f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -16,7 +16,9 @@ mod minimizers; // Public API #[cfg(feature = "cli")] -pub use filter::{FilterRunConfig, FilterSummary, run as run_filter, run_with_index}; +pub use filter::{ + DEFAULT_CBQ_BLOCK_SIZE_MIB, FilterRunConfig, FilterSummary, run as run_filter, run_with_index, +}; pub use filter_kernel::{FilterDecision, FilterKernel, FilterParams}; #[cfg(feature = "fetch")] pub use index::fetch as index_fetch; @@ -331,7 +333,7 @@ pub struct FilterConfig<'a> { /// Path to JSON summary file pub summary_path: Option<&'a PathBuf>, - /// Deplete mode (remove sequences WITH matches, original deacon behavior) + /// Deplete mode (remove sequences WITH matches, original deacon behaviour) pub deplete: bool, /// Replace sequence headers with incrementing numbers (1, 2, 3...) @@ -349,6 +351,9 @@ pub struct FilterConfig<'a> { /// Compression level for output files (1-22 for zst, 1-9 for gz) pub compression_level: u8, + /// cbq output block size in MiB (raised to the cbq input's block size if larger) + pub cbq_block_size: u16, + /// Number of threads for compression (0 = auto-calculate as ceil(total/2)) pub compression_threads: u16, diff --git a/src/main.rs b/src/main.rs index cdb2c94..76cb65f 100644 --- a/src/main.rs +++ b/src/main.rs @@ -3,8 +3,9 @@ use clap::{Parser, Subcommand}; #[cfg(feature = "fetch")] use deacon::index_fetch; use deacon::{ - ComplexityAlgorithm, DEFAULT_KMER_LENGTH, DEFAULT_WINDOW_SIZE, FilterConfig, IndexConfig, - index_diff, index_dump, index_filter, index_freeze, index_info, index_intersect, index_union, + ComplexityAlgorithm, DEFAULT_CBQ_BLOCK_SIZE_MIB, DEFAULT_KMER_LENGTH, DEFAULT_WINDOW_SIZE, + FilterConfig, IndexConfig, index_diff, index_dump, index_filter, index_freeze, index_info, + index_intersect, index_union, }; use serde::{Deserialize, Serialize}; use std::io::{Read, Write}; @@ -29,16 +30,16 @@ enum Command { #[command(subcommand)] command: IndexCommand, }, - /// Retain or deplete sequence records with sufficient minimizer hits to an indexed query + /// Retain or deplete sequence records with sufficient minimizer hits to the index Filter { /// Path to minimizer index file index: PathBuf, - /// Optional path to fastx or CBQ file (or - for stdin) + /// Optional path to fastx or binseq cbq file (or - for stdin) #[arg(default_value = "-")] input: String, - /// Optional path to second paired fastx file (not supported with CBQ input) + /// Optional path to second paired fastx file input2: Option, /// Minimum absolute number of minimizer hits for a match @@ -69,11 +70,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; .cbq suffix writes CBQ) + /// Path to output file (stdout by default; detects fastx with .gz, .zst, .xz, or binseq with .cbq) #[arg(short = 'o', long = "output")] output: Option, - /// Optional path to second paired output fastx file (detects .gz and .zst; not supported with CBQ output) + /// Optional path to second paired output fastx file (detects .gz, .zst, .xz) #[arg(short = 'O', long = "output2")] output2: Option, @@ -89,11 +90,19 @@ enum Command { #[arg(long = "compression-threads", default_value_t = 0)] compression_threads: u16, - /// Output compression level (1-9 for gz & xz; 1-22 for zstd) + /// Output compression level (1-9 for gz & xz; 1-22 for zstd, including cbq) #[arg(long = "compression-level", default_value_t = 2)] compression_level: u8, - /// Treat INPUT as interleaved paired reads from a file or stdin + /// cbq output block size in MiB (or cbq input block size if higher) + #[arg( + long = "cbq-block-size", + default_value_t = DEFAULT_CBQ_BLOCK_SIZE_MIB, + value_parser = clap::value_parser!(u16).range(1..=1024) + )] + cbq_block_size: u16, + + /// Treat INPUT as interleaved paired records from single file or stdin #[arg( long = "interleaved", default_value_t = false, @@ -109,13 +118,13 @@ enum Command { #[arg(long = "check-pairs", default_value_t = false)] check_pairs: bool, - /// Suppress progress reporting - #[arg(short = 'q', long = "quiet", default_value_t = false)] - quiet: bool, - /// Output sequences with minimizer hits to stderr #[arg(long = "debug", default_value_t = false)] debug: bool, + + /// Suppress progress reporting + #[arg(short = 'q', long = "quiet", default_value_t = false)] + quiet: bool, }, /// Start/stop a server process for reduced latency filtering Server { @@ -344,7 +353,7 @@ fn main() -> Result<()> { rayon::ThreadPoolBuilder::new() .num_threads(*threads as usize) .build_global() - .context("Failed to initialize thread pool")?; + .context("Failed to initialise thread pool")?; // Remove existing socket if present let _ = std::fs::remove_file("deacon_server_socket"); @@ -533,6 +542,7 @@ fn process_command(command: &Command) -> Result<(), anyhow::Error> { output_fasta, threads, compression_level, + cbq_block_size, compression_threads, ordered, check_pairs, @@ -565,6 +575,7 @@ fn process_command(command: &Command) -> Result<(), anyhow::Error> { ordered: *ordered, threads: *threads, compression_level: *compression_level, + cbq_block_size: *cbq_block_size, compression_threads: *compression_threads, debug: *debug, quiet: *quiet, diff --git a/tests/filter_tests.rs b/tests/filter_tests.rs index 2fe1c2d..afec2f2 100644 --- a/tests/filter_tests.rs +++ b/tests/filter_tests.rs @@ -177,7 +177,7 @@ fn test_filter_to_file() { assert!(output_path.exists(), "Output file wasn't created"); assert!(summary_path.exists(), "Summary file wasn't created"); - // With new default behavior: sequences without matches are filtered out (sequences too short for k=31) + // With new default behaviour: sequences without matches are filtered out (sequences too short for k=31) let output_content = fs::read_to_string(&output_path).unwrap(); assert!( output_content.is_empty(), @@ -1113,7 +1113,7 @@ fn test_shared_minimizer_counted_once() { // If shared minimizers are counted once (correct): total hits = 1, pair kept (1 < 2) // If shared minimizers are counted twice (bug): total hits = 2+, pair filtered (2+ >= 2) - // Using --deplete to restore original behavior for this bug test + // Using --deplete to restore original behaviour for this bug test let mut cmd = cargo::cargo_bin_cmd!("deacon"); cmd.arg("filter") .arg("--deplete") @@ -2786,14 +2786,17 @@ 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; + use binseq::write::{BinseqWriterBuilder, Format}; 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()] { + for seq in [ + b"ACGTACGTACGTACGTACGT".as_slice(), + b"TGCAACGTACGTACGTACGT".as_slice(), + ] { writer .push( SequencingRecordBuilder::default() @@ -2807,6 +2810,91 @@ fn write_headerless_paired_cbq(path: &Path) { writer.finish().unwrap(); } +/// CBQ output block size is --cbq-block-size raised to the CBQ input's own block +/// size, so ultra-long records always round-trip; FASTX has none to inherit. +#[test] +fn cbq_output_block_size() { + use binseq::SequencingRecordBuilder; + use binseq::write::{BinseqWriterBuilder, Format}; + + 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"); + create_test_fasta(&fasta_path); + create_test_fastq(&fastq_path); + build_index(&fasta_path, &bin_path); + + let input = temp_dir.path().join("in.cbq"); + let mut writer = BinseqWriterBuilder::new(Format::Cbq) + .block_size(32 << 20) + .build(File::create(&input).unwrap()) + .unwrap(); + writer + .push( + SequencingRecordBuilder::default() + .s_seq(b"ACGTACGTACGTACGTACGT") + .build() + .unwrap(), + ) + .unwrap(); + writer.finish().unwrap(); + + let cases = [ + (&fastq_path, None, 16 << 20), // default + (&fastq_path, Some("64"), 64 << 20), // flag honoured + (&input, None, 32 << 20), // input's block size wins over the default + (&input, Some("4"), 32 << 20), // clamped up to the input's + (&input, Some("64"), 64 << 20), // flag wins when larger + ]; + + for (case, (source, block_size, expected)) in cases.iter().enumerate() { + let output = temp_dir.path().join(format!("out{case}.cbq")); + let mut cmd = cargo::cargo_bin_cmd!("deacon"); + cmd.args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(source) + .arg("--output") + .arg(&output); + if let Some(block_size) = block_size { + cmd.args(["--cbq-block-size", block_size]); + } + cmd.assert().success(); + let reader = binseq::cbq::MmapReader::new(&output).unwrap(); + assert_eq!( + reader.header().block_size as usize, + *expected, + "case {case}" + ); + } +} + +/// --cbq-block-size is bounded by clap, before any work starts +#[test] +fn cbq_block_size_out_of_range_rejected() { + 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"); + create_test_fasta(&fasta_path); + create_test_fastq(&fastq_path); + build_index(&fasta_path, &bin_path); + + for block_size in ["0", "2000"] { + let output = temp_dir.path().join(format!("no_create{block_size}.cbq")); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "--cbq-block-size", block_size]) + .arg(&bin_path) + .arg(&fastq_path) + .arg("--output") + .arg(&output) + .assert() + .code(2) + .stderr(predicate::str::contains("not in 1..=1024")); + assert!(!output.exists()); + } +} + /// CBQ is only an I/O concern: FASTX -> CBQ -> FASTX must preserve reads, /// pairing, headers, qualities, and summary counts exactly. #[test] @@ -3007,7 +3095,10 @@ fn cbq_roundtrip_matches_fastx() { .arg(&no_create) .assert() .failure(); - assert!(!no_create.exists(), "CBQ input + INPUT2 must fail before output creation"); + assert!( + !no_create.exists(), + "CBQ input + INPUT2 must fail before output creation" + ); cargo::cargo_bin_cmd!("deacon") .args(["filter", "--interleaved", "-t", "1"]) @@ -3017,7 +3108,10 @@ fn cbq_roundtrip_matches_fastx() { .arg(&no_create) .assert() .failure(); - assert!(!no_create.exists(), "CBQ input + --interleaved must fail before output creation"); + assert!( + !no_create.exists(), + "CBQ input + --interleaved must fail before output creation" + ); cargo::cargo_bin_cmd!("deacon") .args(["filter", "-t", "1"]) @@ -3029,32 +3123,429 @@ fn cbq_roundtrip_matches_fastx() { .arg(temp_dir.path().join("x.fastq")) // CBQ output + OUTPUT2 .assert() .failure(); - assert!(!no_create.exists(), "CBQ output + OUTPUT2 must fail before output creation"); + assert!( + !no_create.exists(), + "CBQ output + OUTPUT2 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" + ); +} + +/// Write `n` distinct single-end FASTQ records that do not match the AAA index +fn create_many_fastq(path: &Path, n: usize) -> String { + let mut content = String::new(); + let seq = "CGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT"; + let qual = "~".repeat(seq.len()); + for i in 0..n { + content.push_str(&format!("@read_{i}\n{seq}\n+\n{qual}\n")); + } + fs::write(path, &content).unwrap(); + content +} + +/// --check-pairs accepts a natively paired CBQ with /1 /2 names +#[test] +fn cbq_check_pairs_paired_input() { + let temp_dir = tempdir().unwrap(); + let fasta_path = temp_dir.path().join("ref.fasta"); + let bin_path = temp_dir.path().join("ref.bin"); + let interleaved_path = temp_dir.path().join("interleaved.fastq"); + create_test_fasta(&fasta_path); + create_test_interleaved_fastq(&interleaved_path); + build_index(&fasta_path, &bin_path); + + // Interleaved FASTQ with read1/1, read1/2 names -> paired CBQ + let paired_cbq = temp_dir.path().join("paired.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1", "--interleaved"]) + .arg(&bin_path) + .arg(&interleaved_path) + .arg("--output") + .arg(&paired_cbq) + .assert() + .success(); + + // Paired CBQ + --check-pairs succeeds despite no INPUT2/--interleaved + let out = temp_dir.path().join("out.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1", "--check-pairs"]) + .arg(&bin_path) + .arg(&paired_cbq) + .arg("--output") + .arg(&out) + .assert() + .success(); + assert_eq!(count_records(&fs::read_to_string(&out).unwrap()), 4); +} + +/// A file bearing the CBQ magic but truncated must error, not panic +#[test] +fn cbq_truncated_input_fails_cleanly() { + 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"); + create_test_fasta(&fasta_path); + create_test_fastq(&fastq_path); + build_index(&fasta_path, &bin_path); + + // Valid CBQ, then truncate it to just past its magic + let cbq_path = temp_dir.path().join("reads.cbq"); + 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) + .assert() + .success(); + let bytes = fs::read(&cbq_path).unwrap(); + let truncated = temp_dir.path().join("truncated.cbq"); + fs::write(&truncated, &bytes[..16]).unwrap(); - let gz_cbq = temp_dir.path().join("bad.cbq.gz"); cargo::cargo_bin_cmd!("deacon") .args(["filter", "-t", "1"]) .arg(&bin_path) + .arg(&truncated) + .assert() + .failure() + .stderr(predicate::str::contains("Truncated or corrupt CBQ")); +} + +/// CBQ per-record flags survive CBQ -> CBQ filtering +#[test] +fn cbq_flags_preserved() { + use binseq::write::{BinseqWriterBuilder, Format}; + use binseq::{ParallelReader, SequencingRecordBuilder}; + use std::sync::{Arc, Mutex}; + + let temp_dir = tempdir().unwrap(); + 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); + + // CBQ input with flags 7 and 9 + let flagged = temp_dir.path().join("flagged.cbq"); + let mut writer = BinseqWriterBuilder::new(Format::Cbq) + .headers(true) + .flags(true) + .build(File::create(&flagged).unwrap()) + .unwrap(); + for (i, flag) in [7u64, 9u64].iter().enumerate() { + writer + .push( + SequencingRecordBuilder::default() + .s_seq(b"CGTACGTACGTACGTACGTACGTACGTACGTACGTACGT".as_slice()) + .s_header(format!("read_{i}").as_bytes()) + .opt_flag(Some(*flag)) + .build() + .unwrap(), + ) + .unwrap(); + } + writer.finish().unwrap(); + + // Deplete against the AAA index keeps everything + let out_cbq = temp_dir.path().join("out.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "-t", "1"]) + .arg(&aaa_bin) + .arg(&flagged) + .arg("--output") + .arg(&out_cbq) + .assert() + .success(); + + // Collect the flags back out of the filtered CBQ + #[derive(Clone)] + struct FlagCollector(Arc>>>); + impl binseq::ParallelProcessor for FlagCollector { + fn process_record(&mut self, record: R) -> binseq::Result<()> { + self.0.lock().unwrap().push(record.flag()); + Ok(()) + } + } + let reader = binseq::cbq::MmapReader::new(&out_cbq).unwrap(); + assert!(reader.header().has_flags(), "flags must survive filtering"); + let collector = FlagCollector(Arc::new(Mutex::new(Vec::new()))); + reader.process_parallel(collector.clone(), 1).unwrap(); + assert_eq!(*collector.0.lock().unwrap(), vec![Some(7), Some(9)]); +} + +/// --ordered CBQ output matches the input order exactly, multithreaded +#[test] +fn cbq_ordered_output_is_input_ordered() { + let temp_dir = tempdir().unwrap(); + 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 fastq_path = temp_dir.path().join("many.fastq"); + let content = create_many_fastq(&fastq_path, 500); + + // Deplete against the AAA index keeps all 500 reads + let cbq_path = temp_dir.path().join("ordered.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "--ordered", "-t", "4"]) + .arg(&aaa_bin) + .arg(&fastq_path) + .arg("--output") + .arg(&cbq_path) + .assert() + .success(); + + // Decoding the CBQ reproduces the input byte-for-byte (CBQ input with + // --ordered falls back to one thread) + let out = temp_dir.path().join("out.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "--ordered", "-t", "4"]) + .arg(&aaa_bin) + .arg(&cbq_path) + .arg("--output") + .arg(&out) + .assert() + .success(); + assert_eq!( + fs::read_to_string(&out).unwrap(), + content, + "ordered CBQ roundtrip must preserve input order" + ); +} + +/// --ordered --rename CBQ output numbers records 1..N in input order +#[test] +fn cbq_ordered_rename_is_sequential() { + let temp_dir = tempdir().unwrap(); + 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 fastq_path = temp_dir.path().join("many.fastq"); + create_many_fastq(&fastq_path, 200); + + let cbq_path = temp_dir.path().join("renamed.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "--ordered", "--rename", "-t", "4"]) + .arg(&aaa_bin) + .arg(&fastq_path) + .arg("--output") + .arg(&cbq_path) + .assert() + .success(); + + let out = temp_dir.path().join("out.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "-t", "1"]) + .arg(&aaa_bin) + .arg(&cbq_path) + .arg("--output") + .arg(&out) + .assert() + .success(); + let content = fs::read_to_string(&out).unwrap(); + let ids: Vec<&str> = content + .lines() + .step_by(4) + .map(|h| h.trim_start_matches('@')) + .collect(); + let expected: Vec = (1..=200).map(|i| i.to_string()).collect(); + assert_eq!( + ids, expected, + "renamed CBQ records must be numbered in order" + ); +} + +/// Multithreaded CBQ roundtrips preserve record counts +#[test] +fn cbq_multithreaded_roundtrip_counts() { + let temp_dir = tempdir().unwrap(); + 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 fastq_path = temp_dir.path().join("many.fastq"); + create_many_fastq(&fastq_path, 500); + + let cbq_path = temp_dir.path().join("many.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "-t", "4"]) + .arg(&aaa_bin) .arg(&fastq_path) .arg("--output") - .arg(&gz_cbq) // compressed CBQ output + .arg(&cbq_path) + .assert() + .success(); + let reader = binseq::cbq::MmapReader::new(&cbq_path).unwrap(); + assert_eq!(reader.num_records(), 500); + + let out = temp_dir.path().join("out.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "-t", "4"]) + .arg(&aaa_bin) + .arg(&cbq_path) + .arg("--output") + .arg(&out) + .assert() + .success(); + assert_eq!(count_records(&fs::read_to_string(&out).unwrap()), 500); +} + +/// OUTPUT2 must not silently write FASTX bytes into a .cbq-named file +#[test] +fn cbq_output2_suffix_rejected() { + let temp_dir = tempdir().unwrap(); + let fasta_path = temp_dir.path().join("ref.fasta"); + let bin_path = temp_dir.path().join("ref.bin"); + 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_paired_fastq(&r1_path, &r2_path); + build_index(&fasta_path, &bin_path); + + let out1 = temp_dir.path().join("out1.fastq"); + let out2 = temp_dir.path().join("out2.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(&out1) + .arg("--output2") + .arg(&out2) .assert() .failure(); - assert!(!gz_cbq.exists(), ".cbq.gz output must fail before output creation"); + assert!(!out1.exists(), "OUTPUT must not be created"); + assert!(!out2.exists(), "OUTPUT2 must not be created"); +} - // --check-pairs + headerless CBQ input +/// Headerless, quality-free paired CBQ filters to CBQ (still headerless) and +/// to FASTX (record-index ids, FASTA since there are no qualities) +#[test] +fn cbq_headerless_roundtrip() { + let temp_dir = tempdir().unwrap(); + 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 headerless = temp_dir.path().join("headerless.cbq"); write_headerless_paired_cbq(&headerless); + + // Deplete against the AAA index keeps both pairs + let out_cbq = temp_dir.path().join("out.cbq"); cargo::cargo_bin_cmd!("deacon") - .args(["filter", "--check-pairs", "-t", "1"]) - .arg(&bin_path) + .args(["filter", "-d", "-t", "1"]) + .arg(&aaa_bin) .arg(&headerless) .arg("--output") - .arg(&no_create) + .arg(&out_cbq) .assert() - .failure(); + .success(); + let reader = binseq::cbq::MmapReader::new(&out_cbq).unwrap(); + assert!(reader.is_paired()); assert!( - !no_create.exists(), - "--check-pairs + headerless CBQ must fail before output creation" + !reader.header().has_headers(), + "output must stay headerless" + ); + assert!(!reader.header().has_qualities()); + assert_eq!(reader.num_records(), 2); + + // Quality-free CBQ decodes to FASTA + let out_fasta = temp_dir.path().join("out.fasta"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "-t", "1"]) + .arg(&aaa_bin) + .arg(&headerless) + .arg("--output") + .arg(&out_fasta) + .assert() + .success(); + let content = fs::read_to_string(&out_fasta).unwrap(); + assert_eq!( + content.matches('>').count(), + 4, + "two pairs, four FASTA records" + ); +} + +/// Ultra-long reads exceed binseq's 1 MiB default block; deacon's larger CBQ +/// block must roundtrip them +#[test] +fn cbq_ultralong_read_roundtrip() { + let temp_dir = tempdir().unwrap(); + 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); + + // 2 Mb read: ~2.5 MB embedded size, over the 1 MiB binseq default + let fastq_path = temp_dir.path().join("ultra.fastq"); + let seq: String = "ACGTGCTA".repeat(250_000); + let qual = "~".repeat(seq.len()); + let content = format!("@ultra1\n{seq}\n+\n{qual}\n"); + fs::write(&fastq_path, &content).unwrap(); + + let cbq_path = temp_dir.path().join("ultra.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "-t", "1"]) + .arg(&aaa_bin) + .arg(&fastq_path) + .arg("--output") + .arg(&cbq_path) + .assert() + .success(); + + let out = temp_dir.path().join("ultra_out.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "-t", "1"]) + .arg(&aaa_bin) + .arg(&cbq_path) + .arg("--output") + .arg(&out) + .assert() + .success(); + assert_eq!( + fs::read_to_string(&out).unwrap(), + content, + "ultra-long read must roundtrip through CBQ" ); } + +/// CBQ output validates the compression level before creating the output +#[test] +fn cbq_invalid_compression_level_rejected() { + 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"); + create_test_fasta(&fasta_path); + create_test_fastq(&fastq_path); + build_index(&fasta_path, &bin_path); + + let no_create = temp_dir.path().join("bad_level.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-t", "1", "--compression-level", "23"]) + .arg(&bin_path) + .arg(&fastq_path) + .arg("--output") + .arg(&no_create) + .assert() + .failure() + .stderr(predicate::str::contains("Invalid CBQ compression level")); + assert!(!no_create.exists()); +} From 73232d1ee58d1d1f05eb71efbd90eea9a1f29e1b Mon Sep 17 00:00:00 2001 From: Bede Constantinides Date: Wed, 12 Aug 2026 16:34:54 +0100 Subject: [PATCH 03/12] Test cbq EF N encoding with fastx to cbq to fastx roundtrips --- tests/filter_tests.rs | 63 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) diff --git a/tests/filter_tests.rs b/tests/filter_tests.rs index afec2f2..df8cddb 100644 --- a/tests/filter_tests.rs +++ b/tests/filter_tests.rs @@ -3145,6 +3145,69 @@ fn cbq_roundtrip_matches_fastx() { ); } +/// CBQ must preserve Ns and filter them like FASTX. +#[test] +fn cbq_ns_match_fastx() { + let temp_dir = tempdir().unwrap(); + let ref_path = temp_dir.path().join("ref.fasta"); + let bin_path = temp_dir.path().join("ref.bin"); + fs::write(&ref_path, format!(">ref\n{}\n", "A".repeat(100))).unwrap(); + build_index(&ref_path, &bin_path); + + let run = |input: &Path, output: &Path, deplete: bool, fasta: bool| { + let mut cmd = cargo::cargo_bin_cmd!("deacon"); + cmd.arg("filter"); + if deplete { + cmd.args(["-d", "-a", "999"]); + } else { + cmd.args(["-a", "1"]); + } + cmd.args(["-r", "0", "-t", "1"]); + if fasta { + cmd.arg("--fasta"); + } + cmd.arg(&bin_path) + .arg(input) + .arg("--output") + .arg(output) + .assert() + .success(); + }; + + let all_n = "N".repeat(64); + let split_n = format!("{}N{}", "A".repeat(64), "A".repeat(64)); + let fastq = format!( + "@all-n\n{all_n}\n+\n{}\n@split-n\n{split_n}\n+\n{}\n", + "I".repeat(all_n.len()), + "I".repeat(split_n.len()), + ); + let fasta = format!(">all-n\n{all_n}\n>split-n\n{split_n}\n"); + + for (ext, content, fasta_output, marker) in [ + ("fastq", fastq.as_str(), false, "@"), + ("fasta", fasta.as_str(), true, ">"), + ] { + let input = temp_dir.path().join(format!("reads.{ext}")); + let cbq = temp_dir.path().join(format!("reads-{ext}.cbq")); + let roundtrip = temp_dir.path().join(format!("roundtrip.{ext}")); + let direct = temp_dir.path().join(format!("direct-filtered.{ext}")); + let via_cbq = temp_dir.path().join(format!("cbq-filtered.{ext}")); + fs::write(&input, content).unwrap(); + + run(&input, &cbq, true, fasta_output); + run(&cbq, &roundtrip, true, fasta_output); + assert_eq!(fs::read(&roundtrip).unwrap(), content.as_bytes()); + + run(&input, &direct, false, fasta_output); + run(&cbq, &via_cbq, false, fasta_output); + let filtered = fs::read(&direct).unwrap(); + assert_eq!(filtered, fs::read(&via_cbq).unwrap()); + let filtered = String::from_utf8_lossy(&filtered); + assert!(!filtered.contains(&format!("{marker}all-n"))); + assert!(filtered.contains(&format!("{marker}split-n"))); + } +} + /// Write `n` distinct single-end FASTQ records that do not match the AAA index fn create_many_fastq(path: &Path, n: usize) -> String { let mut content = String::new(); From 98caa6a0010ba81574c7f32a6a2c802fb5158533 Mon Sep 17 00:00:00 2001 From: Bede Constantinides Date: Wed, 12 Aug 2026 17:12:46 +0100 Subject: [PATCH 04/12] Canonicalise non-N ambig YRWSKM to N before writing cbq, matching fastx behaviour --- src/filter.rs | 31 +++++++++++++++++-- tests/filter_tests.rs | 72 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 100 insertions(+), 3 deletions(-) diff --git a/src/filter.rs b/src/filter.rs index a98bd51..a86b0a4 100644 --- a/src/filter.rs +++ b/src/filter.rs @@ -157,8 +157,9 @@ impl Output { Output::Cbq { local, .. } => { let number = rename.then(|| rename_counter.fetch_add(1, Ordering::Relaxed) + 1); let header = cbq_header(read.id, number, b""); + let seq = sanitise_for_cbq(read.seq); let seq_record = SequencingRecordBuilder::default() - .s_seq(read.seq) + .s_seq(&seq) .s_header(&header) .opt_s_qual(read.qual) .opt_flag(read.flag) @@ -229,11 +230,13 @@ impl Output { let number = rename.then(|| rename_counter.fetch_add(1, Ordering::Relaxed) + 1); let header1 = cbq_header(read1.id, number, b"/1"); let header2 = cbq_header(read2.id, number, b"/2"); + let seq1 = sanitise_for_cbq(read1.seq); + let seq2 = sanitise_for_cbq(read2.seq); let seq_record = SequencingRecordBuilder::default() - .s_seq(read1.seq) + .s_seq(&seq1) .s_header(&header1) .opt_s_qual(read1.qual) - .x_seq(read2.seq) + .x_seq(&seq2) .x_header(&header2) .opt_x_qual(read2.qual) .opt_flag(read1.flag) @@ -343,6 +346,28 @@ fn cbq_header<'a>(id: &'a [u8], number: Option, suffix: &[u8]) -> Cow<'a, [ } } +/// Map unsupported bases to N. +fn sanitise_for_cbq(seq: &[u8]) -> Cow<'_, [u8]> { + #[inline] + fn supported(b: u8) -> bool { + matches!( + b, + b'A' | b'C' | b'G' | b'T' | b'N' | b'a' | b'c' | b'g' | b't' + ) + } + + let Some(first) = seq.iter().position(|&b| !supported(b)) else { + return Cow::Borrowed(seq); + }; + let mut sanitised = seq.to_vec(); + sanitised[first..].iter_mut().for_each(|b| { + if !supported(*b) { + *b = b'N'; + } + }); + Cow::Owned(sanitised) +} + /// Filtering config for an already-loaded index (no index path; see [`FilterConfig`]). pub struct FilterRunConfig { /// Path to input fastx file (or - for stdin) diff --git a/tests/filter_tests.rs b/tests/filter_tests.rs index df8cddb..838ae28 100644 --- a/tests/filter_tests.rs +++ b/tests/filter_tests.rs @@ -3208,6 +3208,78 @@ fn cbq_ns_match_fastx() { } } +/// CBQ encodes IUPAC codes as N across 32bp block boundaries. +#[test] +fn cbq_iupac_codes_become_n() { + let temp_dir = tempdir().unwrap(); + let ref_path = temp_dir.path().join("ref.fasta"); + let bin_path = temp_dir.path().join("ref.bin"); + fs::write(&ref_path, format!(">ref\n{}\n", "A".repeat(100))).unwrap(); + build_index(&ref_path, &bin_path); + + let mut seq = "ACGT".repeat(25).into_bytes(); + for (code, pos) in [ + ('Y', 0), + ('R', 15), + ('W', 31), + ('S', 32), + ('K', 33), + ('M', 63), + ('y', 1), + ('r', 16), + ('w', 30), + ('s', 34), + ('k', 62), + ('m', 64), + ('n', 99), + ] { + seq[pos] = code as u8; + } + let seq = String::from_utf8(seq).unwrap(); + let expected_seq: String = seq + .bytes() + .map(|b| match b { + b'A' | b'C' | b'G' | b'T' => b as char, + _ => 'N', + }) + .collect(); + + for (ext, input, expected, fasta) in [ + ( + "fasta", + format!(">iupac\n{seq}\n"), + format!(">iupac\n{expected_seq}\n"), + true, + ), + ( + "fastq", + format!("@iupac\n{seq}\n+\n{}\n", "I".repeat(seq.len())), + format!("@iupac\n{expected_seq}\n+\n{}\n", "I".repeat(seq.len())), + false, + ), + ] { + let input_path = temp_dir.path().join(format!("reads.{ext}")); + let cbq_path = temp_dir.path().join(format!("reads-{ext}.cbq")); + let output_path = temp_dir.path().join(format!("roundtrip.{ext}")); + fs::write(&input_path, input).unwrap(); + + for (from, to) in [(&input_path, &cbq_path), (&cbq_path, &output_path)] { + let mut cmd = cargo::cargo_bin_cmd!("deacon"); + cmd.args(["filter", "-d", "-a", "65535", "-r", "1", "-t", "1"]); + if fasta { + cmd.arg("--fasta"); + } + cmd.arg(&bin_path) + .arg(from) + .arg("--output") + .arg(to) + .assert() + .success(); + } + assert_eq!(fs::read_to_string(output_path).unwrap(), expected); + } +} + /// Write `n` distinct single-end FASTQ records that do not match the AAA index fn create_many_fastq(path: &Path, n: usize) -> String { let mut content = String::new(); From 42dad2332dced2d1823b58b769b7268722ce1884 Mon Sep 17 00:00:00 2001 From: Bede Constantinides Date: Mon, 17 Aug 2026 11:23:57 +0100 Subject: [PATCH 05/12] Drop non-ACGTN remapping for binseq io now that latest binseq does it by default --- Cargo.lock | 24 ++++++++++++++++++++---- Cargo.toml | 7 ++++--- src/filter.rs | 31 +++---------------------------- 3 files changed, 27 insertions(+), 35 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 3e3759c..00a8713 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -168,13 +168,13 @@ dependencies = [ [[package]] name = "binseq" -version = "0.9.5" +version = "0.9.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e468ca3b52411b87624a45d8bbdefd4e510c9121ae7da2a652fa6b90c90181dc" +checksum = "5fea1e47c3296e0073dc373f24660b1bd8ec36d6c06a7de115e4c75bf5b1753d" dependencies = [ - "anyhow", "auto_impl", - "bitnuc", + "bitnuc 0.4.1", + "bitnuc 0.5.2", "bytemuck", "byteorder", "itoa", @@ -205,6 +205,16 @@ version = "0.4.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d7871516cbb4e097e220623917e1491c995ee58c8018d929d4b1e76c378002bf" +[[package]] +name = "bitnuc" +version = "0.5.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "93612e75f32ab9cca9750b822f77c5dcf213782cfb67e12f6c8cf91679cb112a" +dependencies = [ + "fearless_simd", + "thiserror 2.0.20", +] + [[package]] name = "bstr" version = "1.13.1" @@ -635,6 +645,12 @@ version = "2.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "da7c62ceae207dd37ea5b845da6a0696c799f85e97da1ab5b7910be3c1c80223" +[[package]] +name = "fearless_simd" +version = "0.7.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f4beca3cb2444e3304ac30843cc091f44ed58932353cd492ce740067bfce6b12" + [[package]] name = "find-msvc-tools" version = "0.1.10" diff --git a/Cargo.toml b/Cargo.toml index 8720214..4a08c8b 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -41,8 +41,9 @@ indicatif = { version = "0.18", optional = true } serde_json = { version = "1.0", optional = true } niffler = { version = "3.0.1", default-features = false, optional = true } paraseq = { version = "0.5.0", 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.5", default-features = false, features = ["anyhow"], optional = true } +# CBQ support. Defaults disabled since we don't need binseq's own paraseq integration +# We need binseq >=0.9.6 for safe canonicalisation of non-N ambiguous bases to N +binseq = { version = "0.9.6", default-features = false, optional = true } ensure_simd = { version = "0.1.0", optional = true } # Compression deps (optional) @@ -61,7 +62,7 @@ predicates = "3.0" tempfile = "3.20" rstest = "0.26" nix = { version = "0.31", features = ["fs"] } -binseq = { version = "0.9.5", default-features = false, features = ["anyhow"] } +binseq = { version = "0.9.6", default-features = false } [features] scalar = ["simd-minimizers/scalar", "packed-seq/scalar"] diff --git a/src/filter.rs b/src/filter.rs index a86b0a4..a98bd51 100644 --- a/src/filter.rs +++ b/src/filter.rs @@ -157,9 +157,8 @@ impl Output { Output::Cbq { local, .. } => { let number = rename.then(|| rename_counter.fetch_add(1, Ordering::Relaxed) + 1); let header = cbq_header(read.id, number, b""); - let seq = sanitise_for_cbq(read.seq); let seq_record = SequencingRecordBuilder::default() - .s_seq(&seq) + .s_seq(read.seq) .s_header(&header) .opt_s_qual(read.qual) .opt_flag(read.flag) @@ -230,13 +229,11 @@ impl Output { let number = rename.then(|| rename_counter.fetch_add(1, Ordering::Relaxed) + 1); let header1 = cbq_header(read1.id, number, b"/1"); let header2 = cbq_header(read2.id, number, b"/2"); - let seq1 = sanitise_for_cbq(read1.seq); - let seq2 = sanitise_for_cbq(read2.seq); let seq_record = SequencingRecordBuilder::default() - .s_seq(&seq1) + .s_seq(read1.seq) .s_header(&header1) .opt_s_qual(read1.qual) - .x_seq(&seq2) + .x_seq(read2.seq) .x_header(&header2) .opt_x_qual(read2.qual) .opt_flag(read1.flag) @@ -346,28 +343,6 @@ fn cbq_header<'a>(id: &'a [u8], number: Option, suffix: &[u8]) -> Cow<'a, [ } } -/// Map unsupported bases to N. -fn sanitise_for_cbq(seq: &[u8]) -> Cow<'_, [u8]> { - #[inline] - fn supported(b: u8) -> bool { - matches!( - b, - b'A' | b'C' | b'G' | b'T' | b'N' | b'a' | b'c' | b'g' | b't' - ) - } - - let Some(first) = seq.iter().position(|&b| !supported(b)) else { - return Cow::Borrowed(seq); - }; - let mut sanitised = seq.to_vec(); - sanitised[first..].iter_mut().for_each(|b| { - if !supported(*b) { - *b = b'N'; - } - }); - Cow::Owned(sanitised) -} - /// Filtering config for an already-loaded index (no index path; see [`FilterConfig`]). pub struct FilterRunConfig { /// Path to input fastx file (or - for stdin) From 8f1f2d69e460bba539f042a9057118014734ecdc Mon Sep 17 00:00:00 2001 From: Bede Constantinides Date: Mon, 17 Aug 2026 11:35:59 +0100 Subject: [PATCH 06/12] Test consistent fastx & cbq ambiguity handling --- Cargo.lock | 1 + Cargo.toml | 4 ++-- src/minimizers.rs | 18 ++++++++++++++++++ tests/filter_tests.rs | 8 ++++++++ 4 files changed, 29 insertions(+), 2 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 00a8713..2d72798 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -172,6 +172,7 @@ version = "0.9.6" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "5fea1e47c3296e0073dc373f24660b1bd8ec36d6c06a7de115e4c75bf5b1753d" dependencies = [ + "anyhow", "auto_impl", "bitnuc 0.4.1", "bitnuc 0.5.2", diff --git a/Cargo.toml b/Cargo.toml index 4a08c8b..626f006 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -43,7 +43,7 @@ niffler = { version = "3.0.1", default-features = false, optional = true } paraseq = { version = "0.5.0", default-features = false, features = ["anyhow", "niffler"], optional = true } # CBQ support. Defaults disabled since we don't need binseq's own paraseq integration # We need binseq >=0.9.6 for safe canonicalisation of non-N ambiguous bases to N -binseq = { version = "0.9.6", default-features = false, optional = true } +binseq = { version = "0.9.6", default-features = false, features = ["anyhow"], optional = true } ensure_simd = { version = "0.1.0", optional = true } # Compression deps (optional) @@ -62,7 +62,7 @@ predicates = "3.0" tempfile = "3.20" rstest = "0.26" nix = { version = "0.31", features = ["fs"] } -binseq = { version = "0.9.6", default-features = false } +binseq = { version = "0.9.6", default-features = false, features = ["anyhow"] } [features] scalar = ["simd-minimizers/scalar", "packed-seq/scalar"] diff --git a/src/minimizers.rs b/src/minimizers.rs index 1fbf2ff..d7e498b 100644 --- a/src/minimizers.rs +++ b/src/minimizers.rs @@ -258,6 +258,24 @@ mod tests { assert!(short_minimizers.is_empty()); } + #[test] + fn non_acgt_mapped_to_n() { + let hasher = KmerHasher::new(31); + let minimizers = |base| { + let mut seq = b"ACGT".repeat(40); + seq[64] = base; + match compute_minimizers(&seq, &hasher, 31, 15) { + crate::MinimizerVec::U64(vec) => vec, + crate::MinimizerVec::U128(_) => unreachable!(), + } + }; + let expected = minimizers(b'N'); + + for &base in b"nRrYySsWwKkMmBbDdHhVvUu-" { + assert_eq!(minimizers(base), expected, "base {}", base as char); + } + } + #[test] fn test_calculate_scaled_entropy() { // Test short k-mers (should return 1.0 for k < 10) diff --git a/tests/filter_tests.rs b/tests/filter_tests.rs index 838ae28..ecb7d3a 100644 --- a/tests/filter_tests.rs +++ b/tests/filter_tests.rs @@ -3224,6 +3224,14 @@ fn cbq_iupac_codes_become_n() { ('W', 31), ('S', 32), ('K', 33), + ('B', 40), + ('D', 41), + ('H', 42), + ('V', 43), + ('b', 44), + ('d', 45), + ('h', 46), + ('v', 47), ('M', 63), ('y', 1), ('r', 16), From ab02e0ac951ed9f20c0f0f21a9ca0a5e52045258 Mon Sep 17 00:00:00 2001 From: Bede Constantinides Date: Mon, 17 Aug 2026 12:12:56 +0100 Subject: [PATCH 07/12] Update readme for CBQ --- CHANGELOG.md | 2 +- README.md | 42 ++++++++++++++++++++++++------------------ 2 files changed, 25 insertions(+), 19 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 658028e..ca79fe7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,7 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added - `deacon filter --ordered` preserves output order via paraseq 0.5.0; also exposed as `ordered` in the Python bindings and JSON summary. -- `deacon filter` accepts CBQ (BINSEQ) input, detected automatically, and writes CBQ when given a `.cbq` path; paired reads are stored as native paired records in a single file. Block size is set with `--cbq-block-size` (default 16 MiB). +- `deacon filter` accepts CBQ (BINSEQ) input, detected automatically, and writes CBQ when given a `.cbq` path; paired reads are stored as native paired records in a single file. Output block size is the greater of `--cbq-block-size` (default 16 MiB) and the input block size for CBQ input. ### Changed diff --git a/README.md b/README.md index 58e9046..1920524 100644 --- a/README.md +++ b/README.md @@ -10,7 +10,7 @@ Deacon filters DNA sequences in FASTA/Q files and streams using SIMD-accelerated minimizer comparison with query sequence(s), emitting either matching sequences (**search mode**), or sequences without matches (**deplete mode**). Sequences match when they share enough distinct minimizers with the indexed query to exceed chosen absolute and relative thresholds. Query size has little impact on filtering speed, enabling ultrafast search and depletion with gene-, genome- and pangenome-scale queries using a laptop. Deacon filters uncompressed FASTA/Q at **gigabases per second** on recent AMD, Intel (`x86_64`), and Apple `arm64` systems. Built with panhuman host depletion in mind—yet broadly useful for searching large sequence collections—Deacon delivers [leading classification accuracy](https://doi.org/10.1101/2025.06.09.658732) for host depletion and unrivalled speed using 5GB of RAM. -Default parameters are carefully chosen but easily changed. Classification sensitivity, specificity and memory requirements may be tuned by varying *k*-mer length (`-k`), window size (`-w`), absolute match threshold (`-a`) and relative match threshold (`-r`) . Minimizer `k` and `w` are chosen at query index time, while the match thresholds can be chosen at filter time. Matching sequences are those that share enough distinct minimizers with the indexed query to exceed *both* the absolute threshold (`-a`, default 2 shared minimizers) and the relative threshold (`-r`, default 0.01 [1%] shared minimizers). For paired sequences, hits in either mate counts towards a single match threshold for the pair. Deacon reports filtering performance during execution and optionally writes a JSON `--summary` upon completion. Sequences can optionally be renamed using `--rename` for privacy and smaller file sizes. Deacon fully supports stdin, stdout and natively handles gz, zst and xz compression formats, detected by file extension. +Default parameters are carefully chosen but easily changed. Classification sensitivity, specificity and memory requirements may be tuned by varying *k*-mer length (`-k`), window size (`-w`), absolute match threshold (`-a`) and relative match threshold (`-r`) . Minimizer `k` and `w` are chosen at query index time, while the match thresholds can be chosen at filter time. Matching sequences are those that share enough distinct minimizers with the indexed query to exceed *both* the absolute threshold (`-a`, default 2 shared minimizers) and the relative threshold (`-r`, default 0.01 [1%] shared minimizers). For paired sequences, hits in either mate counts towards a single match threshold for the pair. Deacon reports filtering performance during execution and optionally writes a JSON `--summary` upon completion. Sequences can optionally be renamed using `--rename` for privacy and smaller file sizes. Deacon fully supports stdin/out (uncompressed fastx) and natively handles .gz, .zst and .xz fastx file IO, as well as BINSEQ CBQ (.cbq). Benchmarks for panhuman host depletion of complex microbial metagenomes are described in a [preprint](https://www.biorxiv.org/content/10.1101/2025.06.09.658732v1). Deacon with the `panhuman-1` (*k*=31, w=15) index exhibited the highest balanced accuracy for both long and short simulated reads. Deacon was less specific only than Hostile for short reads. @@ -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`. 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`). With `--ordered`, CBQ input (and CBQ output combined with `--rename`) falls back to one filtering thread to guarantee deterministic output; without `--ordered`, CBQ `--rename` numbers are unique but not sequential across threads. +The main command `deacon filter` accepts an index path followed by up to two sequence 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. Paired read headers can be validated using `--check-pairs`. 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 compressed FASTX formats are detected automatically by file extension. Single and paired [BINSEQ](https://www.biorxiv.org/content/10.1101/2025.04.08.647863v2) CBQ files are natively supported for both input and output. #### Examples @@ -130,9 +130,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 and output -deacon filter panhuman-1.k31w15.idx reads.fq.gz -o filt.cbq -deacon filter panhuman-1.k31w15.idx filt.cbq -o filt.fq.gz +# BINSEQ CBQ input and output; paired records are stored in one file +deacon filter -d panhuman-1.k31w15.idx reads.fq.gz -o filt.cbq +deacon filter -d panhuman-1.k31w15.idx r1.fq.gz r2.fq.gz -o filt12.cbq +deacon filter -d panhuman-1.k31w15.idx filt12.cbq -o filt12.fq.gz # Replace read headers with incrementing integers deacon filter -d -R panhuman-1.k31w15.idx reads.fq.gz > filt.fq @@ -286,9 +287,9 @@ Options: Use `-s summary.json` to save detailed filtering statistics: ```json { - "version": "deacon 0.9.0", + "version": "deacon 0.17.0", "index": "panhuman-1.k31w15.idx", - "input": "HG02334.1m.fastq.gz", + "input": "HG02334.100MB.fastq.gz", "input2": null, "output": "-", "output2": null, @@ -299,18 +300,23 @@ Use `-s summary.json` to save detailed filtering statistics: "prefix_length": 0, "deplete": true, "rename": false, + "ordered": false, "check_pairs": false, - "seqs_in": 1000000, - "seqs_out": 13452, - "seqs_removed": 986548, - "seqs_removed_proportion": 0.986548, - "bp_in": 5477122928, - "bp_out": 5710050, - "bp_removed": 5471412878, - "bp_removed_proportion": 0.9989574727324798, - "time": 125.755103875, - "seqs_per_second": 7951, - "bp_per_second": 43553881 + "seqs_in": 37500, + "seqs_out": 454, + "seqs_out_proportion": 0.012106666666666667, + "seqs_removed": 37046, + "seqs_removed_proportion": 0.9878933333333333, + "bp_in": 141474280, + "bp_out": 227079, + "bp_out_proportion": 0.001605090338682056, + "bp_removed": 141247201, + "bp_removed_proportion": 0.9983949096613179, + "time": 0.446945667, + "seqs_per_second": 84129, + "bp_per_second": 317392822, + "seqs_per_second_total": 83902, + "bp_per_second_total": 316535745 } ``` From 44ae3f22ff7404fa500bec8936f09f5b3a8e1248 Mon Sep 17 00:00:00 2001 From: Bede Constantinides Date: Mon, 17 Aug 2026 13:34:05 +0100 Subject: [PATCH 08/12] Rename --fasta (-f) to --discard-quality, given we can now do cbq --- CHANGELOG.md | 1 + README.md | 12 ++++---- deacon-py/python/README.md | 2 +- deacon-py/python/deacon/_deacon.pyi | 2 +- deacon-py/src/lib.rs | 6 ++-- deacon-wasm/src/lib.rs | 22 +++++++------- deacon-wasm/worker.js | 4 +-- scripts/check-wasm-paired-parity.sh | 14 ++++----- scripts/check-wasm-parity.sh | 18 ++++++------ scripts/wasm_bench.mjs | 2 +- scripts/wasm_filter_to_file.mjs | 4 +-- scripts/wasm_paired_filter_to_file.mjs | 4 +-- scripts/wasm_smoke_test.mjs | 2 +- src/filter.rs | 40 +++++++++++++------------- src/lib.rs | 4 +-- src/main.rs | 20 ++++++------- tests/filter_tests.rs | 16 +++++------ 17 files changed, 87 insertions(+), 86 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ca79fe7..cb1382d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed - Python `Index.filter()` options after `fastq2` are keyword-only. +- `deacon filter --fasta` (`-f`) is now `--discard-quality` (no short form), since it discards quality scores for both fastx output (yielding FASTA) and CBQ output (yielding a quality-free CBQ). The Python binding's `output_fasta` keyword becomes `discard_quality`. ### Removed diff --git a/README.md b/README.md index 1920524..f44d03f 100644 --- a/README.md +++ b/README.md @@ -207,9 +207,7 @@ Options: -d, --deplete Discard matching sequences (invert filtering behaviour) -R, --rename - Replace sequence headers with incrementing numbers (reproducible with --ordered) - -f, --fasta - Output FASTA format regardless of input format + Replace sequence headers with incrementing numbers (deterministic with --ordered) -o, --output Path to output file (stdout by default; detects fastx with .gz, .zst, .xz, or binseq with .cbq) -O, --output2 @@ -221,9 +219,11 @@ Options: --compression-threads Number of threads used for output compression (0 = auto) [default: 0] --compression-level - Output compression level (1-9 for gz & xz; 1-22 for zstd, including cbq) [default: 2] + Output compression level (1-9 for gz & xz; 1-22 for zstd including cbq) [default: 2] --cbq-block-size cbq output block size in MiB (or cbq input block size if higher) [default: 16] + --discard-quality + Emit fasta or quality-free cbq regardless of input format --interleaved Treat INPUT as interleaved paired records from single file or stdin --ordered @@ -231,7 +231,7 @@ Options: --check-pairs Validate paired record names (Illumina CASAVA or /1 /2 suffixes) --debug - Output sequences with minimizer hits to stderr + Emit sequences with minimizer hits to stderr -q, --quiet Suppress progress reporting -h, --help @@ -277,7 +277,7 @@ Options: -w Minimizer window size used for indexing [default: 15] -o, --output Path to output file (stdout if not specified) -t, --threads Number of execution threads (0 = auto) [default: 8] - -q, --quiet Suppress sequence header output + -q, --quiet Suppress progress reporting -h, --help Print help ``` diff --git a/deacon-py/python/README.md b/deacon-py/python/README.md index 74c66c7..624e4f6 100644 --- a/deacon-py/python/README.md +++ b/deacon-py/python/README.md @@ -69,7 +69,7 @@ def filter( abs_threshold=2, # min absolute minimizer hits to call a match rel_threshold=0.01, # min proportion of minimizers hitting to call a match prefix_length=0, # only use the first N bp of each read (0 = whole read) - output_fasta=False, # emit FASTA instead of FASTQ + discard_quality=False, # discard quality scores ordered=False, # preserve input record ordering (deterministic, slightly slower) threads=8, # worker threads for filtering compression_level=2, # output compression level diff --git a/deacon-py/python/deacon/_deacon.pyi b/deacon-py/python/deacon/_deacon.pyi index 59db6f9..e51cb5e 100644 --- a/deacon-py/python/deacon/_deacon.pyi +++ b/deacon-py/python/deacon/_deacon.pyi @@ -12,7 +12,7 @@ class Index: """ Download a prebuilt index, then load and return it. """ - def filter(self, /, fastq: str, fastq2: str |None = None, *, interleaved: bool = False, deplete: bool = False, rename: bool = False, output: str |None = None, output2: str |None = None, abs_threshold: int = 2, rel_threshold: float = 0.01, prefix_length: int = 0, output_fasta: bool = False, ordered: bool = False, threads: int = 8, compression_level: int = 2, compression_threads: int = 0, debug: bool = False, quiet: bool = True) -> Any: ... + def filter(self, /, fastq: str, fastq2: str |None = None, *, interleaved: bool = False, deplete: bool = False, rename: bool = False, output: str |None = None, output2: str |None = None, abs_threshold: int = 2, rel_threshold: float = 0.01, prefix_length: int = 0, discard_quality: bool = False, ordered: bool = False, threads: int = 8, compression_level: int = 2, compression_threads: int = 0, debug: bool = False, quiet: bool = True) -> Any: ... def info(self, /) -> dict: """ Index metadata: k, w, format and minimizer/key count. diff --git a/deacon-py/src/lib.rs b/deacon-py/src/lib.rs index a37fab6..1183ade 100644 --- a/deacon-py/src/lib.rs +++ b/deacon-py/src/lib.rs @@ -100,7 +100,7 @@ impl Index { abs_threshold=2, rel_threshold=0.01, prefix_length=0, - output_fasta=false, + discard_quality=false, ordered=false, threads=8, compression_level=2, @@ -121,7 +121,7 @@ impl Index { abs_threshold: usize, rel_threshold: f64, prefix_length: usize, - output_fasta: bool, + discard_quality: bool, ordered: bool, threads: u16, compression_level: u8, @@ -147,7 +147,7 @@ impl Index { summary_path: None, deplete, rename, - output_fasta, + discard_quality, ordered, threads, compression_level, diff --git a/deacon-wasm/src/lib.rs b/deacon-wasm/src/lib.rs index 9feffd5..31ca209 100644 --- a/deacon-wasm/src/lib.rs +++ b/deacon-wasm/src/lib.rs @@ -70,7 +70,7 @@ pub struct FilterSession { index: Arc, kernel: FilterKernel, rename: bool, - output_fasta: bool, + discard_quality: bool, rename_counter: u64, parser: SeqChunkParser, stats: FilterStats, @@ -91,7 +91,7 @@ impl FilterSession { decompress_input: bool, compress_output: bool, rename: bool, - output_fasta: bool, + discard_quality: bool, ) -> Result { let k = index.inner.header.kmer_length(); let w = index.inner.header.window_size(); @@ -121,7 +121,7 @@ impl FilterSession { ) .map_err(|e| JsValue::from_str(&e.to_string()))?, rename, - output_fasta, + discard_quality, rename_counter: 0, parser: SeqChunkParser::new(), stats: FilterStats::default(), @@ -268,7 +268,7 @@ impl FilterSession { header, seq, qual, - self.output_fasta, + self.discard_quality, self.rename, counter, b"", @@ -291,7 +291,7 @@ pub struct PairedFilterSession { index: Arc, kernel: FilterKernel, rename: bool, - output_fasta: bool, + discard_quality: bool, rename_counter: u64, parser_r1: SeqChunkParser, parser_r2: SeqChunkParser, @@ -323,7 +323,7 @@ impl PairedFilterSession { compress_r1: bool, compress_r2: bool, rename: bool, - output_fasta: bool, + discard_quality: bool, ) -> Result { let k = index.inner.header.kmer_length(); let w = index.inner.header.window_size(); @@ -341,7 +341,7 @@ impl PairedFilterSession { ) .map_err(|e| JsValue::from_str(&e.to_string()))?, rename, - output_fasta, + discard_quality, rename_counter: 0, parser_r1: SeqChunkParser::new(), parser_r2: SeqChunkParser::new(), @@ -469,7 +469,7 @@ impl PairedFilterSession { &record1.header, &record1.seq, record1.qual.as_deref(), - self.output_fasta, + self.discard_quality, self.rename, counter, b"/1", @@ -480,7 +480,7 @@ impl PairedFilterSession { &record2.header, &record2.seq, record2.qual.as_deref(), - self.output_fasta, + self.discard_quality, self.rename, counter, b"/2", @@ -1006,12 +1006,12 @@ fn write_record( header: &[u8], seq: &[u8], qual: Option<&[u8]>, - output_fasta: bool, + discard_quality: bool, rename: bool, counter: u64, read_suffix: &[u8], ) -> std::io::Result<()> { - let is_fasta = output_fasta || qual.is_none(); + let is_fasta = discard_quality || qual.is_none(); output.write_all(if is_fasta { b">" } else { b"@" })?; if rename { output.write_all(counter.to_string().as_bytes())?; diff --git a/deacon-wasm/worker.js b/deacon-wasm/worker.js index 17143fb..975695d 100644 --- a/deacon-wasm/worker.js +++ b/deacon-wasm/worker.js @@ -35,7 +35,7 @@ async function streamFilterFile(file, opts) { isGz, // decompress_input isGz, // compress_output (match input format) false, // rename - false // output_fasta + false // discard_quality ); const reader = file.stream().getReader(); @@ -133,7 +133,7 @@ async function streamFilterPairedFiles(file1, file2, opts) { r1Gz, r2Gz, // decompress_r1, decompress_r2 r1Gz, r2Gz, // compress_r1, compress_r2 false, // rename - false // output_fasta + false // discard_quality ); const readerR1 = file1.stream().getReader(); diff --git a/scripts/check-wasm-paired-parity.sh b/scripts/check-wasm-paired-parity.sh index f5da565..0db97ba 100755 --- a/scripts/check-wasm-paired-parity.sh +++ b/scripts/check-wasm-paired-parity.sh @@ -53,15 +53,15 @@ echo "Native/WASM paired parity check passed" echo " R1 bytes: $(wc -c < "$native_r1" | tr -d ' ')" echo " R2 bytes: $(wc -c < "$native_r2" | tr -d ' ')" -# Second pass: --rename + --fasta (-R -f) must also match byte-for-byte. +# Second pass: --rename + --discard-quality must also match byte-for-byte. native_rf_r1="$work_dir/native.rf.r1.fasta" native_rf_r2="$work_dir/native.rf.r2.fasta" wasm_rf_r1="$work_dir/wasm.rf.r1.fasta" wasm_rf_r2="$work_dir/wasm.rf.r2.fasta" -echo "Running native paired filter (-R -f)..." -cargo run --release -- filter -d -t 1 -R -f "$index_path" "$reads1_path" "$reads2_path" \ +echo "Running native paired filter (--rename --discard-quality)..." +cargo run --release -- filter -d -t 1 -R --discard-quality "$index_path" "$reads1_path" "$reads2_path" \ -o "$native_rf_r1" -O "$native_rf_r2" -echo "Running WASM paired filter (--rename --fasta)..." +echo "Running WASM paired filter (--rename --discard-quality)..." node "$repo_root/scripts/wasm_paired_filter_to_file.mjs" \ --pkg "$pkg_dir" \ --index "$index_path" \ @@ -69,14 +69,14 @@ node "$repo_root/scripts/wasm_paired_filter_to_file.mjs" \ --reads2 "$reads2_path" \ --output1 "$wasm_rf_r1" \ --output2 "$wasm_rf_r2" \ - --deplete --rename --fasta + --deplete --rename --discard-quality if ! cmp -s "$native_rf_r1" "$wasm_rf_r1" || ! cmp -s "$native_rf_r2" "$wasm_rf_r2"; then - echo "Native/WASM paired output mismatch (rename+fasta)" >&2 + echo "Native/WASM paired output mismatch (rename+discard-quality)" >&2 echo " native R1: $(shasum -a 256 "$native_rf_r1" | awk '{print $1}')" >&2 echo " wasm R1: $(shasum -a 256 "$wasm_rf_r1" | awk '{print $1}')" >&2 echo " native R2: $(shasum -a 256 "$native_rf_r2" | awk '{print $1}')" >&2 echo " wasm R2: $(shasum -a 256 "$wasm_rf_r2" | awk '{print $1}')" >&2 exit 1 fi -echo "Native/WASM paired rename+fasta parity check passed" +echo "Native/WASM paired rename+discard-quality parity check passed" diff --git a/scripts/check-wasm-parity.sh b/scripts/check-wasm-parity.sh index 7969977..dced7ec 100755 --- a/scripts/check-wasm-parity.sh +++ b/scripts/check-wasm-parity.sh @@ -48,23 +48,23 @@ echo "Native/WASM parity check passed" echo " bytes: $size" echo " sha256: $sha" -# Second pass: --rename + --fasta (-R -f) must also match byte-for-byte. -native_rf="$work_dir/native.rename-fasta.fasta" -wasm_rf="$work_dir/wasm.rename-fasta.fasta" -echo "Running native filter (-R -f)..." -cargo run --release -- filter -t 1 -R -f "$index_path" "$reads_path" -o "$native_rf" -echo "Running WASM filter (--rename --fasta)..." +# Second pass: --rename + --discard-quality must also match byte-for-byte. +native_rf="$work_dir/native.rename-discard-quality.fasta" +wasm_rf="$work_dir/wasm.rename-discard-quality.fasta" +echo "Running native filter (--rename --discard-quality)..." +cargo run --release -- filter -t 1 -R --discard-quality "$index_path" "$reads_path" -o "$native_rf" +echo "Running WASM filter (--rename --discard-quality)..." node "$repo_root/scripts/wasm_filter_to_file.mjs" \ --pkg "$pkg_dir" \ --index "$index_path" \ --reads "$reads_path" \ --output "$wasm_rf" \ - --rename --fasta + --rename --discard-quality if ! cmp -s "$native_rf" "$wasm_rf"; then - echo "Native/WASM output mismatch (rename+fasta)" >&2 + echo "Native/WASM output mismatch (rename+discard-quality)" >&2 echo " native: $(shasum -a 256 "$native_rf" | awk '{print $1}')" >&2 echo " wasm: $(shasum -a 256 "$wasm_rf" | awk '{print $1}')" >&2 exit 1 fi -echo "Native/WASM rename+fasta parity check passed ($(wc -c < "$native_rf" | tr -d ' ') bytes)" +echo "Native/WASM rename+discard-quality parity check passed ($(wc -c < "$native_rf" | tr -d ' ') bytes)" diff --git a/scripts/wasm_bench.mjs b/scripts/wasm_bench.mjs index 76e4c09..2382c50 100644 --- a/scripts/wasm_bench.mjs +++ b/scripts/wasm_bench.mjs @@ -64,7 +64,7 @@ for (let iter = 0; iter < iters; iter++) { decompressInput, false, // compress_output false, // rename - false, // output_fasta + false, // discard_quality ); const start = process.hrtime.bigint(); diff --git a/scripts/wasm_filter_to_file.mjs b/scripts/wasm_filter_to_file.mjs index 554e610..715af28 100644 --- a/scripts/wasm_filter_to_file.mjs +++ b/scripts/wasm_filter_to_file.mjs @@ -50,7 +50,7 @@ const indexPath = requireArg(args, "index"); const readsPath = requireArg(args, "reads"); const outputPath = requireArg(args, "output"); const rename = args.flags.has("rename"); -const outputFasta = args.flags.has("fasta"); +const discardQuality = args.flags.has("discard-quality"); const wasmModule = await import(pathToFileURL(path.join(pkgDir, "deacon_wasm.js")).href); @@ -69,7 +69,7 @@ const session = new wasmModule.FilterSession( true, // decompress_input false, // compress_output rename, - outputFasta, + discardQuality, ); const input = createReadStream(readsPath, { highWaterMark: 256 * 1024 }); diff --git a/scripts/wasm_paired_filter_to_file.mjs b/scripts/wasm_paired_filter_to_file.mjs index 9e012fb..8b2cce8 100644 --- a/scripts/wasm_paired_filter_to_file.mjs +++ b/scripts/wasm_paired_filter_to_file.mjs @@ -49,7 +49,7 @@ const output1Path = requireArg(args, "output1"); const output2Path = requireArg(args, "output2"); const deplete = args.flags.has("deplete"); const rename = args.flags.has("rename"); -const outputFasta = args.flags.has("fasta"); +const discardQuality = args.flags.has("discard-quality"); const absThreshold = Number(args.get("abs-threshold") ?? 2); const relThreshold = Number(args.get("rel-threshold") ?? 0.01); @@ -72,7 +72,7 @@ const session = new wasmModule.PairedFilterSession( false, // compress_r1 false, // compress_r2 rename, - outputFasta, + discardQuality, ); const input1 = createReadStream(reads1Path, { highWaterMark: 256 * 1024 }); diff --git a/scripts/wasm_smoke_test.mjs b/scripts/wasm_smoke_test.mjs index 13e9de6..101589b 100644 --- a/scripts/wasm_smoke_test.mjs +++ b/scripts/wasm_smoke_test.mjs @@ -55,7 +55,7 @@ const session = new wasmModule.FilterSession( false, // decompress_input false, // compress_output false, // rename - false, // output_fasta + false, // discard_quality ); const outputChunks = []; diff --git a/src/filter.rs b/src/filter.rs index a98bd51..3795d54 100644 --- a/src/filter.rs +++ b/src/filter.rs @@ -119,10 +119,10 @@ impl Output { read: &ReadView, suffix: &'static [u8], rename: bool, - output_fasta: bool, + discard_quality: bool, ) -> Result<()> { let offset = buffer.len(); - let marker = format_record_to_buffer(read, rename, output_fasta, buffer)?; + let marker = format_record_to_buffer(read, rename, discard_quality, buffer)?; if rename { pending.push(PendingRename { offset, @@ -139,7 +139,7 @@ impl Output { &mut self, read: &ReadView, rename: bool, - output_fasta: bool, + discard_quality: bool, rename_counter: &AtomicU64, ) -> Result<()> { match self { @@ -149,7 +149,7 @@ impl Output { batch_kept, .. } => { - Self::push_fastx(local, pending, *batch_kept, read, b"", rename, output_fasta)?; + Self::push_fastx(local, pending, *batch_kept, read, b"", rename, discard_quality)?; if rename { *batch_kept += 1; } @@ -176,7 +176,7 @@ impl Output { read1: &ReadView, read2: &ReadView, rename: bool, - output_fasta: bool, + discard_quality: bool, rename_counter: &AtomicU64, ) -> Result<()> { match self { @@ -196,7 +196,7 @@ impl Output { read1, b"/1", rename, - output_fasta, + discard_quality, )?; if shared2.is_some() { // Separate outputs @@ -207,7 +207,7 @@ impl Output { read2, b"/2", rename, - output_fasta, + discard_quality, )?; } else { // Interleaved output @@ -218,7 +218,7 @@ impl Output { read2, b"/2", rename, - output_fasta, + discard_quality, )?; } if rename { @@ -369,8 +369,8 @@ pub struct FilterRunConfig { pub deplete: bool, /// Replace sequence headers with incrementing numbers pub rename: bool, - /// Force FASTA output (discards quality scores) - pub output_fasta: bool, + /// Emit fasta or quality-free cbq regardless of input format + pub discard_quality: bool, /// Preserve input record ordering (deterministic, slightly slower) pub ordered: bool, /// Number of execution threads (0 = auto) @@ -396,7 +396,7 @@ struct FilterProcessorConfig { prefix_length: usize, deplete: bool, rename: bool, - output_fasta: bool, + discard_quality: bool, debug: bool, check_pairs: bool, ordered: bool, @@ -722,10 +722,10 @@ struct PendingRename { fn format_record_to_buffer( read: &ReadView, rename: bool, - output_fasta: bool, + discard_quality: bool, buffer: &mut Vec, ) -> Result { - let is_fasta = output_fasta || read.qual.is_none(); + let is_fasta = discard_quality || read.qual.is_none(); let marker = if is_fasta { b'>' } else { b'@' }; // Header (omitted when renaming) @@ -921,7 +921,7 @@ struct FilterProcessor { // Minimizer matching parameters minimizers: Arc, rename: bool, - output_fasta: bool, + discard_quality: bool, debug: bool, check_pairs: bool, /// Write batches in input order, not completion order @@ -963,7 +963,7 @@ impl FilterProcessor { Ok(Self { minimizers, rename: config.rename, - output_fasta: config.output_fasta, + discard_quality: config.discard_quality, debug: config.debug, check_pairs: config.check_pairs, ordered: config.ordered, @@ -1052,7 +1052,7 @@ impl FilterProcessor { if decision.keep { self.local_stats.output_bp += read.seq.len() as u64; self.output - .push_read(read, self.rename, self.output_fasta, &self.rename_counter)?; + .push_read(read, self.rename, self.discard_quality, &self.rename_counter)?; } else { self.local_stats.filtered_seqs += 1; self.local_stats.filtered_bp += read.seq.len() as u64; @@ -1095,7 +1095,7 @@ impl FilterProcessor { read1, read2, self.rename, - self.output_fasta, + self.discard_quality, &self.rename_counter, )?; } else { @@ -1270,7 +1270,7 @@ pub fn run(config: &FilterConfig) -> Result { summary_path: config.summary_path.cloned(), deplete: config.deplete, rename: config.rename, - output_fasta: config.output_fasta, + discard_quality: config.discard_quality, ordered: config.ordered, threads: config.threads, compression_level: config.compression_level, @@ -1460,7 +1460,7 @@ pub fn run_with_index( .expect("CBQ output implies a named path"); let cbq_writer = BinseqWriterBuilder::new(BinseqFormat::Cbq) .paired(layout.paired) - .quality(layout.qualities && !config.output_fasta) + .quality(layout.qualities && !config.discard_quality) .headers(layout.headers || config.rename) .flags(layout.flags) .block_size( @@ -1518,7 +1518,7 @@ pub fn run_with_index( prefix_length: config.prefix_length, deplete: config.deplete, rename: config.rename, - output_fasta: config.output_fasta, + discard_quality: config.discard_quality, debug: config.debug, check_pairs: config.check_pairs, ordered: config.ordered, diff --git a/src/lib.rs b/src/lib.rs index e17a13f..1ce04c2 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -339,8 +339,8 @@ pub struct FilterConfig<'a> { /// Replace sequence headers with incrementing numbers (1, 2, 3...) pub rename: bool, - /// Force FASTA output (discards quality scores) - pub output_fasta: bool, + /// Emit fasta or quality-free cbq regardless of input format + pub discard_quality: bool, /// Preserve input record ordering (deterministic, slightly slower) pub ordered: bool, diff --git a/src/main.rs b/src/main.rs index 76cb65f..38e04fb 100644 --- a/src/main.rs +++ b/src/main.rs @@ -62,14 +62,10 @@ enum Command { #[arg(short = 'd', long = "deplete", default_value_t = false)] deplete: bool, - /// Replace sequence headers with incrementing numbers (reproducible with --ordered) + /// Replace sequence headers with incrementing numbers (deterministic with --ordered) #[arg(short = 'R', long = "rename", default_value_t = false)] rename: bool, - /// Output FASTA format regardless of input format - #[arg(short = 'f', long = "fasta", default_value_t = false)] - output_fasta: bool, - /// Path to output file (stdout by default; detects fastx with .gz, .zst, .xz, or binseq with .cbq) #[arg(short = 'o', long = "output")] output: Option, @@ -90,7 +86,7 @@ enum Command { #[arg(long = "compression-threads", default_value_t = 0)] compression_threads: u16, - /// Output compression level (1-9 for gz & xz; 1-22 for zstd, including cbq) + /// Output compression level (1-9 for gz & xz; 1-22 for zstd including cbq) #[arg(long = "compression-level", default_value_t = 2)] compression_level: u8, @@ -102,6 +98,10 @@ enum Command { )] cbq_block_size: u16, + /// Emit fasta or quality-free cbq regardless of input format + #[arg(long = "discard-quality", default_value_t = false)] + discard_quality: bool, + /// Treat INPUT as interleaved paired records from single file or stdin #[arg( long = "interleaved", @@ -118,7 +118,7 @@ enum Command { #[arg(long = "check-pairs", default_value_t = false)] check_pairs: bool, - /// Output sequences with minimizer hits to stderr + /// Emit sequences with minimizer hits to stderr #[arg(long = "debug", default_value_t = false)] debug: bool, @@ -172,7 +172,7 @@ enum IndexCommand { #[arg(short = 't', long = "threads", default_value_t = 8)] threads: u16, - /// Suppress sequence header output + /// Suppress progress reporting #[arg(short = 'q', long = "quiet")] quiet: bool, }, @@ -539,7 +539,7 @@ fn process_command(command: &Command) -> Result<(), anyhow::Error> { summary, deplete, rename, - output_fasta, + discard_quality, threads, compression_level, cbq_block_size, @@ -571,7 +571,7 @@ fn process_command(command: &Command) -> Result<(), anyhow::Error> { summary_path: summary.as_ref(), deplete: *deplete, rename: *rename, - output_fasta: *output_fasta, + discard_quality: *discard_quality, ordered: *ordered, threads: *threads, compression_level: *compression_level, diff --git a/tests/filter_tests.rs b/tests/filter_tests.rs index ecb7d3a..d5e2c86 100644 --- a/tests/filter_tests.rs +++ b/tests/filter_tests.rs @@ -293,7 +293,7 @@ fn test_filter_rename() { } #[test] -fn test_filter_fasta_flag() { +fn test_filter_discard_quality_flag() { let temp_dir = tempdir().unwrap(); let fasta_path = temp_dir.path().join("ref.fasta"); let fastq_path = temp_dir.path().join("reads.fastq"); @@ -306,7 +306,7 @@ fn test_filter_fasta_flag() { let mut cmd = cargo::cargo_bin_cmd!("deacon"); let output = cmd .arg("filter") - .arg("-f") + .arg("--discard-quality") .arg("-a") .arg("1") .arg("-r") @@ -322,7 +322,7 @@ fn test_filter_fasta_flag() { let output_str = std::str::from_utf8(&output).unwrap(); assert!( output_str.starts_with('>'), - "FASTQ in with -f should gen FASTA out" + "FASTQ in with -Q should gen FASTA out" ); } @@ -3028,10 +3028,10 @@ fn cbq_roundtrip_matches_fastx() { "paired CBQ round trip must match the interleaved FASTQ output" ); - // --fasta creates a quality-free CBQ + // --discard-quality 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"]) + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1", "--discard-quality"]) .arg(&bin_path) .arg(&fastq_path) .arg("--output") @@ -3041,7 +3041,7 @@ fn cbq_roundtrip_matches_fastx() { let reader = binseq::cbq::MmapReader::new(&fasta_cbq).unwrap(); assert!( !reader.header().has_qualities(), - "--fasta must produce a quality-free CBQ" + "--discard-quality must produce a quality-free CBQ" ); assert_eq!(reader.num_records(), 2); @@ -3164,7 +3164,7 @@ fn cbq_ns_match_fastx() { } cmd.args(["-r", "0", "-t", "1"]); if fasta { - cmd.arg("--fasta"); + cmd.arg("--discard-quality"); } cmd.arg(&bin_path) .arg(input) @@ -3275,7 +3275,7 @@ fn cbq_iupac_codes_become_n() { let mut cmd = cargo::cargo_bin_cmd!("deacon"); cmd.args(["filter", "-d", "-a", "65535", "-r", "1", "-t", "1"]); if fasta { - cmd.arg("--fasta"); + cmd.arg("--discard-quality"); } cmd.arg(&bin_path) .arg(from) From c0e986d32d2822abfa5864b1e9a6779515e88724 Mon Sep 17 00:00:00 2001 From: Bede Constantinides Date: Tue, 18 Aug 2026 11:40:25 +0100 Subject: [PATCH 09/12] Accept .cba alternative cbq io extension, by convention (no qual) --- CHANGELOG.md | 2 +- README.md | 9 ++- src/filter.rs | 50 +++++++++---- src/lib.rs | 2 +- src/main.rs | 2 +- tests/filter_tests.rs | 169 +++++++++++++++++++++++++++++++++++------- 6 files changed, 188 insertions(+), 46 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index cb1382d..829fe16 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,7 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added - `deacon filter --ordered` preserves output order via paraseq 0.5.0; also exposed as `ordered` in the Python bindings and JSON summary. -- `deacon filter` accepts CBQ (BINSEQ) input, detected automatically, and writes CBQ when given a `.cbq` path; paired reads are stored as native paired records in a single file. Output block size is the greater of `--cbq-block-size` (default 16 MiB) and the input block size for CBQ input. +- `deacon filter` accepts CBQ (BINSEQ) input and writes CBQ when given a `.cbq` path. Paired CBQ records are supported. Output block size is the greater of `--cbq-block-size` (default 16 MiB) and the input block size given CBQ input. By convention a `.cba` output extension also writes CBQ without quality scores, implying `--discard-quality`. ### Changed diff --git a/README.md b/README.md index f44d03f..a7fe052 100644 --- a/README.md +++ b/README.md @@ -10,7 +10,7 @@ Deacon filters DNA sequences in FASTA/Q files and streams using SIMD-accelerated minimizer comparison with query sequence(s), emitting either matching sequences (**search mode**), or sequences without matches (**deplete mode**). Sequences match when they share enough distinct minimizers with the indexed query to exceed chosen absolute and relative thresholds. Query size has little impact on filtering speed, enabling ultrafast search and depletion with gene-, genome- and pangenome-scale queries using a laptop. Deacon filters uncompressed FASTA/Q at **gigabases per second** on recent AMD, Intel (`x86_64`), and Apple `arm64` systems. Built with panhuman host depletion in mind—yet broadly useful for searching large sequence collections—Deacon delivers [leading classification accuracy](https://doi.org/10.1101/2025.06.09.658732) for host depletion and unrivalled speed using 5GB of RAM. -Default parameters are carefully chosen but easily changed. Classification sensitivity, specificity and memory requirements may be tuned by varying *k*-mer length (`-k`), window size (`-w`), absolute match threshold (`-a`) and relative match threshold (`-r`) . Minimizer `k` and `w` are chosen at query index time, while the match thresholds can be chosen at filter time. Matching sequences are those that share enough distinct minimizers with the indexed query to exceed *both* the absolute threshold (`-a`, default 2 shared minimizers) and the relative threshold (`-r`, default 0.01 [1%] shared minimizers). For paired sequences, hits in either mate counts towards a single match threshold for the pair. Deacon reports filtering performance during execution and optionally writes a JSON `--summary` upon completion. Sequences can optionally be renamed using `--rename` for privacy and smaller file sizes. Deacon fully supports stdin/out (uncompressed fastx) and natively handles .gz, .zst and .xz fastx file IO, as well as BINSEQ CBQ (.cbq). +Default parameters are carefully chosen but easily changed. Classification sensitivity, specificity and memory requirements may be tuned by varying *k*-mer length (`-k`), window size (`-w`), absolute match threshold (`-a`) and relative match threshold (`-r`) . Minimizer `k` and `w` are chosen at query index time, while the match thresholds can be chosen at filter time. Matching sequences are those that share enough distinct minimizers with the indexed query to exceed *both* the absolute threshold (`-a`, default 2 shared minimizers) and the relative threshold (`-r`, default 0.01 [1%] shared minimizers). For paired sequences, hits in either mate counts towards a single match threshold for the pair. Deacon reports filtering performance during execution and optionally writes a JSON `--summary` upon completion. Sequences can optionally be renamed using `--rename` for privacy and smaller file sizes. Deacon fully supports stdin/out (uncompressed fastx) and natively handles .gz, .zst and .xz fastx file IO, as well as BINSEQ CBQ (.cbq & .cba). Benchmarks for panhuman host depletion of complex microbial metagenomes are described in a [preprint](https://www.biorxiv.org/content/10.1101/2025.06.09.658732v1). Deacon with the `panhuman-1` (*k*=31, w=15) index exhibited the highest balanced accuracy for both long and short simulated reads. Deacon was less specific only than Hostile for short reads. @@ -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 sequence 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. Paired read headers can be validated using `--check-pairs`. 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 compressed FASTX formats are detected automatically by file extension. Single and paired [BINSEQ](https://www.biorxiv.org/content/10.1101/2025.04.08.647863v2) CBQ files are natively supported for both input and output. +The main command `deacon filter` accepts an index path followed by up to two sequence 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. Paired read headers can be validated using `--check-pairs`. 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 compressed FASTX formats are detected automatically by file extension. Single and paired [BINSEQ](https://www.biorxiv.org/content/10.1101/2025.04.08.647863v2) CBQ files are natively supported for both input and output. CBQ input is detected by content rather than file extension, while a `.cbq` output path writes CBQ, and a `.cba` output path writes CBQ without quality. #### Examples @@ -135,6 +135,9 @@ deacon filter -d panhuman-1.k31w15.idx reads.fq.gz -o filt.cbq deacon filter -d panhuman-1.k31w15.idx r1.fq.gz r2.fq.gz -o filt12.cbq deacon filter -d panhuman-1.k31w15.idx filt12.cbq -o filt12.fq.gz +# A .cba extension writes cbq without quality values, if present +deacon filter -d panhuman-1.k31w15.idx reads.fq.gz -o filt.cba + # Replace read headers with incrementing integers deacon filter -d -R panhuman-1.k31w15.idx reads.fq.gz > filt.fq @@ -209,7 +212,7 @@ Options: -R, --rename Replace sequence headers with incrementing numbers (deterministic with --ordered) -o, --output - Path to output file (stdout by default; detects fastx with .gz, .zst, .xz, or binseq with .cbq) + Path to output file (fastx to stdout by default; detects .gz, .zst, .xz, .cbq, .cba) -O, --output2 Optional path to second paired output fastx file (detects .gz, .zst, .xz) -s, --summary diff --git a/src/filter.rs b/src/filter.rs index 3795d54..5c768ac 100644 --- a/src/filter.rs +++ b/src/filter.rs @@ -149,7 +149,15 @@ impl Output { batch_kept, .. } => { - Self::push_fastx(local, pending, *batch_kept, read, b"", rename, discard_quality)?; + Self::push_fastx( + local, + pending, + *batch_kept, + read, + b"", + rename, + discard_quality, + )?; if rename { *batch_kept += 1; } @@ -369,7 +377,7 @@ pub struct FilterRunConfig { pub deplete: bool, /// Replace sequence headers with incrementing numbers pub rename: bool, - /// Emit fasta or quality-free cbq regardless of input format + /// Emit fasta or quality-free cbq regardless of input format (implied by a .cba output) pub discard_quality: bool, /// Preserve input record ordering (deterministic, slightly slower) pub ordered: bool, @@ -519,10 +527,20 @@ fn create_paraseq_reader(path: Option<&str>) -> Result bool { + path.ends_with(".cbq") || is_cba_path(path) +} + +/// `.cba` names a quality-free CBQ, implying `--discard-quality` for that output +fn is_cba_path(path: &str) -> bool { + path.ends_with(".cba") +} + /// Resolve the output format from the output path suffix fn resolve_output_format(config: &FilterRunConfig) -> Format { match config.output_path.as_deref() { - Some(path) if path.to_string_lossy().ends_with(".cbq") => Format::Cbq, + Some(path) if is_cbq_path(&path.to_string_lossy()) => Format::Cbq, _ => Format::Fastx, } } @@ -682,11 +700,7 @@ fn validate_input_output( if output_format == Format::Cbq && config.output2_path.is_some() { anyhow::bail!("CBQ output does not support OUTPUT2; CBQ pairing is native"); } - if config - .output2_path - .as_deref() - .is_some_and(|output2| output2.ends_with(".cbq")) - { + if config.output2_path.as_deref().is_some_and(is_cbq_path) { anyhow::bail!("OUTPUT2 cannot be CBQ; CBQ pairing is native (use a single --output)"); } if output_format == Format::Cbq && !(1..=22).contains(&config.compression_level) { @@ -1051,8 +1065,12 @@ impl FilterProcessor { if decision.keep { self.local_stats.output_bp += read.seq.len() as u64; - self.output - .push_read(read, self.rename, self.discard_quality, &self.rename_counter)?; + self.output.push_read( + read, + self.rename, + self.discard_quality, + &self.rename_counter, + )?; } else { self.local_stats.filtered_seqs += 1; self.local_stats.filtered_bp += read.seq.len() as u64; @@ -1241,7 +1259,7 @@ pub fn run(config: &FilterConfig) -> Result { if config.check_pairs && !config.interleaved && config.input2_path.is_none() - && !config.input_path.ends_with(".cbq") + && !is_cbq_path(config.input_path) { validate_check_pairs_mode(true, false)?; } @@ -1385,6 +1403,12 @@ pub fn run_with_index( 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); + // `.cba` names a quality-free CBQ, implying --discard-quality for the output + let discard_quality = config.discard_quality + || config + .output_path + .as_deref() + .is_some_and(|path| is_cba_path(&path.to_string_lossy())); let (layout, input) = open_input(config, interleaved_input)?; validate_input_output(&layout, output_format, config)?; @@ -1460,7 +1484,7 @@ pub fn run_with_index( .expect("CBQ output implies a named path"); let cbq_writer = BinseqWriterBuilder::new(BinseqFormat::Cbq) .paired(layout.paired) - .quality(layout.qualities && !config.discard_quality) + .quality(layout.qualities && !discard_quality) .headers(layout.headers || config.rename) .flags(layout.flags) .block_size( @@ -1518,7 +1542,7 @@ pub fn run_with_index( prefix_length: config.prefix_length, deplete: config.deplete, rename: config.rename, - discard_quality: config.discard_quality, + discard_quality, debug: config.debug, check_pairs: config.check_pairs, ordered: config.ordered, diff --git a/src/lib.rs b/src/lib.rs index 1ce04c2..7612ccd 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -339,7 +339,7 @@ pub struct FilterConfig<'a> { /// Replace sequence headers with incrementing numbers (1, 2, 3...) pub rename: bool, - /// Emit fasta or quality-free cbq regardless of input format + /// Emit fasta or quality-free cbq regardless of input format (implied by a .cba output) pub discard_quality: bool, /// Preserve input record ordering (deterministic, slightly slower) diff --git a/src/main.rs b/src/main.rs index 38e04fb..d97c00b 100644 --- a/src/main.rs +++ b/src/main.rs @@ -66,7 +66,7 @@ enum Command { #[arg(short = 'R', long = "rename", default_value_t = false)] rename: bool, - /// Path to output file (stdout by default; detects fastx with .gz, .zst, .xz, or binseq with .cbq) + /// Path to output file (fastx to stdout by default; detects .gz, .zst, .xz, .cbq, .cba) #[arg(short = 'o', long = "output")] output: Option, diff --git a/tests/filter_tests.rs b/tests/filter_tests.rs index d5e2c86..a761d85 100644 --- a/tests/filter_tests.rs +++ b/tests/filter_tests.rs @@ -3031,7 +3031,16 @@ fn cbq_roundtrip_matches_fastx() { // --discard-quality 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", "--discard-quality"]) + .args([ + "filter", + "-a", + "1", + "-r", + "0.0", + "-t", + "1", + "--discard-quality", + ]) .arg(&bin_path) .arg(&fastq_path) .arg("--output") @@ -3311,28 +3320,132 @@ fn cbq_check_pairs_paired_input() { create_test_interleaved_fastq(&interleaved_path); build_index(&fasta_path, &bin_path); - // Interleaved FASTQ with read1/1, read1/2 names -> paired CBQ - let paired_cbq = temp_dir.path().join("paired.cbq"); + // Interleaved FASTQ with read1/1, read1/2 names -> paired CBQ, under both + // CBQ suffixes; --check-pairs is pre-flighted from the input path + for suffix in ["cbq", "cba"] { + let paired_cbq = temp_dir.path().join(format!("paired.{suffix}")); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1", "--interleaved"]) + .arg(&bin_path) + .arg(&interleaved_path) + .arg("--output") + .arg(&paired_cbq) + .assert() + .success(); + + // Paired CBQ + --check-pairs succeeds despite no INPUT2/--interleaved + let out = temp_dir.path().join(format!("out-{suffix}.fastx")); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1", "--check-pairs"]) + .arg(&bin_path) + .arg(&paired_cbq) + .arg("--output") + .arg(&out) + .assert() + .success(); + let content = fs::read_to_string(&out).unwrap(); + if suffix == "cbq" { + assert_eq!(count_records(&content), 4); + } else { + // .cba dropped quality, so it decodes to FASTA + assert_eq!(content.matches('>').count(), 4); + } + } +} + +/// A `.cba` output path writes a CBQ without qualities, implying --discard-quality +#[test] +fn cba_output_is_quality_free_cbq() { + 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); + + // FASTQ -> .cba, without --discard-quality + let cba_path = temp_dir.path().join("out.cba"); cargo::cargo_bin_cmd!("deacon") - .args(["filter", "-a", "1", "-r", "0.0", "-t", "1", "--interleaved"]) + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) .arg(&bin_path) - .arg(&interleaved_path) + .arg(&fastq_path) .arg("--output") - .arg(&paired_cbq) + .arg(&cba_path) .assert() .success(); + assert!( + fs::read(&cba_path).unwrap().starts_with(b"CBQFILE"), + ".cba output must be a CBQ file" + ); + let reader = binseq::cbq::MmapReader::new(&cba_path).unwrap(); + assert!( + !reader.header().has_qualities(), + ".cba must imply --discard-quality" + ); + assert!(reader.header().has_headers(), ".cba must keep headers"); + assert_eq!(reader.num_records(), 2); - // Paired CBQ + --check-pairs succeeds despite no INPUT2/--interleaved - let out = temp_dir.path().join("out.fastq"); + // ... and is byte-identical to the same run written as .cbq with the flag + let cbq_path = temp_dir.path().join("out.cbq"); cargo::cargo_bin_cmd!("deacon") - .args(["filter", "-a", "1", "-r", "0.0", "-t", "1", "--check-pairs"]) + .args([ + "filter", + "-a", + "1", + "-r", + "0.0", + "-t", + "1", + "--discard-quality", + ]) .arg(&bin_path) - .arg(&paired_cbq) + .arg(&fastq_path) .arg("--output") - .arg(&out) + .arg(&cbq_path) .assert() .success(); - assert_eq!(count_records(&fs::read_to_string(&out).unwrap()), 4); + assert_eq!( + fs::read(&cba_path).unwrap(), + fs::read(&cbq_path).unwrap(), + ".cba must match .cbq written with --discard-quality" + ); + + // A .cba reads back as CBQ input (magic sniffed) and decodes to FASTA + let roundtrip = temp_dir.path().join("roundtrip.fastx"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&cba_path) + .arg("--output") + .arg(&roundtrip) + .assert() + .success(); + let content = fs::read_to_string(&roundtrip).unwrap(); + assert!( + content.starts_with('>'), + "quality-free CBQ decodes to FASTA" + ); + assert_eq!(content.matches('>').count(), 2); + + // Paired input stays natively paired in a .cba + let paired_cba = temp_dir.path().join("paired.cba"); + 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_cba) + .assert() + .success(); + let reader = binseq::cbq::MmapReader::new(&paired_cba).unwrap(); + assert!(reader.is_paired(), "paired input must stay paired in .cba"); + assert!(!reader.header().has_qualities()); + assert_eq!(reader.num_records(), 2); } /// A file bearing the CBQ magic but truncated must error, not panic @@ -3560,21 +3673,23 @@ fn cbq_output2_suffix_rejected() { create_test_paired_fastq(&r1_path, &r2_path); build_index(&fasta_path, &bin_path); - let out1 = temp_dir.path().join("out1.fastq"); - let out2 = temp_dir.path().join("out2.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(&out1) - .arg("--output2") - .arg(&out2) - .assert() - .failure(); - assert!(!out1.exists(), "OUTPUT must not be created"); - assert!(!out2.exists(), "OUTPUT2 must not be created"); + for suffix in ["cbq", "cba"] { + let out1 = temp_dir.path().join(format!("out1-{suffix}.fastq")); + let out2 = temp_dir.path().join(format!("out2.{suffix}")); + 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(&out1) + .arg("--output2") + .arg(&out2) + .assert() + .failure(); + assert!(!out1.exists(), "OUTPUT must not be created"); + assert!(!out2.exists(), "OUTPUT2 must not be created"); + } } /// Headerless, quality-free paired CBQ filters to CBQ (still headerless) and From a895081202d512c7bb253daa12fd0b09cb81ba02 Mon Sep 17 00:00:00 2001 From: Bede Constantinides Date: Tue, 18 Aug 2026 11:55:49 +0100 Subject: [PATCH 10/12] Fix problems with --ordered cbq output --- src/filter.rs | 20 ++++++++------- tests/filter_tests.rs | 59 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 70 insertions(+), 9 deletions(-) diff --git a/src/filter.rs b/src/filter.rs index 5c768ac..96a0b96 100644 --- a/src/filter.rs +++ b/src/filter.rs @@ -28,7 +28,7 @@ pub const DEFAULT_CBQ_BLOCK_SIZE_MIB: u16 = 16; type BoxedWriter = Box; /// CBQ output is always a named file, no stdout -type CbqWriter = binseq::BinseqWriter>; +type CbqWriter = binseq::BinseqWriter; /// Sequence file format #[derive(Clone, Copy, PartialEq, Eq)] @@ -263,11 +263,15 @@ impl Output { ) -> Result<()> { match self { Output::Cbq { local, shared } => { + let mut shared = shared.lock(); if ordered { - // Drain the incomplete block too, or records linger across batches - shared.lock().ingest(local)?; + // Keep the shared partial block ahead of the next batch. + if let binseq::BinseqWriter::Cbq(writer) = &mut *shared { + writer.flush()?; + } + shared.ingest(local)?; } else { - shared.lock().ingest_completed(local)?; + shared.ingest_completed(local)?; } } Output::Fastx { @@ -824,11 +828,9 @@ fn count_compressed_outputs(config: &FilterRunConfig) -> u8 { count } -/// Open a CBQ output file directly, bypassing the compressed writer stack -fn open_cbq_output(path: &std::path::Path) -> Result> { - let file = File::create(path) - .with_context(|| format!("Failed to create output file: {}", path.display()))?; - Ok(BufWriter::with_capacity(OUTPUT_BUFFER_SIZE, file)) +/// Open unbuffered so embedded-index write errors surface from `finish`. +fn open_cbq_output(path: &std::path::Path) -> Result { + File::create(path).with_context(|| format!("Failed to create output file: {}", path.display())) } /// Return a suitable writer for the output path extension diff --git a/tests/filter_tests.rs b/tests/filter_tests.rs index a761d85..b6bf680 100644 --- a/tests/filter_tests.rs +++ b/tests/filter_tests.rs @@ -3583,6 +3583,65 @@ fn cbq_ordered_output_is_input_ordered() { ); } +/// Write reads that complete CBQ blocks mid-batch. +fn create_many_long_fastq(path: &Path, n: usize, len: usize) -> String { + let mut content = String::new(); + let seq = "CGTA".repeat(len / 4); + let qual = "~".repeat(seq.len()); + for i in 0..n { + content.push_str(&format!("@read_{i:06}\n{seq}\n+\n{qual}\n")); + } + fs::write(path, &content).unwrap(); + content +} + +/// --ordered CBQ output stays ordered across completed blocks. +#[test] +fn cbq_ordered_block_completes_mid_batch() { + let temp_dir = tempdir().unwrap(); + 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); + + // Force each batch across multiple 1 MiB blocks. + let fastq_path = temp_dir.path().join("long.fastq"); + let content = create_many_long_fastq(&fastq_path, 2048, 1000); + + let cbq_path = temp_dir.path().join("ordered.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args([ + "filter", + "-d", + "--ordered", + "-t", + "4", + "--cbq-block-size", + "1", + ]) + .arg(&aaa_bin) + .arg(&fastq_path) + .arg("--output") + .arg(&cbq_path) + .assert() + .success(); + + let out = temp_dir.path().join("out.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "--ordered", "-t", "1"]) + .arg(&aaa_bin) + .arg(&cbq_path) + .arg("--output") + .arg(&out) + .assert() + .success(); + assert_eq!( + fs::read_to_string(&out).unwrap(), + content, + "multi-block CBQ output must preserve input order" + ); +} + /// --ordered --rename CBQ output numbers records 1..N in input order #[test] fn cbq_ordered_rename_is_sequential() { From 5b2a2e62bf740ea18d158582ad580bed04a8c66c Mon Sep 17 00:00:00 2001 From: Bede Constantinides Date: Tue, 18 Aug 2026 12:04:10 +0100 Subject: [PATCH 11/12] Make --check-pairs sniff cbq rather than rely on extension --- src/filter.rs | 10 ---------- tests/filter_tests.rs | 34 +++++++++++++++++++++++++++------- 2 files changed, 27 insertions(+), 17 deletions(-) diff --git a/src/filter.rs b/src/filter.rs index 96a0b96..d4af53d 100644 --- a/src/filter.rs +++ b/src/filter.rs @@ -1256,16 +1256,6 @@ pub fn run(config: &FilterConfig) -> Result { validate_unit_interval("complexity threshold", threshold)?; } - // Fail --check-pairs before the index loads; CBQ may be natively paired, - // so skip it there and let validate_input_output recheck once the header is read - if config.check_pairs - && !config.interleaved - && config.input2_path.is_none() - && !is_cbq_path(config.input_path) - { - validate_check_pairs_mode(true, false)?; - } - // Validate the index path once here; run_with_index never touches it again. if !config.minimizers_path.exists() { return Err(anyhow::anyhow!( diff --git a/tests/filter_tests.rs b/tests/filter_tests.rs index b6bf680..70e6ae7 100644 --- a/tests/filter_tests.rs +++ b/tests/filter_tests.rs @@ -530,18 +530,25 @@ fn test_check_pairs_rejects_mismatched_names() { } #[test] -fn test_check_pairs_requires_paired_input_before_index_load() { +fn test_check_pairs_requires_paired_input() { + 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"); + create_test_fasta(&fasta_path); + create_test_fastq(&fastq_path); + build_index(&fasta_path, &bin_path); + let mut cmd = cargo::cargo_bin_cmd!("deacon"); cmd.arg("filter") .arg("--check-pairs") - .arg("missing.idx") - .arg("sequences.fastq") + .arg(&bin_path) + .arg(&fastq_path) .assert() .failure() .stderr(predicates::str::contains( "--check-pairs requires paired input", - )) - .stderr(predicates::str::contains("Index file does not exist").not()); + )); } #[test] @@ -3320,8 +3327,7 @@ fn cbq_check_pairs_paired_input() { create_test_interleaved_fastq(&interleaved_path); build_index(&fasta_path, &bin_path); - // Interleaved FASTQ with read1/1, read1/2 names -> paired CBQ, under both - // CBQ suffixes; --check-pairs is pre-flighted from the input path + // Interleaved FASTQ with read1/1, read1/2 names -> paired CBQ. for suffix in ["cbq", "cba"] { let paired_cbq = temp_dir.path().join(format!("paired.{suffix}")); cargo::cargo_bin_cmd!("deacon") @@ -3351,6 +3357,20 @@ fn cbq_check_pairs_paired_input() { assert_eq!(content.matches('>').count(), 4); } } + + // Input format is detected from content, not the suffix. + let disguised = temp_dir.path().join("paired.dat"); + fs::copy(temp_dir.path().join("paired.cbq"), &disguised).unwrap(); + let out = temp_dir.path().join("out-dat.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1", "--check-pairs"]) + .arg(&bin_path) + .arg(&disguised) + .arg("--output") + .arg(&out) + .assert() + .success(); + assert_eq!(count_records(&fs::read_to_string(out).unwrap()), 4); } /// A `.cba` output path writes a CBQ without qualities, implying --discard-quality From fd5708b0f3a0e9a6904885053ca34338ea51b861 Mon Sep 17 00:00:00 2001 From: Bede Constantinides Date: Tue, 18 Aug 2026 13:23:12 +0100 Subject: [PATCH 12/12] Bring python bindings up to date, include in CI --- .github/workflows/release-pypi.yaml | 4 +- .github/workflows/release.yml | 4 +- .github/workflows/test.yml | 26 ++++- CHANGELOG.md | 12 +- deacon-py/python/README.md | 42 ++++--- deacon-py/python/deacon/_deacon.pyi | 9 +- deacon-py/src/lib.rs | 96 +++++++++++----- deacon-py/tests/test_api.py | 169 ++++++++++++++++++++++++++++ 8 files changed, 306 insertions(+), 56 deletions(-) create mode 100644 deacon-py/tests/test_api.py diff --git a/.github/workflows/release-pypi.yaml b/.github/workflows/release-pypi.yaml index 58bf114..ee7c1a3 100644 --- a/.github/workflows/release-pypi.yaml +++ b/.github/workflows/release-pypi.yaml @@ -21,7 +21,7 @@ jobs: runs-on: ${{ matrix.os }} steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v7 # portable SIMD floor (x86-64-v3/AVX2, apple-m1) - never target-cpu=native for wheels - name: Use portable cargo config @@ -45,7 +45,7 @@ jobs: sdist: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v7 - name: Build sdist uses: PyO3/maturin-action@v1 with: diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index 7204905..9034f76 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -26,7 +26,7 @@ jobs: runs-on: ${{ matrix.os }} steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v7 - uses: dtolnay/rust-toolchain@stable - uses: Swatinem/rust-cache@v2 - name: Build & test & rename @@ -65,7 +65,7 @@ jobs: # Skip prerelease tags (SemVer prereleases contain a hyphen, e.g. 0.17.0-rc.1). if: ${{ !contains(github.ref_name, '-') }} steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v7 - uses: dtolnay/rust-toolchain@stable - name: Publish deacon to crates.io run: cargo publish -p deacon diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 8b06e63..5858a15 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -16,7 +16,7 @@ jobs: os: [ubuntu-24.04, ubuntu-24.04-arm, macos-15] steps: - name: Checkout repository - uses: actions/checkout@v4 + uses: actions/checkout@v7 - name: Install Rust toolchain uses: dtolnay/rust-toolchain@stable - name: Cache cargo registry @@ -28,3 +28,27 @@ jobs: # Not a default workspace member, so cargo test skips it - name: Run deacon-wasm unit tests run: cargo test -p deacon-wasm + + python-test: + name: Python bindings + runs-on: ubuntu-24.04 + steps: + - name: Checkout repository + uses: actions/checkout@v7 + - name: Install Rust toolchain + uses: dtolnay/rust-toolchain@stable + - name: Set up Python + uses: actions/setup-python@v7 + with: + python-version: "3.12" + - name: Cache cargo registry + uses: Swatinem/rust-cache@v2 + - name: Use portable config + run: mv .cargo/config-portable.toml .cargo/config.toml + - name: Build bindings + run: | + python -m venv .venv + .venv/bin/python -m pip install "maturin>=1.7,<2" + .venv/bin/maturin develop -m deacon-py/Cargo.toml + - name: Run Python tests + run: .venv/bin/python -m unittest discover -s deacon-py/tests -v diff --git a/CHANGELOG.md b/CHANGELOG.md index 829fe16..399b507 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,17 +5,19 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## [Unreleased] +## [0.17.0] - 2026-08-18 ### Added -- `deacon filter --ordered` preserves output order via paraseq 0.5.0; also exposed as `ordered` in the Python bindings and JSON summary. -- `deacon filter` accepts CBQ (BINSEQ) input and writes CBQ when given a `.cbq` path. Paired CBQ records are supported. Output block size is the greater of `--cbq-block-size` (default 16 MiB) and the input block size given CBQ input. By convention a `.cba` output extension also writes CBQ without quality scores, implying `--discard-quality`. +- `deacon filter --ordered` preserves output order via paraseq 0.5.0. Arg `ordered` added to Python bindings and JSON summary. +- `deacon filter` accepts CBQ (BINSEQ) input and writes CBQ when given a `.cbq` path. Single and paired CBQ is supported. Output block size is the greater of `--cbq-block-size` (default 16 MiB) and the input block size given CBQ input. By convention an alternate `.cba` output extension writes CBQ without quality scores, implying `--discard-quality`. +- Python `Index.filter()` now exposes `check_pairs`, `summary`, and `cbq_block_size`, and accepts path-like objects for all paths. ### Changed -- Python `Index.filter()` options after `fastq2` are keyword-only. -- `deacon filter --fasta` (`-f`) is now `--discard-quality` (no short form), since it discards quality scores for both fastx output (yielding FASTA) and CBQ output (yielding a quality-free CBQ). The Python binding's `output_fasta` keyword becomes `discard_quality`. +- **Breaking (Python):** index/filter input paths are positional-only; options and all `Index.fetch()` arguments are keyword-only. `filter` also renames `fastq`/`fastq2` to `input`/`input2`. +- Python filtering rejects `abs_threshold=0`, matching the CLI. +- Renamed `--fasta` (`-f`)/`output_fasta` to `--discard-quality`/`discard_quality`. ### Removed diff --git a/deacon-py/python/README.md b/deacon-py/python/README.md index 624e4f6..9eb0635 100644 --- a/deacon-py/python/README.md +++ b/deacon-py/python/README.md @@ -5,7 +5,7 @@ Python bindings for [Deacon](https://github.com/bede/deacon), enabling fast mult ## Installation ```bash -uv install deacon +uv pip install deacon ``` ## Quickstart @@ -33,47 +33,61 @@ index = Index("panhuman-1.k31w15.idx", complexity_threshold=None) Pass `complexity_threshold` (0.0–1.0, e.g. `0.9`) to discard low-complexity index minimizers once at load using kdust; the filtered set is then reused across all `filter` calls. Not supported for `bff` (binary fuse filter) indexes. +`path` is positional-only and accepts a string or path-like object. `complexity_threshold` is keyword-only. + ## `Index.fetch()` Download a prebuilt index, then load and return it (a static method, so `Index.fetch(...)` returns an `Index`). `output` is the local path to save to; when omitted it defaults to `"{name}.k{k}w{w}.idx"` in the working directory. The index is downloaded on every call — there is no local cache, so an existing file at that path is overwritten. ```python -index = Index.fetch("panhuman-1", k=31, w=15, output=None, complexity_threshold=None) +index = Index.fetch( + name="panhuman-1", + k=31, + w=15, + output=None, + complexity_threshold=None, +) ``` +All `fetch()` arguments are keyword-only. + ## `Index.info()` -`Index.info(index_path)` Returns a `dict` of index metadata: +`index.info()` returns a `dict` of index metadata: | Key | Meaning | | --- | --- | | `k` | *k*-mer length | | `w` | minimizer window size | | `format` | `exact-u64`, `exact-u128`, or `bff` (binary fuse filter) | -| `count` | number of keys in index (fingerprint slot count for `bff`) | +| `count` | number of minimizers/keys represented by the index | ## `Index.filter()` -Filter a FASTA/FASTQ file or file pair against the index and return a `dict` of summary statistics. Auto-detects `.gz`/`.zst`/`.xz` compression on both input and output based on file extension. The Python GIL is released while filtering meaning calls benefit from multithreading. Refer to the [main Deacon readme](https://github.com/bede/deacon) for more detailed usage examples. +Filter FASTA, FASTQ, or CBQ input against the index and return a `dict` of summary statistics. FASTA/FASTQ compression (`.gz`, `.zst`, or `.xz`) is detected automatically. A `.cbq` output writes CBQ; `.cba` writes quality-free CBQ. The Python GIL is released while filtering, so calls benefit from multithreading. Refer to the [main Deacon readme](https://github.com/bede/deacon) for more detailed usage examples. ```python def filter( - fastq, # input path (FASTA/FASTQ, optionally .gz/.zst/.xz) - fastq2=None, # second mate for paired reads - *, # remaining options are keyword-only - interleaved=False, # treat fastq as interleaved paired reads (cannot combine with fastq2) + input, # positional-only FASTA, FASTQ, or CBQ path + /, + *, # every remaining argument is keyword-only + input2=None, # second FASTA/FASTQ mate + interleaved=False, # treat input as interleaved pairs (cannot combine with input2) + check_pairs=False, # validate paired read names deplete=False, # False = search (keep matches); True = deplete (remove matches) rename=False, # replace read names with sequential integers - output=None, # output path; None writes to stdout + output=None, # output path; None writes FASTA/FASTQ to stdout output2=None, # second output path for paired reads + summary=None, # optional JSON summary output path abs_threshold=2, # min absolute minimizer hits to call a match rel_threshold=0.01, # min proportion of minimizers hitting to call a match prefix_length=0, # only use the first N bp of each read (0 = whole read) - discard_quality=False, # discard quality scores + discard_quality=False, # discard quality scores ordered=False, # preserve input record ordering (deterministic, slightly slower) threads=8, # worker threads for filtering compression_level=2, # output compression level compression_threads=0, # threads for output compression (0 = auto) + cbq_block_size=16, # CBQ output block size in MiB (1-1024) quiet=True, # suppress progress/log output on stderr debug=False, # verbose per-read debug output ) -> dict @@ -81,9 +95,11 @@ def filter( **Modes.** With `deplete=False` (the default, *search* mode) reads that match the index are kept; with `deplete=True` reads that match are removed (host depletion). A read is a match only when it clears **both** thresholds: at least `abs_threshold` minimizer hits **and** at least `rel_threshold` of its minimizers hitting the index. -**Output.** When `output` is `None` the filtered records are written to stdout. To count without keeping the filtered sequences, pass `output="/dev/null"`. Statistics are returned regardless. +**Paired and CBQ input.** Use `input2=` for separate FASTA/FASTQ mates or `interleaved=True` for an interleaved stream. `check_pairs=True` validates Illumina CASAVA or `/1` and `/2` names. CBQ stores pairing internally, so CBQ input cannot be combined with `input2` or `interleaved`; paired CBQ output uses one `output` and cannot use `output2`. + +**Output.** When `output` is `None` the filtered records are written to stdout. To count without keeping the filtered sequences, pass `output="/dev/null"`. Pass `summary=` to write the same statistics returned by the call as JSON. For CBQ output, `cbq_block_size` is a lower bound: a larger input CBQ block size is preserved. -**Return value.** A `dict` including the run configuration (`version`, `index`, `input`/`input2`, `output`/`output2`, `k`, `w`, `abs_threshold`, `rel_threshold`, `prefix_length`, `deplete`, `rename`, `ordered`) and the results: +**Return value.** A `dict` including the run configuration (`version`, `index`, `input`/`input2`, `output`/`output2`, `k`, `w`, `abs_threshold`, `rel_threshold`, `prefix_length`, `deplete`, `rename`, `ordered`, `check_pairs`) and the results: | Key | Meaning | | --- | --- | diff --git a/deacon-py/python/deacon/_deacon.pyi b/deacon-py/python/deacon/_deacon.pyi index e51cb5e..0870901 100644 --- a/deacon-py/python/deacon/_deacon.pyi +++ b/deacon-py/python/deacon/_deacon.pyi @@ -1,18 +1,19 @@ from _typeshed import Incomplete -from typing import Any, final +from os import PathLike +from typing import final @final class Index: """ A loaded minimizer index, reusable across many `filter` calls. """ - def __new__(cls, /, path: str, complexity_threshold: float |None = None) -> Index: ... + def __new__(cls, path: str |PathLike[str], /, *, complexity_threshold: float |None = None) -> Index: ... @staticmethod - def fetch(name: str = "panhuman-1", k: int = 31, w: int = 15, output: str |None = None, complexity_threshold: float |None = None) -> Index: + def fetch(*, name: str = "panhuman-1", k: int = 31, w: int = 15, output: str |PathLike[str] |None = None, complexity_threshold: float |None = None) -> Index: """ Download a prebuilt index, then load and return it. """ - def filter(self, /, fastq: str, fastq2: str |None = None, *, interleaved: bool = False, deplete: bool = False, rename: bool = False, output: str |None = None, output2: str |None = None, abs_threshold: int = 2, rel_threshold: float = 0.01, prefix_length: int = 0, discard_quality: bool = False, ordered: bool = False, threads: int = 8, compression_level: int = 2, compression_threads: int = 0, debug: bool = False, quiet: bool = True) -> Any: ... + def filter(self, input: str |PathLike[str], /, *, input2: str |PathLike[str] |None = None, interleaved: bool = False, check_pairs: bool = False, deplete: bool = False, rename: bool = False, output: str |PathLike[str] |None = None, output2: str |PathLike[str] |None = None, summary: str |PathLike[str] |None = None, abs_threshold: int = 2, rel_threshold: float = 0.01, prefix_length: int = 0, discard_quality: bool = False, ordered: bool = False, threads: int = 8, compression_level: int = 2, compression_threads: int = 0, cbq_block_size: int = 16, quiet: bool = True, debug: bool = False) -> dict: ... def info(self, /) -> dict: """ Index metadata: k, w, format and minimizer/key count. diff --git a/deacon-py/src/lib.rs b/deacon-py/src/lib.rs index 1183ade..e111d9f 100644 --- a/deacon-py/src/lib.rs +++ b/deacon-py/src/lib.rs @@ -12,10 +12,19 @@ use pyo3::exceptions::{PyRuntimeError, PyValueError}; use pyo3::prelude::*; use pyo3::types::PyDict; +// Keep the literal in the Python signature/stub while detecting a changed core default at compile time. +const _: () = assert!(DEFAULT_CBQ_BLOCK_SIZE_MIB == 16); + fn to_pyerr(e: anyhow::Error) -> PyErr { PyRuntimeError::new_err(e.to_string()) } +fn path_to_string(path: PathBuf, argument: &str) -> PyResult { + path.into_os_string().into_string().map_err(|_| { + PyValueError::new_err(format!("{argument} must be representable as valid UTF-8")) + }) +} + /// A loaded minimizer index, reusable across many `filter` calls. #[pyclass(frozen)] struct Index { @@ -25,13 +34,9 @@ struct Index { minimizers: Arc, } -#[pymethods] impl Index { - #[new] - #[pyo3(signature = (path, complexity_threshold=None))] - fn new(path: &str, complexity_threshold: Option) -> PyResult { - let (mut minimizers, header) = - load_index_from_path_auto(Path::new(path)).map_err(to_pyerr)?; + fn load(path: &Path, complexity_threshold: Option) -> PyResult { + let (mut minimizers, header) = load_index_from_path_auto(path).map_err(to_pyerr)?; // Discard low-complexity index minimizers once at load (kdust); reused across filters. if let Some(threshold) = complexity_threshold { if matches!(minimizers, MinimizerSet::Fuse(_)) { @@ -49,28 +54,35 @@ impl Index { .map_err(to_pyerr)?; } Ok(Index { - label: path.to_string(), + label: path.to_string_lossy().into_owned(), k: header.kmer_length(), w: header.window_size(), minimizers: Arc::new(minimizers), }) } +} + +#[pymethods] +impl Index { + #[new] + #[pyo3(signature = (path, /, *, complexity_threshold=None))] + fn new(path: PathBuf, complexity_threshold: Option) -> PyResult { + Self::load(&path, complexity_threshold) + } /// Download a prebuilt index, then load and return it. #[staticmethod] - #[pyo3(signature = (name="panhuman-1", k=31, w=15, output=None, complexity_threshold=None))] + #[pyo3(signature = (*, name="panhuman-1", k=31, w=15, output=None, complexity_threshold=None))] fn fetch( name: &str, k: u8, w: u8, - output: Option, + output: Option, complexity_threshold: Option, ) -> PyResult { - let out_path = output - .map(PathBuf::from) - .unwrap_or_else(|| PathBuf::from(format!("{name}.k{k}w{w}.idx"))); + let out_path = output.unwrap_or_else(|| PathBuf::from(format!("{name}.k{k}w{w}.idx"))); index_fetch(name, k, w, Some(&out_path)).map_err(to_pyerr)?; - Index::new(&out_path.to_string_lossy(), complexity_threshold) + Index::load(&out_path, complexity_threshold) } /// Index metadata: k, w, format and minimizer/key count. @@ -89,14 +101,17 @@ impl Index { } #[pyo3(signature = ( - fastq, - fastq2=None, + input, + /, *, + input2=None, interleaved=false, + check_pairs=false, deplete=false, rename=false, output=None, output2=None, + summary=None, abs_threshold=2, rel_threshold=0.01, prefix_length=0, @@ -105,19 +120,22 @@ impl Index { threads=8, compression_level=2, compression_threads=0, + cbq_block_size=16, quiet=true, debug=false, ))] fn filter( &self, py: Python<'_>, - fastq: String, - fastq2: Option, + input: PathBuf, + input2: Option, interleaved: bool, + check_pairs: bool, deplete: bool, rename: bool, - output: Option, - output2: Option, + output: Option, + output2: Option, + summary: Option, abs_threshold: usize, rel_threshold: f64, prefix_length: usize, @@ -126,32 +144,50 @@ impl Index { threads: u16, compression_level: u8, compression_threads: u16, + cbq_block_size: u16, quiet: bool, debug: bool, - ) -> PyResult> { - if interleaved && fastq2.is_some() { + ) -> PyResult> { + if interleaved && input2.is_some() { + return Err(PyValueError::new_err( + "interleaved cannot be combined with input2 (interleaved input is a single file/stream)", + )); + } + if abs_threshold == 0 { + return Err(PyValueError::new_err("abs_threshold must be at least 1")); + } + if !(1..=1024).contains(&cbq_block_size) { return Err(PyValueError::new_err( - "interleaved cannot be combined with fastq2 (interleaved input is a single file/stream)", + "cbq_block_size must be between 1 and 1024 MiB inclusive", )); } + + let input = path_to_string(input, "input")?; + let input2 = input2 + .map(|path| path_to_string(path, "input2")) + .transpose()?; + let output2 = output2 + .map(|path| path_to_string(path, "output2")) + .transpose()?; + let cfg = FilterRunConfig { - input_path: fastq, - input2_path: fastq2, + input_path: input, + input2_path: input2, interleaved, - check_pairs: false, - output_path: output.map(PathBuf::from), + check_pairs, + output_path: output, output2_path: output2, abs_threshold, rel_threshold, prefix_length, - summary_path: None, + summary_path: summary, deplete, rename, discard_quality, ordered, threads, compression_level, - cbq_block_size: DEFAULT_CBQ_BLOCK_SIZE_MIB, + cbq_block_size, compression_threads, debug, quiet, @@ -163,7 +199,9 @@ impl Index { let summary = py .detach(|| run_with_index(mins, &IndexHeader::new(k, w), &cfg)) .map_err(to_pyerr)?; - Ok(pythonize::pythonize(py, &summary)?.unbind()) + Ok(pythonize::pythonize(py, &summary)? + .cast_into::()? + .unbind()) } } diff --git a/deacon-py/tests/test_api.py b/deacon-py/tests/test_api.py new file mode 100644 index 0000000..20e9675 --- /dev/null +++ b/deacon-py/tests/test_api.py @@ -0,0 +1,169 @@ +import base64 +import inspect +import json +import tempfile +import unittest +from pathlib import Path + +from deacon import Index + + +ROOT = Path(__file__).resolve().parents[2] +READS = ROOT / "tests/data/test_small_1.fastq.gz" +# Exact index built from READS with the default k=31, w=15 parameters. +INDEX_BYTES = base64.b64decode( + "Ax8PQQAAAAAAAACve33W9VXXBESYNVc0k/ANrPWbvtZvcxLn2gd9//W+G89JPdM1xNcVtW1Qym/XCivXZ11fdU0UFa1UtVGcmykSrVPjVL3GpxEJZtxZ51QZB0QNr3t91vUVfff33d+Xax+t1xAXVm6mCL0XtVHXF1ED1V1Xd13ddS2x4wlm3FnnFFB1vRe1UdcXbdbLWX/+vRdUbqZIhFlzBUXkvsPlUP8G1fVe1EZdXwR91TVRVG6mCH0RNbzu9VkXlPLbtcLaMw1d70Vt1PVFFNT3r+8bfb8CpNA3lfA6TQ/1VddEUbmZIm3WkELfVMIrReVmikSYNRe9bdaQQt9UAjmjHEaev7UttVLVRnFupgg94y7CTJzvNBSVmykSYdYcvdP3Xd8zfR/QtJiVr9nKGm3U9UXU8LoXkedvbRuU8htRw+ten3V9FdS/YW/QtJgVKb9dK6w90xR0DfF1ddfVHd3fd39frn0bLWbla7ay1i9Z11ddE0V9FY997/R91/cMXddGXWN1rRTx+27ff33PNb1B02JWvmYrfde1UddYXSs1VtdKVRvFOdC9bdaQQt8U1bVTFFRd7wVEUbmZIhFmDdU1UdRXXes1FVGjUF07RQHXeg1xYeVmCn3X90zfd38fTfVQMNVDQjPVXVd3XZt1PT3T1K1F5L4Dd13ddW3W9R3dWkTuO1wONVHfv75v9P0K" +) +TEST_DIRECTORY = None +INDEX = None + + +def setUpModule(): + global TEST_DIRECTORY, INDEX + TEST_DIRECTORY = tempfile.TemporaryDirectory() + INDEX = Path(TEST_DIRECTORY.name) / "test.idx" + INDEX.write_bytes(INDEX_BYTES) + + +def tearDownModule(): + TEST_DIRECTORY.cleanup() + + +class SignatureTests(unittest.TestCase): + @classmethod + def setUpClass(cls): + cls.index = Index(INDEX) + + def test_signatures_are_future_proof(self): + self.assertEqual( + str(inspect.signature(Index)), + "(path, /, *, complexity_threshold=None)", + ) + self.assertEqual( + str(inspect.signature(Index.fetch)), + "(*, name='panhuman-1', k=31, w=15, output=None, complexity_threshold=None)", + ) + + parameters = inspect.signature(Index.filter).parameters + self.assertEqual(parameters["input"].kind, inspect.Parameter.POSITIONAL_ONLY) + self.assertEqual(parameters["cbq_block_size"].default, 16) + for name, parameter in parameters.items(): + if name not in {"self", "input"}: + self.assertEqual(parameter.kind, inspect.Parameter.KEYWORD_ONLY, name) + + def test_optional_positional_arguments_are_rejected(self): + with self.assertRaises(TypeError): + Index(INDEX, None) + with self.assertRaises(TypeError): + Index.fetch("panhuman-1") + with self.assertRaises(TypeError): + self.index.filter(READS, None) + + def test_primary_paths_are_positional_only(self): + with self.assertRaises(TypeError): + Index(path=INDEX) + with self.assertRaises(TypeError): + self.index.filter(input=READS) + + +class FilteringTests(unittest.TestCase): + @classmethod + def setUpClass(cls): + cls.index = Index(INDEX) + + def test_pathlike_output_and_json_summary(self): + with tempfile.TemporaryDirectory() as directory: + directory = Path(directory) + output = directory / "output.fastq" + summary_path = directory / "summary.json" + + summary = self.index.filter( + READS, + output=output, + summary=summary_path, + threads=1, + ) + + self.assertTrue(output.is_file()) + self.assertEqual(json.loads(summary_path.read_text()), summary) + self.assertEqual(summary["check_pairs"], False) + self.assertGreater(summary["seqs_in"], 0) + + def test_check_pairs(self): + with tempfile.TemporaryDirectory() as directory: + directory = Path(directory) + mate1 = directory / "mate1.fastq" + mate2 = directory / "mate2.fastq" + mate1.write_text("@read/1\nACGT\n+\nIIII\n") + mate2.write_text("@read/2\nACGT\n+\nIIII\n") + summary = self.index.filter( + mate1, + input2=mate2, + output=directory / "output1.fastq", + output2=directory / "output2.fastq", + check_pairs=True, + threads=1, + ) + self.assertTrue(summary["check_pairs"]) + + mate1.write_text("@read-a/1\nACGT\n+\nIIII\n") + mate2.write_text("@read-b/2\nACGT\n+\nIIII\n") + with self.assertRaisesRegex(RuntimeError, "Paired record name mismatch"): + self.index.filter( + mate1, + input2=mate2, + output=directory / "bad1.fastq", + output2=directory / "bad2.fastq", + check_pairs=True, + threads=1, + ) + + def test_cbq_options_and_round_trip(self): + with tempfile.TemporaryDirectory() as directory: + directory = Path(directory) + cbq = directory / "output.cbq" + first = self.index.filter( + READS, + output=cbq, + cbq_block_size=1, + threads=1, + ) + second = self.index.filter( + cbq, + output=directory / "roundtrip.fastq", + threads=1, + ) + self.assertEqual(second["seqs_in"], first["seqs_out"]) + + cba = directory / "output.cba" + self.index.filter(READS, output=cba, cbq_block_size=1, threads=1) + fasta = directory / "quality-free.fasta" + self.index.filter(cba, output=fasta, threads=1) + if fasta.stat().st_size: + self.assertEqual(fasta.read_bytes()[:1], b">") + + with self.assertRaisesRegex(RuntimeError, "CBQ input does not support INPUT2"): + self.index.filter(cbq, input2=READS, threads=1) + with self.assertRaisesRegex(RuntimeError, "CBQ output does not support OUTPUT2"): + self.index.filter( + READS, + output=directory / "invalid.cbq", + output2=directory / "invalid.fastq", + threads=1, + ) + + def test_argument_validation(self): + with self.assertRaisesRegex(ValueError, "abs_threshold"): + self.index.filter(READS, abs_threshold=0) + with self.assertRaisesRegex(ValueError, "cbq_block_size"): + self.index.filter(READS, cbq_block_size=0) + with self.assertRaisesRegex(ValueError, "cbq_block_size"): + self.index.filter(READS, cbq_block_size=1025) + with self.assertRaisesRegex(ValueError, "input2"): + self.index.filter(READS, input2=READS, interleaved=True) + with self.assertRaisesRegex(RuntimeError, "check-pairs requires paired input"): + self.index.filter(READS, check_pairs=True) + + +if __name__ == "__main__": + unittest.main()