From aafabac1ca49e531a3d77186fcdf3197919fa742 Mon Sep 17 00:00:00 2001 From: alejandrogzi Date: Thu, 6 Aug 2026 23:01:50 +0200 Subject: [PATCH 1/3] feat(io): add native CBQ input support Add single-file CBQ read input (single-end or interleaved paired-end), detected by file header rather than extension, with parallel decoding controlled by '--readFilesNthreads'. Decoded reads are bounded by a fixed per-window record budget so peak input memory is independent of thread count and block size. Read types and base helpers shared with the FASTQ path move to a new 'reads' module. --- Cargo.lock | 215 +++- Cargo.toml | 1 + README.md | 21 + .../content/docs/reference/cli-parameters.md | 3 +- src/align/seed.rs | 2 +- src/chimeric/output.rs | 2 +- src/clip/mod.rs | 2 +- src/error.rs | 43 + src/genome/transform.rs | 2 +- src/io/cbq.rs | 1102 +++++++++++++++++ src/io/fastq.rs | 129 +- src/io/input.rs | 468 +++++++ src/io/mod.rs | 3 + src/io/reads.rs | 99 ++ src/io/sam.rs | 81 +- src/lib.rs | 142 ++- src/params/mod.rs | 14 + src/solo/count.rs | 2 +- src/solo/mod.rs | 7 +- src/solo/whitelist.rs | 2 +- src/wasp/mod.rs | 2 +- tests/alignment_features.rs | 278 +++++ 22 files changed, 2402 insertions(+), 218 deletions(-) create mode 100644 src/io/cbq.rs create mode 100644 src/io/input.rs create mode 100644 src/io/reads.rs diff --git a/Cargo.lock b/Cargo.lock index 9f77c86..3b6e39e 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -97,12 +97,44 @@ 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.117", +] + [[package]] name = "autocfg" version = "1.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "c08606f8c3cbf4ce6ec8e28fb0014a2c086708fe954eaa885384a6165172e7e8" +[[package]] +name = "binseq" +version = "0.9.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c03ce47e3c043402933b71d00e26586b91f70d1cc237251097e2eb1ee9621d83" +dependencies = [ + "anyhow", + "auto_impl", + "bitnuc", + "bytemuck", + "byteorder", + "itoa", + "memchr", + "memmap2", + "num_cpus", + "rand", + "sucds", + "thiserror", + "zstd", +] + [[package]] name = "bit-vec" version = "0.10.1" @@ -119,6 +151,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 = "borsh" version = "1.8.0" @@ -160,6 +198,26 @@ version = "3.19.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "5dd9dc738b7a8311c7ade152424974d8115f2cdad61e8dab8dac9f2362298510" +[[package]] +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.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f65693059b6b9c588b9f62fed1cedbf0a8b805631457ea162d68f0de186f3de5" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.117", +] + [[package]] name = "byteorder" version = "1.5.0" @@ -189,6 +247,8 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "556e016178bb5662a08681bbe0f00f8e17631781a4dfc8c45e466e4b185ec27f" dependencies = [ "find-msvc-tools", + "jobserver", + "libc", "shlex", ] @@ -411,6 +471,18 @@ version = "0.1.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d9c4f5dac5e15c24eb999c26181a6ca40b39fe946cbe4c263c7209467bc83af2" +[[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.2" @@ -419,7 +491,7 @@ checksum = "0de51e6874e94e7bf76d726fc5d13ba782deca734ff60d5bb2fb2607c7406555" dependencies = [ "cfg-if", "libc", - "r-efi", + "r-efi 6.0.0", "wasip2", "wasip3", ] @@ -451,6 +523,12 @@ version = "0.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" +[[package]] +name = "hermit-abi" +version = "0.5.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fc0fef456e4baa96da950455cd02c081ca953b141298e41db3fc7e36b1da849c" + [[package]] name = "iana-time-zone" version = "0.1.65" @@ -529,6 +607,16 @@ dependencies = [ "syn 2.0.117", ] +[[package]] +name = "jobserver" +version = "0.1.35" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1c00acbd29eabad4a2392fa0e921c874934dbbf4194312ad20f04a0ed67a3cb3" +dependencies = [ + "getrandom 0.4.2", + "libc", +] + [[package]] name = "js-sys" version = "0.3.85" @@ -659,9 +747,9 @@ checksum = "0ceec5bc11778974d1bcb055b18002eba7f4b3518b6a0081b3af5f21666da9ad" [[package]] name = "memchr" -version = "2.7.6" +version = "2.8.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f52b00d39961fc5b2736ea853c9cc86238e165017a493d1d5c8eac6bdc4cc273" +checksum = "cf8baf1c55e62ffcace7a9f06f4bd9cd3f0c4beb022d3b367256b91b87513d98" [[package]] name = "memmap2" @@ -794,6 +882,16 @@ dependencies = [ "autocfg", ] +[[package]] +name = "num_cpus" +version = "1.17.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91df4bbde75afed763b708b7eee1e8e7651e02d97f6d5dd763e89367e957b23b" +dependencies = [ + "hermit-abi", + "libc", +] + [[package]] name = "once_cell" version = "1.21.4" @@ -819,6 +917,12 @@ dependencies = [ "windows-link", ] +[[package]] +name = "pkg-config" +version = "0.3.33" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "19f132c84eca552bf34cab8ec81f1c1dcc229b811638f9d283dceabe58c5569e" + [[package]] name = "portable-atomic" version = "1.13.1" @@ -834,6 +938,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" @@ -901,12 +1014,47 @@ 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", +] + +[[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", +] + +[[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]] name = "rayon" version = "1.12.0" @@ -971,6 +1119,7 @@ version = "0.2.0" dependencies = [ "anyhow", "assert_cmd", + "binseq", "bitflags", "bstr", "byteorder", @@ -1099,6 +1248,16 @@ version = "0.11.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "7da8b5736845d9f2fcb837ea5d9e2628564b3b043a70948a3f0b778838c5fb4f" +[[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.117" @@ -1128,7 +1287,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "32497e9a4c7b38532efcdebeef879707aa9f794296a4f0244f6f69e9bc8574bd" dependencies = [ "fastrand", - "getrandom", + "getrandom 0.4.2", "once_cell", "rustix", "windows-sys", @@ -1479,6 +1638,26 @@ dependencies = [ "wasmparser", ] +[[package]] +name = "zerocopy" +version = "0.8.55" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b5a105cd7b140f6eeec8acff2ea38135d3cab283ada58540f629fe51e46696eb" +dependencies = [ + "zerocopy-derive", +] + +[[package]] +name = "zerocopy-derive" +version = "0.8.55" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0fe976fb70c78cd64cccfe3a6fc142244e8a77b70959b30faf9d0ac37ee228eb" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.117", +] + [[package]] name = "zlib-rs" version = "0.6.3" @@ -1490,3 +1669,31 @@ name = "zmij" version = "1.0.21" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b8848ee67ecc8aedbaf3e4122217aff892639231befc6a1b58d29fff4c2cabaa" + +[[package]] +name = "zstd" +version = "0.13.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e91ee311a569c327171651566e07972200e76fcfe2242a4fa446149a3881c08a" +dependencies = [ + "zstd-safe", +] + +[[package]] +name = "zstd-safe" +version = "7.2.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8f49c4d5f0abb602a93fb8736af2a4f4dd9512e36f7f570d66e65ff867ed3b9d" +dependencies = [ + "zstd-sys", +] + +[[package]] +name = "zstd-sys" +version = "2.0.16+zstd.1.5.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91e19ebc2adc8f83e43039e79776e3fda8ca919132d68a1fed6a5faca2683748" +dependencies = [ + "cc", + "pkg-config", +] diff --git a/Cargo.toml b/Cargo.toml index 8a4638f..30fa57d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -65,6 +65,7 @@ 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" noodles-bgzf = { version = "0.49", features = ["libdeflate"] } +binseq = { version = "0.9.4", default-features = false, features = ["anyhow"] } [dev-dependencies] assert_cmd = "2" diff --git a/README.md b/README.md index d6536e6..cb60d83 100644 --- a/README.md +++ b/README.md @@ -45,6 +45,26 @@ target/release/rustar-aligner \ --outFileNamePrefix /path/to/output_ ``` +### CBQ input + +CBQ is detected from the file header, so the filename extension is irrelevant. +A paired-end CBQ stores both mates in one file. +CBQ is currently limited to ordinary alignment input; it cannot be combined +with `--readFilesCommand`, STARsolo/SmartSeq workflows, or qualityless +unmapped-FASTQ output. + +`--readFilesNthreads` controls decoder concurrency only (`0` derives a bounded +default from `--runThreadN`); it does not affect how much input is held in +memory, which is fixed by the window budget. + +```bash +target/release/rustar-aligner \ + --genomeDir /path/to/genome_index \ + --readFilesIn reads.cbq \ + --readFilesNthreads 4 \ + --outFileNamePrefix /path/to/output_ +``` + ### BAM output ```bash @@ -207,6 +227,7 @@ resident; the 16 GB sparse index is stable at ~54 s. ## Supported Features - Single-end and paired-end alignment with mate rescue +- Native single-file CBQ input (single-end or interleaved paired-end), with parallel decoding controlled by `--readFilesNthreads`. Decoded reads are bounded by a fixed per-window record budget, so peak input memory is independent of both the thread count and the input file's block size - Read-end alignment mode (`--alignEndsType Local` (default) / `EndToEnd` / `Extend5pOfRead1` / `Extend5pOfReads12` / `Extend3pOfRead1`) - SAM, unsorted BAM, and coordinate-sorted BAM output (`--outSAMtype SAM`, `BAM Unsorted`, or `BAM SortedByCoordinate`) - Multi-threaded parallel alignment (`--runThreadN`) diff --git a/docs/src/content/docs/reference/cli-parameters.md b/docs/src/content/docs/reference/cli-parameters.md index f749039..d4909d4 100644 --- a/docs/src/content/docs/reference/cli-parameters.md +++ b/docs/src/content/docs/reference/cli-parameters.md @@ -29,8 +29,9 @@ Run `rustar-aligner --help` for the full machine-generated listing. | Parameter | Default | Description | |-----------|---------|-------------| -| `--readFilesIn` | — | Input FASTQ file(s); second file is mate 2 for paired-end (required for `alignReads`). | +| `--readFilesIn` | — | Input FASTQ file(s), or one CBQ file. A second FASTQ is mate 2; paired-end CBQ stores both mates in one file (required for `alignReads`). | | `--readFilesCommand` | — | Decompression command, e.g. `zcat` for `.gz`. | +| `--readFilesNthreads` | `0` | CBQ decoder threads (`0` = derive a bounded default from `--runThreadN`; ignored for FASTQ). | | `--readMapNumber` | `-1` | Number of reads to map (`-1` = all). | | `--clip5pNbases` | `0` | Bases to clip from the 5' end of each mate. | | `--clip3pNbases` | `0` | Bases to clip from the 3' end of each mate. | diff --git a/src/align/seed.rs b/src/align/seed.rs index 7c6deb3..6e0be56 100644 --- a/src/align/seed.rs +++ b/src/align/seed.rs @@ -1,6 +1,6 @@ use crate::error::Error; use crate::index::GenomeIndex; -use crate::io::fastq::complement_base; +use crate::io::reads::complement_base; use crate::params::Parameters; /// A seed represents an exact match between a read position and genome location(s). diff --git a/src/chimeric/output.rs b/src/chimeric/output.rs index 627c5a8..452da20 100644 --- a/src/chimeric/output.rs +++ b/src/chimeric/output.rs @@ -193,7 +193,7 @@ fn build_segment_record( is_supplementary: bool, sa_tag: &str, ) -> Result { - use crate::io::fastq::{complement_base, decode_base}; + use crate::io::reads::{complement_base, decode_base}; use noodles::sam::alignment::record::data::field::Tag; let mut record = RecordBuf::default(); diff --git a/src/clip/mod.rs b/src/clip/mod.rs index 20901ec..c7d69ee 100644 --- a/src/clip/mod.rs +++ b/src/clip/mod.rs @@ -21,7 +21,7 @@ //! adapter is not a thing STAR itself supports either (only `CellRanger4` mode //! clips a 5' adapter, the 10x TSO) — that mode is out of scope here. -use crate::io::fastq::encode_base; +use crate::io::reads::encode_base; use crate::params::Parameters; /// One end's clipping parameters (STAR `ClipMate`). `n` = fixed clip, `adapter` = diff --git a/src/error.rs b/src/error.rs index dee967f..5f55c1d 100644 --- a/src/error.rs +++ b/src/error.rs @@ -26,6 +26,49 @@ pub enum Error { #[error("chimeric detection error: {0}")] Chimeric(String), + + #[error("read input error: {0}")] + ReadInput(String), + + #[error("failed to open CBQ input {path}: {source}")] + CbqOpen { + path: PathBuf, + #[source] + source: binseq::Error, + }, + + #[error("failed to decode CBQ input {path} records {start}..{end}: {source}")] + CbqDecode { + path: PathBuf, + start: usize, + end: usize, + #[source] + source: binseq::Error, + }, + + #[error("invalid CBQ record {record_index}{mate_context} in {path}: {message}")] + CbqRecord { + path: PathBuf, + record_index: usize, + mate_context: String, + message: String, + }, + + #[error( + "CBQ ordering error in {path}: expected record {expected_index}, observed {observed_index}" + )] + CbqOrdering { + path: PathBuf, + expected_index: usize, + observed_index: usize, + }, + + #[error("CBQ worker panicked while processing {path} records {start}..{end}")] + CbqWorkerPanic { + path: PathBuf, + start: usize, + end: usize, + }, } impl Error { diff --git a/src/genome/transform.rs b/src/genome/transform.rs index 328e4f0..cc09558 100644 --- a/src/genome/transform.rs +++ b/src/genome/transform.rs @@ -21,7 +21,7 @@ use std::collections::BTreeMap; use std::fmt::Write as _; -use crate::io::fastq::encode_base; +use crate::io::reads::encode_base; use super::compute_chr_starts; use super::fasta::Chromosome; diff --git a/src/io/cbq.rs b/src/io/cbq.rs new file mode 100644 index 0000000..ffe2689 --- /dev/null +++ b/src/io/cbq.rs @@ -0,0 +1,1102 @@ +//! Native CBQ decoding for the existing ordered alignment pipeline. + +use crate::error::Error; +use crate::io::reads::{EncodedRead, PairedRead, encode_base, strip_mate_suffix}; +use binseq::{BinseqRecord, ParallelProcessor, ParallelReader}; +use std::collections::BTreeMap; +use std::path::{Path, PathBuf}; +use std::sync::mpsc::{SyncSender, sync_channel}; + +/// Hard ceiling on decoded records held for one CBQ window. +/// +/// The window is the reorder buffer's memory limit (a smaller pending-map limit +/// could deadlock), so this constant — not the decoder-thread count and not the +/// input file's block size — is what bounds peak decoded-read memory. At ~330 B +/// per 150 bp single-end read that is ~10 MiB (~21 MiB for paired records). +/// +/// A single CBQ block is the unavoidable floor: binseq decodes a whole block per +/// callback, so a file whose blocks each exceed the budget is bounded by one +/// block rather than by this value. +/// +/// Windows must also be large enough to amortize a fixed per-window cost: binseq +/// loads the block that starts exactly at a window's exclusive end +/// (`cbq::MmapReader::process_parallel_range` selects blocks with +/// `iv_start <= range.end`), so every window decompresses one block it discards. +/// +/// Selected by measurement, not assumption — decoder-only median wall time for +/// 200k x 150 bp SE records in 1 MiB quality-bearing blocks (~4651 records per +/// block), median of 5, AMD Ryzen 7 5700X: +/// +/// ```text +/// window 1 thread 2 threads 4 threads 8 threads +/// 9 456 73.4 ms 60.5 ms 62.4 ms 61.1 ms +/// 16 384 68.3 42.9 44.9 44.1 +/// 32 768 61.4 43.0 29.9 29.6 +/// 49 152 58.4 43.5 33.0 28.7 +/// 65 536 57.6 41.7 33.6 27.3 +/// 200 000 53.9 44.1 30.7 28.2 +/// ``` +/// +/// Cost falls steeply to ~32k records and then flattens, so this sits at the knee: +/// best or statistically tied-best for the 1-4 decoder threads the automatic +/// policy can select, at half the peak memory of the next size up. Against the +/// previous `decoder_threads * 2` block policy on the same file: -20.1% at 1 +/// thread, -6.3% at 2, -14.1% at 4, +1.9% at 8 (noise; 8 is not selectable +/// automatically). +const MAX_RECORDS_PER_WINDOW: usize = 32_768; + +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +struct CbqWindow { + record_start: usize, + record_end: usize, + block_count: usize, +} + +#[derive(Debug)] +struct IndexedRecord { + index: usize, + record: T, +} + +#[derive(Debug)] +struct DecodedBlock { + first_index: usize, + end_index: usize, + records: Vec>, +} + +pub trait RecordConverter: Clone + Send + 'static { + type Output: Send + 'static; + + fn convert(&self, record: R, path: &Path) -> Result; +} + +#[derive(Clone, Copy)] +pub struct SingleEndConverter; + +impl RecordConverter for SingleEndConverter { + type Output = EncodedRead; + + fn convert(&self, record: R, path: &Path) -> Result { + let index = record_index(&record, path)?; + if record.is_paired() { + return Err(record_error( + path, + index, + None, + "paired record found in single-end CBQ", + )); + } + convert_mate(&record, path, index, 1, false) + } +} + +#[derive(Clone, Copy)] +pub struct PairedEndConverter; + +impl RecordConverter for PairedEndConverter { + type Output = PairedRead; + + fn convert(&self, record: R, path: &Path) -> Result { + let index = record_index(&record, path)?; + if !record.is_paired() { + return Err(record_error( + path, + index, + None, + "single-end record found in paired CBQ", + )); + } + + let mate1 = convert_mate(&record, path, index, 1, false)?; + let mate2 = convert_mate(&record, path, index, 2, true)?; + let name1 = strip_mate_suffix(&mate1.name); + let name2 = strip_mate_suffix(&mate2.name); + if name1 != name2 { + return Err(record_error( + path, + index, + None, + format!( + "paired read names do not match: '{}' vs '{}'", + mate1.name, mate2.name + ), + )); + } + + Ok(PairedRead { + name: name1, + mate1, + mate2, + }) + } +} + +struct BlockProcessor { + converter: C, + path: PathBuf, + current: Vec>, + sender: SyncSender>, +} + +impl Clone for BlockProcessor { + fn clone(&self) -> Self { + Self { + converter: self.converter.clone(), + path: self.path.clone(), + current: Vec::new(), + sender: self.sender.clone(), + } + } +} + +impl BlockProcessor { + fn flush_block(&mut self) -> binseq::Result<()> { + if self.current.is_empty() { + return Ok(()); + } + + let first_index = self.current[0].index; + let mut expected = first_index; + for indexed in &self.current { + if indexed.index != expected { + return Err(into_binseq_error(Error::CbqOrdering { + path: self.path.clone(), + expected_index: expected, + observed_index: indexed.index, + })); + } + expected = expected.checked_add(1).ok_or_else(|| { + into_binseq_error(Error::ReadInput("CBQ record index overflow".to_string())) + })?; + } + + let records = std::mem::take(&mut self.current); + self.sender + .send(DecodedBlock { + first_index, + end_index: expected, + records, + }) + .map_err(|_| { + std::io::Error::new( + std::io::ErrorKind::BrokenPipe, + "CBQ block coordinator disconnected", + ) + .into() + }) + } +} + +impl ParallelProcessor for BlockProcessor { + fn process_record(&mut self, record: R) -> binseq::Result<()> { + let index = record_index(&record, &self.path).map_err(into_binseq_error)?; + let owned = self + .converter + .convert(record, &self.path) + .map_err(into_binseq_error)?; + self.current.push(IndexedRecord { + index, + record: owned, + }); + Ok(()) + } + + fn on_batch_complete(&mut self) -> binseq::Result<()> { + self.flush_block() + } + + fn on_thread_complete(&mut self) -> binseq::Result<()> { + self.flush_block() + } +} + +pub struct CbqProducer { + path: PathBuf, + converter: C, + decoder_threads: usize, + /// Records per decode window; always `MAX_RECORDS_PER_WINDOW` in production. + /// Tests lower it to reach multi-window decoding on small fixtures. + window_records: usize, +} + +pub fn single_end_producer( + path: PathBuf, + decoder_threads: usize, +) -> CbqProducer { + CbqProducer { + path, + converter: SingleEndConverter, + decoder_threads, + window_records: MAX_RECORDS_PER_WINDOW, + } +} + +pub fn paired_end_producer( + path: PathBuf, + decoder_threads: usize, +) -> CbqProducer { + CbqProducer { + path, + converter: PairedEndConverter, + decoder_threads, + window_records: MAX_RECORDS_PER_WINDOW, + } +} + +impl CbqProducer { + pub fn produce( + self, + batch_size: usize, + max_records: usize, + sender: &SyncSender, Error>>, + ) -> Result<(), Error> { + let reader = open_reader(&self.path)?; + let total_to_decode = reader.num_records().min(max_records); + if total_to_decode == 0 { + let _ = sender.send(Ok(Vec::new())); + return Ok(()); + } + + let block_record_counts = read_block_record_counts(&reader, &self.path)?; + let windows = plan_windows(&block_record_counts, total_to_decode, self.window_records)?; + log::info!( + "CBQ input {}: {} records, {} decoder thread(s), {} bounded window(s) of \u{2264}{} records", + self.path.display(), + total_to_decode, + self.decoder_threads, + windows.len(), + windows + .iter() + .map(|window| window.record_end - window.record_start) + .max() + .unwrap_or(0) + ); + let mut output_batch = Vec::with_capacity(batch_size); + + for window in windows { + if !decode_window( + reader.clone(), + &self.path, + self.converter.clone(), + window, + self.decoder_threads.min(window.block_count).max(1), + batch_size, + &mut output_batch, + sender, + )? { + return Ok(()); + } + } + + if !output_batch.is_empty() + && sender + .send(Ok(std::mem::replace( + &mut output_batch, + Vec::with_capacity(batch_size), + ))) + .is_err() + { + return Ok(()); + } + let _ = sender.send(Ok(Vec::new())); + Ok(()) + } +} + +#[allow(clippy::too_many_arguments)] +fn decode_window( + reader: binseq::cbq::MmapReader, + path: &Path, + converter: C, + window: CbqWindow, + decoder_threads: usize, + batch_size: usize, + output_batch: &mut Vec, + output_sender: &SyncSender, Error>>, +) -> Result { + let (block_sender, block_receiver) = sync_channel(decoder_threads.max(1)); + let processor = BlockProcessor { + converter, + path: path.to_path_buf(), + current: Vec::new(), + sender: block_sender, + }; + let range = window.record_start..window.record_end; + let runner = std::thread::spawn(move || { + std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| { + reader.process_parallel_range(processor, decoder_threads, range) + })) + }); + + let mut coordinator = OrderedBlocks::new(path.to_path_buf(), window.record_start); + let mut coordinator_error = None; + let mut cancelled = false; + + while let Ok(block) = block_receiver.recv() { + if coordinator_error.is_some() || cancelled { + continue; + } + match coordinator.push(block) { + Ok(ready) => { + for record in ready { + output_batch.push(record); + if output_batch.len() == batch_size { + let full = std::mem::replace(output_batch, Vec::with_capacity(batch_size)); + if output_sender.send(Ok(full)).is_err() { + cancelled = true; + break; + } + } + } + } + Err(error) => coordinator_error = Some(error), + } + } + + let runner_result = runner.join().map_err(|_| Error::CbqWorkerPanic { + path: path.to_path_buf(), + start: window.record_start, + end: window.record_end, + })?; + + if cancelled { + return Ok(false); + } + if let Some(error) = coordinator_error { + return Err(error); + } + match runner_result { + Err(_) => Err(Error::CbqWorkerPanic { + path: path.to_path_buf(), + start: window.record_start, + end: window.record_end, + }), + Ok(Err(source)) => Err(Error::CbqDecode { + path: path.to_path_buf(), + start: window.record_start, + end: window.record_end, + source, + }), + Ok(Ok(())) => { + coordinator.finish(window.record_end)?; + Ok(true) + } + } +} + +struct OrderedBlocks { + path: PathBuf, + next_record_index: usize, + pending: BTreeMap>, +} + +impl OrderedBlocks { + fn new(path: PathBuf, next_record_index: usize) -> Self { + Self { + path, + next_record_index, + pending: BTreeMap::new(), + } + } + + fn push(&mut self, block: DecodedBlock) -> Result, Error> { + if block.records.is_empty() || block.first_index >= block.end_index { + return Err(Error::CbqOrdering { + path: self.path.clone(), + expected_index: self.next_record_index, + observed_index: block.first_index, + }); + } + let expected_end = block + .first_index + .checked_add(block.records.len()) + .ok_or_else(|| Error::ReadInput("CBQ block index overflow".to_string()))?; + if expected_end != block.end_index { + return Err(Error::CbqOrdering { + path: self.path.clone(), + expected_index: expected_end, + observed_index: block.end_index, + }); + } + if block.first_index < self.next_record_index + || self.pending.contains_key(&block.first_index) + { + return Err(Error::CbqOrdering { + path: self.path.clone(), + expected_index: self.next_record_index, + observed_index: block.first_index, + }); + } + + self.pending.insert(block.first_index, block); + let mut ready = Vec::new(); + while let Some(block) = self.pending.remove(&self.next_record_index) { + for indexed in block.records { + if indexed.index != self.next_record_index { + return Err(Error::CbqOrdering { + path: self.path.clone(), + expected_index: self.next_record_index, + observed_index: indexed.index, + }); + } + ready.push(indexed.record); + self.next_record_index += 1; + } + } + Ok(ready) + } + + fn finish(&self, expected_end: usize) -> Result<(), Error> { + if self.next_record_index != expected_end || !self.pending.is_empty() { + let observed = self + .pending + .first_key_value() + .map_or(self.next_record_index, |(index, _)| *index); + return Err(Error::CbqOrdering { + path: self.path.clone(), + expected_index: expected_end, + observed_index: observed, + }); + } + Ok(()) + } +} + +fn open_reader(path: &Path) -> Result { + std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| { + binseq::cbq::MmapReader::new(path) + })) + .map_err(|_| Error::CbqWorkerPanic { + path: path.to_path_buf(), + start: 0, + end: 0, + })? + .map_err(|source| Error::CbqOpen { + path: path.to_path_buf(), + source, + }) +} + +fn read_block_record_counts( + reader: &binseq::cbq::MmapReader, + path: &Path, +) -> Result, Error> { + let headers = std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| { + reader + .iter_block_headers() + .collect::>>() + })) + .map_err(|_| Error::CbqWorkerPanic { + path: path.to_path_buf(), + start: 0, + end: reader.num_records(), + })? + .map_err(|source| Error::CbqOpen { + path: path.to_path_buf(), + source, + })?; + + let counts: Vec = headers + .into_iter() + .map(|header| { + usize::try_from(header.num_records).map_err(|_| { + Error::ReadInput(format!( + "CBQ block record count does not fit usize in {}", + path.display() + )) + }) + }) + .collect::>()?; + let indexed_records = counts.iter().try_fold(0usize, |total, count| { + total + .checked_add(*count) + .ok_or_else(|| Error::ReadInput("CBQ record-count overflow".to_string())) + })?; + if indexed_records != reader.num_records() { + return Err(Error::ReadInput(format!( + "CBQ block headers contain {indexed_records} records but the index advertises {} in {}", + reader.num_records(), + path.display() + ))); + } + Ok(counts) +} + +/// Split the CBQ block index into bounded, block-aligned decode windows. +/// +/// `max_records_per_window` is a hard ceiling with one exception: the first block +/// of a window is always admitted, so an oversized block still makes progress +/// instead of stalling the plan. +fn plan_windows( + block_record_counts: &[usize], + total_records: usize, + max_records_per_window: usize, +) -> Result, Error> { + if total_records == 0 { + return Ok(Vec::new()); + } + let budget = max_records_per_window.max(1); + let mut windows = Vec::new(); + let mut block_index = 0usize; + let mut record_start = 0usize; + + while record_start < total_records { + if block_index >= block_record_counts.len() { + return Err(Error::ReadInput(format!( + "CBQ index contains fewer records than advertised ({record_start} < {total_records})" + ))); + } + let mut block_count = 0usize; + let mut boundary_end = record_start; + + while block_index < block_record_counts.len() { + let block_records = block_record_counts[block_index]; + if block_count > 0 + && boundary_end + .saturating_sub(record_start) + .saturating_add(block_records) + > budget + { + break; + } + boundary_end = boundary_end + .checked_add(block_records) + .ok_or_else(|| Error::ReadInput("CBQ record-count overflow".to_string()))?; + block_index += 1; + block_count += 1; + if boundary_end >= total_records { + break; + } + } + + let record_end = boundary_end.min(total_records); + if record_end <= record_start { + return Err(Error::ReadInput( + "CBQ index contains an empty or non-advancing block window".to_string(), + )); + } + windows.push(CbqWindow { + record_start, + record_end, + block_count, + }); + record_start = record_end; + if record_end < boundary_end { + break; + } + } + + Ok(windows) +} + +fn record_index(record: &R, path: &Path) -> Result { + usize::try_from(record.index()) + .map_err(|_| record_error(path, usize::MAX, None, "record index does not fit usize")) +} + +fn convert_mate( + record: &R, + path: &Path, + index: usize, + mate: u8, + extended: bool, +) -> Result { + let (header, sequence, quality) = if extended { + (record.xheader(), record.xseq(), record.xqual()) + } else { + (record.sheader(), record.sseq(), record.squal()) + }; + + let name = std::str::from_utf8(header) + .map_err(|source| { + record_error( + path, + index, + Some(mate), + format!("invalid UTF-8 in read name: {source}"), + ) + })? + .split(|c: char| c.is_ascii_whitespace()) + .next() + .unwrap_or_default() + .to_string(); + let encoded_sequence = sequence.iter().map(|&base| encode_base(base)).collect(); + let encoded_quality = if !record.has_quality() { + // No placeholder: qualities are never read by alignment, clipping or + // scoring, and `clip_read` preserves an empty quality through clipping. + Vec::new() + } else { + if quality.len() != sequence.len() { + return Err(record_error( + path, + index, + Some(mate), + format!( + "sequence/quality length mismatch: {} bases, {} qualities", + sequence.len(), + quality.len() + ), + )); + } + quality.to_vec() + }; + + Ok(EncodedRead { + name, + sequence: encoded_sequence, + quality: encoded_quality, + }) +} + +fn record_error( + path: &Path, + record_index: usize, + mate: Option, + message: impl Into, +) -> Error { + Error::CbqRecord { + path: path.to_path_buf(), + record_index, + mate_context: mate.map_or_else(String::new, |mate| format!(" mate {mate}")), + message: message.into(), + } +} + +fn into_binseq_error(error: Error) -> binseq::Error { + binseq::Error::GenericError(Box::new(error)) +} + +/// Resolve CBQ decoder concurrency from CLI values. +pub fn decoder_threads(params: &crate::params::Parameters) -> usize { + let requested = if params.read_files_n_threads > 0 { + params.read_files_n_threads + } else { + (usize::from(params.run_thread_n) / 4).clamp(1, 4) + }; + let available = std::thread::available_parallelism().map_or(1, usize::from); + requested.clamp(1, available) +} + +#[cfg(test)] +mod tests { + use super::*; + use binseq::SequencingRecordBuilder; + use binseq::write::{BinseqWriterBuilder, Format}; + use std::fs::{File, OpenOptions}; + use std::io::{Read, Seek, SeekFrom, Write}; + use tempfile::{NamedTempFile, TempPath}; + + fn write_cbq(paired: bool, qualities: bool, headers: bool, n_records: usize) -> TempPath { + let temp = NamedTempFile::new().unwrap(); + let path = temp.into_temp_path(); + let file = File::create(&path).unwrap(); + let mut writer = BinseqWriterBuilder::new(Format::Cbq) + .paired(paired) + .quality(qualities) + .headers(headers) + .block_size(256) + .build(file) + .unwrap(); + + for i in 0..n_records { + let h1 = format!("read{i}/1"); + let h2 = format!("read{i}/2"); + let mut builder = SequencingRecordBuilder::default().s_seq(b"ACGTNACGTN"); + if headers { + builder = builder.s_header(h1.as_bytes()); + } + if qualities { + builder = builder.s_qual(b"ABCDEFGHIJ"); + } + if paired { + builder = builder.x_seq(b"TGCANTGCAN"); + if headers { + builder = builder.x_header(h2.as_bytes()); + } + if qualities { + builder = builder.x_qual(b"JKLMNOPQRS"); + } + } + writer.push(builder.build().unwrap()).unwrap(); + } + writer.finish().unwrap(); + path + } + + fn collect_single(path: &Path, threads: usize, max_records: usize) -> Vec { + collect_single_windowed(path, threads, max_records, MAX_RECORDS_PER_WINDOW) + } + + fn collect_single_windowed( + path: &Path, + threads: usize, + max_records: usize, + window_records: usize, + ) -> Vec { + let (sender, receiver) = sync_channel(128); + CbqProducer { + path: path.to_path_buf(), + converter: SingleEndConverter, + decoder_threads: threads, + window_records, + } + .produce(7, max_records, &sender) + .unwrap(); + drop(sender); + receiver + .into_iter() + .map(Result::unwrap) + .take_while(|batch| !batch.is_empty()) + .flatten() + .collect() + } + + fn collect_paired(path: &Path, threads: usize) -> Vec { + collect_paired_windowed(path, threads, MAX_RECORDS_PER_WINDOW) + } + + fn collect_paired_windowed( + path: &Path, + threads: usize, + window_records: usize, + ) -> Vec { + let (sender, receiver) = sync_channel(128); + CbqProducer { + path: path.to_path_buf(), + converter: PairedEndConverter, + decoder_threads: threads, + window_records, + } + .produce(9, usize::MAX, &sender) + .unwrap(); + drop(sender); + receiver + .into_iter() + .map(Result::unwrap) + .take_while(|batch| !batch.is_empty()) + .flatten() + .collect() + } + + fn write_custom_headers(header1: &[u8], header2: Option<&[u8]>) -> TempPath { + let temp = NamedTempFile::new().unwrap(); + let path = temp.into_temp_path(); + let paired = header2.is_some(); + let mut writer = BinseqWriterBuilder::new(Format::Cbq) + .paired(paired) + .quality(true) + .headers(true) + .block_size(256) + .build(File::create(&path).unwrap()) + .unwrap(); + let mut record = SequencingRecordBuilder::default() + .s_seq(b"ACGT") + .s_qual(b"IIII") + .s_header(header1); + if let Some(header2) = header2 { + record = record.x_seq(b"TGCA").x_qual(b"IIII").x_header(header2); + } + writer.push(record.build().unwrap()).unwrap(); + writer.finish().unwrap(); + path + } + + /// Every window must tile `0..total` exactly once, in order, at block + /// boundaries, and honor the record budget (one oversized block excepted). + fn assert_windows_tile(windows: &[CbqWindow], counts: &[usize], total: usize, budget: usize) { + let mut expected_start = 0usize; + let mut consumed_blocks = 0usize; + for window in windows { + assert_eq!( + window.record_start, expected_start, + "windows must be gapless" + ); + assert!(window.record_end > window.record_start, "empty window"); + let records = window.record_end - window.record_start; + let largest = counts[consumed_blocks..consumed_blocks + window.block_count] + .iter() + .copied() + .max() + .unwrap_or(0); + assert!( + records <= budget.max(largest), + "window of {records} records exceeds max(budget {budget}, largest block {largest})" + ); + expected_start = window.record_end; + consumed_blocks += window.block_count; + } + assert_eq!(expected_start, total, "windows must cover every record"); + } + + #[test] + fn window_planner_fills_windows_up_to_the_record_budget() { + // 20 blocks of 10 records, budget 45 -> 4 whole blocks per window. + let windows = plan_windows(&[10; 20], 200, 45).unwrap(); + assert_windows_tile(&windows, &[10; 20], 200, 45); + assert!(windows.iter().all(|window| window.block_count == 4)); + assert_eq!(windows.len(), 5); + } + + #[test] + fn window_planner_is_independent_of_decoder_threads() { + // The plan is a pure function of the block index and the budget: no + // decoder-thread input, so peak memory cannot scale with --readFilesNthreads. + // Counts are the measured shape of a 1 MiB-block, quality-bearing, + // 150 bp single-end CBQ (~4651 records per block). + let counts = [4_651usize; 43]; + let total: usize = counts.iter().sum(); + let windows = plan_windows(&counts, total, MAX_RECORDS_PER_WINDOW).unwrap(); + assert_windows_tile(&windows, &counts, total, MAX_RECORDS_PER_WINDOW); + assert_eq!(windows.first().unwrap().block_count, 7); + assert!( + windows.len() <= 7, + "windows must stay large enough to amortize the extra trailing block \ + decompression (one per window): {windows:?}" + ); + } + + #[test] + fn window_planner_caps_records_regardless_of_block_size() { + // A file whose blocks dwarf the budget: each window is exactly one block, + // so peak memory is one block rather than a multiple of it. + let counts = [800_000usize; 4]; + let windows = plan_windows(&counts, 3_200_000, MAX_RECORDS_PER_WINDOW).unwrap(); + assert_windows_tile(&windows, &counts, 3_200_000, MAX_RECORDS_PER_WINDOW); + assert!(windows.iter().all(|window| window.block_count == 1)); + + // Mixed sizes: an oversized block never drags neighbours into its window. + let counts = [10usize, 500, 10, 10]; + let windows = plan_windows(&counts, 530, 100).unwrap(); + assert_windows_tile(&windows, &counts, 530, 100); + assert_eq!( + windows, + vec![ + CbqWindow { + record_start: 0, + record_end: 10, + block_count: 1 + }, + CbqWindow { + record_start: 10, + record_end: 510, + block_count: 1 + }, + CbqWindow { + record_start: 510, + record_end: 530, + block_count: 2 + }, + ] + ); + } + + #[test] + fn window_planner_stops_inside_final_block() { + let windows = plan_windows(&[10, 10, 10], 17, 25).unwrap(); + assert_eq!( + windows, + vec![CbqWindow { + record_start: 0, + record_end: 17, + block_count: 2, + }] + ); + } + + #[test] + fn window_planner_rejects_a_short_block_index() { + // The index advertises more records than its block headers account for. + let error = plan_windows(&[10, 10], 30, 25).unwrap_err(); + assert!(error.to_string().contains("fewer records than advertised")); + } + + #[test] + fn single_end_decoding_is_ordered_and_honors_limit() { + let path = write_cbq(false, true, true, 75); + let one_thread = collect_single(&path, 1, 53); + let four_threads = collect_single(&path, 4, 53); + + assert_eq!(one_thread.len(), 53); + assert_eq!(four_threads.len(), 53); + for (i, (one, four)) in one_thread.iter().zip(&four_threads).enumerate() { + assert_eq!(one.name, format!("read{i}/1")); + assert_eq!(one.name, four.name); + assert_eq!(one.sequence, vec![0, 1, 2, 3, 4, 0, 1, 2, 3, 4]); + assert_eq!(one.sequence, four.sequence); + assert_eq!(one.quality, b"ABCDEFGHIJ"); + } + } + + /// The record budget makes windows thread-independent, so multi-window + /// decoding must be exercised explicitly: the alignment-batch buffer is + /// carried across window boundaries while blocks still arrive out of order + /// inside each window. Output must be identical for every combination. + #[test] + fn decoding_is_identical_across_decoder_threads_and_window_sizes() { + let path = write_cbq(false, true, true, 300); + let counts: Vec = open_reader(&path) + .unwrap() + .iter_block_headers() + .map(|header| header.unwrap().num_records as usize) + .collect(); + assert!(counts.len() >= 8, "fixture needs many blocks: {counts:?}"); + + let reference = collect_single(&path, 1, usize::MAX); + assert_eq!(reference.len(), 300); + for window_records in [1usize, 17, 64, 250, MAX_RECORDS_PER_WINDOW] { + let windows = plan_windows(&counts, 300, window_records).unwrap(); + for threads in [1usize, 2, 3, 4, 8] { + let decoded = collect_single_windowed(&path, threads, usize::MAX, window_records); + assert_eq!( + decoded.len(), + reference.len(), + "threads={threads} window_records={window_records}" + ); + for (expected, actual) in reference.iter().zip(&decoded) { + assert_eq!( + (&expected.name, &expected.sequence, &expected.quality), + (&actual.name, &actual.sequence, &actual.quality), + "threads={threads} window_records={window_records}" + ); + } + // A truncating limit must also not leak later records. + let limited = collect_single_windowed(&path, threads, 137, window_records); + assert_eq!(limited.len(), 137); + assert_eq!(limited.last().unwrap().name, reference[136].name); + } + if window_records <= 64 { + assert!( + windows.len() >= 3, + "expected multi-window decoding for window_records={window_records}" + ); + } + } + + // Same for paired records, where each window carries two mates per index. + let paired_path = write_cbq(true, true, true, 300); + let paired_reference = collect_paired(&paired_path, 1); + assert_eq!(paired_reference.len(), 300); + for window_records in [1usize, 23, 100] { + for threads in [1usize, 3, 8] { + let decoded = collect_paired_windowed(&paired_path, threads, window_records); + let names: Vec<_> = decoded.iter().map(|read| read.name.clone()).collect(); + let expected: Vec<_> = paired_reference + .iter() + .map(|read| read.name.clone()) + .collect(); + assert_eq!( + names, expected, + "threads={threads} window_records={window_records}" + ); + } + } + } + + #[test] + fn paired_qualityless_decoding_uses_generated_names_without_exporting_placeholders() { + let path = write_cbq(true, false, false, 41); + let reads = collect_paired(&path, 4); + assert_eq!(reads.len(), 41); + for (i, read) in reads.iter().enumerate() { + assert_eq!(read.name, i.to_string()); + assert_eq!(read.mate1.name, i.to_string()); + assert_eq!(read.mate2.name, i.to_string()); + // No placeholder buffer is allocated for a quality-less record. + assert!(read.mate1.quality.is_empty()); + assert!(read.mate2.quality.is_empty()); + } + } + + #[test] + fn paired_headers_are_normalized() { + let path = write_cbq(true, true, true, 17); + let reads = collect_paired(&path, 3); + assert_eq!(reads.len(), 17); + for (i, read) in reads.iter().enumerate() { + assert_eq!(read.name, format!("read{i}")); + assert_eq!(read.mate1.name, format!("read{i}/1")); + assert_eq!(read.mate2.name, format!("read{i}/2")); + } + } + + #[test] + fn invalid_utf8_and_mismatched_pair_names_are_contextual_errors() { + let invalid = write_custom_headers(&[0xff], None); + let (sender, _receiver) = sync_channel(8); + let error = single_end_producer(invalid.to_path_buf(), 1) + .produce(10, usize::MAX, &sender) + .unwrap_err(); + let message = error.to_string(); + assert!(message.contains("invalid UTF-8")); + assert!(message.contains("record 0 mate 1")); + + let mismatched = write_custom_headers(b"left/1", Some(b"right/2")); + let (sender, _receiver) = sync_channel(8); + let error = paired_end_producer(mismatched.to_path_buf(), 2) + .produce(10, usize::MAX, &sender) + .unwrap_err(); + assert!(error.to_string().contains("paired read names do not match")); + } + + #[test] + fn dropped_alignment_receiver_cancels_without_deadlock() { + let path = write_cbq(false, true, true, 75); + let (sender, receiver) = sync_channel(0); + drop(receiver); + single_end_producer(path.to_path_buf(), 4) + .produce(1, usize::MAX, &sender) + .unwrap(); + } + + #[test] + fn ordered_coordinator_reassembles_blocks_and_rejects_duplicates() { + fn block(first: usize, values: &[usize]) -> DecodedBlock { + DecodedBlock { + first_index: first, + end_index: first + values.len(), + records: values + .iter() + .enumerate() + .map(|(offset, value)| IndexedRecord { + index: first + offset, + record: *value, + }) + .collect(), + } + } + + let mut coordinator = OrderedBlocks::new(PathBuf::from("test.cbq"), 0); + assert!(coordinator.push(block(2, &[2, 3])).unwrap().is_empty()); + assert_eq!( + coordinator.push(block(0, &[0, 1])).unwrap(), + vec![0, 1, 2, 3] + ); + assert!(coordinator.push(block(4, &[4])).is_ok()); + assert!(coordinator.push(block(4, &[4])).is_err()); + } + + #[test] + fn corrupt_compressed_block_returns_error_without_hanging() { + let path = write_cbq(false, true, true, 30); + let data_offset = (std::mem::size_of::() + + std::mem::size_of::()) as u64; + let mut file = OpenOptions::new() + .read(true) + .write(true) + .open(&path) + .unwrap(); + file.seek(SeekFrom::Start(data_offset)).unwrap(); + let mut byte = [0u8; 1]; + file.read_exact(&mut byte).unwrap(); + byte[0] ^= 0xff; + file.seek(SeekFrom::Start(data_offset)).unwrap(); + file.write_all(&byte).unwrap(); + file.flush().unwrap(); + + let (sender, _receiver) = sync_channel(64); + let result = + single_end_producer(path.to_path_buf(), 3).produce(10, usize::MAX, &sender); + assert!(result.is_err()); + } +} diff --git a/src/io/fastq.rs b/src/io/fastq.rs index 71b1f81..7a8641d 100644 --- a/src/io/fastq.rs +++ b/src/io/fastq.rs @@ -1,5 +1,6 @@ /// FASTQ reader with base encoding and decompression support use crate::error::Error; +use crate::io::reads::{EncodedRead, PairedRead, decode_base, encode_base, strip_mate_suffix}; use flate2::read::GzDecoder; use noodles::fastq; use std::fs::File; @@ -46,29 +47,6 @@ impl UnmappedFastqWriter { } } -/// A read from a FASTQ file with encoded bases -#[derive(Debug, Clone)] -pub struct EncodedRead { - /// Read identifier - pub name: String, - /// Base sequence encoded as 0=A, 1=C, 2=G, 3=T, 4=N - pub sequence: Vec, - /// FASTQ ASCII quality bytes (Phred+33 encoded) - subtract 33 before - /// writing to BAM binary QUAL. - pub quality: Vec, -} - -/// A paired-end read from two FASTQ files -#[derive(Debug, Clone)] -pub struct PairedRead { - /// Base read name (without /1 or /2 suffix) - pub name: String, - /// First mate in pair - pub mate1: EncodedRead, - /// Second mate in pair - pub mate2: EncodedRead, -} - /// FASTQ reader that handles decompression and base encoding pub struct FastqReader { inner: fastq::io::Reader>, @@ -268,113 +246,10 @@ impl PairedFastqReader { } } -/// Strip mate suffix from read name for pairing -/// -/// Removes common paired-end suffixes: -/// - /1 or /2 (Illumina convention) -/// - .R1 or .R2 (alternative convention) -/// - _1 or _2 (another convention) -/// - space and everything after (e.g., "READ_NAME 1:N:0:0" -> "READ_NAME") -/// -/// # Arguments -/// * `name` - Original read name from FASTQ -/// -/// # Returns -/// Base name with mate suffix removed -#[allow(clippy::case_sensitive_file_extension_comparisons)] // false positive -pub fn strip_mate_suffix(name: &str) -> String { - // First, strip space and everything after (Illumina format) - let name = if let Some(pos) = name.find(' ') { - &name[..pos] - } else { - name - }; - - // Strip common mate suffixes - if name.ends_with("/1") || name.ends_with("/2") { - name[..name.len() - 2].to_string() - } else if name.ends_with(".R1") || name.ends_with(".R2") { - name[..name.len() - 3].to_string() - } else if name.ends_with("_1") || name.ends_with("_2") { - name[..name.len() - 2].to_string() - } else { - name.to_string() - } -} - -/// Convert FASTQ base character to genome encoding -/// -/// # Arguments -/// * `base` - ASCII base character (A, C, G, T, N, or lowercase variants) -/// -/// # Returns -/// Encoded base: 0=A, 1=C, 2=G, 3=T, 4=N (or any ambiguous base) -pub fn encode_base(base: u8) -> u8 { - match base.to_ascii_uppercase() { - b'A' => 0, - b'C' => 1, - b'G' => 2, - b'T' => 3, - _ => 4, // N or any ambiguous base (R, Y, S, W, K, M, etc.) - } -} - -/// Decode genome encoding to ASCII base character -/// -/// # Arguments -/// * `encoded` - Encoded base (0-4) -/// -/// # Returns -/// ASCII base character (A, C, G, T, or N) -pub fn decode_base(encoded: u8) -> u8 { - match encoded { - 0 => b'A', - 1 => b'C', - 2 => b'G', - 3 => b'T', - _ => b'N', - } -} - -/// Complement an encoded base (A=0↔T=3, C=1↔G=2, N=4→N=4). -pub fn complement_base(encoded: u8) -> u8 { - match encoded { - 0 => 3, // A -> T - 1 => 2, // C -> G - 2 => 1, // G -> C - 3 => 0, // T -> A - _ => encoded, // N -> N - } -} - -/// Apply read clipping from 5' and 3' ends -/// -/// # Arguments -/// * `seq` - Original sequence -/// * `qual` - Original quality scores -/// * `clip5p` - Number of bases to clip from 5' end -/// * `clip3p` - Number of bases to clip from 3' end -/// -/// # Returns -/// Tuple of (clipped_sequence, clipped_quality) -pub fn clip_read(seq: &[u8], qual: &[u8], clip5p: usize, clip3p: usize) -> (Vec, Vec) { - let len = seq.len(); - - // Handle edge cases - if clip5p + clip3p >= len { - // Clipping removes entire read - return (Vec::new(), Vec::new()); - } - - let start = clip5p; - let end = len - clip3p; - - (seq[start..end].to_vec(), qual[start..end].to_vec()) -} - #[cfg(test)] mod tests { use super::*; + use crate::io::reads::clip_read; use std::io::Write; use tempfile::NamedTempFile; diff --git a/src/io/input.rs b/src/io/input.rs new file mode 100644 index 0000000..897e580 --- /dev/null +++ b/src/io/input.rs @@ -0,0 +1,468 @@ +//! Input format detection and layout planning. + +use crate::error::Error; +use crate::io::cbq::{ + CbqProducer, PairedEndConverter, SingleEndConverter, decoder_threads, paired_end_producer, + single_end_producer, +}; +use crate::io::fastq::{FastqReader, PairedFastqReader}; +use crate::io::reads::{EncodedRead, PairedRead}; +use crate::params::{OutReadsUnmapped, Parameters}; +use binseq::write::Format; +use std::fs::File; +use std::io::Read; +use std::path::{Path, PathBuf}; +use std::sync::mpsc::SyncSender; + +const MAGIC_PEEK_LEN: usize = 7; + +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum ReadLayout { + SingleEnd, + PairedEnd, +} + +/// A validated ordinary-alignment input configuration. +/// +/// Enum variants encode the only supported path/layout combinations, so callers +/// cannot accidentally construct (for example) paired CBQ from two paths. +#[derive(Debug, Clone, PartialEq, Eq)] +pub enum ReadInputPlan { + FastqSingle { + path: PathBuf, + }, + FastqPaired { + mate1: PathBuf, + mate2: PathBuf, + }, + Cbq { + path: PathBuf, + layout: ReadLayout, + }, +} + +#[derive(Debug)] +enum DetectedInput { + Fastq, + Cbq { + layout: ReadLayout, + has_qualities: bool, + }, + Unsupported(Format), +} + +pub enum SingleEndProducer { + Fastq(FastqReader), + Cbq(CbqProducer), +} + +pub enum PairedEndProducer { + Fastq(PairedFastqReader), + Cbq(CbqProducer), +} + +impl ReadInputPlan { + pub fn resolve(params: &Parameters) -> Result { + if params.read_files_in.is_empty() { + return Err(Error::ReadInput( + "no read files specified (--readFilesIn)".to_string(), + )); + } + + let detected = params + .read_files_in + .iter() + .map(|path| detect_input(path)) + .collect::, _>>()?; + + for (path, format) in params.read_files_in.iter().zip(&detected) { + if let DetectedInput::Unsupported(format) = format { + return Err(Error::ReadInput(format!( + "unsupported BINSEQ format {} for {}; only CBQ is supported", + format_name(*format), + path.display() + ))); + } + } + + if (params.solo_enabled() || params.read_files_manifest.is_some()) + && detected + .iter() + .any(|input| !matches!(input, DetectedInput::Fastq)) + { + return Err(Error::ReadInput( + "CBQ is not supported with STARsolo or --readFilesManifest".to_string(), + )); + } + + match (params.read_files_in.as_slice(), detected.as_slice()) { + ([path], [DetectedInput::Fastq]) => Ok(Self::FastqSingle { path: path.clone() }), + ( + [path], + [DetectedInput::Cbq { + layout, + has_qualities, + .. + }], + ) => { + if params.read_files_command.is_some() { + return Err(Error::ReadInput( + "CBQ is incompatible with --readFilesCommand".to_string(), + )); + } + if !has_qualities && params.out_reads_unmapped == OutReadsUnmapped::Fastx { + return Err(Error::ReadInput( + "quality-less CBQ cannot be used with --outReadsUnmapped Fastx" + .to_string(), + )); + } + Ok(Self::Cbq { + path: path.clone(), + layout: *layout, + }) + } + ([mate1, mate2], [DetectedInput::Fastq, DetectedInput::Fastq]) => { + Ok(Self::FastqPaired { + mate1: mate1.clone(), + mate2: mate2.clone(), + }) + } + ([_, _], _) => Err(Error::ReadInput( + "paired CBQ is stored in one file; two CBQ paths and mixed FASTQ/CBQ inputs are unsupported" + .to_string(), + )), + (paths, _) => Err(Error::ReadInput(format!( + "invalid number of read files: {} (expected 1 or 2)", + paths.len() + ))), + } + } + + pub fn label(&self) -> &'static str { + match self { + Self::FastqSingle { .. } | Self::FastqPaired { .. } => "Fastq", + Self::Cbq { .. } => "Cbq", + } + } + + pub fn layout(&self) -> ReadLayout { + match self { + Self::FastqSingle { .. } + | Self::Cbq { + layout: ReadLayout::SingleEnd, + .. + } => ReadLayout::SingleEnd, + Self::FastqPaired { .. } + | Self::Cbq { + layout: ReadLayout::PairedEnd, + .. + } => ReadLayout::PairedEnd, + } + } + + pub fn open_single_end_producer( + &self, + params: &Parameters, + ) -> Result { + match self { + Self::FastqSingle { path } => Ok(SingleEndProducer::Fastq(FastqReader::open( + path, + params.read_files_command.as_deref(), + )?)), + Self::Cbq { + path, + layout: ReadLayout::SingleEnd, + .. + } => Ok(SingleEndProducer::Cbq(single_end_producer( + path.clone(), + decoder_threads(params), + ))), + Self::FastqPaired { .. } + | Self::Cbq { + layout: ReadLayout::PairedEnd, + .. + } => Err(Error::ReadInput( + "paired input cannot be opened as single-end".to_string(), + )), + } + } + + pub fn open_paired_end_producer( + &self, + params: &Parameters, + ) -> Result { + match self { + Self::FastqPaired { mate1, mate2 } => Ok(PairedEndProducer::Fastq( + PairedFastqReader::open(mate1, mate2, params.read_files_command.as_deref())?, + )), + Self::Cbq { + path, + layout: ReadLayout::PairedEnd, + .. + } => Ok(PairedEndProducer::Cbq(paired_end_producer( + path.clone(), + decoder_threads(params), + ))), + Self::FastqSingle { .. } + | Self::Cbq { + layout: ReadLayout::SingleEnd, + .. + } => Err(Error::ReadInput( + "single-end input cannot be opened as paired-end".to_string(), + )), + } + } +} + +impl SingleEndProducer { + pub fn produce( + self, + batch_size: usize, + max_records: usize, + sender: &SyncSender, Error>>, + ) -> Result<(), Error> { + match self { + Self::Cbq(producer) => producer.produce(batch_size, max_records, sender), + Self::Fastq(mut reader) => produce_fastq( + |n| reader.read_batch(n), + batch_size, + max_records, + sender, + ), + } + } +} + +impl PairedEndProducer { + pub fn produce( + self, + batch_size: usize, + max_records: usize, + sender: &SyncSender, Error>>, + ) -> Result<(), Error> { + match self { + Self::Cbq(producer) => producer.produce(batch_size, max_records, sender), + Self::Fastq(mut reader) => produce_fastq( + |n| reader.read_paired_batch(n), + batch_size, + max_records, + sender, + ), + } + } +} + +/// Shared FASTQ producer loop: bounded batches, stopping at `max_records`. +fn produce_fastq( + mut read_batch: impl FnMut(usize) -> Result, Error>, + batch_size: usize, + max_records: usize, + sender: &SyncSender, Error>>, +) -> Result<(), Error> { + let mut produced = 0usize; + loop { + let remaining = max_records.saturating_sub(produced); + if remaining == 0 { + let _ = sender.send(Ok(Vec::new())); + return Ok(()); + } + let batch = read_batch(batch_size.min(remaining))?; + let finished = batch.is_empty(); + produced += batch.len(); + if sender.send(Ok(batch)).is_err() || finished { + return Ok(()); + } + } +} + +fn detect_input(path: &Path) -> Result { + let metadata = std::fs::metadata(path).map_err(|source| Error::io(source, path))?; + if !metadata.is_file() { + // Preserve the existing FASTQ behavior for named pipes and other streams. + // CBQ itself is only recognized from regular files because it requires mmap. + return Ok(DetectedInput::Fastq); + } + let mut file = File::open(path).map_err(|source| Error::io(source, path))?; + let mut magic = [0u8; MAGIC_PEEK_LEN]; + let n = file + .read(&mut magic) + .map_err(|source| Error::io(source, path))?; + + match Format::sniff(&magic[..n]) { + None => Ok(DetectedInput::Fastq), + Some(Format::Bq) => Ok(DetectedInput::Unsupported(Format::Bq)), + Some(Format::Vbq) => Ok(DetectedInput::Unsupported(Format::Vbq)), + Some(Format::Cbq) => { + let reader = std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| { + binseq::cbq::MmapReader::new(path) + })) + .map_err(|_| Error::CbqWorkerPanic { + path: path.to_path_buf(), + start: 0, + end: 0, + })? + .map_err(|source| Error::CbqOpen { + path: path.to_path_buf(), + source, + })?; + let header = reader.header(); + Ok(DetectedInput::Cbq { + layout: if reader.is_paired() { + ReadLayout::PairedEnd + } else { + ReadLayout::SingleEnd + }, + has_qualities: header.has_qualities(), + }) + } + } +} + +fn format_name(format: Format) -> &'static str { + match format { + Format::Bq => "BQ", + Format::Vbq => "VBQ", + Format::Cbq => "CBQ", + } +} + +#[cfg(test)] +mod tests { + use super::*; + use binseq::SequencingRecordBuilder; + use binseq::write::BinseqWriterBuilder; + use std::ffi::OsString; + use std::io::Write; + use tempfile::{NamedTempFile, TempPath}; + + fn params_for(paths: &[&Path]) -> Parameters { + let mut args = vec![ + OsString::from("rustar-aligner"), + OsString::from("--readFilesIn"), + ]; + args.extend(paths.iter().map(|path| path.as_os_str().to_os_string())); + Parameters::parse_from(args) + } + + fn write_cbq(paired: bool, qualities: bool, suffix: &str) -> TempPath { + let temp = tempfile::Builder::new().suffix(suffix).tempfile().unwrap(); + let path = temp.into_temp_path(); + let mut writer = BinseqWriterBuilder::new(Format::Cbq) + .paired(paired) + .quality(qualities) + .headers(true) + .block_size(256) + .build(File::create(&path).unwrap()) + .unwrap(); + let mut record = SequencingRecordBuilder::default() + .s_seq(b"ACGTN") + .s_header(b"read/1"); + if qualities { + record = record.s_qual(b"IIIII"); + } + if paired { + record = record.x_seq(b"TGCAN").x_header(b"read/2"); + if qualities { + record = record.x_qual(b"IIIII"); + } + } + writer.push(record.build().unwrap()).unwrap(); + writer.finish().unwrap(); + path + } + + #[test] + fn detects_single_and_paired_cbq_by_magic_not_extension() { + let single = write_cbq(false, true, ".data"); + let paired = write_cbq(true, true, ""); + + assert!(matches!( + ReadInputPlan::resolve(¶ms_for(&[&single])).unwrap(), + ReadInputPlan::Cbq { + layout: ReadLayout::SingleEnd, + .. + } + )); + assert!(matches!( + ReadInputPlan::resolve(¶ms_for(&[&paired])).unwrap(), + ReadInputPlan::Cbq { + layout: ReadLayout::PairedEnd, + .. + } + )); + } + + #[test] + fn fastq_named_cbq_remains_fastq() { + let mut fastq = tempfile::Builder::new().suffix(".cbq").tempfile().unwrap(); + fastq.write_all(b"@r\nACGT\n+\nIIII\n").unwrap(); + let plan = ReadInputPlan::resolve(¶ms_for(&[fastq.path()])).unwrap(); + assert!(matches!(plan, ReadInputPlan::FastqSingle { .. })); + } + + #[test] + fn rejects_bq_and_vbq_magic_explicitly() { + for (magic, expected) in [ + (binseq::bq::FILE_MAGIC.as_slice(), "BQ"), + (binseq::vbq::FILE_MAGIC.as_slice(), "VBQ"), + ] { + let mut input = NamedTempFile::new().unwrap(); + input.write_all(magic).unwrap(); + let error = ReadInputPlan::resolve(¶ms_for(&[input.path()])).unwrap_err(); + assert!(error.to_string().contains(expected)); + } + } + + #[test] + fn recognized_corrupt_cbq_is_not_treated_as_fastq() { + let mut input = NamedTempFile::new().unwrap(); + input.write_all(binseq::cbq::FILE_MAGIC).unwrap(); + let error = ReadInputPlan::resolve(¶ms_for(&[input.path()])).unwrap_err(); + assert!(error.to_string().contains("CBQ")); + } + + #[test] + fn rejects_mixed_and_two_path_cbq_inputs() { + let cbq1 = write_cbq(false, true, ".cbq"); + let cbq2 = write_cbq(false, true, ".cbq"); + let mut fastq = NamedTempFile::new().unwrap(); + fastq.write_all(b"@r\nACGT\n+\nIIII\n").unwrap(); + + let mixed = ReadInputPlan::resolve(¶ms_for(&[&cbq1, fastq.path()])).unwrap_err(); + assert!(mixed.to_string().contains("mixed FASTQ/CBQ")); + let two = ReadInputPlan::resolve(¶ms_for(&[&cbq1, &cbq2])).unwrap_err(); + assert!(two.to_string().contains("one file")); + } + + #[test] + fn rejects_incompatible_cbq_options_upfront() { + let qualityless = write_cbq(false, false, ".cbq"); + let mut params = params_for(&[&qualityless]); + params.out_reads_unmapped = OutReadsUnmapped::Fastx; + assert!( + ReadInputPlan::resolve(¶ms) + .unwrap_err() + .to_string() + .contains("quality-less") + ); + + let quality = write_cbq(false, true, ".cbq"); + let mut params = params_for(&[&quality]); + params.read_files_command = Some("cat".to_string()); + assert!( + ReadInputPlan::resolve(¶ms) + .unwrap_err() + .to_string() + .contains("readFilesCommand") + ); + + let mut params = params_for(&[&quality]); + params.solo_type = crate::params::SoloType::CbUmiSimple; + assert!( + ReadInputPlan::resolve(¶ms) + .unwrap_err() + .to_string() + .contains("STARsolo") + ); + } +} diff --git a/src/io/mod.rs b/src/io/mod.rs index a458f83..6106435 100644 --- a/src/io/mod.rs +++ b/src/io/mod.rs @@ -1,6 +1,9 @@ // Phase 6+: FASTQ reader, SAM/BAM output, SJ.out.tab pub mod bam; +pub mod cbq; pub mod fastq; +pub mod input; pub mod log; +pub mod reads; pub mod sam; diff --git a/src/io/reads.rs b/src/io/reads.rs new file mode 100644 index 0000000..8e81e1b --- /dev/null +++ b/src/io/reads.rs @@ -0,0 +1,99 @@ +//! Format-independent owned read types and transformations. + +/// An owned read with bases in ruSTAR's internal encoding. +#[derive(Debug, Clone)] +pub struct EncodedRead { + /// Read identifier. + pub name: String, + /// Base sequence encoded as 0=A, 1=C, 2=G, 3=T, 4=N. + pub sequence: Vec, + /// FASTQ ASCII quality bytes (Phred+33 encoded). Empty when the input format + /// omitted qualities — nothing in the alignment, clipping or scoring path + /// reads qualities, so no placeholder is materialized. SAM/BAM builders treat + /// the empty slice as the format-level missing value (`*` / 0xff bytes). + pub quality: Vec, +} + +/// An owned paired-end read. +#[derive(Debug, Clone)] +pub struct PairedRead { + /// Base read name (without a mate suffix). + pub name: String, + /// First mate in the pair. + pub mate1: EncodedRead, + /// Second mate in the pair. + pub mate2: EncodedRead, +} + +/// Strip a common mate suffix from a read name. +#[allow(clippy::case_sensitive_file_extension_comparisons)] // false positive +pub fn strip_mate_suffix(name: &str) -> String { + let name = if let Some(pos) = name.find(' ') { + &name[..pos] + } else { + name + }; + + if name.ends_with("/1") || name.ends_with("/2") { + name[..name.len() - 2].to_string() + } else if name.ends_with(".R1") || name.ends_with(".R2") { + name[..name.len() - 3].to_string() + } else if name.ends_with("_1") || name.ends_with("_2") { + name[..name.len() - 2].to_string() + } else { + name.to_string() + } +} + +/// Convert an ASCII nucleotide to ruSTAR's genome encoding. +pub fn encode_base(base: u8) -> u8 { + match base.to_ascii_uppercase() { + b'A' => 0, + b'C' => 1, + b'G' => 2, + b'T' => 3, + _ => 4, + } +} + +/// Decode a ruSTAR-encoded nucleotide to ASCII. +pub fn decode_base(encoded: u8) -> u8 { + match encoded { + 0 => b'A', + 1 => b'C', + 2 => b'G', + 3 => b'T', + _ => b'N', + } +} + +/// Complement an encoded base (A=0↔T=3, C=1↔G=2, N=4→N=4). +pub fn complement_base(encoded: u8) -> u8 { + match encoded { + 0 => 3, + 1 => 2, + 2 => 1, + 3 => 0, + _ => encoded, + } +} + +/// Apply fixed 5' and 3' clipping to a sequence and its processing qualities. +/// +/// An empty `qual` (an input format that omits qualities) clips to an empty +/// quality, preserving the SAM/BAM missing-quality sentinel through clipping. +pub fn clip_read(seq: &[u8], qual: &[u8], clip5p: usize, clip3p: usize) -> (Vec, Vec) { + let len = seq.len(); + if clip5p + clip3p >= len { + return (Vec::new(), Vec::new()); + } + + let start = clip5p; + let end = len - clip3p; + let clipped_qual = if qual.is_empty() { + Vec::new() + } else { + qual[start..end].to_vec() + }; + (seq[start..end].to_vec(), clipped_qual) +} diff --git a/src/io/sam.rs b/src/io/sam.rs index 343e9e0..a1ecf7d 100644 --- a/src/io/sam.rs +++ b/src/io/sam.rs @@ -3,7 +3,7 @@ use crate::align::read_align::PairedAlignment; use crate::align::transcript::{Transcript, cigar_to_string}; use crate::error::Error; use crate::genome::Genome; -use crate::io::fastq::{complement_base, decode_base}; +use crate::io::reads::{complement_base, decode_base}; use crate::junction::encode_motif; use crate::mapq::calculate_mapq; use crate::params::{Parameters, SamAttributes}; @@ -263,7 +263,7 @@ impl SamWriter { // (CIGAR core, MD, NM) is built from that aligned slice, then apply_read_clips // restores the full read + soft-clips. clip5p/clip3p == 0 => aligned == full. let aligned_seq = &read_seq[clip5p..read_seq.len() - clip3p]; - let aligned_qual = &read_qual[clip5p..read_qual.len() - clip3p]; + let aligned_qual = clipped_quality(read_qual, clip5p, clip3p); let n_alignments = transcripts.len(); let max_output = if params.out_sam_mult_nmax < 0 { @@ -505,7 +505,7 @@ impl SamWriter { }; // Aligned slice of the mapped mate (core SEQ/CIGAR/MD built from this). let aligned_mapped_seq = &mapped_seq[mapped_clip5p..mapped_seq.len() - mapped_clip3p]; - let aligned_mapped_qual = &mapped_qual[mapped_clip5p..mapped_qual.len() - mapped_clip3p]; + let aligned_mapped_qual = clipped_quality(mapped_qual, mapped_clip5p, mapped_clip3p); // --- Build mapped mate record --- let mut mapped_rec = RecordBuf::default(); @@ -1093,6 +1093,16 @@ fn fastq_qual_to_phred(qual: &[u8]) -> Vec { qual.iter().map(|&b| b.saturating_sub(33)).collect() } +/// Return the clipped stored qualities, preserving an empty slice as the +/// SAM/BAM missing-quality sentinel for quality-less CBQ input. +fn clipped_quality(qual: &[u8], clip5p: usize, clip3p: usize) -> &[u8] { + if qual.is_empty() { + qual + } else { + &qual[clip5p..qual.len() - clip3p] + } +} + /// Re-attach clipped bases to a mapped record as soft-clips, matching STAR (and /// STAR-rs): the record is first built from the *clipped* read (so CIGAR core, MD /// and NM cover only the aligned bases), then this restores the full original read @@ -1576,7 +1586,7 @@ fn build_paired_mate_record( // Core SEQ/QUAL/MD are built from the aligned (clipped) slice; apply_read_clips // below restores the full mate read + soft-clips (clip==0 => aligned == mate_seq). let aligned_seq = &mate_seq[clip5p..mate_seq.len() - clip3p]; - let aligned_qual = &mate_qual[clip5p..mate_qual.len() - clip3p]; + let aligned_qual = clipped_quality(mate_qual, clip5p, clip3p); // Sequence and quality scores (reverse complement for reverse strand) if transcript.is_reverse { @@ -1911,6 +1921,69 @@ mod tests { assert_eq!(stored, &[0u8, 0, 1]); } + #[test] + fn missing_quality_survives_clipping_and_is_emitted_as_absent() { + use cigar::op::{Kind, Op}; + + let genome = make_test_genome(); + let params = Parameters::parse_from(["rustar-aligner", "--readFilesIn", "test.cbq"]); + let read_seq = vec![0, 1, 2, 3]; + let transcript = Transcript { + chr_idx: 0, + genome_start: 1, + genome_end: 3, + is_reverse: false, + exons: vec![], + cigar: vec![Op::new(Kind::Match, 2)], + score: 2, + n_mismatch: 0, + n_gap: 0, + n_junction: 0, + junction_motifs: vec![], + junction_annotated: vec![], + }; + + let records = SamWriter::build_alignment_records( + "read1", + &read_seq, + &[], + 1, + 1, + &[transcript], + &genome, + ¶ms, + 1, + ) + .unwrap(); + assert_eq!(records.len(), 1); + assert_eq!(records[0].sequence().len(), read_seq.len()); + assert!(records[0].quality_scores().is_empty()); + + let unmapped = SamWriter::build_unmapped_record( + "read2", + &read_seq, + &[], + ¶ms, + UnmappedReason::Other, + ) + .unwrap(); + assert!(unmapped.quality_scores().is_empty()); + + let output = NamedTempFile::new().unwrap(); + { + let mut writer = SamWriter::create(output.path(), &genome, ¶ms).unwrap(); + writer.write_batch(&[unmapped]).unwrap(); + } + let text = std::fs::read_to_string(output.path()).unwrap(); + let fields: Vec<_> = text + .lines() + .find(|line| !line.starts_with('@')) + .unwrap() + .split('\t') + .collect(); + assert_eq!(fields[10], "*"); + } + #[test] fn test_unmapped_record_tags_emitted() { let params = Parameters::parse_from(["rustar-aligner", "--readFilesIn", "t.fq"]); diff --git a/src/lib.rs b/src/lib.rs index e5f4f4e..3a6c9fb 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -236,6 +236,7 @@ impl AlignmentWriter for crate::io::bam::SortedBamStdoutWriter { fn align_reads(params: &Parameters) -> anyhow::Result<()> { use crate::index::GenomeIndex; + use crate::io::input::ReadInputPlan; use crate::params::TwopassMode; @@ -259,6 +260,14 @@ fn align_reads(params: &Parameters) -> anyhow::Result<()> { anyhow::bail!("No read files specified (--readFilesIn)"); } + // Resolve content-based input format and logical layout once. SmartSeq reads + // come from its manifest and do not participate in ordinary input planning. + let input_plan = if params.solo_type == params::SoloType::SmartSeq { + None + } else { + Some(ReadInputPlan::resolve(params)?) + }; + // 1. Load genome index info!("Loading genome index from {}", params.genome_dir.display()); let index = Arc::new(GenomeIndex::load(¶ms.genome_dir, params)?); @@ -355,6 +364,7 @@ fn align_reads(params: &Parameters) -> anyhow::Result<()> { run_single_pass( &index, ¶ms, + input_plan.as_ref(), quant_ctx.as_ref(), tr_idx.as_ref(), solo_ctx.as_ref(), @@ -365,6 +375,7 @@ fn align_reads(params: &Parameters) -> anyhow::Result<()> { run_two_pass( &index, ¶ms, + input_plan.as_ref(), quant_ctx.as_ref(), tr_idx.as_ref(), solo_ctx.as_ref(), @@ -617,6 +628,7 @@ fn run_smartseq( fn run_single_pass( index: &std::sync::Arc, params: &Parameters, + input_plan: Option<&crate::io::input::ReadInputPlan>, quant_ctx: Option<&std::sync::Arc>, tr_idx: Option<&std::sync::Arc>, solo_ctx: Option<&std::sync::Arc>, @@ -650,7 +662,9 @@ fn run_single_pass( use crate::io::fastq::UnmappedFastqWriter; use crate::params::OutReadsUnmapped; - let is_paired = params.read_files_in.len() == 2 && !params.solo_enabled(); + let is_paired = solo_ctx.is_none() + && input_plan + .is_some_and(|input| input.layout() == crate::io::input::ReadLayout::PairedEnd); let mut unmapped_w1: Option = if params.out_reads_unmapped == OutReadsUnmapped::Fastx { let path = params.output_path("Unmapped.out.mate1"); @@ -762,10 +776,11 @@ fn run_single_pass( return Ok(stats); } - let n_align_files = params.read_files_in.len(); - match n_align_files { - 1 => align_reads_single_end( + let input = input_plan.ok_or_else(|| anyhow::anyhow!("ordinary alignment input is missing"))?; + match input.layout() { + crate::io::input::ReadLayout::SingleEnd => align_reads_single_end( params, + input, index, writer.as_mut(), &stats, @@ -775,8 +790,9 @@ fn run_single_pass( tr_writer.as_mut(), unmapped_w1.as_mut(), ), - 2 => align_reads_paired_end( + crate::io::input::ReadLayout::PairedEnd => align_reads_paired_end( params, + input, index, writer.as_mut(), &stats, @@ -787,7 +803,6 @@ fn run_single_pass( unmapped_w1.as_mut(), unmapped_w2.as_mut(), ), - n => anyhow::bail!("Invalid number of read files: {n} (expected 1 or 2)"), }?; writer.finish()?; @@ -817,6 +832,7 @@ fn run_single_pass( fn run_two_pass( index: &std::sync::Arc, params: &Parameters, + input_plan: Option<&crate::io::input::ReadInputPlan>, quant_ctx: Option<&std::sync::Arc>, tr_idx: Option<&std::sync::Arc>, solo_ctx: Option<&std::sync::Arc>, @@ -825,7 +841,7 @@ fn run_two_pass( // PASS 1: Junction discovery (no quant counting in pass 1) info!("Two-pass mode: Pass 1 - Junction discovery"); - let (sj_stats_pass1, novel_junctions) = run_pass1(index, params)?; + let (sj_stats_pass1, novel_junctions) = run_pass1(index, params, input_plan)?; let pass1_dir = params.output_path("_STARpass1"); std::fs::create_dir_all(&pass1_dir)?; @@ -850,7 +866,14 @@ fn run_two_pass( // PASS 2: Re-alignment with merged DB (quant counts happen here) info!("Two-pass mode: Pass 2 - Re-alignment"); - let stats = run_single_pass(&Arc::new(merged_index), params, quant_ctx, tr_idx, solo_ctx)?; + let stats = run_single_pass( + &Arc::new(merged_index), + params, + input_plan, + quant_ctx, + tr_idx, + solo_ctx, + )?; Ok(stats) } @@ -859,6 +882,7 @@ fn run_two_pass( fn run_pass1( index: &std::sync::Arc, params: &Parameters, + input_plan: Option<&crate::io::input::ReadInputPlan>, ) -> anyhow::Result<( crate::junction::SpliceJunctionStats, Vec<( @@ -885,14 +909,19 @@ fn run_pass1( // Align reads (single-end or paired-end); no quant counting in pass 1. // Solo runs align only the cDNA read (file 0) — route to the SE path. - let n_align_files = if params.solo_enabled() { - 1 + let solo_input; + let input = if params.solo_enabled() { + solo_input = crate::io::input::ReadInputPlan::FastqSingle { + path: params.read_files_in[0].clone(), + }; + &solo_input } else { - params.read_files_in.len() + input_plan.ok_or_else(|| anyhow::anyhow!("ordinary alignment input is missing"))? }; - match n_align_files { - 1 => align_reads_single_end( + match input.layout() { + crate::io::input::ReadLayout::SingleEnd => align_reads_single_end( ¶ms_pass1, + input, index, &mut null_writer, &stats, @@ -902,8 +931,9 @@ fn run_pass1( None, None, )?, - 2 => align_reads_paired_end( + crate::io::input::ReadLayout::PairedEnd => align_reads_paired_end( ¶ms_pass1, + input, index, &mut null_writer, &stats, @@ -914,7 +944,6 @@ fn run_pass1( None, None, )?, - n => anyhow::bail!("Invalid number of read files: {n} (expected 1 or 2)"), } info!("Pass 1 aligned {} reads", stats.total_reads()); @@ -933,7 +962,7 @@ fn run_pass1( fn rc_encode(seq: &[u8]) -> Vec { seq.iter() .rev() - .map(|&b| crate::io::fastq::complement_base(b)) + .map(|&b| crate::io::reads::complement_base(b)) .collect() } @@ -1031,7 +1060,7 @@ where let take = ((max_reads - read_count) as usize).min(batch.len()); batch.truncate(take); - // Batch base offset in input (FASTQ) order — threaded to `align` so per-read + // Batch base offset in input order — threaded to `align` so per-read // indices (e.g. --outSAMreadID Number) are deterministic regardless of pool // scheduling: dispatch is sequential here, so `read_count` is this batch's base. let base = read_count; @@ -1350,6 +1379,7 @@ fn write_signal_tracks( #[allow(clippy::too_many_arguments)] fn align_reads_single_end( params: &Parameters, + input: &crate::io::input::ReadInputPlan, index: &std::sync::Arc, writer: &mut W, stats: &std::sync::Arc, @@ -1360,7 +1390,7 @@ fn align_reads_single_end( unmapped_writer: Option<&mut crate::io::fastq::UnmappedFastqWriter>, ) -> anyhow::Result<()> { use crate::align::read_align::align_read; - use crate::io::fastq::{FastqReader, clip_read}; + use crate::io::reads::clip_read; use crate::io::sam::{BufferedSamRecords, SamWriter}; use crate::params::OutFilterType; use rayon::prelude::*; @@ -1369,10 +1399,8 @@ fn align_reads_single_end( let quant = quant_ctx.map(Arc::clone); let tr = tr_idx.map(Arc::clone); - let read_file = ¶ms.read_files_in[0]; - info!("Reading single-end from {}", read_file.display()); - - let reader = FastqReader::open(read_file, params.read_files_command.as_deref())?; + info!("Reading single-end {} input", input.label()); + let producer = input.open_single_end_producer(params)?; // Create chimeric output writer if enabled let chimeric_writer = if params.chim_segment_min > 0 && params.chim_out_junctions() { @@ -1393,6 +1421,7 @@ fn align_reads_single_end( } else { params.read_map_number as u64 }; + let producer_max_records = usize::try_from(max_reads).unwrap_or(usize::MAX); let batch_size = 10000; let max_multimaps = params.out_filter_multimap_nmax as usize; @@ -1424,9 +1453,9 @@ fn align_reads_single_end( let bysj_meta: Vec = Vec::new(); info!("Aligning reads..."); - // Three-stage pipeline: a producer thread decodes the next FASTQ batch, rayon + // Three-stage pipeline: a producer thread decodes the next input batch, rayon // aligns the current batch in parallel, and a dedicated writer thread serializes - // SAM/BAM output — so gzip inflate, alignment, and record encoding all overlap + // SAM/BAM output — so input decoding, alignment, and record encoding all overlap // instead of running one-after-another per batch. Bounded channels (depth 2) give // backpressure. Output order is preserved: aligned batches flow through the // channel in input order and the writer consumes them in that order. @@ -1437,7 +1466,7 @@ fn align_reads_single_end( let params_arc = Arc::new(params.clone()); std::thread::scope(|scope| -> anyhow::Result<()> { let (read_tx, read_rx) = std::sync::mpsc::sync_channel::< - Result, error::Error>, + Result, error::Error>, >(2); #[allow(clippy::type_complexity)] let (res_tx, res_rx) = @@ -1445,20 +1474,8 @@ fn align_reads_single_end( // Stage 1: decode. scope.spawn(move || { - let mut reader = reader; - loop { - match reader.read_batch(batch_size) { - Ok(batch) => { - let last = batch.is_empty(); - if read_tx.send(Ok(batch)).is_err() || last { - break; - } - } - Err(e) => { - let _ = read_tx.send(Err(e)); - break; - } - } + if let Err(error) = producer.produce(batch_size, producer_max_records, &read_tx) { + let _ = read_tx.send(Err(error)); } }); @@ -1725,7 +1742,7 @@ fn align_reads_single_end( let params_arc = Arc::clone(¶ms_arc); let wasp_ctx = Arc::clone(&wasp_ctx); let align = move |base: u64, - batch: Vec| + batch: Vec| -> BatchOut { let params: &Parameters = ¶ms_arc; // Adapter-aware clip params (fixed 5'/3' Nbases + 3' adapter Hamming @@ -1736,7 +1753,7 @@ fn align_reads_single_end( .enumerate() .map(|(read_idx, read)| { // --outSAMreadID Number: replace the output QNAME with the - // read's 1-based input index (deterministic — from the FASTQ + // read's 1-based input index (deterministic — from the input // order via `base`, not parallel execution order). The seed // name passed to align_read is left as the real read name. let out_read_name = @@ -2002,7 +2019,7 @@ fn align_reads_solo( solo_ctx: &std::sync::Arc, ) -> anyhow::Result<()> { use crate::align::read_align::align_read; - use crate::io::fastq::clip_read; + use crate::io::reads::clip_read; use crate::io::sam::{BufferedSamRecords, SamWriter}; use crate::solo::{SoloCountRecord, SoloMultiRecord}; use rayon::prelude::*; @@ -2316,7 +2333,7 @@ fn align_reads_solo_pe( solo_ctx: &std::sync::Arc, ) -> anyhow::Result<()> { use crate::align::read_align::{PairedAlignment, PairedAlignmentResult, align_paired_read}; - use crate::io::fastq::clip_read; + use crate::io::reads::clip_read; use crate::io::sam::{BufferedSamRecords, SamWriter}; use crate::solo::{SoloCountRecord, SoloMultiRecord}; use rayon::prelude::*; @@ -2681,6 +2698,7 @@ fn align_reads_solo_pe( #[allow(clippy::too_many_arguments)] fn align_reads_paired_end( params: &Parameters, + input: &crate::io::input::ReadInputPlan, index: &std::sync::Arc, writer: &mut W, stats: &std::sync::Arc, @@ -2692,7 +2710,7 @@ fn align_reads_paired_end( unmapped_writer2: Option<&mut crate::io::fastq::UnmappedFastqWriter>, ) -> anyhow::Result<()> { use crate::align::read_align::{PairedAlignment, PairedAlignmentResult, align_paired_read}; - use crate::io::fastq::{PairedFastqReader, clip_read}; + use crate::io::reads::clip_read; use crate::io::sam::{BufferedSamRecords, SamWriter}; use crate::params::OutFilterType; use rayon::prelude::*; @@ -2701,17 +2719,8 @@ fn align_reads_paired_end( let quant = quant_ctx.map(Arc::clone); let tr = tr_idx.map(Arc::clone); - info!( - "Reading paired-end from {} and {}", - params.read_files_in[0].display(), - params.read_files_in[1].display() - ); - - let reader = PairedFastqReader::open( - ¶ms.read_files_in[0], - ¶ms.read_files_in[1], - params.read_files_command.as_deref(), - )?; + info!("Reading paired-end {} input", input.label()); + let producer = input.open_paired_end_producer(params)?; // Create chimeric output writer if enabled let chimeric_writer = if params.chim_segment_min > 0 && params.chim_out_junctions() { @@ -2732,6 +2741,7 @@ fn align_reads_paired_end( } else { params.read_map_number as u64 }; + let producer_max_records = usize::try_from(max_reads).unwrap_or(usize::MAX); let batch_size = 10000; let max_multimaps = params.out_filter_multimap_nmax as usize; @@ -2764,7 +2774,7 @@ fn align_reads_paired_end( info!("Aligning paired-end reads..."); // Three-stage pipeline: producer decodes the next pair batch, rayon aligns the // current batch, and a dedicated writer thread serializes output — overlapping - // gzip inflate of both mate files, alignment, and record encoding. Bounded + // input decoding, alignment, and record encoding. Bounded // channels (depth 2) give backpressure; output order is preserved. let stats_writer = Arc::clone(&stats); let sj_stats_writer = Arc::clone(&sj_stats); @@ -2773,7 +2783,7 @@ fn align_reads_paired_end( let params_arc = Arc::new(params.clone()); std::thread::scope(|scope| -> anyhow::Result<()> { let (read_tx, read_rx) = std::sync::mpsc::sync_channel::< - Result, error::Error>, + Result, error::Error>, >(2); #[allow(clippy::type_complexity)] let (res_tx, res_rx) = @@ -2781,20 +2791,8 @@ fn align_reads_paired_end( // Stage 1: decode. scope.spawn(move || { - let mut reader = reader; - loop { - match reader.read_paired_batch(batch_size) { - Ok(batch) => { - let last = batch.is_empty(); - if read_tx.send(Ok(batch)).is_err() || last { - break; - } - } - Err(e) => { - let _ = read_tx.send(Err(e)); - break; - } - } + if let Err(error) = producer.produce(batch_size, producer_max_records, &read_tx) { + let _ = read_tx.send(Err(error)); } }); @@ -3042,7 +3040,7 @@ fn align_reads_paired_end( let params_arc = Arc::clone(¶ms_arc); let wasp_ctx = Arc::clone(&wasp_ctx); let align = move |base: u64, - batch: Vec| + batch: Vec| -> BatchOut { let params: &Parameters = ¶ms_arc; // Adapter-aware clip params per mate (fixed Nbases are per-mate; the diff --git a/src/params/mod.rs b/src/params/mod.rs index 8114210..2a4c2a7 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -534,6 +534,10 @@ pub struct Parameters { #[arg(long = "readFilesCommand")] pub read_files_command: Option, + /// Number of CBQ decoder threads; 0 selects a conservative automatic value + #[arg(long = "readFilesNthreads", default_value_t = 0)] + pub read_files_n_threads: usize, + /// `--soloType SmartSeq` manifest: a TSV with `read1 read2 cellID` /// per line (`read2` = `-` for single-end). Each line is one plate-well cell; /// reads are counted per gene with no UMI. @@ -2127,6 +2131,16 @@ mod tests { assert_eq!(p.twopass_mode, TwopassMode::Basic); } + #[test] + fn cbq_decoder_thread_parameter() { + let default = try_parse(&["--readFilesIn", "reads.cbq"]).unwrap(); + assert_eq!(default.read_files_n_threads, 0); + + let explicit = + try_parse(&["--readFilesIn", "reads.cbq", "--readFilesNthreads", "6"]).unwrap(); + assert_eq!(explicit.read_files_n_threads, 6); + } + #[test] fn clip_nbases_per_mate() { // A single value applies to both mates. diff --git a/src/solo/count.rs b/src/solo/count.rs index 4ddb414..21801ea 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -2329,7 +2329,7 @@ fn open_maybe_gz(path: &Path) -> Result, Error> { #[cfg(test)] mod tests { use super::*; - use crate::io::fastq::encode_base; + use crate::io::reads::encode_base; use crate::solo::whitelist::pack_barcode; #[test] diff --git a/src/solo/mod.rs b/src/solo/mod.rs index 22d63e0..d603e24 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -29,7 +29,8 @@ pub use whitelist::{ use crate::align::transcript::Transcript; use crate::error::Error; -use crate::io::fastq::{EncodedRead, FastqReader, decode_base}; +use crate::io::fastq::FastqReader; +use crate::io::reads::{EncodedRead, decode_base}; use crate::params::{Parameters, SoloType}; use crate::quant::GeneAnnotation; use std::path::Path; @@ -408,7 +409,7 @@ pub fn clip_adapter_cr4(seq: &[u8], qual: &[u8]) -> (Vec, Vec, usize, us if seq.len() >= TSO_SEQ.len() { let tso: Vec = TSO_SEQ .iter() - .map(|&b| crate::io::fastq::encode_base(b)) + .map(|&b| crate::io::reads::encode_base(b)) .collect(); let mismatches = seq[..tso.len()] .iter() @@ -1020,7 +1021,7 @@ impl SoloContext { #[cfg(test)] mod tests { use super::*; - use crate::io::fastq::encode_base; + use crate::io::reads::encode_base; fn encoded_read(name: &str, seq: &str, qual: &str) -> EncodedRead { EncodedRead { diff --git a/src/solo/whitelist.rs b/src/solo/whitelist.rs index 1d882ae..278ad84 100644 --- a/src/solo/whitelist.rs +++ b/src/solo/whitelist.rs @@ -12,7 +12,7 @@ //! quality, exactly as STAR's `cbMatchString`. use crate::error::Error; -use crate::io::fastq::{decode_base, encode_base}; +use crate::io::reads::{decode_base, encode_base}; use flate2::read::GzDecoder; use std::fs::File; use std::io::{BufRead, BufReader}; diff --git a/src/wasp/mod.rs b/src/wasp/mod.rs index 78f5051..77a9f14 100644 --- a/src/wasp/mod.rs +++ b/src/wasp/mod.rs @@ -34,7 +34,7 @@ use crate::align::read_align::{ use crate::align::transcript::Transcript; use crate::error::Error; use crate::index::GenomeIndex; -use crate::io::fastq::complement_base; +use crate::io::reads::complement_base; use crate::params::{Parameters, SamAttributes}; /// One heterozygous SNV: absolute 0-based genomic `loci`, and `nt = [ref, allele0, diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index 62ee32a..a777136 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -4,6 +4,8 @@ //! GT-AG intron structure for splice tests. use assert_cmd::cargo::cargo_bin_cmd; +use binseq::SequencingRecordBuilder; +use binseq::write::{BinseqWriterBuilder, Format}; use noodles::bam; use std::fs; use std::io::Write; @@ -132,6 +134,50 @@ fn count_sam_records(sam_path: &Path) -> usize { content.lines().filter(|l| !l.starts_with('@')).count() } +fn sam_records(sam_path: &Path) -> Vec { + fs::read_to_string(sam_path) + .unwrap() + .lines() + .filter(|line| !line.starts_with('@')) + .map(str::to_string) + .collect() +} + +type CbqMate = (String, Vec, Option>); +type CbqRecord = (String, Vec, Option>, Option); + +fn write_cbq_records(path: &Path, records: &[CbqRecord]) { + let paired = records.iter().any(|record| record.3.is_some()); + let qualities = records.iter().all(|record| { + record.2.is_some() + && record + .3 + .as_ref() + .is_none_or(|(_, _, quality)| quality.is_some()) + }); + let mut writer = BinseqWriterBuilder::new(Format::Cbq) + .paired(paired) + .quality(qualities) + .headers(true) + .block_size(512) + .build(fs::File::create(path).unwrap()) + .unwrap(); + for (name1, seq1, qual1, mate2) in records { + let mut record = SequencingRecordBuilder::default() + .s_header(name1.as_bytes()) + .s_seq(seq1) + .opt_s_qual(qual1.as_deref()); + if let Some((name2, seq2, qual2)) = mate2 { + record = record + .x_header(name2.as_bytes()) + .x_seq(seq2) + .opt_x_qual(qual2.as_deref()); + } + writer.push(record.build().unwrap()).unwrap(); + } + writer.finish().unwrap(); +} + // --------------------------------------------------------------------------- // Test 1 — BAM unsorted output // --------------------------------------------------------------------------- @@ -2025,6 +2071,238 @@ fn test_wasp_samtag() { ); } +// --------------------------------------------------------------------------- +// Native CBQ input +// --------------------------------------------------------------------------- + +#[test] +fn test_cbq_single_end_matches_fastq_and_honors_read_map_number() { + let tmpdir = TempDir::new().unwrap(); + let genome = build_genome(); + let fasta = write_fasta(&tmpdir, &genome); + let genome_dir = tmpdir.path().join("genome"); + build_index(&fasta, &genome_dir, "7", None); + + let fastq = tmpdir.path().join("reads.fq"); + let cbq = tmpdir.path().join("reads.binary-data"); + let mut logical_records = Vec::new(); + { + let mut writer = fs::File::create(&fastq).unwrap(); + for i in 0..31usize { + let start = 200 + i * 200; + let sequence = genome[start..start + 50].to_vec(); + let name = format!("cbq-se-{i}"); + let quality = vec![b'I'; sequence.len()]; + writeln!(writer, "@{name}").unwrap(); + writer.write_all(&sequence).unwrap(); + writeln!(writer, "\n+\n{}", "I".repeat(sequence.len())).unwrap(); + logical_records.push((name, sequence, Some(quality), None)); + } + } + write_cbq_records(&cbq, &logical_records); + + let fastq_out = tmpdir.path().join("fastq-out"); + let cbq_out = tmpdir.path().join("cbq-out"); + fs::create_dir_all(&fastq_out).unwrap(); + fs::create_dir_all(&cbq_out).unwrap(); + for (input, output, decoder_threads) in + [(&fastq, &fastq_out, None), (&cbq, &cbq_out, Some("4"))] + { + let prefix = format!("{}/", output.display()); + let mut command = cargo_bin_cmd!("rustar-aligner"); + command.args([ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + input.to_str().unwrap(), + "--readMapNumber", + "19", + "--outSAMtype", + "SAM", + "--outFileNamePrefix", + &prefix, + ]); + if let Some(threads) = decoder_threads { + command.args(["--readFilesNthreads", threads]); + } + command.assert().success(); + } + + let fastq_records = sam_records(&fastq_out.join("Aligned.out.sam")); + let cbq_records = sam_records(&cbq_out.join("Aligned.out.sam")); + assert_eq!(fastq_records.len(), 19); + assert_eq!(cbq_records, fastq_records); +} + +#[test] +fn test_qualityless_paired_cbq_two_pass_output_is_ordered_and_has_missing_qual() { + let tmpdir = TempDir::new().unwrap(); + let genome = build_genome(); + let fasta = write_fasta(&tmpdir, &genome); + let genome_dir = tmpdir.path().join("genome"); + build_index(&fasta, &genome_dir, "7", None); + + let cbq = tmpdir.path().join("paired.cbq"); + let mut logical_records = Vec::new(); + for i in 0..13usize { + let mate1_start = 500 + i * 300; + let mate2_start = mate1_start + 100; + let name = format!("cbq-pe-{i}"); + logical_records.push(( + format!("{name}/1"), + genome[mate1_start..mate1_start + 50].to_vec(), + None, + Some(( + format!("{name}/2"), + rc(&genome[mate2_start..mate2_start + 50]), + None, + )), + )); + } + write_cbq_records(&cbq, &logical_records); + + let output = tmpdir.path().join("paired-out"); + fs::create_dir_all(&output).unwrap(); + let prefix = format!("{}/", output.display()); + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + cbq.to_str().unwrap(), + "--readFilesNthreads", + "3", + "--twopassMode", + "Basic", + "--outSAMtype", + "SAM", + "--outFileNamePrefix", + &prefix, + ]) + .assert() + .success(); + + let records = sam_records(&output.join("Aligned.out.sam")); + assert_eq!(records.len(), 26); + for (pair_index, pair) in records.chunks_exact(2).enumerate() { + let first: Vec<_> = pair[0].split('\t').collect(); + let second: Vec<_> = pair[1].split('\t').collect(); + assert_eq!(first[0], format!("cbq-pe-{pair_index}")); + assert_eq!(second[0], first[0]); + assert_eq!(first[10], "*"); + assert_eq!(second[10], "*"); + } +} + +// --------------------------------------------------------------------------- +// Coordinate-sort spilling (--limitBAMsortRAM) +// --------------------------------------------------------------------------- + +/// A low `--limitBAMsortRAM` must spill sorted runs to disk and merge them into +/// a BAM identical to the unbounded in-memory sort — same records, same order. +#[test] +fn test_sorted_bam_spills_to_disk_and_matches_unbounded_sort() { + let tmpdir = TempDir::new().unwrap(); + let genome = build_genome(); + let fasta = write_fasta(&tmpdir, &genome); + let genome_dir = tmpdir.path().join("genome"); + build_index(&fasta, &genome_dir, "7", None); + + // 2000 reads across the genome, emitted in an order that does not match + // coordinate order so the sort is actually doing work. + let fastq_path = tmpdir.path().join("reads.fq"); + { + let mut f = fs::File::create(&fastq_path).unwrap(); + for i in 0..2000usize { + let start = (i * 7919) % (genome.len() - 60); + let seq = &genome[start..start + 50]; + writeln!(f, "@read{i}").unwrap(); + f.write_all(seq).unwrap(); + writeln!(f, "\n+\n{}", "I".repeat(50)).unwrap(); + } + } + + let run = |label: &str, limit: &str| -> (PathBuf, PathBuf) { + let output_dir = tmpdir.path().join(label); + fs::create_dir_all(&output_dir).unwrap(); + let prefix = format!("{}/", output_dir.display()); + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + fastq_path.to_str().unwrap(), + "--outSAMtype", + "BAM", + "SortedByCoordinate", + "--limitBAMsortRAM", + limit, + "--outFileNamePrefix", + &prefix, + ]) + .assert() + .success(); + (output_dir.join("Aligned.sortedByCoord.out.bam"), output_dir) + }; + + // 64 KiB forces many spill runs; 1 GiB holds everything in memory. + let (spilled_bam, spilled_dir) = run("out_sort_spill", "65536"); + let (memory_bam, _) = run("out_sort_memory", "1G"); + + let read_bam = |path: &PathBuf| -> Vec<(String, Option, Option, String)> { + let mut reader = bam::io::Reader::new(fs::File::open(path).unwrap()); + reader.read_header().unwrap(); + reader + .records() + .map(|record| { + let record = record.unwrap(); + ( + String::from_utf8(record.name().unwrap().to_vec()).unwrap(), + record.reference_sequence_id().transpose().unwrap(), + record + .alignment_start() + .transpose() + .unwrap() + .map(|p| p.get()), + format!("{:?}", record.cigar().iter().collect::>()), + ) + }) + .collect() + }; + + let spilled = read_bam(&spilled_bam); + let in_memory = read_bam(&memory_bam); + + assert!(!spilled.is_empty(), "expected alignments"); + assert_eq!( + spilled, in_memory, + "spill+merge output must equal the unbounded in-memory sort" + ); + + // Coordinate-sorted, and no spill scratch left behind. + let keys: Vec<_> = spilled + .iter() + .map(|(_, chr, pos, _)| (chr.unwrap_or(usize::MAX), pos.unwrap_or(0))) + .collect(); + assert!( + keys.windows(2).all(|w| w[0] <= w[1]), + "output must be non-decreasing by (chr, pos)" + ); + let leftover: Vec<_> = fs::read_dir(&spilled_dir) + .unwrap() + .filter_map(Result::ok) + .map(|e| e.file_name().to_string_lossy().to_string()) + .filter(|name| name.starts_with("rustar-bamsort-")) + .collect(); + assert!(leftover.is_empty(), "spill files left behind: {leftover:?}"); +} + // --------------------------------------------------------------------------- // --runMode soloCellFiltering // --------------------------------------------------------------------------- From e9bd3f6c8cd5d8f62e301b2798f6f978ab4b2211 Mon Sep 17 00:00:00 2001 From: alejandrogzi Date: Thu, 6 Aug 2026 23:17:34 +0200 Subject: [PATCH 2/3] test: drop coordinate-sort spill test from CBQ feature The spill+merge test exercises the bounded-memory coordinate sorter, not CBQ input; it ships with the BAM sorter change instead. --- tests/alignment_features.rs | 105 ------------------------------------ 1 file changed, 105 deletions(-) diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index a777136..b1bf40b 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -2198,111 +2198,6 @@ fn test_qualityless_paired_cbq_two_pass_output_is_ordered_and_has_missing_qual() } } -// --------------------------------------------------------------------------- -// Coordinate-sort spilling (--limitBAMsortRAM) -// --------------------------------------------------------------------------- - -/// A low `--limitBAMsortRAM` must spill sorted runs to disk and merge them into -/// a BAM identical to the unbounded in-memory sort — same records, same order. -#[test] -fn test_sorted_bam_spills_to_disk_and_matches_unbounded_sort() { - let tmpdir = TempDir::new().unwrap(); - let genome = build_genome(); - let fasta = write_fasta(&tmpdir, &genome); - let genome_dir = tmpdir.path().join("genome"); - build_index(&fasta, &genome_dir, "7", None); - - // 2000 reads across the genome, emitted in an order that does not match - // coordinate order so the sort is actually doing work. - let fastq_path = tmpdir.path().join("reads.fq"); - { - let mut f = fs::File::create(&fastq_path).unwrap(); - for i in 0..2000usize { - let start = (i * 7919) % (genome.len() - 60); - let seq = &genome[start..start + 50]; - writeln!(f, "@read{i}").unwrap(); - f.write_all(seq).unwrap(); - writeln!(f, "\n+\n{}", "I".repeat(50)).unwrap(); - } - } - - let run = |label: &str, limit: &str| -> (PathBuf, PathBuf) { - let output_dir = tmpdir.path().join(label); - fs::create_dir_all(&output_dir).unwrap(); - let prefix = format!("{}/", output_dir.display()); - cargo_bin_cmd!("rustar-aligner") - .args([ - "--runMode", - "alignReads", - "--genomeDir", - genome_dir.to_str().unwrap(), - "--readFilesIn", - fastq_path.to_str().unwrap(), - "--outSAMtype", - "BAM", - "SortedByCoordinate", - "--limitBAMsortRAM", - limit, - "--outFileNamePrefix", - &prefix, - ]) - .assert() - .success(); - (output_dir.join("Aligned.sortedByCoord.out.bam"), output_dir) - }; - - // 64 KiB forces many spill runs; 1 GiB holds everything in memory. - let (spilled_bam, spilled_dir) = run("out_sort_spill", "65536"); - let (memory_bam, _) = run("out_sort_memory", "1G"); - - let read_bam = |path: &PathBuf| -> Vec<(String, Option, Option, String)> { - let mut reader = bam::io::Reader::new(fs::File::open(path).unwrap()); - reader.read_header().unwrap(); - reader - .records() - .map(|record| { - let record = record.unwrap(); - ( - String::from_utf8(record.name().unwrap().to_vec()).unwrap(), - record.reference_sequence_id().transpose().unwrap(), - record - .alignment_start() - .transpose() - .unwrap() - .map(|p| p.get()), - format!("{:?}", record.cigar().iter().collect::>()), - ) - }) - .collect() - }; - - let spilled = read_bam(&spilled_bam); - let in_memory = read_bam(&memory_bam); - - assert!(!spilled.is_empty(), "expected alignments"); - assert_eq!( - spilled, in_memory, - "spill+merge output must equal the unbounded in-memory sort" - ); - - // Coordinate-sorted, and no spill scratch left behind. - let keys: Vec<_> = spilled - .iter() - .map(|(_, chr, pos, _)| (chr.unwrap_or(usize::MAX), pos.unwrap_or(0))) - .collect(); - assert!( - keys.windows(2).all(|w| w[0] <= w[1]), - "output must be non-decreasing by (chr, pos)" - ); - let leftover: Vec<_> = fs::read_dir(&spilled_dir) - .unwrap() - .filter_map(Result::ok) - .map(|e| e.file_name().to_string_lossy().to_string()) - .filter(|name| name.starts_with("rustar-bamsort-")) - .collect(); - assert!(leftover.is_empty(), "spill files left behind: {leftover:?}"); -} - // --------------------------------------------------------------------------- // --runMode soloCellFiltering // --------------------------------------------------------------------------- From cb0ebb1bd8f143feb808656a4218d5af48c08104 Mon Sep 17 00:00:00 2001 From: alejandrogzi Date: Thu, 6 Aug 2026 23:18:13 +0200 Subject: [PATCH 3/3] fix(solo): use reads::encode_base now that the helper moved out of fastq Upstream added a call to fastq::encode_base inside index_of_barcode after the fork diverged; the move of base helpers to io::reads requires the updated path for the CBQ-era import layout to keep compiling. --- src/solo/whitelist.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solo/whitelist.rs b/src/solo/whitelist.rs index 278ad84..1b7fec1 100644 --- a/src/solo/whitelist.rs +++ b/src/solo/whitelist.rs @@ -412,7 +412,7 @@ impl CbWhitelist { pub fn index_of_barcode(&self, ascii: &[u8]) -> Option { let codes: Vec = ascii .iter() - .map(|&b| crate::io::fastq::encode_base(b)) + .map(|&b| crate::io::reads::encode_base(b)) .collect(); match pack_barcode(&codes) { PackResult::NoN(packed) => self.search(packed),