diff --git a/.github/workflows/release-pypi.yaml b/.github/workflows/release-pypi.yaml index 58bf114..ee7c1a3 100644 --- a/.github/workflows/release-pypi.yaml +++ b/.github/workflows/release-pypi.yaml @@ -21,7 +21,7 @@ jobs: runs-on: ${{ matrix.os }} steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v7 # portable SIMD floor (x86-64-v3/AVX2, apple-m1) - never target-cpu=native for wheels - name: Use portable cargo config @@ -45,7 +45,7 @@ jobs: sdist: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v7 - name: Build sdist uses: PyO3/maturin-action@v1 with: diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index 7204905..9034f76 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -26,7 +26,7 @@ jobs: runs-on: ${{ matrix.os }} steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v7 - uses: dtolnay/rust-toolchain@stable - uses: Swatinem/rust-cache@v2 - name: Build & test & rename @@ -65,7 +65,7 @@ jobs: # Skip prerelease tags (SemVer prereleases contain a hyphen, e.g. 0.17.0-rc.1). if: ${{ !contains(github.ref_name, '-') }} steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v7 - uses: dtolnay/rust-toolchain@stable - name: Publish deacon to crates.io run: cargo publish -p deacon diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 8b06e63..5858a15 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -16,7 +16,7 @@ jobs: os: [ubuntu-24.04, ubuntu-24.04-arm, macos-15] steps: - name: Checkout repository - uses: actions/checkout@v4 + uses: actions/checkout@v7 - name: Install Rust toolchain uses: dtolnay/rust-toolchain@stable - name: Cache cargo registry @@ -28,3 +28,27 @@ jobs: # Not a default workspace member, so cargo test skips it - name: Run deacon-wasm unit tests run: cargo test -p deacon-wasm + + python-test: + name: Python bindings + runs-on: ubuntu-24.04 + steps: + - name: Checkout repository + uses: actions/checkout@v7 + - name: Install Rust toolchain + uses: dtolnay/rust-toolchain@stable + - name: Set up Python + uses: actions/setup-python@v7 + with: + python-version: "3.12" + - name: Cache cargo registry + uses: Swatinem/rust-cache@v2 + - name: Use portable config + run: mv .cargo/config-portable.toml .cargo/config.toml + - name: Build bindings + run: | + python -m venv .venv + .venv/bin/python -m pip install "maturin>=1.7,<2" + .venv/bin/maturin develop -m deacon-py/Cargo.toml + - name: Run Python tests + run: .venv/bin/python -m unittest discover -s deacon-py/tests -v diff --git a/CHANGELOG.md b/CHANGELOG.md index ea5ef0e..399b507 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,19 +5,23 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## [0.17.0] - 2026-08-11 +## [0.17.0] - 2026-08-18 ### Added -- `deacon filter --ordered` uses paraseq 0.5.0 order preservation for deterministic output. Also exposed as `ordered` in the Python bindings and JSON summary. +- `deacon filter --ordered` preserves output order via paraseq 0.5.0. Arg `ordered` added to Python bindings and JSON summary. +- `deacon filter` accepts CBQ (BINSEQ) input and writes CBQ when given a `.cbq` path. Single and paired CBQ is supported. Output block size is the greater of `--cbq-block-size` (default 16 MiB) and the input block size given CBQ input. By convention an alternate `.cba` output extension writes CBQ without quality scores, implying `--discard-quality`. +- Python `Index.filter()` now exposes `check_pairs`, `summary`, and `cbq_block_size`, and accepts path-like objects for all paths. ### Changed -- Python `Index.filter()` options after `fastq2` are now keyword-only. +- **Breaking (Python):** index/filter input paths are positional-only; options and all `Index.fetch()` arguments are keyword-only. `filter` also renames `fastq`/`fastq2` to `input`/`input2`. +- Python filtering rejects `abs_threshold=0`, matching the CLI. +- Renamed `--fasta` (`-f`)/`output_fasta` to `--discard-quality`/`discard_quality`. ### Removed -- `deacon filter --rename-random`, and `rename_random` from the JSON summary and Python bindings. +- Removed `--rename-random` CLI arg and `rename_random` from the JSON summary and Python bindings. ## [0.16.0] - 2026-08-09 diff --git a/Cargo.lock b/Cargo.lock index 2d5dd3f..2d72798 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -94,6 +94,17 @@ dependencies = [ "wait-timeout", ] +[[package]] +name = "auto_impl" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ffdcb70bdbc4d478427380519163274ac86e52916e10f0a8889adf0f96d3fee7" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.119", +] + [[package]] name = "autocfg" version = "1.5.1" @@ -155,6 +166,28 @@ dependencies = [ "virtue", ] +[[package]] +name = "binseq" +version = "0.9.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5fea1e47c3296e0073dc373f24660b1bd8ec36d6c06a7de115e4c75bf5b1753d" +dependencies = [ + "anyhow", + "auto_impl", + "bitnuc 0.4.1", + "bitnuc 0.5.2", + "bytemuck", + "byteorder", + "itoa", + "memchr", + "memmap2", + "num_cpus", + "rand 0.9.5", + "sucds", + "thiserror 2.0.20", + "zstd", +] + [[package]] name = "bitflags" version = "1.3.2" @@ -167,6 +200,22 @@ 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 = "bitnuc" +version = "0.5.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "93612e75f32ab9cca9750b822f77c5dcf213782cfb67e12f6c8cf91679cb112a" +dependencies = [ + "fearless_simd", + "thiserror 2.0.20", +] + [[package]] name = "bstr" version = "1.13.1" @@ -189,6 +238,20 @@ name = "bytemuck" version = "1.25.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "95832e849adfb21180ccb6826a99da14e5d266ae5c2e668e1602cf234f153797" +dependencies = [ + "bytemuck_derive", +] + +[[package]] +name = "bytemuck_derive" +version = "1.12.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fc0e56a716f1e132ff6bf4bdac1c944a3fcdc1cae65f70a4a2a1ac3b401d2d1f" +dependencies = [ + "proc-macro2", + "quote", + "syn 3.0.3", +] [[package]] name = "byteorder" @@ -249,7 +312,7 @@ checksum = "d524456ba66e72eb8b115ff89e01e497f8e6d11d78b70b1aa13c0fbd97540a81" dependencies = [ "cfg-if", "cpufeatures", - "rand_core", + "rand_core 0.10.1", ] [[package]] @@ -425,6 +488,7 @@ dependencies = [ "anyhow", "assert_cmd", "bincode", + "binseq", "clap", "ensure_simd", "flate2", @@ -582,6 +646,12 @@ version = "2.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "da7c62ceae207dd37ea5b845da6a0696c799f85e97da1ab5b7910be3c1c80223" +[[package]] +name = "fearless_simd" +version = "0.7.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f4beca3cb2444e3304ac30843cc091f44ed58932353cd492ce740067bfce6b12" + [[package]] name = "find-msvc-tools" version = "0.1.10" @@ -693,6 +763,18 @@ dependencies = [ "wasm-bindgen", ] +[[package]] +name = "getrandom" +version = "0.3.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "899def5c37c4fd7b2664648c28120ecec138e4d395b459e5ca34f9cce2dd77fd" +dependencies = [ + "cfg-if", + "libc", + "r-efi 5.3.0", + "wasip2", +] + [[package]] name = "getrandom" version = "0.4.3" @@ -702,8 +784,8 @@ dependencies = [ "cfg-if", "js-sys", "libc", - "r-efi", - "rand_core", + "r-efi 6.0.0", + "rand_core 0.10.1", "wasm-bindgen", ] @@ -1010,6 +1092,15 @@ version = "2.8.3" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "cf8baf1c55e62ffcace7a9f06f4bd9cd3f0c4beb022d3b367256b91b87513d98" +[[package]] +name = "memmap2" +version = "0.9.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d1219ed1b7f229ee7104d281dd01d6802fe28bb6e95d292942c4daacdeb798c0" +dependencies = [ + "libc", +] + [[package]] name = "miniz_oxide" version = "0.8.9" @@ -1118,7 +1209,7 @@ dependencies = [ "cfg-if", "ensure_simd", "mem_dbg", - "rand", + "rand 0.10.2", "wide", ] @@ -1198,6 +1289,15 @@ dependencies = [ "portable-atomic", ] +[[package]] +name = "ppv-lite86" +version = "0.2.21" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "85eae3c4ed2f50dcfe72643da4befc30deadb458a9b590d720cde2f2b1e97da9" +dependencies = [ + "zerocopy", +] + [[package]] name = "predicates" version = "3.1.4" @@ -1334,12 +1434,28 @@ dependencies = [ "proc-macro2", ] +[[package]] +name = "r-efi" +version = "5.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "69cdb34c158ceb288df11e18b4bd39de994f6657d83847bdffdbd7f346754b0f" + [[package]] name = "r-efi" version = "6.0.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "f8dcc9c7d52a811697d2151c701e0d08956f92b0e24136cf4cf27b57a6a0d9bf" +[[package]] +name = "rand" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b9ef1d0d795eb7d84685bca4f72f3649f064e6641543d3a8c415898726a57b41" +dependencies = [ + "rand_chacha", + "rand_core 0.9.5", +] + [[package]] name = "rand" version = "0.10.2" @@ -1348,7 +1464,26 @@ checksum = "c7f5fa3a058cd35567ef9bfa5e75732bee0f9e4c55fa90477bef2dfcdbc4be80" dependencies = [ "chacha20", "getrandom 0.4.3", - "rand_core", + "rand_core 0.10.1", +] + +[[package]] +name = "rand_chacha" +version = "0.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d3022b5f1df60f26e1ffddd6c66e8aa15de382ae63b3a0c1bfc0e4d3e3f325cb" +dependencies = [ + "ppv-lite86", + "rand_core 0.9.5", +] + +[[package]] +name = "rand_core" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "76afc826de14238e6e8c374ddcc1fa19e374fd8dd986b0d2af0d02377261d83c" +dependencies = [ + "getrandom 0.3.4", ] [[package]] @@ -1766,6 +1901,16 @@ version = "2.6.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "13c2bddecc57b384dee18652358fb23172facb8a2c51ccc10d74c157bdea3292" +[[package]] +name = "sucds" +version = "0.8.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cd324eaa05be64f105ea5269bb8aabd70e5dd57fa5c673b167f451b07d6c0dcd" +dependencies = [ + "anyhow", + "num-traits", +] + [[package]] name = "syn" version = "2.0.119" @@ -1950,6 +2095,15 @@ version = "0.11.1+wasi-snapshot-preview1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "ccf3ec651a847eb01de73ccad15eb7d99f80485de043efb2f370cd654f4ea44b" +[[package]] +name = "wasip2" +version = "1.0.4+wasi-0.2.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b67efb37e106e55ce722a510d6b5f9c17f083e5fc79afc2badeb12cc313d9487" +dependencies = [ + "wit-bindgen", +] + [[package]] name = "wasm-bindgen" version = "0.2.127" @@ -2228,6 +2382,12 @@ dependencies = [ "memchr", ] +[[package]] +name = "wit-bindgen" +version = "0.57.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1ebf944e87a7c253233ad6766e082e3cd714b5d03812acc24c318f549614536e" + [[package]] name = "xorf" version = "0.12.0" @@ -2238,6 +2398,26 @@ dependencies = [ "libm", ] +[[package]] +name = "zerocopy" +version = "0.8.56" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "556764e583adb45a9f8d413c2a147fa7e8d821e48e12b14fd560b607998b75eb" +dependencies = [ + "zerocopy-derive", +] + +[[package]] +name = "zerocopy-derive" +version = "0.8.56" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f2ab42fc20575779bd240faa45f94a74256f755c0fa9e89f0ede20d91d0cdfc1" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.119", +] + [[package]] name = "zeroize" version = "1.9.0" diff --git a/Cargo.toml b/Cargo.toml index 22c6498..626f006 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -41,6 +41,9 @@ indicatif = { version = "0.18", optional = true } serde_json = { version = "1.0", optional = true } niffler = { version = "3.0.1", default-features = false, optional = true } paraseq = { version = "0.5.0", default-features = false, features = ["anyhow", "niffler"], optional = true } +# CBQ support. Defaults disabled since we don't need binseq's own paraseq integration +# We need binseq >=0.9.6 for safe canonicalisation of non-N ambiguous bases to N +binseq = { version = "0.9.6", default-features = false, features = ["anyhow"], optional = true } ensure_simd = { version = "0.1.0", optional = true } # Compression deps (optional) @@ -59,10 +62,11 @@ predicates = "3.0" tempfile = "3.20" rstest = "0.26" nix = { version = "0.31", features = ["fs"] } +binseq = { version = "0.9.6", default-features = false, features = ["anyhow"] } [features] scalar = ["simd-minimizers/scalar", "packed-seq/scalar"] -cli = ["rayon", "parking_lot", "indicatif", "paraseq", "ensure_simd", "clap", "niffler", "serde_json"] +cli = ["rayon", "parking_lot", "indicatif", "paraseq", "binseq", "ensure_simd", "clap", "niffler", "serde_json"] # Disable for faster builds. compression = ["zstd", "liblzma", "flate2", "gzp", "paraseq/default"] # Use to still handle .gz files when "compression" is not enabled. diff --git a/README.md b/README.md index 8d36384..a7fe052 100644 --- a/README.md +++ b/README.md @@ -10,7 +10,7 @@ Deacon filters DNA sequences in FASTA/Q files and streams using SIMD-accelerated minimizer comparison with query sequence(s), emitting either matching sequences (**search mode**), or sequences without matches (**deplete mode**). Sequences match when they share enough distinct minimizers with the indexed query to exceed chosen absolute and relative thresholds. Query size has little impact on filtering speed, enabling ultrafast search and depletion with gene-, genome- and pangenome-scale queries using a laptop. Deacon filters uncompressed FASTA/Q at **gigabases per second** on recent AMD, Intel (`x86_64`), and Apple `arm64` systems. Built with panhuman host depletion in mind—yet broadly useful for searching large sequence collections—Deacon delivers [leading classification accuracy](https://doi.org/10.1101/2025.06.09.658732) for host depletion and unrivalled speed using 5GB of RAM. -Default parameters are carefully chosen but easily changed. Classification sensitivity, specificity and memory requirements may be tuned by varying *k*-mer length (`-k`), window size (`-w`), absolute match threshold (`-a`) and relative match threshold (`-r`) . Minimizer `k` and `w` are chosen at query index time, while the match thresholds can be chosen at filter time. Matching sequences are those that share enough distinct minimizers with the indexed query to exceed *both* the absolute threshold (`-a`, default 2 shared minimizers) and the relative threshold (`-r`, default 0.01 [1%] shared minimizers). For paired sequences, hits in either mate counts towards a single match threshold for the pair. Deacon reports filtering performance during execution and optionally writes a JSON `--summary` upon completion. Sequences can optionally be renamed using `--rename` for privacy and smaller file sizes. Deacon fully supports stdin, stdout and natively handles gz, zst and xz compression formats, detected by file extension. +Default parameters are carefully chosen but easily changed. Classification sensitivity, specificity and memory requirements may be tuned by varying *k*-mer length (`-k`), window size (`-w`), absolute match threshold (`-a`) and relative match threshold (`-r`) . Minimizer `k` and `w` are chosen at query index time, while the match thresholds can be chosen at filter time. Matching sequences are those that share enough distinct minimizers with the indexed query to exceed *both* the absolute threshold (`-a`, default 2 shared minimizers) and the relative threshold (`-r`, default 0.01 [1%] shared minimizers). For paired sequences, hits in either mate counts towards a single match threshold for the pair. Deacon reports filtering performance during execution and optionally writes a JSON `--summary` upon completion. Sequences can optionally be renamed using `--rename` for privacy and smaller file sizes. Deacon fully supports stdin/out (uncompressed fastx) and natively handles .gz, .zst and .xz fastx file IO, as well as BINSEQ CBQ (.cbq & .cba). Benchmarks for panhuman host depletion of complex microbial metagenomes are described in a [preprint](https://www.biorxiv.org/content/10.1101/2025.06.09.658732v1). Deacon with the `panhuman-1` (*k*=31, w=15) index exhibited the highest balanced accuracy for both long and short simulated reads. Deacon was less specific only than Hostile for short reads. @@ -89,7 +89,7 @@ Prebuilt pangenome indexes are provided. These can be downloaded using the links ### Filtering -The main command `deacon filter` accepts an index path followed by up to two FASTA/FASTQ file paths, depending on whether input sequences originate from stdin, a single file, or paired input files. Indexes are built with `deacon index build`. Paired inputs are supported as either two separate or one interleaved file/stream when using `--interleaved`, and may be written either to separate paired output files or one interleaved file. For paired sequences, distinct minimizer hits originating from either mate are counted. By default, input sequences must meet both an absolute threshold of 2 minimizer hits (`-a 2`) and a relative threshold of 1% of minimizers (`-r 0.01`) to pass the filter. Filtering can be inverted for e.g. host depletion using the `--deplete` (`-d`) flag. Gzip, Zstandard, and xz compression formats are detected automatically by file extension. Paired read headers can be validated using `--check-pairs`. +The main command `deacon filter` accepts an index path followed by up to two sequence file paths, depending on whether input sequences originate from stdin, a single file, or paired input files. Indexes are built with `deacon index build`. Paired inputs are supported as either two separate or one interleaved file/stream when using `--interleaved`, and may be written either to separate paired output files or one interleaved file. For paired sequences, distinct minimizer hits originating from either mate are counted. Paired read headers can be validated using `--check-pairs`. By default, input sequences must meet both an absolute threshold of 2 minimizer hits (`-a 2`) and a relative threshold of 1% of minimizers (`-r 0.01`) to pass the filter. Filtering can be inverted for e.g. host depletion using the `--deplete` (`-d`) flag. Gzip, Zstandard, and xz compressed FASTX formats are detected automatically by file extension. Single and paired [BINSEQ](https://www.biorxiv.org/content/10.1101/2025.04.08.647863v2) CBQ files are natively supported for both input and output. CBQ input is detected by content rather than file extension, while a `.cbq` output path writes CBQ, and a `.cba` output path writes CBQ without quality. #### Examples @@ -130,6 +130,14 @@ zcat r12.fq.gz | deacon filter -d panhuman-1.k31w15.idx - - > filt12.fq # Save summary JSON deacon filter -d panhuman-1.k31w15.idx reads.fq.gz -o filt.fq.gz -s summary.json +# BINSEQ CBQ input and output; paired records are stored in one file +deacon filter -d panhuman-1.k31w15.idx reads.fq.gz -o filt.cbq +deacon filter -d panhuman-1.k31w15.idx r1.fq.gz r2.fq.gz -o filt12.cbq +deacon filter -d panhuman-1.k31w15.idx filt12.cbq -o filt12.fq.gz + +# A .cba extension writes cbq without quality values, if present +deacon filter -d panhuman-1.k31w15.idx reads.fq.gz -o filt.cba + # Replace read headers with incrementing integers deacon filter -d -R panhuman-1.k31w15.idx reads.fq.gz > filt.fq @@ -181,13 +189,13 @@ A differentiating feature of Deacon is the ease of combining, subtracting and in ```bash $ deacon filter -h -Retain or deplete sequence records with sufficient minimizer hits to an indexed query +Retain or deplete sequence records with sufficient minimizer hits to the index Usage: deacon filter [OPTIONS] [INPUT] [INPUT2] Arguments: Path to minimizer index file - [INPUT] Optional path to fastx file (or - for stdin) [default: -] + [INPUT] Optional path to fastx or binseq cbq file (or - for stdin) [default: -] [INPUT2] Optional path to second paired fastx file Options: @@ -202,13 +210,11 @@ Options: -d, --deplete Discard matching sequences (invert filtering behaviour) -R, --rename - Replace sequence headers with incrementing numbers (reproducible with --ordered) - -f, --fasta - Output FASTA format regardless of input format + Replace sequence headers with incrementing numbers (deterministic with --ordered) -o, --output - Path to output fastx file (stdout if not specified; detects .gz and .zst) + Path to output file (fastx to stdout by default; detects .gz, .zst, .xz, .cbq, .cba) -O, --output2 - Optional path to second paired output fastx file (detects .gz and .zst) + Optional path to second paired output fastx file (detects .gz, .zst, .xz) -s, --summary Path to JSON summary output file -t, --threads @@ -216,17 +222,21 @@ Options: --compression-threads Number of threads used for output compression (0 = auto) [default: 0] --compression-level - Output compression level (1-9 for gz & xz; 1-22 for zstd) [default: 2] + Output compression level (1-9 for gz & xz; 1-22 for zstd including cbq) [default: 2] + --cbq-block-size + cbq output block size in MiB (or cbq input block size if higher) [default: 16] + --discard-quality + Emit fasta or quality-free cbq regardless of input format --interleaved - Treat INPUT as interleaved paired reads from a file or stdin + Treat INPUT as interleaved paired records from single file or stdin --ordered Preserve input record ordering (deterministic, slightly slower) --check-pairs Validate paired record names (Illumina CASAVA or /1 /2 suffixes) + --debug + Emit sequences with minimizer hits to stderr -q, --quiet Suppress progress reporting - --debug - Output sequences with minimizer hits to stderr -h, --help Print help ``` @@ -270,7 +280,7 @@ Options: -w Minimizer window size used for indexing [default: 15] -o, --output Path to output file (stdout if not specified) -t, --threads Number of execution threads (0 = auto) [default: 8] - -q, --quiet Suppress sequence header output + -q, --quiet Suppress progress reporting -h, --help Print help ``` @@ -280,9 +290,9 @@ Options: Use `-s summary.json` to save detailed filtering statistics: ```json { - "version": "deacon 0.9.0", + "version": "deacon 0.17.0", "index": "panhuman-1.k31w15.idx", - "input": "HG02334.1m.fastq.gz", + "input": "HG02334.100MB.fastq.gz", "input2": null, "output": "-", "output2": null, @@ -293,18 +303,23 @@ Use `-s summary.json` to save detailed filtering statistics: "prefix_length": 0, "deplete": true, "rename": false, + "ordered": false, "check_pairs": false, - "seqs_in": 1000000, - "seqs_out": 13452, - "seqs_removed": 986548, - "seqs_removed_proportion": 0.986548, - "bp_in": 5477122928, - "bp_out": 5710050, - "bp_removed": 5471412878, - "bp_removed_proportion": 0.9989574727324798, - "time": 125.755103875, - "seqs_per_second": 7951, - "bp_per_second": 43553881 + "seqs_in": 37500, + "seqs_out": 454, + "seqs_out_proportion": 0.012106666666666667, + "seqs_removed": 37046, + "seqs_removed_proportion": 0.9878933333333333, + "bp_in": 141474280, + "bp_out": 227079, + "bp_out_proportion": 0.001605090338682056, + "bp_removed": 141247201, + "bp_removed_proportion": 0.9983949096613179, + "time": 0.446945667, + "seqs_per_second": 84129, + "bp_per_second": 317392822, + "seqs_per_second_total": 83902, + "bp_per_second_total": 316535745 } ``` diff --git a/deacon-py/python/README.md b/deacon-py/python/README.md index 74c66c7..9eb0635 100644 --- a/deacon-py/python/README.md +++ b/deacon-py/python/README.md @@ -5,7 +5,7 @@ Python bindings for [Deacon](https://github.com/bede/deacon), enabling fast mult ## Installation ```bash -uv install deacon +uv pip install deacon ``` ## Quickstart @@ -33,47 +33,61 @@ index = Index("panhuman-1.k31w15.idx", complexity_threshold=None) Pass `complexity_threshold` (0.0–1.0, e.g. `0.9`) to discard low-complexity index minimizers once at load using kdust; the filtered set is then reused across all `filter` calls. Not supported for `bff` (binary fuse filter) indexes. +`path` is positional-only and accepts a string or path-like object. `complexity_threshold` is keyword-only. + ## `Index.fetch()` Download a prebuilt index, then load and return it (a static method, so `Index.fetch(...)` returns an `Index`). `output` is the local path to save to; when omitted it defaults to `"{name}.k{k}w{w}.idx"` in the working directory. The index is downloaded on every call — there is no local cache, so an existing file at that path is overwritten. ```python -index = Index.fetch("panhuman-1", k=31, w=15, output=None, complexity_threshold=None) +index = Index.fetch( + name="panhuman-1", + k=31, + w=15, + output=None, + complexity_threshold=None, +) ``` +All `fetch()` arguments are keyword-only. + ## `Index.info()` -`Index.info(index_path)` Returns a `dict` of index metadata: +`index.info()` returns a `dict` of index metadata: | Key | Meaning | | --- | --- | | `k` | *k*-mer length | | `w` | minimizer window size | | `format` | `exact-u64`, `exact-u128`, or `bff` (binary fuse filter) | -| `count` | number of keys in index (fingerprint slot count for `bff`) | +| `count` | number of minimizers/keys represented by the index | ## `Index.filter()` -Filter a FASTA/FASTQ file or file pair against the index and return a `dict` of summary statistics. Auto-detects `.gz`/`.zst`/`.xz` compression on both input and output based on file extension. The Python GIL is released while filtering meaning calls benefit from multithreading. Refer to the [main Deacon readme](https://github.com/bede/deacon) for more detailed usage examples. +Filter FASTA, FASTQ, or CBQ input against the index and return a `dict` of summary statistics. FASTA/FASTQ compression (`.gz`, `.zst`, or `.xz`) is detected automatically. A `.cbq` output writes CBQ; `.cba` writes quality-free CBQ. The Python GIL is released while filtering, so calls benefit from multithreading. Refer to the [main Deacon readme](https://github.com/bede/deacon) for more detailed usage examples. ```python def filter( - fastq, # input path (FASTA/FASTQ, optionally .gz/.zst/.xz) - fastq2=None, # second mate for paired reads - *, # remaining options are keyword-only - interleaved=False, # treat fastq as interleaved paired reads (cannot combine with fastq2) + input, # positional-only FASTA, FASTQ, or CBQ path + /, + *, # every remaining argument is keyword-only + input2=None, # second FASTA/FASTQ mate + interleaved=False, # treat input as interleaved pairs (cannot combine with input2) + check_pairs=False, # validate paired read names deplete=False, # False = search (keep matches); True = deplete (remove matches) rename=False, # replace read names with sequential integers - output=None, # output path; None writes to stdout + output=None, # output path; None writes FASTA/FASTQ to stdout output2=None, # second output path for paired reads + summary=None, # optional JSON summary output path abs_threshold=2, # min absolute minimizer hits to call a match rel_threshold=0.01, # min proportion of minimizers hitting to call a match prefix_length=0, # only use the first N bp of each read (0 = whole read) - output_fasta=False, # emit FASTA instead of FASTQ + discard_quality=False, # discard quality scores ordered=False, # preserve input record ordering (deterministic, slightly slower) threads=8, # worker threads for filtering compression_level=2, # output compression level compression_threads=0, # threads for output compression (0 = auto) + cbq_block_size=16, # CBQ output block size in MiB (1-1024) quiet=True, # suppress progress/log output on stderr debug=False, # verbose per-read debug output ) -> dict @@ -81,9 +95,11 @@ def filter( **Modes.** With `deplete=False` (the default, *search* mode) reads that match the index are kept; with `deplete=True` reads that match are removed (host depletion). A read is a match only when it clears **both** thresholds: at least `abs_threshold` minimizer hits **and** at least `rel_threshold` of its minimizers hitting the index. -**Output.** When `output` is `None` the filtered records are written to stdout. To count without keeping the filtered sequences, pass `output="/dev/null"`. Statistics are returned regardless. +**Paired and CBQ input.** Use `input2=` for separate FASTA/FASTQ mates or `interleaved=True` for an interleaved stream. `check_pairs=True` validates Illumina CASAVA or `/1` and `/2` names. CBQ stores pairing internally, so CBQ input cannot be combined with `input2` or `interleaved`; paired CBQ output uses one `output` and cannot use `output2`. + +**Output.** When `output` is `None` the filtered records are written to stdout. To count without keeping the filtered sequences, pass `output="/dev/null"`. Pass `summary=` to write the same statistics returned by the call as JSON. For CBQ output, `cbq_block_size` is a lower bound: a larger input CBQ block size is preserved. -**Return value.** A `dict` including the run configuration (`version`, `index`, `input`/`input2`, `output`/`output2`, `k`, `w`, `abs_threshold`, `rel_threshold`, `prefix_length`, `deplete`, `rename`, `ordered`) and the results: +**Return value.** A `dict` including the run configuration (`version`, `index`, `input`/`input2`, `output`/`output2`, `k`, `w`, `abs_threshold`, `rel_threshold`, `prefix_length`, `deplete`, `rename`, `ordered`, `check_pairs`) and the results: | Key | Meaning | | --- | --- | diff --git a/deacon-py/python/deacon/_deacon.pyi b/deacon-py/python/deacon/_deacon.pyi index 59db6f9..0870901 100644 --- a/deacon-py/python/deacon/_deacon.pyi +++ b/deacon-py/python/deacon/_deacon.pyi @@ -1,18 +1,19 @@ from _typeshed import Incomplete -from typing import Any, final +from os import PathLike +from typing import final @final class Index: """ A loaded minimizer index, reusable across many `filter` calls. """ - def __new__(cls, /, path: str, complexity_threshold: float |None = None) -> Index: ... + def __new__(cls, path: str |PathLike[str], /, *, complexity_threshold: float |None = None) -> Index: ... @staticmethod - def fetch(name: str = "panhuman-1", k: int = 31, w: int = 15, output: str |None = None, complexity_threshold: float |None = None) -> Index: + def fetch(*, name: str = "panhuman-1", k: int = 31, w: int = 15, output: str |PathLike[str] |None = None, complexity_threshold: float |None = None) -> Index: """ Download a prebuilt index, then load and return it. """ - def filter(self, /, fastq: str, fastq2: str |None = None, *, interleaved: bool = False, deplete: bool = False, rename: bool = False, output: str |None = None, output2: str |None = None, abs_threshold: int = 2, rel_threshold: float = 0.01, prefix_length: int = 0, output_fasta: bool = False, ordered: bool = False, threads: int = 8, compression_level: int = 2, compression_threads: int = 0, debug: bool = False, quiet: bool = True) -> Any: ... + def filter(self, input: str |PathLike[str], /, *, input2: str |PathLike[str] |None = None, interleaved: bool = False, check_pairs: bool = False, deplete: bool = False, rename: bool = False, output: str |PathLike[str] |None = None, output2: str |PathLike[str] |None = None, summary: str |PathLike[str] |None = None, abs_threshold: int = 2, rel_threshold: float = 0.01, prefix_length: int = 0, discard_quality: bool = False, ordered: bool = False, threads: int = 8, compression_level: int = 2, compression_threads: int = 0, cbq_block_size: int = 16, quiet: bool = True, debug: bool = False) -> dict: ... def info(self, /) -> dict: """ Index metadata: k, w, format and minimizer/key count. diff --git a/deacon-py/src/lib.rs b/deacon-py/src/lib.rs index a5d5043..e111d9f 100644 --- a/deacon-py/src/lib.rs +++ b/deacon-py/src/lib.rs @@ -5,17 +5,26 @@ use std::path::{Path, PathBuf}; use std::sync::Arc; use ::deacon::{ - ComplexityAlgorithm, FilterRunConfig, IndexHeader, MinimizerSet, index_fetch, - load_index_from_path_auto, run_with_index, + ComplexityAlgorithm, DEFAULT_CBQ_BLOCK_SIZE_MIB, FilterRunConfig, IndexHeader, MinimizerSet, + index_fetch, load_index_from_path_auto, run_with_index, }; use pyo3::exceptions::{PyRuntimeError, PyValueError}; use pyo3::prelude::*; use pyo3::types::PyDict; +// Keep the literal in the Python signature/stub while detecting a changed core default at compile time. +const _: () = assert!(DEFAULT_CBQ_BLOCK_SIZE_MIB == 16); + fn to_pyerr(e: anyhow::Error) -> PyErr { PyRuntimeError::new_err(e.to_string()) } +fn path_to_string(path: PathBuf, argument: &str) -> PyResult { + path.into_os_string().into_string().map_err(|_| { + PyValueError::new_err(format!("{argument} must be representable as valid UTF-8")) + }) +} + /// A loaded minimizer index, reusable across many `filter` calls. #[pyclass(frozen)] struct Index { @@ -25,13 +34,9 @@ struct Index { minimizers: Arc, } -#[pymethods] impl Index { - #[new] - #[pyo3(signature = (path, complexity_threshold=None))] - fn new(path: &str, complexity_threshold: Option) -> PyResult { - let (mut minimizers, header) = - load_index_from_path_auto(Path::new(path)).map_err(to_pyerr)?; + fn load(path: &Path, complexity_threshold: Option) -> PyResult { + let (mut minimizers, header) = load_index_from_path_auto(path).map_err(to_pyerr)?; // Discard low-complexity index minimizers once at load (kdust); reused across filters. if let Some(threshold) = complexity_threshold { if matches!(minimizers, MinimizerSet::Fuse(_)) { @@ -49,28 +54,35 @@ impl Index { .map_err(to_pyerr)?; } Ok(Index { - label: path.to_string(), + label: path.to_string_lossy().into_owned(), k: header.kmer_length(), w: header.window_size(), minimizers: Arc::new(minimizers), }) } +} + +#[pymethods] +impl Index { + #[new] + #[pyo3(signature = (path, /, *, complexity_threshold=None))] + fn new(path: PathBuf, complexity_threshold: Option) -> PyResult { + Self::load(&path, complexity_threshold) + } /// Download a prebuilt index, then load and return it. #[staticmethod] - #[pyo3(signature = (name="panhuman-1", k=31, w=15, output=None, complexity_threshold=None))] + #[pyo3(signature = (*, name="panhuman-1", k=31, w=15, output=None, complexity_threshold=None))] fn fetch( name: &str, k: u8, w: u8, - output: Option, + output: Option, complexity_threshold: Option, ) -> PyResult { - let out_path = output - .map(PathBuf::from) - .unwrap_or_else(|| PathBuf::from(format!("{name}.k{k}w{w}.idx"))); + let out_path = output.unwrap_or_else(|| PathBuf::from(format!("{name}.k{k}w{w}.idx"))); index_fetch(name, k, w, Some(&out_path)).map_err(to_pyerr)?; - Index::new(&out_path.to_string_lossy(), complexity_threshold) + Index::load(&out_path, complexity_threshold) } /// Index metadata: k, w, format and minimizer/key count. @@ -89,68 +101,93 @@ impl Index { } #[pyo3(signature = ( - fastq, - fastq2=None, + input, + /, *, + input2=None, interleaved=false, + check_pairs=false, deplete=false, rename=false, output=None, output2=None, + summary=None, abs_threshold=2, rel_threshold=0.01, prefix_length=0, - output_fasta=false, + discard_quality=false, ordered=false, threads=8, compression_level=2, compression_threads=0, + cbq_block_size=16, quiet=true, debug=false, ))] fn filter( &self, py: Python<'_>, - fastq: String, - fastq2: Option, + input: PathBuf, + input2: Option, interleaved: bool, + check_pairs: bool, deplete: bool, rename: bool, - output: Option, - output2: Option, + output: Option, + output2: Option, + summary: Option, abs_threshold: usize, rel_threshold: f64, prefix_length: usize, - output_fasta: bool, + discard_quality: bool, ordered: bool, threads: u16, compression_level: u8, compression_threads: u16, + cbq_block_size: u16, quiet: bool, debug: bool, - ) -> PyResult> { - if interleaved && fastq2.is_some() { + ) -> PyResult> { + if interleaved && input2.is_some() { + return Err(PyValueError::new_err( + "interleaved cannot be combined with input2 (interleaved input is a single file/stream)", + )); + } + if abs_threshold == 0 { + return Err(PyValueError::new_err("abs_threshold must be at least 1")); + } + if !(1..=1024).contains(&cbq_block_size) { return Err(PyValueError::new_err( - "interleaved cannot be combined with fastq2 (interleaved input is a single file/stream)", + "cbq_block_size must be between 1 and 1024 MiB inclusive", )); } + + let input = path_to_string(input, "input")?; + let input2 = input2 + .map(|path| path_to_string(path, "input2")) + .transpose()?; + let output2 = output2 + .map(|path| path_to_string(path, "output2")) + .transpose()?; + let cfg = FilterRunConfig { - input_path: fastq, - input2_path: fastq2, + input_path: input, + input2_path: input2, interleaved, - check_pairs: false, - output_path: output.map(PathBuf::from), + check_pairs, + output_path: output, output2_path: output2, abs_threshold, rel_threshold, prefix_length, - summary_path: None, + summary_path: summary, deplete, rename, - output_fasta, + discard_quality, ordered, threads, compression_level, + cbq_block_size, compression_threads, debug, quiet, @@ -160,9 +197,11 @@ impl Index { let mins = Arc::clone(&self.minimizers); let (k, w) = (self.k, self.w); let summary = py - .detach(|| run_with_index(&mins, &IndexHeader::new(k, w), &cfg)) + .detach(|| run_with_index(mins, &IndexHeader::new(k, w), &cfg)) .map_err(to_pyerr)?; - Ok(pythonize::pythonize(py, &summary)?.unbind()) + Ok(pythonize::pythonize(py, &summary)? + .cast_into::()? + .unbind()) } } diff --git a/deacon-py/tests/test_api.py b/deacon-py/tests/test_api.py new file mode 100644 index 0000000..20e9675 --- /dev/null +++ b/deacon-py/tests/test_api.py @@ -0,0 +1,169 @@ +import base64 +import inspect +import json +import tempfile +import unittest +from pathlib import Path + +from deacon import Index + + +ROOT = Path(__file__).resolve().parents[2] +READS = ROOT / "tests/data/test_small_1.fastq.gz" +# Exact index built from READS with the default k=31, w=15 parameters. +INDEX_BYTES = base64.b64decode( + "Ax8PQQAAAAAAAACve33W9VXXBESYNVc0k/ANrPWbvtZvcxLn2gd9//W+G89JPdM1xNcVtW1Qym/XCivXZ11fdU0UFa1UtVGcmykSrVPjVL3GpxEJZtxZ51QZB0QNr3t91vUVfff33d+Xax+t1xAXVm6mCL0XtVHXF1ED1V1Xd13ddS2x4wlm3FnnFFB1vRe1UdcXbdbLWX/+vRdUbqZIhFlzBUXkvsPlUP8G1fVe1EZdXwR91TVRVG6mCH0RNbzu9VkXlPLbtcLaMw1d70Vt1PVFFNT3r+8bfb8CpNA3lfA6TQ/1VddEUbmZIm3WkELfVMIrReVmikSYNRe9bdaQQt9UAjmjHEaev7UttVLVRnFupgg94y7CTJzvNBSVmykSYdYcvdP3Xd8zfR/QtJiVr9nKGm3U9UXU8LoXkedvbRuU8htRw+ten3V9FdS/YW/QtJgVKb9dK6w90xR0DfF1ddfVHd3fd39frn0bLWbla7ay1i9Z11ddE0V9FY997/R91/cMXddGXWN1rRTx+27ff33PNb1B02JWvmYrfde1UddYXSs1VtdKVRvFOdC9bdaQQt8U1bVTFFRd7wVEUbmZIhFmDdU1UdRXXes1FVGjUF07RQHXeg1xYeVmCn3X90zfd38fTfVQMNVDQjPVXVd3XZt1PT3T1K1F5L4Dd13ddW3W9R3dWkTuO1wONVHfv75v9P0K" +) +TEST_DIRECTORY = None +INDEX = None + + +def setUpModule(): + global TEST_DIRECTORY, INDEX + TEST_DIRECTORY = tempfile.TemporaryDirectory() + INDEX = Path(TEST_DIRECTORY.name) / "test.idx" + INDEX.write_bytes(INDEX_BYTES) + + +def tearDownModule(): + TEST_DIRECTORY.cleanup() + + +class SignatureTests(unittest.TestCase): + @classmethod + def setUpClass(cls): + cls.index = Index(INDEX) + + def test_signatures_are_future_proof(self): + self.assertEqual( + str(inspect.signature(Index)), + "(path, /, *, complexity_threshold=None)", + ) + self.assertEqual( + str(inspect.signature(Index.fetch)), + "(*, name='panhuman-1', k=31, w=15, output=None, complexity_threshold=None)", + ) + + parameters = inspect.signature(Index.filter).parameters + self.assertEqual(parameters["input"].kind, inspect.Parameter.POSITIONAL_ONLY) + self.assertEqual(parameters["cbq_block_size"].default, 16) + for name, parameter in parameters.items(): + if name not in {"self", "input"}: + self.assertEqual(parameter.kind, inspect.Parameter.KEYWORD_ONLY, name) + + def test_optional_positional_arguments_are_rejected(self): + with self.assertRaises(TypeError): + Index(INDEX, None) + with self.assertRaises(TypeError): + Index.fetch("panhuman-1") + with self.assertRaises(TypeError): + self.index.filter(READS, None) + + def test_primary_paths_are_positional_only(self): + with self.assertRaises(TypeError): + Index(path=INDEX) + with self.assertRaises(TypeError): + self.index.filter(input=READS) + + +class FilteringTests(unittest.TestCase): + @classmethod + def setUpClass(cls): + cls.index = Index(INDEX) + + def test_pathlike_output_and_json_summary(self): + with tempfile.TemporaryDirectory() as directory: + directory = Path(directory) + output = directory / "output.fastq" + summary_path = directory / "summary.json" + + summary = self.index.filter( + READS, + output=output, + summary=summary_path, + threads=1, + ) + + self.assertTrue(output.is_file()) + self.assertEqual(json.loads(summary_path.read_text()), summary) + self.assertEqual(summary["check_pairs"], False) + self.assertGreater(summary["seqs_in"], 0) + + def test_check_pairs(self): + with tempfile.TemporaryDirectory() as directory: + directory = Path(directory) + mate1 = directory / "mate1.fastq" + mate2 = directory / "mate2.fastq" + mate1.write_text("@read/1\nACGT\n+\nIIII\n") + mate2.write_text("@read/2\nACGT\n+\nIIII\n") + summary = self.index.filter( + mate1, + input2=mate2, + output=directory / "output1.fastq", + output2=directory / "output2.fastq", + check_pairs=True, + threads=1, + ) + self.assertTrue(summary["check_pairs"]) + + mate1.write_text("@read-a/1\nACGT\n+\nIIII\n") + mate2.write_text("@read-b/2\nACGT\n+\nIIII\n") + with self.assertRaisesRegex(RuntimeError, "Paired record name mismatch"): + self.index.filter( + mate1, + input2=mate2, + output=directory / "bad1.fastq", + output2=directory / "bad2.fastq", + check_pairs=True, + threads=1, + ) + + def test_cbq_options_and_round_trip(self): + with tempfile.TemporaryDirectory() as directory: + directory = Path(directory) + cbq = directory / "output.cbq" + first = self.index.filter( + READS, + output=cbq, + cbq_block_size=1, + threads=1, + ) + second = self.index.filter( + cbq, + output=directory / "roundtrip.fastq", + threads=1, + ) + self.assertEqual(second["seqs_in"], first["seqs_out"]) + + cba = directory / "output.cba" + self.index.filter(READS, output=cba, cbq_block_size=1, threads=1) + fasta = directory / "quality-free.fasta" + self.index.filter(cba, output=fasta, threads=1) + if fasta.stat().st_size: + self.assertEqual(fasta.read_bytes()[:1], b">") + + with self.assertRaisesRegex(RuntimeError, "CBQ input does not support INPUT2"): + self.index.filter(cbq, input2=READS, threads=1) + with self.assertRaisesRegex(RuntimeError, "CBQ output does not support OUTPUT2"): + self.index.filter( + READS, + output=directory / "invalid.cbq", + output2=directory / "invalid.fastq", + threads=1, + ) + + def test_argument_validation(self): + with self.assertRaisesRegex(ValueError, "abs_threshold"): + self.index.filter(READS, abs_threshold=0) + with self.assertRaisesRegex(ValueError, "cbq_block_size"): + self.index.filter(READS, cbq_block_size=0) + with self.assertRaisesRegex(ValueError, "cbq_block_size"): + self.index.filter(READS, cbq_block_size=1025) + with self.assertRaisesRegex(ValueError, "input2"): + self.index.filter(READS, input2=READS, interleaved=True) + with self.assertRaisesRegex(RuntimeError, "check-pairs requires paired input"): + self.index.filter(READS, check_pairs=True) + + +if __name__ == "__main__": + unittest.main() diff --git a/deacon-wasm/src/lib.rs b/deacon-wasm/src/lib.rs index 9feffd5..31ca209 100644 --- a/deacon-wasm/src/lib.rs +++ b/deacon-wasm/src/lib.rs @@ -70,7 +70,7 @@ pub struct FilterSession { index: Arc, kernel: FilterKernel, rename: bool, - output_fasta: bool, + discard_quality: bool, rename_counter: u64, parser: SeqChunkParser, stats: FilterStats, @@ -91,7 +91,7 @@ impl FilterSession { decompress_input: bool, compress_output: bool, rename: bool, - output_fasta: bool, + discard_quality: bool, ) -> Result { let k = index.inner.header.kmer_length(); let w = index.inner.header.window_size(); @@ -121,7 +121,7 @@ impl FilterSession { ) .map_err(|e| JsValue::from_str(&e.to_string()))?, rename, - output_fasta, + discard_quality, rename_counter: 0, parser: SeqChunkParser::new(), stats: FilterStats::default(), @@ -268,7 +268,7 @@ impl FilterSession { header, seq, qual, - self.output_fasta, + self.discard_quality, self.rename, counter, b"", @@ -291,7 +291,7 @@ pub struct PairedFilterSession { index: Arc, kernel: FilterKernel, rename: bool, - output_fasta: bool, + discard_quality: bool, rename_counter: u64, parser_r1: SeqChunkParser, parser_r2: SeqChunkParser, @@ -323,7 +323,7 @@ impl PairedFilterSession { compress_r1: bool, compress_r2: bool, rename: bool, - output_fasta: bool, + discard_quality: bool, ) -> Result { let k = index.inner.header.kmer_length(); let w = index.inner.header.window_size(); @@ -341,7 +341,7 @@ impl PairedFilterSession { ) .map_err(|e| JsValue::from_str(&e.to_string()))?, rename, - output_fasta, + discard_quality, rename_counter: 0, parser_r1: SeqChunkParser::new(), parser_r2: SeqChunkParser::new(), @@ -469,7 +469,7 @@ impl PairedFilterSession { &record1.header, &record1.seq, record1.qual.as_deref(), - self.output_fasta, + self.discard_quality, self.rename, counter, b"/1", @@ -480,7 +480,7 @@ impl PairedFilterSession { &record2.header, &record2.seq, record2.qual.as_deref(), - self.output_fasta, + self.discard_quality, self.rename, counter, b"/2", @@ -1006,12 +1006,12 @@ fn write_record( header: &[u8], seq: &[u8], qual: Option<&[u8]>, - output_fasta: bool, + discard_quality: bool, rename: bool, counter: u64, read_suffix: &[u8], ) -> std::io::Result<()> { - let is_fasta = output_fasta || qual.is_none(); + let is_fasta = discard_quality || qual.is_none(); output.write_all(if is_fasta { b">" } else { b"@" })?; if rename { output.write_all(counter.to_string().as_bytes())?; diff --git a/deacon-wasm/worker.js b/deacon-wasm/worker.js index 17143fb..975695d 100644 --- a/deacon-wasm/worker.js +++ b/deacon-wasm/worker.js @@ -35,7 +35,7 @@ async function streamFilterFile(file, opts) { isGz, // decompress_input isGz, // compress_output (match input format) false, // rename - false // output_fasta + false // discard_quality ); const reader = file.stream().getReader(); @@ -133,7 +133,7 @@ async function streamFilterPairedFiles(file1, file2, opts) { r1Gz, r2Gz, // decompress_r1, decompress_r2 r1Gz, r2Gz, // compress_r1, compress_r2 false, // rename - false // output_fasta + false // discard_quality ); const readerR1 = file1.stream().getReader(); diff --git a/scripts/check-wasm-paired-parity.sh b/scripts/check-wasm-paired-parity.sh index f5da565..0db97ba 100755 --- a/scripts/check-wasm-paired-parity.sh +++ b/scripts/check-wasm-paired-parity.sh @@ -53,15 +53,15 @@ echo "Native/WASM paired parity check passed" echo " R1 bytes: $(wc -c < "$native_r1" | tr -d ' ')" echo " R2 bytes: $(wc -c < "$native_r2" | tr -d ' ')" -# Second pass: --rename + --fasta (-R -f) must also match byte-for-byte. +# Second pass: --rename + --discard-quality must also match byte-for-byte. native_rf_r1="$work_dir/native.rf.r1.fasta" native_rf_r2="$work_dir/native.rf.r2.fasta" wasm_rf_r1="$work_dir/wasm.rf.r1.fasta" wasm_rf_r2="$work_dir/wasm.rf.r2.fasta" -echo "Running native paired filter (-R -f)..." -cargo run --release -- filter -d -t 1 -R -f "$index_path" "$reads1_path" "$reads2_path" \ +echo "Running native paired filter (--rename --discard-quality)..." +cargo run --release -- filter -d -t 1 -R --discard-quality "$index_path" "$reads1_path" "$reads2_path" \ -o "$native_rf_r1" -O "$native_rf_r2" -echo "Running WASM paired filter (--rename --fasta)..." +echo "Running WASM paired filter (--rename --discard-quality)..." node "$repo_root/scripts/wasm_paired_filter_to_file.mjs" \ --pkg "$pkg_dir" \ --index "$index_path" \ @@ -69,14 +69,14 @@ node "$repo_root/scripts/wasm_paired_filter_to_file.mjs" \ --reads2 "$reads2_path" \ --output1 "$wasm_rf_r1" \ --output2 "$wasm_rf_r2" \ - --deplete --rename --fasta + --deplete --rename --discard-quality if ! cmp -s "$native_rf_r1" "$wasm_rf_r1" || ! cmp -s "$native_rf_r2" "$wasm_rf_r2"; then - echo "Native/WASM paired output mismatch (rename+fasta)" >&2 + echo "Native/WASM paired output mismatch (rename+discard-quality)" >&2 echo " native R1: $(shasum -a 256 "$native_rf_r1" | awk '{print $1}')" >&2 echo " wasm R1: $(shasum -a 256 "$wasm_rf_r1" | awk '{print $1}')" >&2 echo " native R2: $(shasum -a 256 "$native_rf_r2" | awk '{print $1}')" >&2 echo " wasm R2: $(shasum -a 256 "$wasm_rf_r2" | awk '{print $1}')" >&2 exit 1 fi -echo "Native/WASM paired rename+fasta parity check passed" +echo "Native/WASM paired rename+discard-quality parity check passed" diff --git a/scripts/check-wasm-parity.sh b/scripts/check-wasm-parity.sh index 7969977..dced7ec 100755 --- a/scripts/check-wasm-parity.sh +++ b/scripts/check-wasm-parity.sh @@ -48,23 +48,23 @@ echo "Native/WASM parity check passed" echo " bytes: $size" echo " sha256: $sha" -# Second pass: --rename + --fasta (-R -f) must also match byte-for-byte. -native_rf="$work_dir/native.rename-fasta.fasta" -wasm_rf="$work_dir/wasm.rename-fasta.fasta" -echo "Running native filter (-R -f)..." -cargo run --release -- filter -t 1 -R -f "$index_path" "$reads_path" -o "$native_rf" -echo "Running WASM filter (--rename --fasta)..." +# Second pass: --rename + --discard-quality must also match byte-for-byte. +native_rf="$work_dir/native.rename-discard-quality.fasta" +wasm_rf="$work_dir/wasm.rename-discard-quality.fasta" +echo "Running native filter (--rename --discard-quality)..." +cargo run --release -- filter -t 1 -R --discard-quality "$index_path" "$reads_path" -o "$native_rf" +echo "Running WASM filter (--rename --discard-quality)..." node "$repo_root/scripts/wasm_filter_to_file.mjs" \ --pkg "$pkg_dir" \ --index "$index_path" \ --reads "$reads_path" \ --output "$wasm_rf" \ - --rename --fasta + --rename --discard-quality if ! cmp -s "$native_rf" "$wasm_rf"; then - echo "Native/WASM output mismatch (rename+fasta)" >&2 + echo "Native/WASM output mismatch (rename+discard-quality)" >&2 echo " native: $(shasum -a 256 "$native_rf" | awk '{print $1}')" >&2 echo " wasm: $(shasum -a 256 "$wasm_rf" | awk '{print $1}')" >&2 exit 1 fi -echo "Native/WASM rename+fasta parity check passed ($(wc -c < "$native_rf" | tr -d ' ') bytes)" +echo "Native/WASM rename+discard-quality parity check passed ($(wc -c < "$native_rf" | tr -d ' ') bytes)" diff --git a/scripts/wasm_bench.mjs b/scripts/wasm_bench.mjs index 76e4c09..2382c50 100644 --- a/scripts/wasm_bench.mjs +++ b/scripts/wasm_bench.mjs @@ -64,7 +64,7 @@ for (let iter = 0; iter < iters; iter++) { decompressInput, false, // compress_output false, // rename - false, // output_fasta + false, // discard_quality ); const start = process.hrtime.bigint(); diff --git a/scripts/wasm_filter_to_file.mjs b/scripts/wasm_filter_to_file.mjs index 554e610..715af28 100644 --- a/scripts/wasm_filter_to_file.mjs +++ b/scripts/wasm_filter_to_file.mjs @@ -50,7 +50,7 @@ const indexPath = requireArg(args, "index"); const readsPath = requireArg(args, "reads"); const outputPath = requireArg(args, "output"); const rename = args.flags.has("rename"); -const outputFasta = args.flags.has("fasta"); +const discardQuality = args.flags.has("discard-quality"); const wasmModule = await import(pathToFileURL(path.join(pkgDir, "deacon_wasm.js")).href); @@ -69,7 +69,7 @@ const session = new wasmModule.FilterSession( true, // decompress_input false, // compress_output rename, - outputFasta, + discardQuality, ); const input = createReadStream(readsPath, { highWaterMark: 256 * 1024 }); diff --git a/scripts/wasm_paired_filter_to_file.mjs b/scripts/wasm_paired_filter_to_file.mjs index 9e012fb..8b2cce8 100644 --- a/scripts/wasm_paired_filter_to_file.mjs +++ b/scripts/wasm_paired_filter_to_file.mjs @@ -49,7 +49,7 @@ const output1Path = requireArg(args, "output1"); const output2Path = requireArg(args, "output2"); const deplete = args.flags.has("deplete"); const rename = args.flags.has("rename"); -const outputFasta = args.flags.has("fasta"); +const discardQuality = args.flags.has("discard-quality"); const absThreshold = Number(args.get("abs-threshold") ?? 2); const relThreshold = Number(args.get("rel-threshold") ?? 0.01); @@ -72,7 +72,7 @@ const session = new wasmModule.PairedFilterSession( false, // compress_r1 false, // compress_r2 rename, - outputFasta, + discardQuality, ); const input1 = createReadStream(reads1Path, { highWaterMark: 256 * 1024 }); diff --git a/scripts/wasm_smoke_test.mjs b/scripts/wasm_smoke_test.mjs index 13e9de6..101589b 100644 --- a/scripts/wasm_smoke_test.mjs +++ b/scripts/wasm_smoke_test.mjs @@ -55,7 +55,7 @@ const session = new wasmModule.FilterSession( false, // decompress_input false, // compress_output false, // rename - false, // output_fasta + false, // discard_quality ); const outputChunks = []; diff --git a/src/filter.rs b/src/filter.rs index b2c38b3..d4af53d 100644 --- a/src/filter.rs +++ b/src/filter.rs @@ -4,14 +4,18 @@ use crate::{ MinimizerSet, validate_unit_interval, }; use anyhow::{Context, Result}; +use binseq::cbq; +use binseq::write::{BinseqWriterBuilder, Format as BinseqFormat}; +use binseq::{BinseqRecord, ParallelReader as BinseqParallelReader, SequencingRecordBuilder}; use indicatif::{ProgressBar, ProgressDrawTarget, ProgressStyle}; use paraseq::Record; use paraseq::fastx::Reader; use paraseq::parallel::{PairedParallelProcessor, ParallelProcessor, ParallelReader}; use parking_lot::Mutex; use serde::{Deserialize, Serialize}; +use std::borrow::Cow; use std::fs::{File, OpenOptions}; -use std::io::{self, BufWriter, Write}; +use std::io::{self, BufWriter, Read, Write}; use std::path::PathBuf; use std::sync::Arc; use std::sync::atomic::{AtomicU64, Ordering}; @@ -19,8 +23,337 @@ use std::time::Instant; const OUTPUT_BUFFER_SIZE: usize = 8 * 1024 * 1024; // Opt: 8MB output buffer const DEFAULT_BUFFER_SIZE: usize = 64 * 1024; +/// Default CBQ block size in MiB; a block must hold the largest single record +pub const DEFAULT_CBQ_BLOCK_SIZE_MIB: u16 = 16; type BoxedWriter = Box; +/// CBQ output is always a named file, no stdout +type CbqWriter = binseq::BinseqWriter; + +/// Sequence file format +#[derive(Clone, Copy, PartialEq, Eq)] +enum Format { + Fastx, + Cbq, +} + +/// Input metadata, resolved before any output is opened +struct InputLayout { + format: Format, + paired: bool, + qualities: bool, + headers: bool, + /// CBQ per-record flags present (always false for FASTX) + flags: bool, + /// CBQ input block size, a lower bound for CBQ output (None for FASTX) + block_size: Option, +} + +#[allow(clippy::large_enum_variant)] +enum Input { + Cbq(cbq::MmapReader), + FastxSingle(Reader>), + FastxInterleaved(Reader>), + FastxPaired( + Reader>, + Reader>, + ), + Empty, +} + +/// Borrowed view of one read, shared by every reader and writer +struct ReadView<'a> { + id: &'a [u8], + seq: &'a [u8], + qual: Option<&'a [u8]>, + /// CBQ per-record flag, preserved on CBQ output (None elsewhere) + flag: Option, +} + +/// Format-specific output state: thread-local buffers plus the shared +/// writers each processor clone merges into +#[allow(clippy::large_enum_variant)] // mirrors binseq's BinseqWriter, which is not boxed upstream +#[derive(Clone)] +enum Output { + Fastx { + local: Vec, + local2: Vec, + pending: Vec, + pending2: Vec, + batch_kept: u64, + shared: Arc>, + shared2: Option>>, + }, + Cbq { + /// Thread-local block writer; completed blocks merge into the shared writer + local: binseq::BinseqWriter>, + shared: Arc>, + }, +} + +impl Output { + fn fastx(writer: BoxedWriter, writer2: Option) -> Self { + Output::Fastx { + local: Vec::with_capacity(DEFAULT_BUFFER_SIZE), + local2: Vec::with_capacity(DEFAULT_BUFFER_SIZE), + pending: Vec::new(), + pending2: Vec::new(), + batch_kept: 0, + shared: Arc::new(Mutex::new(writer)), + shared2: writer2.map(|w| Arc::new(Mutex::new(w))), + } + } + + fn cbq(writer: CbqWriter) -> Result { + let shared = Arc::new(Mutex::new(writer)); + let local = shared.lock().new_headless_buffer()?; + Ok(Output::Cbq { local, shared }) + } + + /// Format one read into a FASTX local buffer, deferring the header to + /// [`write_renamed`] when renaming + fn push_fastx( + buffer: &mut Vec, + pending: &mut Vec, + batch_kept: u64, + read: &ReadView, + suffix: &'static [u8], + rename: bool, + discard_quality: bool, + ) -> Result<()> { + let offset = buffer.len(); + let marker = format_record_to_buffer(read, rename, discard_quality, buffer)?; + if rename { + pending.push(PendingRename { + offset, + ordinal: batch_kept, + marker, + suffix, + }); + } + Ok(()) + } + + /// Write one unpaired read + fn push_read( + &mut self, + read: &ReadView, + rename: bool, + discard_quality: bool, + rename_counter: &AtomicU64, + ) -> Result<()> { + match self { + Output::Fastx { + local, + pending, + batch_kept, + .. + } => { + Self::push_fastx( + local, + pending, + *batch_kept, + read, + b"", + rename, + discard_quality, + )?; + if rename { + *batch_kept += 1; + } + } + Output::Cbq { local, .. } => { + let number = rename.then(|| rename_counter.fetch_add(1, Ordering::Relaxed) + 1); + let header = cbq_header(read.id, number, b""); + let seq_record = SequencingRecordBuilder::default() + .s_seq(read.seq) + .s_header(&header) + .opt_s_qual(read.qual) + .opt_flag(read.flag) + .build()?; + local.push(seq_record)?; + } + } + Ok(()) + } + + /// Write a kept pair: one native CBQ record, two FASTX files, or + /// interleaved FASTX. Both mates take the pair's number when renaming. + fn push_pair( + &mut self, + read1: &ReadView, + read2: &ReadView, + rename: bool, + discard_quality: bool, + rename_counter: &AtomicU64, + ) -> Result<()> { + match self { + Output::Fastx { + local, + local2, + pending, + pending2, + batch_kept, + shared2, + .. + } => { + Self::push_fastx( + local, + pending, + *batch_kept, + read1, + b"/1", + rename, + discard_quality, + )?; + if shared2.is_some() { + // Separate outputs + Self::push_fastx( + local2, + pending2, + *batch_kept, + read2, + b"/2", + rename, + discard_quality, + )?; + } else { + // Interleaved output + Self::push_fastx( + local, + pending, + *batch_kept, + read2, + b"/2", + rename, + discard_quality, + )?; + } + if rename { + *batch_kept += 1; + } + } + Output::Cbq { local, .. } => { + let number = rename.then(|| rename_counter.fetch_add(1, Ordering::Relaxed) + 1); + let header1 = cbq_header(read1.id, number, b"/1"); + let header2 = cbq_header(read2.id, number, b"/2"); + let seq_record = SequencingRecordBuilder::default() + .s_seq(read1.seq) + .s_header(&header1) + .opt_s_qual(read1.qual) + .x_seq(read2.seq) + .x_header(&header2) + .opt_x_qual(read2.qual) + .opt_flag(read1.flag) + .build()?; + local.push(seq_record)?; + } + } + Ok(()) + } + + /// Merge thread-local buffers into the shared writer(s) (per batch). Rename + /// numbers are claimed under the writer lock: monotone in output order, but + /// deterministic only with `--ordered` or `-t 1`. + fn flush_batch( + &mut self, + rename: bool, + rename_counter: &AtomicU64, + ordered: bool, + ) -> Result<()> { + match self { + Output::Cbq { local, shared } => { + let mut shared = shared.lock(); + if ordered { + // Keep the shared partial block ahead of the next batch. + if let binseq::BinseqWriter::Cbq(writer) = &mut *shared { + writer.flush()?; + } + shared.ingest(local)?; + } else { + shared.ingest_completed(local)?; + } + } + Output::Fastx { + local, + local2, + pending, + pending2, + batch_kept, + shared, + shared2, + } => { + if let Some(shared2) = shared2 { + // Atomic paired batch writing + if !local.is_empty() || !local2.is_empty() { + let mut writer1 = shared.lock(); + let mut writer2 = shared2.lock(); + + if rename { + // Both mates share one block of numbers + let base = rename_counter.fetch_add(*batch_kept, Ordering::Relaxed) + 1; + write_renamed(&mut writer1, local, pending, base)?; + write_renamed(&mut writer2, local2, pending2, base)?; + } else { + writer1.write_all(local)?; + writer2.write_all(local2)?; + } + writer1.flush()?; + writer2.flush()?; + } + } else if !local.is_empty() { + let mut writer = shared.lock(); + if rename { + let base = rename_counter.fetch_add(*batch_kept, Ordering::Relaxed) + 1; + write_renamed(&mut writer, local, pending, base)?; + } else { + writer.write_all(local)?; + } + writer.flush()?; + } + local.clear(); + local2.clear(); + pending.clear(); + pending2.clear(); + *batch_kept = 0; + } + } + Ok(()) + } + + /// Merge any remaining local CBQ blocks into the shared writer (per thread) + fn flush_thread(&mut self) -> Result<()> { + if let Output::Cbq { local, shared } = self { + shared.lock().ingest(local)?; + } + Ok(()) + } + + /// Finish a CBQ stream: flush remaining blocks and write the embedded index + fn finish(&self) -> Result<()> { + if let Output::Cbq { shared, .. } = self { + shared.lock().finish()?; + } + Ok(()) + } +} + +/// Header for a CBQ record: the original id, or `number[ suffix]` when +/// renaming. Numbers are claimed per record: unique, but sequential only +/// with `--ordered` or `-t 1`. +fn cbq_header<'a>(id: &'a [u8], number: Option, suffix: &[u8]) -> Cow<'a, [u8]> { + match number { + Some(n) => { + let mut header = n.to_string().into_bytes(); + if !suffix.is_empty() { + header.push(b' '); + header.extend_from_slice(suffix); + } + Cow::Owned(header) + } + None => Cow::Borrowed(id), + } +} /// Filtering config for an already-loaded index (no index path; see [`FilterConfig`]). pub struct FilterRunConfig { @@ -48,14 +381,16 @@ pub struct FilterRunConfig { pub deplete: bool, /// Replace sequence headers with incrementing numbers pub rename: bool, - /// Force FASTA output (discards quality scores) - pub output_fasta: bool, + /// Emit fasta or quality-free cbq regardless of input format (implied by a .cba output) + pub discard_quality: bool, /// Preserve input record ordering (deterministic, slightly slower) pub ordered: bool, /// Number of execution threads (0 = auto) pub threads: u16, /// Compression level for output files (1-22 for zst, 1-9 for gz) pub compression_level: u8, + /// cbq output block size in MiB (raised to the cbq input's block size if larger) + pub cbq_block_size: u16, /// Number of threads for compression (0 = auto) pub compression_threads: u16, /// Debug mode: output sequences with minimizer hits to stderr @@ -73,7 +408,7 @@ struct FilterProcessorConfig { prefix_length: usize, deplete: bool, rename: bool, - output_fasta: bool, + discard_quality: bool, debug: bool, check_pairs: bool, ordered: bool, @@ -196,6 +531,195 @@ fn create_paraseq_reader(path: Option<&str>) -> Result bool { + path.ends_with(".cbq") || is_cba_path(path) +} + +/// `.cba` names a quality-free CBQ, implying `--discard-quality` for that output +fn is_cba_path(path: &str) -> bool { + path.ends_with(".cba") +} + +/// Resolve the output format from the output path suffix +fn resolve_output_format(config: &FilterRunConfig) -> Format { + match config.output_path.as_deref() { + Some(path) if is_cbq_path(&path.to_string_lossy()) => Format::Cbq, + _ => Format::Fastx, + } +} + +/// Resolve the input format and layout, opening all readers up front so input +/// errors surface before any output file is created. +fn open_input(config: &FilterRunConfig, interleaved_input: bool) -> Result<(InputLayout, Input)> { + // CBQ input is file-only (mmap reader); stdin would need binseq's streaming reader + if config.input_path == "-" || is_special_input_path(&config.input_path) { + return open_fastx(config, interleaved_input); + } + + // Regular files: sniff the CBQ magic; everything else is FASTX + let file = File::open(&config.input_path) + .map_err(|e| anyhow::anyhow!("Failed to open file {}: {}", config.input_path, e))?; + let file_len = file.metadata()?.len(); + let mut magic = Vec::with_capacity(64); + file.take(64).read_to_end(&mut magic)?; + match BinseqFormat::sniff(&magic) { + Some(BinseqFormat::Cbq) => { + // The mmap reader slices header and footer unchecked, so guard the length + let min_len = (std::mem::size_of::() + + std::mem::size_of::()) as u64; + if file_len < min_len { + anyhow::bail!("Truncated or corrupt CBQ input: {}", config.input_path); + } + let reader = + cbq::MmapReader::new(&config.input_path).context("Failed to open CBQ input")?; + let header = reader.header(); + let layout = InputLayout { + format: Format::Cbq, + paired: header.is_paired(), + qualities: header.has_qualities(), + headers: header.has_headers(), + flags: header.has_flags(), + block_size: Some(header.block_size as usize), + }; + // binseq's parallel reader rejects an empty record range + let input = if reader.num_records() == 0 { + Input::Empty + } else { + Input::Cbq(reader) + }; + Ok((layout, input)) + } + Some(f) => anyhow::bail!("{f:?} input is not supported"), + None => open_fastx(config, interleaved_input), + } +} + +/// Open a FASTX input and resolve its layout +fn open_fastx(config: &FilterRunConfig, interleaved_input: bool) -> Result<(InputLayout, Input)> { + let layout = InputLayout { + format: Format::Fastx, + paired: interleaved_input || config.input2_path.is_some(), + qualities: false, + headers: true, + flags: false, + block_size: None, + }; + + let input1_empty = is_empty_file(&config.input_path)?; + let input2_empty = config + .input2_path + .as_deref() + .map(is_empty_file) + .transpose()? + .unwrap_or(false); + + if interleaved_input { + if input1_empty { + return Ok((layout, Input::Empty)); + } + return match create_paraseq_reader(Some(config.input_path.as_str())) { + Ok(reader) => { + let qualities = reader.format() == paraseq::fastx::Format::Fastq; + Ok(( + InputLayout { + qualities, + ..layout + }, + Input::FastxInterleaved(reader), + )) + } + Err(e) if is_empty_input_error(&e) => Ok((layout, Input::Empty)), + Err(e) => Err(e), + }; + } + + if let Some(input2_path) = config.input2_path.as_deref() { + if input1_empty && input2_empty { + return Ok((layout, Input::Empty)); + } + if input1_empty || input2_empty { + return Err(anyhow::anyhow!( + "One paired file is empty but the other is not" + )); + } + let r1 = create_paraseq_reader(Some(config.input_path.as_str())); + let r2 = create_paraseq_reader(Some(input2_path)); + return match (r1, r2) { + (Ok(reader1), Ok(reader2)) => { + let qualities = reader1.format() == paraseq::fastx::Format::Fastq; + Ok(( + InputLayout { + qualities, + ..layout + }, + Input::FastxPaired(reader1, reader2), + )) + } + (Err(e1), Err(e2)) if is_empty_input_error(&e1) && is_empty_input_error(&e2) => { + Ok((layout, Input::Empty)) + } + (Err(e), _) if is_empty_input_error(&e) => Err(anyhow::anyhow!( + "First paired file appears empty while second is not" + )), + (_, Err(e)) if is_empty_input_error(&e) => Err(anyhow::anyhow!( + "Second paired file appears empty while first is not" + )), + (Err(e), _) => Err(e), + (_, Err(e)) => Err(e), + }; + } + + if input1_empty { + return Ok((layout, Input::Empty)); + } + match create_paraseq_reader(Some(config.input_path.as_str())) { + Ok(reader) => { + let qualities = reader.format() == paraseq::fastx::Format::Fastq; + Ok(( + InputLayout { + qualities, + ..layout + }, + Input::FastxSingle(reader), + )) + } + Err(e) if is_empty_input_error(&e) => Ok((layout, Input::Empty)), + Err(e) => Err(e), + } +} + +/// Validate format combinations before any output file is opened or truncated +fn validate_input_output( + layout: &InputLayout, + output_format: Format, + config: &FilterRunConfig, +) -> Result<()> { + if layout.format == Format::Cbq && config.input2_path.is_some() { + anyhow::bail!("CBQ input does not support INPUT2"); + } + if layout.format == Format::Cbq && config.interleaved { + anyhow::bail!("CBQ input does not support --interleaved"); + } + if output_format == Format::Cbq && config.output2_path.is_some() { + anyhow::bail!("CBQ output does not support OUTPUT2; CBQ pairing is native"); + } + if config.output2_path.as_deref().is_some_and(is_cbq_path) { + anyhow::bail!("OUTPUT2 cannot be CBQ; CBQ pairing is native (use a single --output)"); + } + if output_format == Format::Cbq && !(1..=22).contains(&config.compression_level) { + anyhow::bail!( + "Invalid CBQ compression level {} (must be 1-22)", + config.compression_level + ); + } + if config.check_pairs && layout.format == Format::Cbq && !layout.headers { + anyhow::bail!("--check-pairs requires CBQ input with headers"); + } + validate_check_pairs_mode(config.check_pairs, layout.paired)?; + Ok(()) +} + /// A record buffered without its header, which is written later by [`write_renamed`] /// once its final number is known (numbering depends on earlier batches finishing). #[derive(Clone)] @@ -213,32 +737,31 @@ struct PendingRename { /// Format a record into a buffer (FASTA/FASTQ), returning its line prefix /// `seq` is the newline-stripped sequence from `record.seq()`. /// When renaming, the header is left to [`write_renamed`]. -fn format_record_to_buffer( - record: &R, - seq: &[u8], +fn format_record_to_buffer( + read: &ReadView, rename: bool, - output_fasta: bool, + discard_quality: bool, buffer: &mut Vec, ) -> Result { - let is_fasta = output_fasta || record.qual().is_none(); + let is_fasta = discard_quality || read.qual.is_none(); let marker = if is_fasta { b'>' } else { b'@' }; // Header (omitted when renaming) if !rename { buffer.push(marker); - buffer.extend_from_slice(record.id()); + buffer.extend_from_slice(read.id); buffer.write_all(b"\n")?; } // Sequence - buffer.extend_from_slice(seq); + buffer.extend_from_slice(read.seq); if is_fasta { buffer.write_all(b"\n")?; } else { // FASTQ: plus and qual lines buffer.write_all(b"\n+\n")?; - if let Some(qual) = record.qual() { + if let Some(qual) = read.qual { buffer.extend_from_slice(qual); } buffer.write_all(b"\n")?; @@ -305,7 +828,12 @@ fn count_compressed_outputs(config: &FilterRunConfig) -> u8 { count } -// Return a suitable writer for the output path extension +/// Open unbuffered so embedded-index write errors surface from `finish`. +fn open_cbq_output(path: &std::path::Path) -> Result { + File::create(path).with_context(|| format!("Failed to create output file: {}", path.display())) +} + +/// Return a suitable writer for the output path extension #[cfg_attr(not(feature = "compression"), allow(unused_variables))] fn get_writer( output_path: Option<&std::path::Path>, @@ -405,32 +933,24 @@ pub struct FilterSummary { } #[derive(Clone)] -struct FilterProcessor<'a> { +struct FilterProcessor { // Minimizer matching parameters - minimizers: &'a MinimizerSet, + minimizers: Arc, rename: bool, - output_fasta: bool, + discard_quality: bool, debug: bool, check_pairs: bool, /// Write batches in input order, not completion order ordered: bool, kernel: FilterKernel, - // Local buffers - local_buffer: Vec, - local_buffer2: Vec, // Second buffer for paired output + output: Output, local_stats: ProcessingStats, - // Headers awaiting a number, and how many records/pairs this batch kept - pending_renames: Vec, - pending_renames2: Vec, // Parallel to local_buffer2 - batch_kept: u64, - /// Shared across workers, handing each batch the first number of its block + /// Shared across workers, handing out rename numbers rename_counter: Arc, // Global state - global_writer: Arc>, - global_writer2: Option>>, global_stats: Arc>, spinner: Option>>, filtering_start_time: Instant, @@ -446,21 +966,20 @@ pub(crate) struct ProcessingStats { pub last_reported: u64, } -impl<'a> FilterProcessor<'a> { +impl FilterProcessor { fn new( - minimizers: &'a MinimizerSet, + minimizers: Arc, kmer_length: u8, window_size: u8, config: &FilterProcessorConfig, - writer: BoxedWriter, - writer2: Option, + output: Output, spinner: Option>>, filtering_start_time: Instant, ) -> Result { Ok(Self { minimizers, rename: config.rename, - output_fasta: config.output_fasta, + discard_quality: config.discard_quality, debug: config.debug, check_pairs: config.check_pairs, ordered: config.ordered, @@ -474,15 +993,9 @@ impl<'a> FilterProcessor<'a> { prefix_length: config.prefix_length, }, )?, - local_buffer: Vec::with_capacity(DEFAULT_BUFFER_SIZE), - local_buffer2: Vec::with_capacity(DEFAULT_BUFFER_SIZE), - pending_renames: Vec::new(), - pending_renames2: Vec::new(), - batch_kept: 0, + output, local_stats: ProcessingStats::default(), rename_counter: Arc::new(AtomicU64::new(0)), - global_writer: Arc::new(Mutex::new(writer)), - global_writer2: writer2.map(|w| Arc::new(Mutex::new(w))), global_stats: Arc::new(Mutex::new(ProcessingStats::default())), spinner, filtering_start_time, @@ -490,72 +1003,12 @@ impl<'a> FilterProcessor<'a> { } fn should_keep_sequence(&mut self, seq: &[u8]) -> FilterDecision { - self.kernel.classify_read(self.minimizers, seq, self.debug) + self.kernel.classify_read(&self.minimizers, seq, self.debug) } fn should_keep_pair(&mut self, seq1: &[u8], seq2: &[u8]) -> FilterDecision { self.kernel - .classify_pair(self.minimizers, seq1, seq2, self.debug) - } - - fn write_record( - &mut self, - record: &Rf, - seq: &[u8], - read_suffix: &'static [u8], - ) -> Result<()> { - let offset = self.local_buffer.len(); - let marker = format_record_to_buffer( - record, - seq, - self.rename, - self.output_fasta, - &mut self.local_buffer, - )?; - if self.rename { - self.pending_renames.push(PendingRename { - offset, - ordinal: self.batch_kept, - marker, - suffix: read_suffix, - }); - } - Ok(()) - } - - fn write_record_to_buffer2( - &mut self, - record: &Rf, - seq: &[u8], - read_suffix: &'static [u8], - ) -> Result<()> { - let offset = self.local_buffer2.len(); - let marker = format_record_to_buffer( - record, - seq, - self.rename, - self.output_fasta, - &mut self.local_buffer2, - )?; - if self.rename { - self.pending_renames2.push(PendingRename { - offset, - ordinal: self.batch_kept, - marker, - suffix: read_suffix, - }); - } - Ok(()) - } - - /// Claim this batch's block of numbers, returning the first - /// - /// Called under the writer lock: numbers are monotone in output order, - /// but deterministic only with `--ordered` or `-t 1`. - fn reserve_rename_ids(&self) -> u64 { - self.rename_counter - .fetch_add(self.batch_kept, Ordering::Relaxed) - + 1 + .classify_pair(&self.minimizers, seq1, seq2, self.debug) } fn update_spinner(&self) { @@ -592,25 +1045,19 @@ impl<'a> FilterProcessor<'a> { )); } } -} - -impl<'a, Rf: Record> ParallelProcessor for FilterProcessor<'a> { - fn requires_ordering(&self) -> bool { - self.ordered - } - fn process_record(&mut self, record: Rf) -> paraseq::parallel::Result<()> { - let seq = record.seq(); + /// Shared per-read logic for every reader (paraseq FASTX, CBQ mmap) + fn handle_read(&mut self, read: &ReadView) -> Result<()> { self.local_stats.total_seqs += 1; - self.local_stats.total_bp += seq.len() as u64; + self.local_stats.total_bp += read.seq.len() as u64; - let decision = self.should_keep_sequence(&seq); + let decision = self.should_keep_sequence(read.seq); // Show debug info for sequences with hits if self.debug { eprintln!( "DEBUG: {} hits={}/{} keep={} kmers=[{}]", - String::from_utf8_lossy(record.id()), + String::from_utf8_lossy(read.id), decision.hit_count, decision.total_minimizers, decision.keep, @@ -619,91 +1066,42 @@ impl<'a, Rf: Record> ParallelProcessor for FilterProcessor<'a> { } if decision.keep { - self.local_stats.output_bp += seq.len() as u64; - self.write_record(&record, &seq, b"")?; - if self.rename { - self.batch_kept += 1; - } + self.local_stats.output_bp += read.seq.len() as u64; + self.output.push_read( + read, + self.rename, + self.discard_quality, + &self.rename_counter, + )?; } else { self.local_stats.filtered_seqs += 1; - self.local_stats.filtered_bp += seq.len() as u64; + self.local_stats.filtered_bp += read.seq.len() as u64; } Ok(()) } - fn on_batch_complete(&mut self) -> paraseq::parallel::Result<()> { - // Write buffer to output - if !self.local_buffer.is_empty() { - let mut global_writer = self.global_writer.lock(); - if self.rename { - let base = self.reserve_rename_ids(); - write_renamed( - &mut global_writer, - &self.local_buffer, - &self.pending_renames, - base, - )?; - } else { - global_writer.write_all(&self.local_buffer)?; - } - global_writer.flush()?; - } - - // Clear buffer after releasing the lock - self.local_buffer.clear(); - self.pending_renames.clear(); - self.batch_kept = 0; - - // Update global stats - { - let mut stats = self.global_stats.lock(); - stats.total_seqs += self.local_stats.total_seqs; - stats.filtered_seqs += self.local_stats.filtered_seqs; - stats.total_bp += self.local_stats.total_bp; - stats.output_bp += self.local_stats.output_bp; - stats.filtered_bp += self.local_stats.filtered_bp; - } - - // Update spinner - self.update_spinner(); - - // Reset local stats - self.local_stats = ProcessingStats::default(); - - Ok(()) - } -} - -impl<'a, Rf: Record> PairedParallelProcessor for FilterProcessor<'a> { - fn requires_ordering(&self) -> bool { - self.ordered - } - - fn process_record_pair(&mut self, record1: Rf, record2: Rf) -> paraseq::parallel::Result<()> { - if self.check_pairs && !paired_record_names_match(record1.id(), record2.id()) { + /// Shared per-pair logic for every reader (paraseq FASTX, CBQ mmap) + fn handle_pair(&mut self, read1: &ReadView, read2: &ReadView) -> Result<()> { + if self.check_pairs && !paired_record_names_match(read1.id, read2.id) { return Err(anyhow::anyhow!( "Paired record name mismatch: R1='{}', R2='{}'. Expected matching Illumina CASAVA 1: and 2: fields or names suffixed with /1 and /2", - String::from_utf8_lossy(record1.id()), - String::from_utf8_lossy(record2.id()) - ) - .into()); + String::from_utf8_lossy(read1.id), + String::from_utf8_lossy(read2.id) + )); } - let seq1 = record1.seq(); - let seq2 = record2.seq(); - self.local_stats.total_seqs += 2; - self.local_stats.total_bp += (seq1.len() + seq2.len()) as u64; + self.local_stats.total_bp += (read1.seq.len() + read2.seq.len()) as u64; - let decision = self.should_keep_pair(&seq1, &seq2); + let decision = self.should_keep_pair(read1.seq, read2.seq); // Debug info for interleaved pairs if self.debug && decision.hit_count > 0 { eprintln!( "DEBUG: {}/{} hits={}/{} keep={} kmers=[{}]", - String::from_utf8_lossy(record1.id()), - String::from_utf8_lossy(record2.id()), + String::from_utf8_lossy(read1.id), + String::from_utf8_lossy(read2.id), decision.hit_count, decision.total_minimizers, decision.keep, @@ -712,77 +1110,25 @@ impl<'a, Rf: Record> PairedParallelProcessor for FilterProcessor<'a> { } if decision.keep { - self.local_stats.output_bp += (seq1.len() + seq2.len()) as u64; - // Both mates take the pair's number - if self.global_writer2.is_some() { - // Separate outputs - self.write_record(&record1, &seq1, b"/1")?; - self.write_record_to_buffer2(&record2, &seq2, b"/2")?; - } else { - // Interleaved output - self.write_record(&record1, &seq1, b"/1")?; - self.write_record(&record2, &seq2, b"/2")?; - } - if self.rename { - self.batch_kept += 1; - } + self.local_stats.output_bp += (read1.seq.len() + read2.seq.len()) as u64; + self.output.push_pair( + read1, + read2, + self.rename, + self.discard_quality, + &self.rename_counter, + )?; } else { self.local_stats.filtered_seqs += 2; - self.local_stats.filtered_bp += (seq1.len() + seq2.len()) as u64; + self.local_stats.filtered_bp += (read1.seq.len() + read2.seq.len()) as u64; } Ok(()) } - fn on_batch_complete(&mut self) -> paraseq::parallel::Result<()> { - if let Some(ref writer2) = self.global_writer2 { - // Atomic paired batch writing - if !self.local_buffer.is_empty() || !self.local_buffer2.is_empty() { - let mut writer1 = self.global_writer.lock(); - let mut writer2 = writer2.lock(); - - if self.rename { - // Both mates share one block of numbers - let base = self.reserve_rename_ids(); - write_renamed( - &mut writer1, - &self.local_buffer, - &self.pending_renames, - base, - )?; - write_renamed( - &mut writer2, - &self.local_buffer2, - &self.pending_renames2, - base, - )?; - } else { - writer1.write_all(&self.local_buffer)?; - writer2.write_all(&self.local_buffer2)?; - } - writer1.flush()?; - writer2.flush()?; - } - } else { - // Interleaved output - if !self.local_buffer.is_empty() { - let mut writer = self.global_writer.lock(); - if self.rename { - let base = self.reserve_rename_ids(); - write_renamed(&mut writer, &self.local_buffer, &self.pending_renames, base)?; - } else { - writer.write_all(&self.local_buffer)?; - } - writer.flush()?; - } - } - - // Clear buffer after releasing the lock for better performance - self.local_buffer.clear(); - self.local_buffer2.clear(); - self.pending_renames.clear(); - self.pending_renames2.clear(); - self.batch_kept = 0; + fn flush_batch(&mut self) -> Result<()> { + self.output + .flush_batch(self.rename, &self.rename_counter, self.ordered)?; // Update global stats { @@ -802,6 +1148,106 @@ impl<'a, Rf: Record> PairedParallelProcessor for FilterProcessor<'a> { Ok(()) } + + fn flush_thread(&mut self) -> Result<()> { + self.output.flush_thread() + } +} + +impl ParallelProcessor for FilterProcessor { + fn requires_ordering(&self) -> bool { + self.ordered + } + + fn process_record(&mut self, record: Rf) -> paraseq::parallel::Result<()> { + let seq = record.seq(); + self.handle_read(&ReadView { + id: record.id(), + seq: &seq, + qual: record.qual(), + flag: None, + })?; + Ok(()) + } + + fn on_batch_complete(&mut self) -> paraseq::parallel::Result<()> { + self.flush_batch()?; + Ok(()) + } + + fn on_thread_complete(&mut self) -> paraseq::parallel::Result<()> { + self.flush_thread()?; + Ok(()) + } +} + +impl PairedParallelProcessor for FilterProcessor { + fn requires_ordering(&self) -> bool { + self.ordered + } + + fn process_record_pair(&mut self, record1: Rf, record2: Rf) -> paraseq::parallel::Result<()> { + let seq1 = record1.seq(); + let seq2 = record2.seq(); + self.handle_pair( + &ReadView { + id: record1.id(), + seq: &seq1, + qual: record1.qual(), + flag: None, + }, + &ReadView { + id: record2.id(), + seq: &seq2, + qual: record2.qual(), + flag: None, + }, + )?; + Ok(()) + } + + fn on_batch_complete(&mut self) -> paraseq::parallel::Result<()> { + self.flush_batch()?; + Ok(()) + } + + fn on_thread_complete(&mut self) -> paraseq::parallel::Result<()> { + self.flush_thread()?; + Ok(()) + } +} + +impl binseq::ParallelProcessor for FilterProcessor { + fn process_record(&mut self, record: R) -> binseq::Result<()> { + let read1 = ReadView { + id: record.sheader(), + seq: record.sseq(), + qual: record.has_quality().then(|| record.squal()), + flag: record.flag(), + }; + if record.is_paired() { + let read2 = ReadView { + id: record.xheader(), + seq: record.xseq(), + qual: record.has_quality().then(|| record.xqual()), + flag: record.flag(), + }; + self.handle_pair(&read1, &read2)?; + } else { + self.handle_read(&read1)?; + } + Ok(()) + } + + fn on_batch_complete(&mut self) -> binseq::Result<()> { + self.flush_batch()?; + Ok(()) + } + + fn on_thread_complete(&mut self) -> binseq::Result<()> { + self.flush_thread()?; + Ok(()) + } } pub fn run(config: &FilterConfig) -> Result { @@ -809,10 +1255,6 @@ pub fn run(config: &FilterConfig) -> Result { if let Some(threshold) = config.complexity_threshold { validate_unit_interval("complexity threshold", threshold)?; } - validate_check_pairs_mode( - config.check_pairs, - config.interleaved || config.input2_path.is_some(), - )?; // Validate the index path once here; run_with_index never touches it again. if !config.minimizers_path.exists() { @@ -838,10 +1280,11 @@ pub fn run(config: &FilterConfig) -> Result { summary_path: config.summary_path.cloned(), deplete: config.deplete, rename: config.rename, - output_fasta: config.output_fasta, + discard_quality: config.discard_quality, ordered: config.ordered, threads: config.threads, compression_level: config.compression_level, + cbq_block_size: config.cbq_block_size, compression_threads: config.compression_threads, debug: config.debug, quiet: config.quiet, @@ -874,7 +1317,7 @@ pub fn run(config: &FilterConfig) -> Result { threshold ); } - return run_with_index(&minimizers, &header, &run_config); + return run_with_index(Arc::new(minimizers), &header, &run_config); } let (minimizers, header) = load_minimizers_cached(config.minimizers_path)?; @@ -895,15 +1338,11 @@ pub fn run(config: &FilterConfig) -> Result { /// Reusable entry point behind the Python bindings... load the index once, call repeatedly. /// Does no index-path validation; the index is already in memory. pub fn run_with_index( - minimizers: &MinimizerSet, + minimizers: Arc, header: &IndexHeader, config: &FilterRunConfig, ) -> Result { validate_unit_interval("relative threshold", config.rel_threshold)?; - validate_check_pairs_mode( - config.check_pairs, - config.interleaved || config.input2_path.is_some(), - )?; let start_time = Instant::now(); let version: String = env!("CARGO_PKG_VERSION").to_string(); @@ -925,7 +1364,7 @@ pub fn run_with_index( let compressed_output_count = count_compressed_outputs(config); // Allocate threads between filtering (rayon) and compression (gzp). - // Rayon pool can only be initialized once, so calculate before build_global(). + // Rayon pool can only be initialised once, so calculate before build_global(). let (filtering_threads, compression_threads_per_output) = if compressed_output_count > 0 { let compression_threads_total = if config.compression_threads > 0 { config.compression_threads as usize @@ -943,22 +1382,48 @@ pub fn run_with_index( }; if filtering_threads > 0 { - // error is OK here when we initialize a 2nd time in server mode. + // error is OK here when we initialise a 2nd time in server mode. let _ = rayon::ThreadPoolBuilder::new() .num_threads(filtering_threads) .build_global() - .context("Failed to initialize thread pool"); + .context("Failed to initialise thread pool"); } + check_input_paths(config)?; + + // Resolve formats and input metadata before opening or truncating outputs + let interleaved_stdin = config.input_path == "-" && config.input2_path.as_deref() == Some("-"); + let interleaved_input = config.interleaved || interleaved_stdin; + let output_format = resolve_output_format(config); + // `.cba` names a quality-free CBQ, implying --discard-quality for the output + let discard_quality = config.discard_quality + || config + .output_path + .as_deref() + .is_some_and(|path| is_cba_path(&path.to_string_lossy())); + let (layout, input) = open_input(config, interleaved_input)?; + validate_input_output(&layout, output_format, config)?; + + // The binseq reader and CBQ rename numbering cannot honour --ordered + // multithreaded + let ordered_cbq = config.ordered + && (layout.format == Format::Cbq || (output_format == Format::Cbq && config.rename)); + let filtering_threads = if ordered_cbq && filtering_threads > 1 { + if !quiet { + eprintln!("Using 1 filtering thread: --ordered with CBQ input or renamed CBQ output"); + } + 1 + } else { + filtering_threads + }; + let mode = if config.deplete { "deplete" } else { "search" }; let mut input_type = String::new(); let mut options = Vec::::new(); - let interleaved_stdin = config.input_path == "-" && config.input2_path.as_deref() == Some("-"); - let interleaved_input = config.interleaved || interleaved_stdin; if interleaved_input { input_type.push_str("interleaved"); - } else if config.input2_path.is_some() { + } else if layout.paired { input_type.push_str("paired"); } else { input_type.push_str("single"); @@ -1003,25 +1468,47 @@ pub fn run_with_index( ); } - check_input_paths(config)?; - - let writer = get_writer( - config.output_path.as_deref(), - config.compression_level, - compression_threads_per_output, - )?; - let writer2 = if let Some(output2) = config.output2_path.as_deref() { - if config.input2_path.is_some() || interleaved_input { - Some(get_writer( - Some(std::path::Path::new(output2)), + let output = match output_format { + Format::Cbq => { + let path = config + .output_path + .as_deref() + .expect("CBQ output implies a named path"); + let cbq_writer = BinseqWriterBuilder::new(BinseqFormat::Cbq) + .paired(layout.paired) + .quality(layout.qualities && !discard_quality) + .headers(layout.headers || config.rename) + .flags(layout.flags) + .block_size( + (config.cbq_block_size as usize * 1024 * 1024) + .max(layout.block_size.unwrap_or(0)), + ) + .compression_level(i32::from(config.compression_level)) + .build(open_cbq_output(path)?) + .context("Failed to create CBQ writer")?; + Output::cbq(cbq_writer)? + } + Format::Fastx => { + let writer = get_writer( + config.output_path.as_deref(), config.compression_level, compression_threads_per_output, - )?) - } else { - None + )?; + let writer2 = if let Some(output2) = config.output2_path.as_deref() { + if layout.paired { + Some(get_writer( + Some(std::path::Path::new(output2)), + config.compression_level, + compression_threads_per_output, + )?) + } else { + None + } + } else { + None + }; + Output::fastx(writer, writer2) } - } else { - None }; // Progress bar setup if not quiet @@ -1047,7 +1534,7 @@ pub fn run_with_index( prefix_length: config.prefix_length, deplete: config.deplete, rename: config.rename, - output_fasta: config.output_fasta, + discard_quality, debug: config.debug, check_pairs: config.check_pairs, ordered: config.ordered, @@ -1057,100 +1544,31 @@ pub fn run_with_index( kmer_length, window_size, &processor_config, - writer, - writer2, + output, spinner.clone(), filtering_start_time, )?; - // Check for empty files via metadata (fast path for uncompressed files <5 bytes) - let input1_empty = is_empty_file(&config.input_path)?; - let input2_empty = config - .input2_path - .as_deref() - .map(is_empty_file) - .transpose()? - .unwrap_or(false); - // Process based on input type - use filtering threads (already calculated above) let num_threads = filtering_threads; - if interleaved_input { - // Interleaved paired input from stdin or a file - if input1_empty { - if !quiet { - eprintln!("Empty input file(s) detected"); - } - } else { - match create_paraseq_reader(Some(config.input_path.as_str())) { - Ok(reader) => { - reader.process_parallel_interleaved(&mut processor, num_threads)?; - } - Err(e) if is_empty_input_error(&e) => { - if !quiet { - eprintln!("Empty input file(s) detected"); - } - } - Err(e) => return Err(e), - } + match input { + Input::Cbq(reader) => { + reader.process_parallel(processor.clone(), num_threads)?; } - } else if let Some(input2_path) = config.input2_path.as_deref() { - // Paired files - both must be empty or both non-empty - if input1_empty && input2_empty { - if !quiet { - eprintln!("Empty input file(s) detected"); - } - } else if input1_empty || input2_empty { - return Err(anyhow::anyhow!( - "One paired file is empty but the other is not" - )); - } else { - // Try to create readers, catching empty compressed files - let r1_result = create_paraseq_reader(Some(config.input_path.as_str())); - let r2_result = create_paraseq_reader(Some(input2_path)); - - match (r1_result, r2_result) { - (Ok(r1), Ok(r2)) => { - r1.process_parallel_paired(r2, &mut processor, num_threads)?; - } - (Err(e1), Err(e2)) if is_empty_input_error(&e1) && is_empty_input_error(&e2) => { - if !quiet { - eprintln!("Empty input file(s) detected"); - } - } - (Err(e), _) if is_empty_input_error(&e) => { - return Err(anyhow::anyhow!( - "First paired file appears empty while second is not" - )); - } - (_, Err(e)) if is_empty_input_error(&e) => { - return Err(anyhow::anyhow!( - "Second paired file appears empty while first is not" - )); - } - (Err(e), _) => return Err(e), - (_, Err(e)) => return Err(e), - } + Input::FastxSingle(reader) => { + reader.process_parallel(&mut processor, num_threads)?; } - } else { - // Single file or stdin - if input1_empty { + Input::FastxInterleaved(reader) => { + reader.process_parallel_interleaved(&mut processor, num_threads)?; + } + Input::FastxPaired(reader1, reader2) => { + reader1.process_parallel_paired(reader2, &mut processor, num_threads)?; + } + Input::Empty => { if !quiet { eprintln!("Empty input file(s) detected"); } - } else { - // Try to create reader, catching empty compressed files - match create_paraseq_reader(Some(config.input_path.as_str())) { - Ok(reader) => { - reader.process_parallel(&mut processor, num_threads)?; - } - Err(e) if is_empty_input_error(&e) => { - if !quiet { - eprintln!("Empty input file(s) detected"); - } - } - Err(e) => return Err(e), - } } } @@ -1163,11 +1581,10 @@ pub fn run_with_index( drop(final_stats); // Release lock - // Flush writers - they should auto-flush on drop - drop(processor.global_writer); - if let Some(w2) = processor.global_writer2 { - drop(w2); - } + // Finish any CBQ stream (writes the embedded index), then drop the + // processor so the writers flush + processor.output.finish()?; + drop(processor); let total_time = start_time.elapsed(); let filtering_time = filtering_start_time.elapsed(); diff --git a/src/index.rs b/src/index.rs index 7e4426b..5fbf49b 100644 --- a/src/index.rs +++ b/src/index.rs @@ -167,7 +167,7 @@ pub fn load_header_and_count>(path: &P) -> Result<(IndexHeader, u } #[cfg(feature = "cli")] -static INDEX: OnceLock<(PathBuf, crate::MinimizerSet, IndexHeader)> = OnceLock::new(); +static INDEX: OnceLock<(PathBuf, Arc, IndexHeader)> = OnceLock::new(); #[cfg(feature = "cli")] pub fn current_index_path() -> Option { @@ -177,18 +177,18 @@ pub fn current_index_path() -> Option { #[cfg(feature = "cli")] pub fn load_minimizers_cached( path: &Path, -) -> Result<(&'static crate::MinimizerSet, &'static IndexHeader)> { +) -> Result<(Arc, &'static IndexHeader)> { let (p, minimizers, header) = INDEX.get_or_init(|| { // Auto-detect exact vs BFF format let (m, h) = load_index_from_path_auto(path).unwrap(); - (path.to_owned(), m, h) + (path.to_owned(), Arc::new(m), h) }); assert_eq!( p, path, "Currently, the server can only have one index loaded." ); - Ok((minimizers, header)) + Ok((Arc::clone(minimizers), header)) } /// Load minimizers from a reader (generic over any Read impl) diff --git a/src/lib.rs b/src/lib.rs index dc12d4f..7612ccd 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -16,7 +16,9 @@ mod minimizers; // Public API #[cfg(feature = "cli")] -pub use filter::{FilterRunConfig, FilterSummary, run as run_filter, run_with_index}; +pub use filter::{ + DEFAULT_CBQ_BLOCK_SIZE_MIB, FilterRunConfig, FilterSummary, run as run_filter, run_with_index, +}; pub use filter_kernel::{FilterDecision, FilterKernel, FilterParams}; #[cfg(feature = "fetch")] pub use index::fetch as index_fetch; @@ -331,14 +333,14 @@ pub struct FilterConfig<'a> { /// Path to JSON summary file pub summary_path: Option<&'a PathBuf>, - /// Deplete mode (remove sequences WITH matches, original deacon behavior) + /// Deplete mode (remove sequences WITH matches, original deacon behaviour) pub deplete: bool, /// Replace sequence headers with incrementing numbers (1, 2, 3...) pub rename: bool, - /// Force FASTA output (discards quality scores) - pub output_fasta: bool, + /// Emit fasta or quality-free cbq regardless of input format (implied by a .cba output) + pub discard_quality: bool, /// Preserve input record ordering (deterministic, slightly slower) pub ordered: bool, @@ -349,6 +351,9 @@ pub struct FilterConfig<'a> { /// Compression level for output files (1-22 for zst, 1-9 for gz) pub compression_level: u8, + /// cbq output block size in MiB (raised to the cbq input's block size if larger) + pub cbq_block_size: u16, + /// Number of threads for compression (0 = auto-calculate as ceil(total/2)) pub compression_threads: u16, diff --git a/src/main.rs b/src/main.rs index 449e79e..d97c00b 100644 --- a/src/main.rs +++ b/src/main.rs @@ -3,8 +3,9 @@ use clap::{Parser, Subcommand}; #[cfg(feature = "fetch")] use deacon::index_fetch; use deacon::{ - ComplexityAlgorithm, DEFAULT_KMER_LENGTH, DEFAULT_WINDOW_SIZE, FilterConfig, IndexConfig, - index_diff, index_dump, index_filter, index_freeze, index_info, index_intersect, index_union, + ComplexityAlgorithm, DEFAULT_CBQ_BLOCK_SIZE_MIB, DEFAULT_KMER_LENGTH, DEFAULT_WINDOW_SIZE, + FilterConfig, IndexConfig, index_diff, index_dump, index_filter, index_freeze, index_info, + index_intersect, index_union, }; use serde::{Deserialize, Serialize}; use std::io::{Read, Write}; @@ -29,12 +30,12 @@ enum Command { #[command(subcommand)] command: IndexCommand, }, - /// Retain or deplete sequence records with sufficient minimizer hits to an indexed query + /// Retain or deplete sequence records with sufficient minimizer hits to the index Filter { /// Path to minimizer index file index: PathBuf, - /// Optional path to fastx file (or - for stdin) + /// Optional path to fastx or binseq cbq file (or - for stdin) #[arg(default_value = "-")] input: String, @@ -61,19 +62,15 @@ enum Command { #[arg(short = 'd', long = "deplete", default_value_t = false)] deplete: bool, - /// Replace sequence headers with incrementing numbers (reproducible with --ordered) + /// Replace sequence headers with incrementing numbers (deterministic with --ordered) #[arg(short = 'R', long = "rename", default_value_t = false)] rename: bool, - /// Output FASTA format regardless of input format - #[arg(short = 'f', long = "fasta", default_value_t = false)] - output_fasta: bool, - - /// Path to output fastx file (stdout if not specified; detects .gz and .zst) + /// Path to output file (fastx to stdout by default; detects .gz, .zst, .xz, .cbq, .cba) #[arg(short = 'o', long = "output")] output: Option, - /// Optional path to second paired output fastx file (detects .gz and .zst) + /// Optional path to second paired output fastx file (detects .gz, .zst, .xz) #[arg(short = 'O', long = "output2")] output2: Option, @@ -89,11 +86,23 @@ enum Command { #[arg(long = "compression-threads", default_value_t = 0)] compression_threads: u16, - /// Output compression level (1-9 for gz & xz; 1-22 for zstd) + /// Output compression level (1-9 for gz & xz; 1-22 for zstd including cbq) #[arg(long = "compression-level", default_value_t = 2)] compression_level: u8, - /// Treat INPUT as interleaved paired reads from a file or stdin + /// cbq output block size in MiB (or cbq input block size if higher) + #[arg( + long = "cbq-block-size", + default_value_t = DEFAULT_CBQ_BLOCK_SIZE_MIB, + value_parser = clap::value_parser!(u16).range(1..=1024) + )] + cbq_block_size: u16, + + /// Emit fasta or quality-free cbq regardless of input format + #[arg(long = "discard-quality", default_value_t = false)] + discard_quality: bool, + + /// Treat INPUT as interleaved paired records from single file or stdin #[arg( long = "interleaved", default_value_t = false, @@ -109,13 +118,13 @@ enum Command { #[arg(long = "check-pairs", default_value_t = false)] check_pairs: bool, + /// Emit sequences with minimizer hits to stderr + #[arg(long = "debug", default_value_t = false)] + debug: bool, + /// Suppress progress reporting #[arg(short = 'q', long = "quiet", default_value_t = false)] quiet: bool, - - /// Output sequences with minimizer hits to stderr - #[arg(long = "debug", default_value_t = false)] - debug: bool, }, /// Start/stop a server process for reduced latency filtering Server { @@ -163,7 +172,7 @@ enum IndexCommand { #[arg(short = 't', long = "threads", default_value_t = 8)] threads: u16, - /// Suppress sequence header output + /// Suppress progress reporting #[arg(short = 'q', long = "quiet")] quiet: bool, }, @@ -344,7 +353,7 @@ fn main() -> Result<()> { rayon::ThreadPoolBuilder::new() .num_threads(*threads as usize) .build_global() - .context("Failed to initialize thread pool")?; + .context("Failed to initialise thread pool")?; // Remove existing socket if present let _ = std::fs::remove_file("deacon_server_socket"); @@ -530,9 +539,10 @@ fn process_command(command: &Command) -> Result<(), anyhow::Error> { summary, deplete, rename, - output_fasta, + discard_quality, threads, compression_level, + cbq_block_size, compression_threads, ordered, check_pairs, @@ -561,10 +571,11 @@ fn process_command(command: &Command) -> Result<(), anyhow::Error> { summary_path: summary.as_ref(), deplete: *deplete, rename: *rename, - output_fasta: *output_fasta, + discard_quality: *discard_quality, ordered: *ordered, threads: *threads, compression_level: *compression_level, + cbq_block_size: *cbq_block_size, compression_threads: *compression_threads, debug: *debug, quiet: *quiet, diff --git a/src/minimizers.rs b/src/minimizers.rs index 1fbf2ff..d7e498b 100644 --- a/src/minimizers.rs +++ b/src/minimizers.rs @@ -258,6 +258,24 @@ mod tests { assert!(short_minimizers.is_empty()); } + #[test] + fn non_acgt_mapped_to_n() { + let hasher = KmerHasher::new(31); + let minimizers = |base| { + let mut seq = b"ACGT".repeat(40); + seq[64] = base; + match compute_minimizers(&seq, &hasher, 31, 15) { + crate::MinimizerVec::U64(vec) => vec, + crate::MinimizerVec::U128(_) => unreachable!(), + } + }; + let expected = minimizers(b'N'); + + for &base in b"nRrYySsWwKkMmBbDdHhVvUu-" { + assert_eq!(minimizers(base), expected, "base {}", base as char); + } + } + #[test] fn test_calculate_scaled_entropy() { // Test short k-mers (should return 1.0 for k < 10) diff --git a/tests/filter_tests.rs b/tests/filter_tests.rs index cb2f2b6..70e6ae7 100644 --- a/tests/filter_tests.rs +++ b/tests/filter_tests.rs @@ -177,7 +177,7 @@ fn test_filter_to_file() { assert!(output_path.exists(), "Output file wasn't created"); assert!(summary_path.exists(), "Summary file wasn't created"); - // With new default behavior: sequences without matches are filtered out (sequences too short for k=31) + // With new default behaviour: sequences without matches are filtered out (sequences too short for k=31) let output_content = fs::read_to_string(&output_path).unwrap(); assert!( output_content.is_empty(), @@ -293,7 +293,7 @@ fn test_filter_rename() { } #[test] -fn test_filter_fasta_flag() { +fn test_filter_discard_quality_flag() { let temp_dir = tempdir().unwrap(); let fasta_path = temp_dir.path().join("ref.fasta"); let fastq_path = temp_dir.path().join("reads.fastq"); @@ -306,7 +306,7 @@ fn test_filter_fasta_flag() { let mut cmd = cargo::cargo_bin_cmd!("deacon"); let output = cmd .arg("filter") - .arg("-f") + .arg("--discard-quality") .arg("-a") .arg("1") .arg("-r") @@ -322,7 +322,7 @@ fn test_filter_fasta_flag() { let output_str = std::str::from_utf8(&output).unwrap(); assert!( output_str.starts_with('>'), - "FASTQ in with -f should gen FASTA out" + "FASTQ in with -Q should gen FASTA out" ); } @@ -530,18 +530,25 @@ fn test_check_pairs_rejects_mismatched_names() { } #[test] -fn test_check_pairs_requires_paired_input_before_index_load() { +fn test_check_pairs_requires_paired_input() { + let temp_dir = tempdir().unwrap(); + let fasta_path = temp_dir.path().join("ref.fasta"); + let bin_path = temp_dir.path().join("ref.bin"); + let fastq_path = temp_dir.path().join("reads.fastq"); + create_test_fasta(&fasta_path); + create_test_fastq(&fastq_path); + build_index(&fasta_path, &bin_path); + let mut cmd = cargo::cargo_bin_cmd!("deacon"); cmd.arg("filter") .arg("--check-pairs") - .arg("missing.idx") - .arg("sequences.fastq") + .arg(&bin_path) + .arg(&fastq_path) .assert() .failure() .stderr(predicates::str::contains( "--check-pairs requires paired input", - )) - .stderr(predicates::str::contains("Index file does not exist").not()); + )); } #[test] @@ -1113,7 +1120,7 @@ fn test_shared_minimizer_counted_once() { // If shared minimizers are counted once (correct): total hits = 1, pair kept (1 < 2) // If shared minimizers are counted twice (bug): total hits = 2+, pair filtered (2+ >= 2) - // Using --deplete to restore original behavior for this bug test + // Using --deplete to restore original behaviour for this bug test let mut cmd = cargo::cargo_bin_cmd!("deacon"); cmd.arg("filter") .arg("--deplete") @@ -2783,3 +2790,1099 @@ fn test_filter_paired_deplete_with_rename() { ); } } + +/// Write a small paired, headerless, quality-free CBQ file for validation tests +fn write_headerless_paired_cbq(path: &Path) { + use binseq::SequencingRecordBuilder; + use binseq::write::{BinseqWriterBuilder, Format}; + let file = File::create(path).unwrap(); + let mut writer = BinseqWriterBuilder::new(Format::Cbq) + .paired(true) + .build(file) + .unwrap(); + for seq in [ + b"ACGTACGTACGTACGTACGT".as_slice(), + b"TGCAACGTACGTACGTACGT".as_slice(), + ] { + writer + .push( + SequencingRecordBuilder::default() + .s_seq(seq) + .x_seq(seq) + .build() + .unwrap(), + ) + .unwrap(); + } + writer.finish().unwrap(); +} + +/// CBQ output block size is --cbq-block-size raised to the CBQ input's own block +/// size, so ultra-long records always round-trip; FASTX has none to inherit. +#[test] +fn cbq_output_block_size() { + use binseq::SequencingRecordBuilder; + use binseq::write::{BinseqWriterBuilder, Format}; + + let temp_dir = tempdir().unwrap(); + let fasta_path = temp_dir.path().join("ref.fasta"); + let bin_path = temp_dir.path().join("ref.bin"); + let fastq_path = temp_dir.path().join("reads.fastq"); + create_test_fasta(&fasta_path); + create_test_fastq(&fastq_path); + build_index(&fasta_path, &bin_path); + + let input = temp_dir.path().join("in.cbq"); + let mut writer = BinseqWriterBuilder::new(Format::Cbq) + .block_size(32 << 20) + .build(File::create(&input).unwrap()) + .unwrap(); + writer + .push( + SequencingRecordBuilder::default() + .s_seq(b"ACGTACGTACGTACGTACGT") + .build() + .unwrap(), + ) + .unwrap(); + writer.finish().unwrap(); + + let cases = [ + (&fastq_path, None, 16 << 20), // default + (&fastq_path, Some("64"), 64 << 20), // flag honoured + (&input, None, 32 << 20), // input's block size wins over the default + (&input, Some("4"), 32 << 20), // clamped up to the input's + (&input, Some("64"), 64 << 20), // flag wins when larger + ]; + + for (case, (source, block_size, expected)) in cases.iter().enumerate() { + let output = temp_dir.path().join(format!("out{case}.cbq")); + let mut cmd = cargo::cargo_bin_cmd!("deacon"); + cmd.args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(source) + .arg("--output") + .arg(&output); + if let Some(block_size) = block_size { + cmd.args(["--cbq-block-size", block_size]); + } + cmd.assert().success(); + let reader = binseq::cbq::MmapReader::new(&output).unwrap(); + assert_eq!( + reader.header().block_size as usize, + *expected, + "case {case}" + ); + } +} + +/// --cbq-block-size is bounded by clap, before any work starts +#[test] +fn cbq_block_size_out_of_range_rejected() { + let temp_dir = tempdir().unwrap(); + let fasta_path = temp_dir.path().join("ref.fasta"); + let bin_path = temp_dir.path().join("ref.bin"); + let fastq_path = temp_dir.path().join("reads.fastq"); + create_test_fasta(&fasta_path); + create_test_fastq(&fastq_path); + build_index(&fasta_path, &bin_path); + + for block_size in ["0", "2000"] { + let output = temp_dir.path().join(format!("no_create{block_size}.cbq")); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "--cbq-block-size", block_size]) + .arg(&bin_path) + .arg(&fastq_path) + .arg("--output") + .arg(&output) + .assert() + .code(2) + .stderr(predicate::str::contains("not in 1..=1024")); + assert!(!output.exists()); + } +} + +/// CBQ is only an I/O concern: FASTX -> CBQ -> FASTX must preserve reads, +/// pairing, headers, qualities, and summary counts exactly. +#[test] +fn cbq_roundtrip_matches_fastx() { + let temp_dir = tempdir().unwrap(); + let fasta_path = temp_dir.path().join("ref.fasta"); + let bin_path = temp_dir.path().join("ref.bin"); + let fastq_path = temp_dir.path().join("reads.fastq"); + let r1_path = temp_dir.path().join("reads_1.fastq"); + let r2_path = temp_dir.path().join("reads_2.fastq"); + + create_test_fasta(&fasta_path); + create_test_fastq(&fastq_path); + create_test_paired_fastq(&r1_path, &r2_path); + build_index(&fasta_path, &bin_path); + + // Single-end: FASTQ -> FASTQ baseline vs FASTQ -> CBQ + let baseline = temp_dir.path().join("baseline.fastq"); + let cbq_path = temp_dir.path().join("reads.cbq"); + let summary1 = temp_dir.path().join("summary1.json"); + let summary2 = temp_dir.path().join("summary2.json"); + + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&fastq_path) + .arg("--output") + .arg(&baseline) + .arg("--summary") + .arg(&summary1) + .assert() + .success(); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&fastq_path) + .arg("--output") + .arg(&cbq_path) + .arg("--summary") + .arg(&summary2) + .assert() + .success(); + + let baseline_content = fs::read_to_string(&baseline).unwrap(); + assert_eq!(count_records(&baseline_content), 2); + + // The CBQ itself: single, with headers and qualities, same record count + let reader = binseq::cbq::MmapReader::new(&cbq_path).unwrap(); + assert!(!reader.is_paired(), "single-end CBQ must be unpaired"); + assert!(reader.header().has_headers(), "CBQ must keep headers"); + assert!(reader.header().has_qualities(), "CBQ must keep qualities"); + assert_eq!(reader.num_records(), 2); + + // CBQ -> FASTQ matches the FASTQ baseline byte-for-byte + let roundtrip = temp_dir.path().join("roundtrip.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&cbq_path) + .arg("--output") + .arg(&roundtrip) + .assert() + .success(); + assert_eq!( + fs::read(&roundtrip).unwrap(), + fs::read(&baseline).unwrap(), + "CBQ round trip must match the FASTQ output" + ); + + // Summary sequence/base counts agree across formats + let s1: serde_json::Value = + serde_json::from_str(&fs::read_to_string(&summary1).unwrap()).unwrap(); + let s2: serde_json::Value = + serde_json::from_str(&fs::read_to_string(&summary2).unwrap()).unwrap(); + for field in ["seqs_in", "seqs_out", "bp_in", "bp_out"] { + assert_eq!(s1[field], s2[field], "summary field {field} differs"); + } + + // Paired: FASTQ -> one paired CBQ (2 native paired records) + let paired_cbq = temp_dir.path().join("paired.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&r1_path) + .arg(&r2_path) + .arg("--output") + .arg(&paired_cbq) + .assert() + .success(); + let reader = binseq::cbq::MmapReader::new(&paired_cbq).unwrap(); + assert!(reader.is_paired(), "paired FASTQ must produce a paired CBQ"); + assert_eq!(reader.num_records(), 2); + + // CBQ -> CBQ preserves pairing, headers, sequences, qualities + let cbq2 = temp_dir.path().join("paired2.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&paired_cbq) + .arg("--output") + .arg(&cbq2) + .assert() + .success(); + let reader = binseq::cbq::MmapReader::new(&cbq2).unwrap(); + assert!(reader.is_paired(), "CBQ -> CBQ must preserve pairing"); + assert_eq!(reader.num_records(), 2); + + // CBQ -> FASTQ matches the paired FASTQ -> FASTQ interleaved baseline + let paired_baseline = temp_dir.path().join("paired_baseline.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&r1_path) + .arg(&r2_path) + .arg("--output") + .arg(&paired_baseline) + .assert() + .success(); + let paired_roundtrip = temp_dir.path().join("paired_roundtrip.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&cbq2) + .arg("--output") + .arg(&paired_roundtrip) + .assert() + .success(); + assert_eq!( + fs::read(&paired_roundtrip).unwrap(), + fs::read(&paired_baseline).unwrap(), + "paired CBQ round trip must match the interleaved FASTQ output" + ); + + // --discard-quality creates a quality-free CBQ + let fasta_cbq = temp_dir.path().join("fasta.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args([ + "filter", + "-a", + "1", + "-r", + "0.0", + "-t", + "1", + "--discard-quality", + ]) + .arg(&bin_path) + .arg(&fastq_path) + .arg("--output") + .arg(&fasta_cbq) + .assert() + .success(); + let reader = binseq::cbq::MmapReader::new(&fasta_cbq).unwrap(); + assert!( + !reader.header().has_qualities(), + "--discard-quality must produce a quality-free CBQ" + ); + assert_eq!(reader.num_records(), 2); + + // Empty retained output is written as a structurally valid zero-record CBQ + let aaa_path = temp_dir.path().join("aaa.fasta"); + let aaa_bin = temp_dir.path().join("aaa.bin"); + create_test_fasta_aaa(&aaa_path); + build_index(&aaa_path, &aaa_bin); + let empty_cbq = temp_dir.path().join("empty.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-t", "1"]) + .arg(&aaa_bin) + .arg(&fastq_path) + .arg("--output") + .arg(&empty_cbq) + .assert() + .success(); + let empty_bytes = fs::read(&empty_cbq).unwrap(); + assert!( + empty_bytes.starts_with(b"CBQFILE"), + "empty retained output must be a CBQ file" + ); + let reader = binseq::cbq::MmapReader::new(&empty_cbq).unwrap(); + assert_eq!(reader.num_records(), 0, "empty CBQ must hold zero records"); + assert!( + reader.header().has_qualities(), + "empty CBQ must keep the writer quality flag" + ); + + // ... and it reads back as empty FASTX + let empty_out = temp_dir.path().join("empty_out.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&empty_cbq) + .arg("--output") + .arg(&empty_out) + .assert() + .success(); + assert!(fs::read_to_string(&empty_out).unwrap().is_empty()); + + // Invalid CBQ argument combinations fail before output creation + let no_create = temp_dir.path().join("should_not_exist.cbq"); + + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-t", "1"]) + .arg(&bin_path) + .arg(&cbq_path) + .arg(&r2_path) // CBQ input + INPUT2 + .arg("--output") + .arg(&no_create) + .assert() + .failure(); + assert!( + !no_create.exists(), + "CBQ input + INPUT2 must fail before output creation" + ); + + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "--interleaved", "-t", "1"]) + .arg(&bin_path) + .arg(&cbq_path) // CBQ input + --interleaved + .arg("--output") + .arg(&no_create) + .assert() + .failure(); + assert!( + !no_create.exists(), + "CBQ input + --interleaved must fail before output creation" + ); + + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-t", "1"]) + .arg(&bin_path) + .arg(&fastq_path) + .arg("--output") + .arg(&no_create) + .arg("--output2") + .arg(temp_dir.path().join("x.fastq")) // CBQ output + OUTPUT2 + .assert() + .failure(); + assert!( + !no_create.exists(), + "CBQ output + OUTPUT2 must fail before output creation" + ); + + // --check-pairs + headerless CBQ input + let headerless = temp_dir.path().join("headerless.cbq"); + write_headerless_paired_cbq(&headerless); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "--check-pairs", "-t", "1"]) + .arg(&bin_path) + .arg(&headerless) + .arg("--output") + .arg(&no_create) + .assert() + .failure(); + assert!( + !no_create.exists(), + "--check-pairs + headerless CBQ must fail before output creation" + ); +} + +/// CBQ must preserve Ns and filter them like FASTX. +#[test] +fn cbq_ns_match_fastx() { + let temp_dir = tempdir().unwrap(); + let ref_path = temp_dir.path().join("ref.fasta"); + let bin_path = temp_dir.path().join("ref.bin"); + fs::write(&ref_path, format!(">ref\n{}\n", "A".repeat(100))).unwrap(); + build_index(&ref_path, &bin_path); + + let run = |input: &Path, output: &Path, deplete: bool, fasta: bool| { + let mut cmd = cargo::cargo_bin_cmd!("deacon"); + cmd.arg("filter"); + if deplete { + cmd.args(["-d", "-a", "999"]); + } else { + cmd.args(["-a", "1"]); + } + cmd.args(["-r", "0", "-t", "1"]); + if fasta { + cmd.arg("--discard-quality"); + } + cmd.arg(&bin_path) + .arg(input) + .arg("--output") + .arg(output) + .assert() + .success(); + }; + + let all_n = "N".repeat(64); + let split_n = format!("{}N{}", "A".repeat(64), "A".repeat(64)); + let fastq = format!( + "@all-n\n{all_n}\n+\n{}\n@split-n\n{split_n}\n+\n{}\n", + "I".repeat(all_n.len()), + "I".repeat(split_n.len()), + ); + let fasta = format!(">all-n\n{all_n}\n>split-n\n{split_n}\n"); + + for (ext, content, fasta_output, marker) in [ + ("fastq", fastq.as_str(), false, "@"), + ("fasta", fasta.as_str(), true, ">"), + ] { + let input = temp_dir.path().join(format!("reads.{ext}")); + let cbq = temp_dir.path().join(format!("reads-{ext}.cbq")); + let roundtrip = temp_dir.path().join(format!("roundtrip.{ext}")); + let direct = temp_dir.path().join(format!("direct-filtered.{ext}")); + let via_cbq = temp_dir.path().join(format!("cbq-filtered.{ext}")); + fs::write(&input, content).unwrap(); + + run(&input, &cbq, true, fasta_output); + run(&cbq, &roundtrip, true, fasta_output); + assert_eq!(fs::read(&roundtrip).unwrap(), content.as_bytes()); + + run(&input, &direct, false, fasta_output); + run(&cbq, &via_cbq, false, fasta_output); + let filtered = fs::read(&direct).unwrap(); + assert_eq!(filtered, fs::read(&via_cbq).unwrap()); + let filtered = String::from_utf8_lossy(&filtered); + assert!(!filtered.contains(&format!("{marker}all-n"))); + assert!(filtered.contains(&format!("{marker}split-n"))); + } +} + +/// CBQ encodes IUPAC codes as N across 32bp block boundaries. +#[test] +fn cbq_iupac_codes_become_n() { + let temp_dir = tempdir().unwrap(); + let ref_path = temp_dir.path().join("ref.fasta"); + let bin_path = temp_dir.path().join("ref.bin"); + fs::write(&ref_path, format!(">ref\n{}\n", "A".repeat(100))).unwrap(); + build_index(&ref_path, &bin_path); + + let mut seq = "ACGT".repeat(25).into_bytes(); + for (code, pos) in [ + ('Y', 0), + ('R', 15), + ('W', 31), + ('S', 32), + ('K', 33), + ('B', 40), + ('D', 41), + ('H', 42), + ('V', 43), + ('b', 44), + ('d', 45), + ('h', 46), + ('v', 47), + ('M', 63), + ('y', 1), + ('r', 16), + ('w', 30), + ('s', 34), + ('k', 62), + ('m', 64), + ('n', 99), + ] { + seq[pos] = code as u8; + } + let seq = String::from_utf8(seq).unwrap(); + let expected_seq: String = seq + .bytes() + .map(|b| match b { + b'A' | b'C' | b'G' | b'T' => b as char, + _ => 'N', + }) + .collect(); + + for (ext, input, expected, fasta) in [ + ( + "fasta", + format!(">iupac\n{seq}\n"), + format!(">iupac\n{expected_seq}\n"), + true, + ), + ( + "fastq", + format!("@iupac\n{seq}\n+\n{}\n", "I".repeat(seq.len())), + format!("@iupac\n{expected_seq}\n+\n{}\n", "I".repeat(seq.len())), + false, + ), + ] { + let input_path = temp_dir.path().join(format!("reads.{ext}")); + let cbq_path = temp_dir.path().join(format!("reads-{ext}.cbq")); + let output_path = temp_dir.path().join(format!("roundtrip.{ext}")); + fs::write(&input_path, input).unwrap(); + + for (from, to) in [(&input_path, &cbq_path), (&cbq_path, &output_path)] { + let mut cmd = cargo::cargo_bin_cmd!("deacon"); + cmd.args(["filter", "-d", "-a", "65535", "-r", "1", "-t", "1"]); + if fasta { + cmd.arg("--discard-quality"); + } + cmd.arg(&bin_path) + .arg(from) + .arg("--output") + .arg(to) + .assert() + .success(); + } + assert_eq!(fs::read_to_string(output_path).unwrap(), expected); + } +} + +/// Write `n` distinct single-end FASTQ records that do not match the AAA index +fn create_many_fastq(path: &Path, n: usize) -> String { + let mut content = String::new(); + let seq = "CGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT"; + let qual = "~".repeat(seq.len()); + for i in 0..n { + content.push_str(&format!("@read_{i}\n{seq}\n+\n{qual}\n")); + } + fs::write(path, &content).unwrap(); + content +} + +/// --check-pairs accepts a natively paired CBQ with /1 /2 names +#[test] +fn cbq_check_pairs_paired_input() { + let temp_dir = tempdir().unwrap(); + let fasta_path = temp_dir.path().join("ref.fasta"); + let bin_path = temp_dir.path().join("ref.bin"); + let interleaved_path = temp_dir.path().join("interleaved.fastq"); + create_test_fasta(&fasta_path); + create_test_interleaved_fastq(&interleaved_path); + build_index(&fasta_path, &bin_path); + + // Interleaved FASTQ with read1/1, read1/2 names -> paired CBQ. + for suffix in ["cbq", "cba"] { + let paired_cbq = temp_dir.path().join(format!("paired.{suffix}")); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1", "--interleaved"]) + .arg(&bin_path) + .arg(&interleaved_path) + .arg("--output") + .arg(&paired_cbq) + .assert() + .success(); + + // Paired CBQ + --check-pairs succeeds despite no INPUT2/--interleaved + let out = temp_dir.path().join(format!("out-{suffix}.fastx")); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1", "--check-pairs"]) + .arg(&bin_path) + .arg(&paired_cbq) + .arg("--output") + .arg(&out) + .assert() + .success(); + let content = fs::read_to_string(&out).unwrap(); + if suffix == "cbq" { + assert_eq!(count_records(&content), 4); + } else { + // .cba dropped quality, so it decodes to FASTA + assert_eq!(content.matches('>').count(), 4); + } + } + + // Input format is detected from content, not the suffix. + let disguised = temp_dir.path().join("paired.dat"); + fs::copy(temp_dir.path().join("paired.cbq"), &disguised).unwrap(); + let out = temp_dir.path().join("out-dat.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1", "--check-pairs"]) + .arg(&bin_path) + .arg(&disguised) + .arg("--output") + .arg(&out) + .assert() + .success(); + assert_eq!(count_records(&fs::read_to_string(out).unwrap()), 4); +} + +/// A `.cba` output path writes a CBQ without qualities, implying --discard-quality +#[test] +fn cba_output_is_quality_free_cbq() { + let temp_dir = tempdir().unwrap(); + let fasta_path = temp_dir.path().join("ref.fasta"); + let bin_path = temp_dir.path().join("ref.bin"); + let fastq_path = temp_dir.path().join("reads.fastq"); + let r1_path = temp_dir.path().join("reads_1.fastq"); + let r2_path = temp_dir.path().join("reads_2.fastq"); + create_test_fasta(&fasta_path); + create_test_fastq(&fastq_path); + create_test_paired_fastq(&r1_path, &r2_path); + build_index(&fasta_path, &bin_path); + + // FASTQ -> .cba, without --discard-quality + let cba_path = temp_dir.path().join("out.cba"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&fastq_path) + .arg("--output") + .arg(&cba_path) + .assert() + .success(); + assert!( + fs::read(&cba_path).unwrap().starts_with(b"CBQFILE"), + ".cba output must be a CBQ file" + ); + let reader = binseq::cbq::MmapReader::new(&cba_path).unwrap(); + assert!( + !reader.header().has_qualities(), + ".cba must imply --discard-quality" + ); + assert!(reader.header().has_headers(), ".cba must keep headers"); + assert_eq!(reader.num_records(), 2); + + // ... and is byte-identical to the same run written as .cbq with the flag + let cbq_path = temp_dir.path().join("out.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args([ + "filter", + "-a", + "1", + "-r", + "0.0", + "-t", + "1", + "--discard-quality", + ]) + .arg(&bin_path) + .arg(&fastq_path) + .arg("--output") + .arg(&cbq_path) + .assert() + .success(); + assert_eq!( + fs::read(&cba_path).unwrap(), + fs::read(&cbq_path).unwrap(), + ".cba must match .cbq written with --discard-quality" + ); + + // A .cba reads back as CBQ input (magic sniffed) and decodes to FASTA + let roundtrip = temp_dir.path().join("roundtrip.fastx"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&cba_path) + .arg("--output") + .arg(&roundtrip) + .assert() + .success(); + let content = fs::read_to_string(&roundtrip).unwrap(); + assert!( + content.starts_with('>'), + "quality-free CBQ decodes to FASTA" + ); + assert_eq!(content.matches('>').count(), 2); + + // Paired input stays natively paired in a .cba + let paired_cba = temp_dir.path().join("paired.cba"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&r1_path) + .arg(&r2_path) + .arg("--output") + .arg(&paired_cba) + .assert() + .success(); + let reader = binseq::cbq::MmapReader::new(&paired_cba).unwrap(); + assert!(reader.is_paired(), "paired input must stay paired in .cba"); + assert!(!reader.header().has_qualities()); + assert_eq!(reader.num_records(), 2); +} + +/// A file bearing the CBQ magic but truncated must error, not panic +#[test] +fn cbq_truncated_input_fails_cleanly() { + let temp_dir = tempdir().unwrap(); + let fasta_path = temp_dir.path().join("ref.fasta"); + let bin_path = temp_dir.path().join("ref.bin"); + let fastq_path = temp_dir.path().join("reads.fastq"); + create_test_fasta(&fasta_path); + create_test_fastq(&fastq_path); + build_index(&fasta_path, &bin_path); + + // Valid CBQ, then truncate it to just past its magic + let cbq_path = temp_dir.path().join("reads.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&fastq_path) + .arg("--output") + .arg(&cbq_path) + .assert() + .success(); + let bytes = fs::read(&cbq_path).unwrap(); + let truncated = temp_dir.path().join("truncated.cbq"); + fs::write(&truncated, &bytes[..16]).unwrap(); + + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-t", "1"]) + .arg(&bin_path) + .arg(&truncated) + .assert() + .failure() + .stderr(predicate::str::contains("Truncated or corrupt CBQ")); +} + +/// CBQ per-record flags survive CBQ -> CBQ filtering +#[test] +fn cbq_flags_preserved() { + use binseq::write::{BinseqWriterBuilder, Format}; + use binseq::{ParallelReader, SequencingRecordBuilder}; + use std::sync::{Arc, Mutex}; + + let temp_dir = tempdir().unwrap(); + let aaa_path = temp_dir.path().join("aaa.fasta"); + let aaa_bin = temp_dir.path().join("aaa.bin"); + create_test_fasta_aaa(&aaa_path); + build_index(&aaa_path, &aaa_bin); + + // CBQ input with flags 7 and 9 + let flagged = temp_dir.path().join("flagged.cbq"); + let mut writer = BinseqWriterBuilder::new(Format::Cbq) + .headers(true) + .flags(true) + .build(File::create(&flagged).unwrap()) + .unwrap(); + for (i, flag) in [7u64, 9u64].iter().enumerate() { + writer + .push( + SequencingRecordBuilder::default() + .s_seq(b"CGTACGTACGTACGTACGTACGTACGTACGTACGTACGT".as_slice()) + .s_header(format!("read_{i}").as_bytes()) + .opt_flag(Some(*flag)) + .build() + .unwrap(), + ) + .unwrap(); + } + writer.finish().unwrap(); + + // Deplete against the AAA index keeps everything + let out_cbq = temp_dir.path().join("out.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "-t", "1"]) + .arg(&aaa_bin) + .arg(&flagged) + .arg("--output") + .arg(&out_cbq) + .assert() + .success(); + + // Collect the flags back out of the filtered CBQ + #[derive(Clone)] + struct FlagCollector(Arc>>>); + impl binseq::ParallelProcessor for FlagCollector { + fn process_record(&mut self, record: R) -> binseq::Result<()> { + self.0.lock().unwrap().push(record.flag()); + Ok(()) + } + } + let reader = binseq::cbq::MmapReader::new(&out_cbq).unwrap(); + assert!(reader.header().has_flags(), "flags must survive filtering"); + let collector = FlagCollector(Arc::new(Mutex::new(Vec::new()))); + reader.process_parallel(collector.clone(), 1).unwrap(); + assert_eq!(*collector.0.lock().unwrap(), vec![Some(7), Some(9)]); +} + +/// --ordered CBQ output matches the input order exactly, multithreaded +#[test] +fn cbq_ordered_output_is_input_ordered() { + let temp_dir = tempdir().unwrap(); + let aaa_path = temp_dir.path().join("aaa.fasta"); + let aaa_bin = temp_dir.path().join("aaa.bin"); + create_test_fasta_aaa(&aaa_path); + build_index(&aaa_path, &aaa_bin); + let fastq_path = temp_dir.path().join("many.fastq"); + let content = create_many_fastq(&fastq_path, 500); + + // Deplete against the AAA index keeps all 500 reads + let cbq_path = temp_dir.path().join("ordered.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "--ordered", "-t", "4"]) + .arg(&aaa_bin) + .arg(&fastq_path) + .arg("--output") + .arg(&cbq_path) + .assert() + .success(); + + // Decoding the CBQ reproduces the input byte-for-byte (CBQ input with + // --ordered falls back to one thread) + let out = temp_dir.path().join("out.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "--ordered", "-t", "4"]) + .arg(&aaa_bin) + .arg(&cbq_path) + .arg("--output") + .arg(&out) + .assert() + .success(); + assert_eq!( + fs::read_to_string(&out).unwrap(), + content, + "ordered CBQ roundtrip must preserve input order" + ); +} + +/// Write reads that complete CBQ blocks mid-batch. +fn create_many_long_fastq(path: &Path, n: usize, len: usize) -> String { + let mut content = String::new(); + let seq = "CGTA".repeat(len / 4); + let qual = "~".repeat(seq.len()); + for i in 0..n { + content.push_str(&format!("@read_{i:06}\n{seq}\n+\n{qual}\n")); + } + fs::write(path, &content).unwrap(); + content +} + +/// --ordered CBQ output stays ordered across completed blocks. +#[test] +fn cbq_ordered_block_completes_mid_batch() { + let temp_dir = tempdir().unwrap(); + let aaa_path = temp_dir.path().join("aaa.fasta"); + let aaa_bin = temp_dir.path().join("aaa.bin"); + create_test_fasta_aaa(&aaa_path); + build_index(&aaa_path, &aaa_bin); + + // Force each batch across multiple 1 MiB blocks. + let fastq_path = temp_dir.path().join("long.fastq"); + let content = create_many_long_fastq(&fastq_path, 2048, 1000); + + let cbq_path = temp_dir.path().join("ordered.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args([ + "filter", + "-d", + "--ordered", + "-t", + "4", + "--cbq-block-size", + "1", + ]) + .arg(&aaa_bin) + .arg(&fastq_path) + .arg("--output") + .arg(&cbq_path) + .assert() + .success(); + + let out = temp_dir.path().join("out.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "--ordered", "-t", "1"]) + .arg(&aaa_bin) + .arg(&cbq_path) + .arg("--output") + .arg(&out) + .assert() + .success(); + assert_eq!( + fs::read_to_string(&out).unwrap(), + content, + "multi-block CBQ output must preserve input order" + ); +} + +/// --ordered --rename CBQ output numbers records 1..N in input order +#[test] +fn cbq_ordered_rename_is_sequential() { + let temp_dir = tempdir().unwrap(); + let aaa_path = temp_dir.path().join("aaa.fasta"); + let aaa_bin = temp_dir.path().join("aaa.bin"); + create_test_fasta_aaa(&aaa_path); + build_index(&aaa_path, &aaa_bin); + let fastq_path = temp_dir.path().join("many.fastq"); + create_many_fastq(&fastq_path, 200); + + let cbq_path = temp_dir.path().join("renamed.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "--ordered", "--rename", "-t", "4"]) + .arg(&aaa_bin) + .arg(&fastq_path) + .arg("--output") + .arg(&cbq_path) + .assert() + .success(); + + let out = temp_dir.path().join("out.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "-t", "1"]) + .arg(&aaa_bin) + .arg(&cbq_path) + .arg("--output") + .arg(&out) + .assert() + .success(); + let content = fs::read_to_string(&out).unwrap(); + let ids: Vec<&str> = content + .lines() + .step_by(4) + .map(|h| h.trim_start_matches('@')) + .collect(); + let expected: Vec = (1..=200).map(|i| i.to_string()).collect(); + assert_eq!( + ids, expected, + "renamed CBQ records must be numbered in order" + ); +} + +/// Multithreaded CBQ roundtrips preserve record counts +#[test] +fn cbq_multithreaded_roundtrip_counts() { + let temp_dir = tempdir().unwrap(); + let aaa_path = temp_dir.path().join("aaa.fasta"); + let aaa_bin = temp_dir.path().join("aaa.bin"); + create_test_fasta_aaa(&aaa_path); + build_index(&aaa_path, &aaa_bin); + let fastq_path = temp_dir.path().join("many.fastq"); + create_many_fastq(&fastq_path, 500); + + let cbq_path = temp_dir.path().join("many.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "-t", "4"]) + .arg(&aaa_bin) + .arg(&fastq_path) + .arg("--output") + .arg(&cbq_path) + .assert() + .success(); + let reader = binseq::cbq::MmapReader::new(&cbq_path).unwrap(); + assert_eq!(reader.num_records(), 500); + + let out = temp_dir.path().join("out.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "-t", "4"]) + .arg(&aaa_bin) + .arg(&cbq_path) + .arg("--output") + .arg(&out) + .assert() + .success(); + assert_eq!(count_records(&fs::read_to_string(&out).unwrap()), 500); +} + +/// OUTPUT2 must not silently write FASTX bytes into a .cbq-named file +#[test] +fn cbq_output2_suffix_rejected() { + let temp_dir = tempdir().unwrap(); + let fasta_path = temp_dir.path().join("ref.fasta"); + let bin_path = temp_dir.path().join("ref.bin"); + let r1_path = temp_dir.path().join("reads_1.fastq"); + let r2_path = temp_dir.path().join("reads_2.fastq"); + create_test_fasta(&fasta_path); + create_test_paired_fastq(&r1_path, &r2_path); + build_index(&fasta_path, &bin_path); + + for suffix in ["cbq", "cba"] { + let out1 = temp_dir.path().join(format!("out1-{suffix}.fastq")); + let out2 = temp_dir.path().join(format!("out2.{suffix}")); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-a", "1", "-r", "0.0", "-t", "1"]) + .arg(&bin_path) + .arg(&r1_path) + .arg(&r2_path) + .arg("--output") + .arg(&out1) + .arg("--output2") + .arg(&out2) + .assert() + .failure(); + assert!(!out1.exists(), "OUTPUT must not be created"); + assert!(!out2.exists(), "OUTPUT2 must not be created"); + } +} + +/// Headerless, quality-free paired CBQ filters to CBQ (still headerless) and +/// to FASTX (record-index ids, FASTA since there are no qualities) +#[test] +fn cbq_headerless_roundtrip() { + let temp_dir = tempdir().unwrap(); + let aaa_path = temp_dir.path().join("aaa.fasta"); + let aaa_bin = temp_dir.path().join("aaa.bin"); + create_test_fasta_aaa(&aaa_path); + build_index(&aaa_path, &aaa_bin); + let headerless = temp_dir.path().join("headerless.cbq"); + write_headerless_paired_cbq(&headerless); + + // Deplete against the AAA index keeps both pairs + let out_cbq = temp_dir.path().join("out.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "-t", "1"]) + .arg(&aaa_bin) + .arg(&headerless) + .arg("--output") + .arg(&out_cbq) + .assert() + .success(); + let reader = binseq::cbq::MmapReader::new(&out_cbq).unwrap(); + assert!(reader.is_paired()); + assert!( + !reader.header().has_headers(), + "output must stay headerless" + ); + assert!(!reader.header().has_qualities()); + assert_eq!(reader.num_records(), 2); + + // Quality-free CBQ decodes to FASTA + let out_fasta = temp_dir.path().join("out.fasta"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "-t", "1"]) + .arg(&aaa_bin) + .arg(&headerless) + .arg("--output") + .arg(&out_fasta) + .assert() + .success(); + let content = fs::read_to_string(&out_fasta).unwrap(); + assert_eq!( + content.matches('>').count(), + 4, + "two pairs, four FASTA records" + ); +} + +/// Ultra-long reads exceed binseq's 1 MiB default block; deacon's larger CBQ +/// block must roundtrip them +#[test] +fn cbq_ultralong_read_roundtrip() { + let temp_dir = tempdir().unwrap(); + let aaa_path = temp_dir.path().join("aaa.fasta"); + let aaa_bin = temp_dir.path().join("aaa.bin"); + create_test_fasta_aaa(&aaa_path); + build_index(&aaa_path, &aaa_bin); + + // 2 Mb read: ~2.5 MB embedded size, over the 1 MiB binseq default + let fastq_path = temp_dir.path().join("ultra.fastq"); + let seq: String = "ACGTGCTA".repeat(250_000); + let qual = "~".repeat(seq.len()); + let content = format!("@ultra1\n{seq}\n+\n{qual}\n"); + fs::write(&fastq_path, &content).unwrap(); + + let cbq_path = temp_dir.path().join("ultra.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "-t", "1"]) + .arg(&aaa_bin) + .arg(&fastq_path) + .arg("--output") + .arg(&cbq_path) + .assert() + .success(); + + let out = temp_dir.path().join("ultra_out.fastq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-d", "-t", "1"]) + .arg(&aaa_bin) + .arg(&cbq_path) + .arg("--output") + .arg(&out) + .assert() + .success(); + assert_eq!( + fs::read_to_string(&out).unwrap(), + content, + "ultra-long read must roundtrip through CBQ" + ); +} + +/// CBQ output validates the compression level before creating the output +#[test] +fn cbq_invalid_compression_level_rejected() { + let temp_dir = tempdir().unwrap(); + let fasta_path = temp_dir.path().join("ref.fasta"); + let bin_path = temp_dir.path().join("ref.bin"); + let fastq_path = temp_dir.path().join("reads.fastq"); + create_test_fasta(&fasta_path); + create_test_fastq(&fastq_path); + build_index(&fasta_path, &bin_path); + + let no_create = temp_dir.path().join("bad_level.cbq"); + cargo::cargo_bin_cmd!("deacon") + .args(["filter", "-t", "1", "--compression-level", "23"]) + .arg(&bin_path) + .arg(&fastq_path) + .arg("--output") + .arg(&no_create) + .assert() + .failure() + .stderr(predicate::str::contains("Invalid CBQ compression level")); + assert!(!no_create.exists()); +}