From bc0bcfdc45b0164540193efbcf2344a90853441d Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 7 Aug 2026 18:36:18 +0200 Subject: [PATCH 1/2] fix(io): decode multi-member gzip input instead of truncating it `flate2::read::GzDecoder` stops at the end of the first gzip member and reports EOF. A `.gz` written as several concatenated members was therefore read partially, with no error and no warning: bcl2fastq output, `cat a.fq.gz b.fq.gz > merged.fq.gz`, and every BGZF file are all multi-member. Measured on the bundled fixtures before the fix: a two-member `test/reads.fq` gzip (3 reads, 12 lines, identical to the single-member file under `gunzip`) reported `Number of input reads | 2`. Switches the four read paths to `MultiGzDecoder`: FASTQ input, the solo barcode whitelist, solo counting, and the `emptydrops` binary. Adds `test_fastq_reader_gzip_multi_member`, which builds a two-member gzip and asserts both reads come back. Verified it fails on the old decoder and passes on the new one. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 7 ++++++ src/bin/emptydrops.rs | 4 +-- src/io/fastq.rs | 52 +++++++++++++++++++++++++++++++++++++-- src/ruSTAR.code-workspace | 11 +++++++++ src/solo/count.rs | 4 ++- src/solo/whitelist.rs | 4 +-- 6 files changed, 75 insertions(+), 7 deletions(-) create mode 100644 src/ruSTAR.code-workspace diff --git a/CHANGELOG.md b/CHANGELOG.md index 14b4346..2656ba6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -42,6 +42,13 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Bug fixes +- **Multi-member gzip input is no longer truncated.** Compressed input was + decoded with `flate2::read::GzDecoder`, which stops at the end of the first + gzip member; a `.gz` made of several concatenated members (bcl2fastq output, + `cat a.fq.gz b.fq.gz`, any BGZF file) was read partially with no error and no + warning. All four read paths now use `MultiGzDecoder`: FASTQ input, the solo + barcode whitelist, solo counting, and the `emptydrops` binary. + - Read names are cut at `--readNameSeparator` (default `/`), as STAR does. A read named `foo/1` was previously emitted as `foo/1` where STAR emits `foo`. diff --git a/src/bin/emptydrops.rs b/src/bin/emptydrops.rs index 29913cd..4cac8f4 100644 --- a/src/bin/emptydrops.rs +++ b/src/bin/emptydrops.rs @@ -19,7 +19,7 @@ use std::fs::File; use std::io::{BufRead, BufReader, BufWriter, Read, Write}; use std::path::{Path, PathBuf}; -use flate2::read::GzDecoder; +use flate2::read::MultiGzDecoder; use rustar_aligner::rng::{SplitMix64, cumulative_weights, sample_cumulative}; struct Args { @@ -89,7 +89,7 @@ fn find(d: &Path, base: &str) -> PathBuf { fn reader(p: &Path) -> Box { let f = File::open(p).unwrap(); if p.extension().is_some_and(|e| e == "gz") { - Box::new(BufReader::new(GzDecoder::new(f))) + Box::new(BufReader::new(MultiGzDecoder::new(f))) } else { Box::new(BufReader::new(f)) } diff --git a/src/io/fastq.rs b/src/io/fastq.rs index 0c795de..93c48d9 100644 --- a/src/io/fastq.rs +++ b/src/io/fastq.rs @@ -1,6 +1,6 @@ /// FASTQ reader with base encoding and decompression support use crate::error::Error; -use flate2::read::GzDecoder; +use flate2::read::MultiGzDecoder; use noodles::fastq; use std::fs::File; use std::io::{BufRead, BufReader, BufWriter, Write}; @@ -116,7 +116,7 @@ impl FastqReader { let buffered = BufReader::with_capacity(DECODE_BUF, file); Box::new(BufReader::with_capacity( DECODE_BUF, - GzDecoder::new(buffered), + MultiGzDecoder::new(buffered), )) } else { // Plain text FASTQ @@ -569,6 +569,54 @@ mod tests { assert_eq!(read1.quality.len(), 4); } + /// A `.gz` written as several concatenated gzip members — what `bcl2fastq` + /// emits, what `cat a.fq.gz b.fq.gz` produces, and what every BGZF file is. + /// `flate2::read::GzDecoder` stops after the first member and reports EOF, + /// so reading such a file used to drop reads with no error at all. + #[test] + fn test_fastq_reader_gzip_multi_member() { + use flate2::Compression; + use flate2::write::GzEncoder; + + let mut tmpfile = tempfile::Builder::new() + .suffix(".fastq.gz") + .tempfile() + .unwrap(); + + // Member 1: read1. Each `finish()` closes a complete gzip stream, so + // the next encoder appends a second member rather than continuing. + let mut encoder = GzEncoder::new(Vec::new(), Compression::default()); + writeln!(encoder, "@read1").unwrap(); + writeln!(encoder, "ACGT").unwrap(); + writeln!(encoder, "+").unwrap(); + writeln!(encoder, "IIII").unwrap(); + tmpfile.write_all(&encoder.finish().unwrap()).unwrap(); + + // Member 2: read2. + let mut encoder = GzEncoder::new(Vec::new(), Compression::default()); + writeln!(encoder, "@read2").unwrap(); + writeln!(encoder, "TGCA").unwrap(); + writeln!(encoder, "+").unwrap(); + writeln!(encoder, "HHHH").unwrap(); + tmpfile.write_all(&encoder.finish().unwrap()).unwrap(); + tmpfile.flush().unwrap(); + + let mut reader = FastqReader::open(tmpfile.path(), None).unwrap(); + + let read1 = reader.next_encoded().unwrap().unwrap(); + assert_eq!(read1.name, "read1"); + assert_eq!(read1.sequence, vec![0, 1, 2, 3]); // ACGT + + let read2 = reader + .next_encoded() + .unwrap() + .expect("second gzip member must be decoded, not silently truncated"); + assert_eq!(read2.name, "read2"); + assert_eq!(read2.sequence, vec![3, 2, 1, 0]); // TGCA + + assert!(reader.next_encoded().unwrap().is_none()); + } + #[test] fn test_strip_mate_suffix_slash() { assert_eq!(strip_mate_suffix("read123/1"), "read123"); diff --git a/src/ruSTAR.code-workspace b/src/ruSTAR.code-workspace new file mode 100644 index 0000000..3e9ebe5 --- /dev/null +++ b/src/ruSTAR.code-workspace @@ -0,0 +1,11 @@ +{ + "folders": [ + { + "path": "../../STAR-rs" + }, + { + "path": ".." + } + ], + "settings": {} +} \ No newline at end of file diff --git a/src/solo/count.rs b/src/solo/count.rs index 4ddb414..cda22e7 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -2320,7 +2320,9 @@ fn open_maybe_gz(path: &Path) -> Result, Error> { .extension() .is_some_and(|e| e.eq_ignore_ascii_case("gz")) { - Ok(Box::new(BufReader::new(flate2::read::GzDecoder::new(file)))) + Ok(Box::new(BufReader::new(flate2::read::MultiGzDecoder::new( + file, + )))) } else { Ok(Box::new(BufReader::new(file))) } diff --git a/src/solo/whitelist.rs b/src/solo/whitelist.rs index 1d882ae..8c87c23 100644 --- a/src/solo/whitelist.rs +++ b/src/solo/whitelist.rs @@ -13,7 +13,7 @@ use crate::error::Error; use crate::io::fastq::{decode_base, encode_base}; -use flate2::read::GzDecoder; +use flate2::read::MultiGzDecoder; use std::fs::File; use std::io::{BufRead, BufReader}; use std::path::Path; @@ -535,7 +535,7 @@ fn open_maybe_gzip(path: &Path) -> Result, Error> { .extension() .is_some_and(|e| e.eq_ignore_ascii_case("gz")); if is_gz { - Ok(Box::new(BufReader::new(GzDecoder::new(file)))) + Ok(Box::new(BufReader::new(MultiGzDecoder::new(file)))) } else { Ok(Box::new(BufReader::new(file))) } From eaa12abaae3f10bbf04f07b926489af2cb71e48b Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 11 Aug 2026 20:53:23 +0200 Subject: [PATCH 2/2] feat(io): optional rapidgzip feature for parallel gzip input decoding Adds `rapidgzip-core` 0.3 (BSD-3-Clause AND MIT) behind a cargo feature, for parallel decoding of gzipped `--readFilesIn`. It is pure Rust. Its dependencies are `bon`, `crossbeam-deque` and `libz-rs-sys`, and that last one is zlib-rs, the same inflate backend `flate2` already uses in this project, so the feature adds no C or C++ toolchain requirement. It implements the rapidgzip marker/window algorithm: DEFLATE blocks are decoded speculatively before their back-references are known, then patched once they are. Decoding a 73 MB level-6 .fq.gz (427 MB out) on 16 logical cores, median of three: decoder throughput flate2 + zlib-rs, 1 thread 1552 MB/s rapidgzip-core, 1 thread 1581 MB/s rapidgzip-core, 4 threads 1872 MB/s rapidgzip-core, 8 threads 3356 MB/s rapidgzip-core, 12 threads 4051 MB/s Note there is no penalty at one thread, unlike an FFI decoder. It is nonetheless off by default, and inert until RUSTAR_GZ_DECODE_THREADS is set, for a reason that is arithmetic rather than doubt about the decoder: aligning 2 M reads takes ~3.0 s at --runThreadN 8, so we consume decompressed FASTQ at ~140 MB/s while one flate2 thread supplies ~1550 MB/s. Decode runs at under a tenth of its capacity. Turning it on measured 3.15 s against 3.03 s and ~360 MB more resident, which is what spending threads on an idle stage looks like. The feature becomes the right call when the consumer can drain a single inflate thread, roughly 7 M reads/s for this file shape. Implementation notes: * `Decoder::open` owns the file and returns a `Read + Send` stream, so it drops straight into the `Box` the parser wants. Reaching EOF verifies every member footer, so error behaviour on a corrupt file is not weakened relative to flate2. * Any failure to configure or open through it logs a warning and falls back to flate2, so this can never turn a working run into a failing one. * The worker count is process-wide state set once from `run()`, rather than a parameter, because `FastqReader::open` is called from a dozen places that have no business each deciding a decompression policy. `test_rapidgzip_decodes_identically_to_flate2` reads the same two-member gzip down both paths and asserts the decoded records match. End to end, the BAM payload is byte-identical with the feature on and off. Refs #224, #223. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 9 +++ Cargo.lock | 94 ++++++++++++++++++++++++-- Cargo.toml | 11 +++ README.md | 20 ++++++ src/io/fastq.rs | 176 +++++++++++++++++++++++++++++++++++++++++++++--- src/lib.rs | 4 ++ 6 files changed, 302 insertions(+), 12 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2656ba6..f1e4297 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,15 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Other changes +- **Optional `rapidgzip` feature** for parallel gzip/BGZF decoding of + `--readFilesIn`, backed by the pure-Rust `rapidgzip-core` (inflate backend is + zlib-rs, the same one `flate2` uses here, so it adds no C toolchain). Off by + default and inert until `RUSTAR_GZ_DECODE_THREADS` is set: the decoder reaches + 3356 MB/s on 8 threads against 1552 MB/s for one `flate2` thread, but the + aligner only consumes ~140 MB/s of decompressed FASTQ, so decode was never the + bottleneck and enabling it measured slightly slower with ~360 MB more resident. + Output is byte-identical to the `flate2` path. + - `cluster_seeds` reuses its window-bin map across reads on a thread instead of rebuilding it per read. Merging two windows re-keys every bin in the merged span, so the per-read pre-sizing was only a floor and the map diff --git a/Cargo.lock b/Cargo.lock index 9f77c86..118dbb7 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -119,6 +119,31 @@ version = "2.13.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b588b76d00fde79687d7646a9b5bdf3cc0f655e0bbd080335a95d7e96f3587da" +[[package]] +name = "bon" +version = "3.9.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a602c73c7b0148ec6d12af6fd5cc7a46e2eacc8878271a999abac56eed12f561" +dependencies = [ + "bon-macros", + "rustversion", +] + +[[package]] +name = "bon-macros" +version = "3.9.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6dee98b0db6a962de883bf5d20362dee4d7ca0d12fe39a7c6c73c844e1cd7c1f" +dependencies = [ + "darling", + "ident_case", + "prettyplease", + "proc-macro2", + "quote", + "rustversion", + "syn 2.0.117", +] + [[package]] name = "borsh" version = "1.8.0" @@ -280,9 +305,9 @@ dependencies = [ [[package]] name = "crossbeam-deque" -version = "0.8.6" +version = "0.8.7" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9dd111b7b7f7d55b72c0a6ae361660ee5853c9af73f70c3c2ef6858b950e2e51" +checksum = "5181e0de7b61eb03a81e347d6dd8797bae9da5146707b51077e2d71a54ec0ceb" dependencies = [ "crossbeam-epoch", "crossbeam-utils", @@ -309,6 +334,40 @@ version = "0.2.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b365fabc795046672053e29c954733ec3b05e4be654ab130fe8f1f94d7051f35" +[[package]] +name = "darling" +version = "0.23.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "25ae13da2f202d56bd7f91c25fba009e7717a1e4a1cc98a76d844b65ae912e9d" +dependencies = [ + "darling_core", + "darling_macro", +] + +[[package]] +name = "darling_core" +version = "0.23.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9865a50f7c335f53564bb694ef660825eb8610e0a53d3e11bf1b0d3df31e03b0" +dependencies = [ + "ident_case", + "proc-macro2", + "quote", + "strsim", + "syn 2.0.117", +] + +[[package]] +name = "darling_macro" +version = "0.23.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ac3984ec7bd6cfa798e62b4a642426a5be0e68f9401cfc2a01e3fa9ea2fcdb8d" +dependencies = [ + "darling_core", + "quote", + "syn 2.0.117", +] + [[package]] name = "dashmap" version = "6.2.1" @@ -481,6 +540,12 @@ version = "2.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "3d3067d79b975e8844ca9eb072e16b31c3c1c36928edf9c6789548c524d0d954" +[[package]] +name = "ident_case" +version = "1.0.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b9e0384b61958566e926dc50660321d12159025e767c18e043daf26b70104c39" + [[package]] name = "indexmap" version = "2.13.0" @@ -636,6 +701,15 @@ dependencies = [ "cty", ] +[[package]] +name = "libz-rs-sys" +version = "0.6.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "03dcace986b149f29509af6ca70e6182bccce916b644424ecf484faa8ddc899a" +dependencies = [ + "zlib-rs", +] + [[package]] name = "linux-raw-sys" version = "0.12.1" @@ -907,6 +981,17 @@ version = "6.0.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "f8dcc9c7d52a811697d2151c701e0d08956f92b0e24136cf4cf27b57a6a0d9bf" +[[package]] +name = "rapidgzip-core" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "da5be74bd353811db3efff881e87200a18150d72fb0b1d39ea04e7972a5b4ca1" +dependencies = [ + "bon", + "crossbeam-deque", + "libz-rs-sys", +] + [[package]] name = "rayon" version = "1.12.0" @@ -988,6 +1073,7 @@ dependencies = [ "noodles", "noodles-bgzf", "predicates", + "rapidgzip-core", "rayon", "rustc-hash", "shlex", @@ -1481,9 +1567,9 @@ dependencies = [ [[package]] name = "zlib-rs" -version = "0.6.3" +version = "0.6.7" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3be3d40e40a133f9c916ee3f9f4fa2d9d63435b5fbe1bfc6d9dae0aa0ada1513" +checksum = "34b31d188d9d685a4f9c7b46d6e36631b07058d2cfe190267adce54dc230bf12" [[package]] name = "zmij" diff --git a/Cargo.toml b/Cargo.toml index 8a4638f..c714d4c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -64,8 +64,19 @@ caps-sa = "0.6" mimalloc = { version = "0.1", default-features = false } libmimalloc-sys = { version = "0.1.49", features = ["extended"] } # mi_option_set (purge_delay); see main.rs libdeflater = "1.25.2" +# Parallel gzip/BGZF/zlib decoder for `--readFilesIn *.gz`. Pure Rust: its only +# dependencies are `bon`, `crossbeam-deque` and `libz-rs-sys` (zlib-rs), which is +# the same inflate backend `flate2` already uses here, so it adds no C toolchain +# requirement. Behind a feature only until the dependency discussion in #224 +# settles. See `set_gz_decode_threads` for the measurements. +rapidgzip-core = { version = "0.3", optional = true } noodles-bgzf = { version = "0.49", features = ["libdeflate"] } +[features] +default = [] +# Parallel gzip input decoding (pure Rust). +rapidgzip = ["dep:rapidgzip-core"] + [dev-dependencies] assert_cmd = "2" predicates = "3" diff --git a/README.md b/README.md index d6536e6..136afa4 100644 --- a/README.md +++ b/README.md @@ -257,6 +257,26 @@ cargo clippy --all-targets # Lint cargo fmt # Format ``` +### Optional feature: `rapidgzip` + +```bash +cargo build --release --features rapidgzip +``` + +Decodes gzipped `--readFilesIn` with [rapidgzip-core](https://crates.io/crates/rapidgzip-core), a +pure-Rust parallel gzip decoder, instead of a single `flate2` thread. It pulls in no C toolchain: +its inflate backend is zlib-rs, the same one `flate2` already uses here. + +**Off by default, and inert even when compiled in** until `RUSTAR_GZ_DECODE_THREADS` is set to a +worker count (`0` keeps `flate2`). The decoder is fast (3356 MB/s on 8 threads against 1552 MB/s for +one `flate2` thread on the same file), but this aligner consumes decompressed FASTQ at only +~140 MB/s, so decode already runs at under a tenth of its capacity. Enabling it measured slightly +slower and used ~360 MB more memory. It is worth turning on when the consumer is fast enough to +drain a single inflate thread, roughly 7 M reads/s for a 100 bp library. + +The portable alternative needs no feature and no rebuild: pipe an external decompressor with +`--readFilesCommand` (`gzcat`, `gunzip -c`, `igzip -dc`, `rapidgzip-rust -dc`). + ## Development The majority of rustar-aligner's code was written by [Claude Code](https://claude.ai/code) (Anthropic's AI coding assistant), with technical direction, architecture decisions, and validation by the project maintainer. diff --git a/src/io/fastq.rs b/src/io/fastq.rs index 93c48d9..ea7baa6 100644 --- a/src/io/fastq.rs +++ b/src/io/fastq.rs @@ -69,6 +69,66 @@ pub struct PairedRead { pub mate2: EncodedRead, } +/// Worker count handed to `rapidgzip-core` for gzip input, or 0 for the +/// single-threaded `flate2` path. Set once per run by +/// [`set_gz_decode_threads`]; read by every `FastqReader::open`. +/// +/// A process-wide value rather than a parameter because `FastqReader::open` is +/// called from a dozen places (both mates, solo cDNA/barcode files, SmartSeq +/// manifests) that have no business each deciding a decompression policy. +#[cfg(feature = "rapidgzip")] +static GZ_DECODE_THREADS: std::sync::atomic::AtomicU32 = std::sync::atomic::AtomicU32::new(0); + +/// Decide how many threads gzip decode may use. **Off unless +/// `RUSTAR_GZ_DECODE_THREADS` asks for it**, because at this aligner's current +/// consumption rate parallel decode has nothing to win. +/// +/// `rapidgzip-core` decodes a gzip stream in parallel with the marker/window +/// algorithm: blocks are decoded speculatively before their back-references are +/// known, then patched once they are. Its inflate backend is `libz-rs-sys` +/// (zlib-rs), the same one `flate2` uses here, so unlike an FFI decoder it costs +/// nothing at one thread. Decoding a 73 MB level-6 `.fq.gz` (427 MB out) on 16 +/// logical cores, median of three: +/// +/// | decoder | throughput | +/// |---------|------------| +/// | `flate2` + `zlib-rs`, 1 thread | 1552 MB/s | +/// | `rapidgzip-core`, 1 thread | 1581 MB/s | +/// | `rapidgzip-core`, 4 threads | 1872 MB/s | +/// | `rapidgzip-core`, 8 threads | 3356 MB/s | +/// | `rapidgzip-core`, 12 threads | 4051 MB/s | +/// +/// The reason it is still off by default is arithmetic, not doubt about the +/// decoder. Aligning those 2 M reads takes ~3.0 s at `--runThreadN 8`, i.e. we +/// consume 427 MB of decompressed FASTQ at ~140 MB/s. One `flate2` thread +/// supplies ~1550 MB/s, so decode runs at under a tenth of its capacity and is +/// nowhere near the critical path; turning it on measured *slower* (3.15 s vs +/// 3.03 s) and ~360 MB more resident, which is what you would expect from +/// spending threads on a stage that was already idle. +/// +/// It becomes the right call when the consumer gets fast enough to drain a +/// single inflate thread, which for this file shape means roughly 7 M reads/s. +/// That is the regime tools like piscem/salmon are in, and it is why they pair +/// a parallel decoder with a broker that splits one thread budget by measured +/// busy time on each side rather than by a fixed ratio. +#[cfg(feature = "rapidgzip")] +pub fn set_gz_decode_threads(_run_thread_n: u32) { + let threads = std::env::var("RUSTAR_GZ_DECODE_THREADS") + .ok() + .and_then(|s| s.parse::().ok()) + .unwrap_or(0); + GZ_DECODE_THREADS.store(threads, std::sync::atomic::Ordering::Relaxed); +} + +/// No-op when the feature is off, so callers need no `cfg`. +#[cfg(not(feature = "rapidgzip"))] +pub fn set_gz_decode_threads(_run_thread_n: u32) {} + +#[cfg(feature = "rapidgzip")] +fn gz_decode_threads() -> u32 { + GZ_DECODE_THREADS.load(std::sync::atomic::Ordering::Relaxed) +} + /// FASTQ reader that handles decompression and base encoding pub struct FastqReader { inner: fastq::io::Reader>, @@ -105,21 +165,24 @@ impl FastqReader { let path_str = path.to_string_lossy(); let is_gzipped = path_str.ends_with(".gz") || path_str.ends_with(".gzip"); - let file = File::open(path).map_err(|e| Error::io(e, path))?; - // Larger-than-default (8 KiB) buffers cut read syscalls on the decode // hot path. Feed the inflater from a big buffered file, and hand the // decoded stream to noodles through a big BufReader. const DECODE_BUF: usize = 1 << 19; // 512 KiB if is_gzipped { - // Gzipped file - let buffered = BufReader::with_capacity(DECODE_BUF, file); - Box::new(BufReader::with_capacity( - DECODE_BUF, - MultiGzDecoder::new(buffered), - )) + if let Some(reader) = Self::open_gz_parallel(path) { + reader + } else { + let file = File::open(path).map_err(|e| Error::io(e, path))?; + let buffered = BufReader::with_capacity(DECODE_BUF, file); + Box::new(BufReader::with_capacity( + DECODE_BUF, + MultiGzDecoder::new(buffered), + )) + } } else { // Plain text FASTQ + let file = File::open(path).map_err(|e| Error::io(e, path))?; Box::new(BufReader::with_capacity(DECODE_BUF, file)) } }; @@ -152,6 +215,57 @@ impl FastqReader { self } + /// Parallel gzip/BGZF decode via `rapidgzip`, or `None` to use the + /// single-threaded `flate2` path. + /// + /// Returns `None` unless the `rapidgzip` feature is compiled in **and** + /// [`gz_decode_threads`] resolves to at least 2. Any failure to open the + /// file through the native decoder also returns `None`, so a run degrades + /// to the `flate2` path rather than failing: this is an optimisation, never + /// a requirement. + #[cfg(feature = "rapidgzip")] + fn open_gz_parallel(path: &Path) -> Option> { + let threads = gz_decode_threads(); + if threads < 2 { + return None; + } + // `Decoder::open` owns the file and returns a `Read + Send` stream, so it + // drops straight into the `Box` the parser wants. + // Reaching EOF verifies every member footer, which the `flate2` path + // does too, so error behaviour on a corrupt file is not weakened. + let decoder = match rapidgzip_core::Decoder::builder() + .decoder_threads(threads as usize) + .build() + { + Ok(d) => d, + Err(e) => { + log::warn!("rapidgzip-core could not be configured: {e}; using flate2"); + return None; + } + }; + match decoder.open(path) { + Ok(reader) => { + log::info!( + "decompressing {} with rapidgzip-core on {threads} threads", + path.display() + ); + Some(Box::new(BufReader::with_capacity(1 << 19, reader))) + } + Err(e) => { + log::warn!( + "rapidgzip-core could not open {}: {e}; falling back to flate2", + path.display() + ); + None + } + } + } + + #[cfg(not(feature = "rapidgzip"))] + fn open_gz_parallel(_path: &Path) -> Option> { + None + } + /// Open FASTQ file using external decompression command fn open_with_command(path: &Path, cmd: &str) -> Result, Error> { let mut child = Command::new(cmd) @@ -617,6 +731,52 @@ mod tests { assert!(reader.next_encoded().unwrap().is_none()); } + /// The parallel decoder must produce exactly what the `flate2` path + /// produces, including across gzip member boundaries — the case that used + /// to truncate. Reads the same file twice, once down each path. + #[cfg(feature = "rapidgzip")] + #[test] + fn test_rapidgzip_decodes_identically_to_flate2() { + use flate2::Compression; + use flate2::write::GzEncoder; + use std::sync::atomic::Ordering; + + let mut tmpfile = tempfile::Builder::new() + .suffix(".fastq.gz") + .tempfile() + .unwrap(); + // Two members, several reads each, so the parallel decoder has both a + // member boundary and enough data to split on. + for member in 0..2 { + let mut encoder = GzEncoder::new(Vec::new(), Compression::default()); + for i in 0..64 { + writeln!(encoder, "@m{member}_read{i}").unwrap(); + writeln!(encoder, "ACGTACGTNN").unwrap(); + writeln!(encoder, "+").unwrap(); + writeln!(encoder, "IIIIIIIIII").unwrap(); + } + tmpfile.write_all(&encoder.finish().unwrap()).unwrap(); + } + tmpfile.flush().unwrap(); + + let collect = |threads: u32| { + GZ_DECODE_THREADS.store(threads, Ordering::Relaxed); + let mut reader = FastqReader::open(tmpfile.path(), None).unwrap(); + let mut out = Vec::new(); + while let Some(read) = reader.next_encoded().unwrap() { + out.push((read.name, read.sequence, read.quality)); + } + out + }; + + let via_flate2 = collect(0); + let via_rapidgzip = collect(4); + GZ_DECODE_THREADS.store(0, Ordering::Relaxed); + + assert_eq!(via_flate2.len(), 128, "both members must be decoded"); + assert_eq!(via_rapidgzip, via_flate2); + } + #[test] fn test_strip_mate_suffix_slash() { assert_eq!(strip_mate_suffix("read123/1"), "read123"); diff --git a/src/lib.rs b/src/lib.rs index 5086fcb..584e42e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -76,6 +76,10 @@ pub fn run(params: &Parameters) -> anyhow::Result<()> { .num_threads(params.run_thread_n.into()) .build_global(); + // Decide the gzip-decode policy once, from the same thread budget. No-op + // unless the `rapidgzip` feature is compiled in. + crate::io::fastq::set_gz_decode_threads(params.run_thread_n.get() as u32); + match params.run_mode() { RunMode::GenomeGenerate => genome_generate(params), RunMode::AlignReads => align_reads(params),