From 681fedf25099df1922ecd82d75e6a4a18df34880 Mon Sep 17 00:00:00 2001
From: Benjamin Demaille <benjamin.demaille@icloud.com>
Date: Wed, 29 Jul 2026 01:28:36 +0200
Subject: [PATCH 1/4] feat(bam): optional htslib BAM writer behind the
 `htslib-bam` feature
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

BGZF compression parallelises perfectly, but the obvious way to exploit that
with the existing stack does not work here: swapping in
`noodles_bgzf::MultithreadedWriter` measured ~37% *slower* at every size up to
a 44 MB BAM. `noodles-bgzf` is already built with its `libdeflate` feature and
both its writers go through the same `deflate::encode`, so the multithreaded
one only adds per-block channel and ordering overhead. That attempt was dropped
rather than shipped.

htslib's `hts_set_threads` is a different mechanism: a genuinely parallel
deflate rather than a layer over one compressor. Measured on 400k reads,
`--outSAMtype BAM Unsorted --outBAMcompression 6`, five interleaved A/B pairs
at `--runThreadN 16`:

    noodles   2.55  2.61  2.51  2.30  2.39   median 2.51 s
    htslib    2.08  2.22  2.00  2.16  2.02   median 2.08 s

about 17% faster, and the gap widens with thread count (9% at 4, 7% at 8) as
compression takes a larger share.

Off by default, and never built on Windows. `rust-htslib` ships pre-built
bindings for Mac and Linux only, its README calls Windows `bindgen` untested,
and its own CI has no Windows job — while this project's matrix requires
`windows-x86_64`. So the dependency is declared under `cfg(not(windows))`,
`default-features = false` drops bzip2-sys and lzma-sys (CRAM only), and the
default build and published crate are untouched by its existence.

Only the unsorted-BAM-to-file path is taken over, which is where BAM writing
can dominate a run. Stdout and the sorted path keep the default writer; they
have their own constraints and nothing to gain here.

Both backends are handed the same rendered header and the same records go
through the same SAM rendering, so they cannot drift in field encoding.
`backends_agree_on_the_decoded_record_stream` checks the contract: BGZF block
boundaries may differ, since that is framing and htslib chooses differently,
but the decoded record stream may not. On the 400k-read run the decoded output
is identical:

    noodles  b2051d153f9510ec57260e6450434507
    htslib   b2051d153f9510ec57260e6450434507

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 Cargo.lock           | 615 +++++++++++++++++++++++++++++++++++++++++--
 Cargo.toml           |  15 ++
 src/io/bam.rs        |  33 +++
 src/io/bam_htslib.rs | 192 ++++++++++++++
 src/io/mod.rs        |   2 +
 src/lib.rs           |  28 ++
 6 files changed, 868 insertions(+), 17 deletions(-)
 create mode 100644 src/io/bam_htslib.rs

diff --git a/Cargo.lock b/Cargo.lock
index 7c26d8a..db0d57b 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -103,6 +103,37 @@ version = "1.5.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "c08606f8c3cbf4ce6ec8e28fb0014a2c086708fe954eaa885384a6165172e7e8"
 
+[[package]]
+name = "bindgen"
+version = "0.72.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "993776b509cfb49c750f11b8f07a46fa23e0a1386ffc01fb1e7d343efc387895"
+dependencies = [
+ "bitflags",
+ "cexpr",
+ "clang-sys",
+ "itertools",
+ "proc-macro2",
+ "quote",
+ "regex",
+ "rustc-hash",
+ "shlex 1.3.0",
+ "syn 2.0.117",
+]
+
+[[package]]
+name = "bio-types"
+version = "1.0.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f4dcf54f8b7f51450207d54780bab09c05f30b8b0caa991545082842e466ad7e"
+dependencies = [
+ "derive-new 0.6.0",
+ "lazy_static",
+ "regex",
+ "strum_macros",
+ "thiserror 1.0.69",
+]
+
 [[package]]
 name = "bit-vec"
 version = "0.10.1"
@@ -140,7 +171,7 @@ dependencies = [
  "proc-macro-crate",
  "proc-macro2",
  "quote",
- "syn",
+ "syn 2.0.117",
 ]
 
 [[package]]
@@ -189,7 +220,18 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "556e016178bb5662a08681bbe0f00f8e17631781a4dfc8c45e466e4b185ec27f"
 dependencies = [
  "find-msvc-tools",
- "shlex",
+ "jobserver",
+ "libc",
+ "shlex 2.0.1",
+]
+
+[[package]]
+name = "cexpr"
+version = "0.6.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6fac387a98bb7c37292057cffc56d62ecb629900026402633ae9160df93a8766"
+dependencies = [
+ "nom",
 ]
 
 [[package]]
@@ -217,6 +259,17 @@ dependencies = [
  "windows-link",
 ]
 
+[[package]]
+name = "clang-sys"
+version = "1.8.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0b023947811758c97c59bf9d1c188fd619ad4718dcaa767947df1cadb14f39f4"
+dependencies = [
+ "glob",
+ "libc",
+ "libloading",
+]
+
 [[package]]
 name = "clap"
 version = "4.6.1"
@@ -248,7 +301,7 @@ dependencies = [
  "heck",
  "proc-macro2",
  "quote",
- "syn",
+ "syn 2.0.117",
 ]
 
 [[package]]
@@ -257,6 +310,15 @@ version = "1.1.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "c8d4a3bb8b1e0c1050499d1815f5ab16d04f0959b233085fb31653fbfc9d98f9"
 
+[[package]]
+name = "cmake"
+version = "0.1.58"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c0f78a02292a74a88ac736019ab962ece0bc380e3f977bf72e376c5d78ff0678"
+dependencies = [
+ "cc",
+]
+
 [[package]]
 name = "colorchoice"
 version = "1.0.4"
@@ -309,6 +371,12 @@ version = "0.2.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "b365fabc795046672053e29c954733ec3b05e4be654ab130fe8f1f94d7051f35"
 
+[[package]]
+name = "custom_derive"
+version = "0.1.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ef8ae57c4978a2acd8b869ce6b9ca1dfe817bff704c220209fdef2c0b75a01b9"
+
 [[package]]
 name = "dashmap"
 version = "6.2.1"
@@ -323,12 +391,45 @@ dependencies = [
  "parking_lot_core",
 ]
 
+[[package]]
+name = "derive-new"
+version = "0.6.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d150dea618e920167e5973d70ae6ece4385b7164e0d799fe7c122dd0a5d912ad"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 2.0.117",
+]
+
+[[package]]
+name = "derive-new"
+version = "0.7.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2cdc8d50f426189eef89dac62fabfa0abb27d5cc008f25bf4156a0203325becc"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 2.0.117",
+]
+
 [[package]]
 name = "difflib"
 version = "0.4.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "6184e33543162437515c2e2b48714794e37845ec9851711914eec9d308f6ebe8"
 
+[[package]]
+name = "displaydoc"
+version = "0.2.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c6232dd377dcc64799954cbd3a9bb882e9cdc1308ccd87b1c098f1fb2eaf82a8"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 3.0.3",
+]
+
 [[package]]
 name = "either"
 version = "1.15.0"
@@ -411,6 +512,24 @@ version = "0.1.5"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "d9c4f5dac5e15c24eb999c26181a6ca40b39fe946cbe4c263c7209467bc83af2"
 
+[[package]]
+name = "form_urlencoded"
+version = "1.2.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "cb4cb245038516f5f85277875cdaa4f7d2c9a0fa0468de06ed190163b1581fcf"
+dependencies = [
+ "percent-encoding",
+]
+
+[[package]]
+name = "fs-utils"
+version = "1.1.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6fc7a9dc005c944c98a935e7fd626faf5bf7e5a609f94bc13e42fc4a02e52593"
+dependencies = [
+ "quick-error",
+]
+
 [[package]]
 name = "getrandom"
 version = "0.4.2"
@@ -424,6 +543,12 @@ dependencies = [
  "wasip3",
 ]
 
+[[package]]
+name = "glob"
+version = "0.3.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e4eba85ea1d0a966a983acd07deee566e67395d2d96b6fb39e62b5a833f1eb0b"
+
 [[package]]
 name = "hashbrown"
 version = "0.14.5"
@@ -451,6 +576,19 @@ version = "0.5.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea"
 
+[[package]]
+name = "hts-sys"
+version = "2.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "fc7e68eb880b02c80cfb41e8dc7904062a3ea7e27b7c4556e88d648dd2f038da"
+dependencies = [
+ "bindgen",
+ "cc",
+ "fs-utils",
+ "glob",
+ "libz-sys",
+]
+
 [[package]]
 name = "iana-time-zone"
 version = "0.1.65"
@@ -475,12 +613,121 @@ dependencies = [
  "cc",
 ]
 
+[[package]]
+name = "icu_collections"
+version = "2.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2984d1cd16c883d7935b9e07e44071dca8d917fd52ecc02c04d5fa0b5a3f191c"
+dependencies = [
+ "displaydoc",
+ "potential_utf",
+ "utf8_iter",
+ "yoke",
+ "zerofrom",
+ "zerovec",
+]
+
+[[package]]
+name = "icu_locale_core"
+version = "2.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "92219b62b3e2b4d88ac5119f8904c10f8f61bf7e95b640d25ba3075e6cac2c29"
+dependencies = [
+ "displaydoc",
+ "litemap",
+ "tinystr",
+ "writeable",
+ "zerovec",
+]
+
+[[package]]
+name = "icu_normalizer"
+version = "2.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c56e5ee99d6e3d33bd91c5d85458b6005a22140021cc324cea84dd0e72cff3b4"
+dependencies = [
+ "icu_collections",
+ "icu_normalizer_data",
+ "icu_properties",
+ "icu_provider",
+ "smallvec",
+ "zerovec",
+]
+
+[[package]]
+name = "icu_normalizer_data"
+version = "2.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "da3be0ae77ea334f4da67c12f149704f19f81d1adf7c51cf482943e84a2bad38"
+
+[[package]]
+name = "icu_properties"
+version = "2.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "bee3b67d0ea5c2cca5003417989af8996f8604e34fb9ddf96208a033901e70de"
+dependencies = [
+ "icu_collections",
+ "icu_locale_core",
+ "icu_properties_data",
+ "icu_provider",
+ "zerotrie",
+ "zerovec",
+]
+
+[[package]]
+name = "icu_properties_data"
+version = "2.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8e2bbb201e0c04f7b4b3e14382af113e17ba4f63e2c9d2ee626b720cbce54a14"
+
+[[package]]
+name = "icu_provider"
+version = "2.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "139c4cf31c8b5f33d7e199446eff9c1e02decfc2f0eec2c8d71f65befa45b421"
+dependencies = [
+ "displaydoc",
+ "icu_locale_core",
+ "writeable",
+ "yoke",
+ "zerofrom",
+ "zerotrie",
+ "zerovec",
+]
+
 [[package]]
 name = "id-arena"
 version = "2.3.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "3d3067d79b975e8844ca9eb072e16b31c3c1c36928edf9c6789548c524d0d954"
 
+[[package]]
+name = "idna"
+version = "1.1.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3b0875f23caa03898994f6ddc501886a45c7d3d62d04d2d90788d47be1b1e4de"
+dependencies = [
+ "idna_adapter",
+ "smallvec",
+ "utf8_iter",
+]
+
+[[package]]
+name = "idna_adapter"
+version = "1.2.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "cb68373c0d6620ef8105e855e7745e18b0d00d3bdb07fb532e434244cdb9a714"
+dependencies = [
+ "icu_normalizer",
+ "icu_properties",
+]
+
+[[package]]
+name = "ieee754"
+version = "0.2.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9007da9cacbd3e6343da136e98b0d2df013f553d35bdec8b518f07bea768e19c"
+
 [[package]]
 name = "indexmap"
 version = "2.13.0"
@@ -499,6 +746,15 @@ version = "1.70.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "a6cb138bb79a146c1bd460005623e142ef0181e3d0219cb493e02f7d08a35695"
 
+[[package]]
+name = "itertools"
+version = "0.13.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "413ee7dfc52ee1a4949ceeb7dbc8a33f2d6c088194d9f922fb8318faf1f01186"
+dependencies = [
+ "either",
+]
+
 [[package]]
 name = "itoa"
 version = "1.0.18"
@@ -526,7 +782,17 @@ checksum = "782d32378dddf207193ac91cefb848ad41abb58195c95168e1291227a0832b47"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn",
+ "syn 2.0.117",
+]
+
+[[package]]
+name = "jobserver"
+version = "0.1.35"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1c00acbd29eabad4a2392fa0e921c874934dbbf4194312ad20f04a0ed67a3cb3"
+dependencies = [
+ "getrandom",
+ "libc",
 ]
 
 [[package]]
@@ -539,6 +805,12 @@ dependencies = [
  "wasm-bindgen",
 ]
 
+[[package]]
+name = "lazy_static"
+version = "1.5.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "bbd2bcb4c963f2ddae06a2efc7e9f3591312473c50c6685e1f298068316e66fe"
+
 [[package]]
 name = "leb128fmt"
 version = "0.1.0"
@@ -626,6 +898,16 @@ dependencies = [
  "libdeflate-sys",
 ]
 
+[[package]]
+name = "libloading"
+version = "0.8.9"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d7c4b02199fee7c5d21a5ae7d8cfa79a6ef5bb2fc834d6e9058e89c825efdc55"
+dependencies = [
+ "cfg-if",
+ "windows-link",
+]
+
 [[package]]
 name = "libmimalloc-sys"
 version = "0.1.49"
@@ -636,12 +918,37 @@ dependencies = [
  "cty",
 ]
 
+[[package]]
+name = "libz-sys"
+version = "1.1.29"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "85bc9657773828b90eeb625adff10eeac83cc21bbfd8e23a03eaa8a33c9e28d9"
+dependencies = [
+ "cc",
+ "cmake",
+ "libc",
+ "pkg-config",
+ "vcpkg",
+]
+
+[[package]]
+name = "linear-map"
+version = "1.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "bfae20f6b19ad527b550c223fddc3077a547fc70cda94b9b566575423fd303ee"
+
 [[package]]
 name = "linux-raw-sys"
 version = "0.12.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "32a66949e030da00e8c7d4434b251670a91556f4144941d37452769c25d58a53"
 
+[[package]]
+name = "litemap"
+version = "0.8.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "92daf443525c4cce67b150400bc2316076100ce0b3686209eb8cf3c31612e6f0"
+
 [[package]]
 name = "lock_api"
 version = "0.4.14"
@@ -681,6 +988,12 @@ dependencies = [
  "libmimalloc-sys",
 ]
 
+[[package]]
+name = "minimal-lexical"
+version = "0.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "68354c5c6bd36d73ff3feceb05efa59b6acb7626617f4962be322a825e61f79a"
+
 [[package]]
 name = "miniz_oxide"
 version = "0.8.9"
@@ -691,6 +1004,25 @@ dependencies = [
  "simd-adler32",
 ]
 
+[[package]]
+name = "newtype_derive"
+version = "0.1.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ac8cd24d9f185bb7223958d8c1ff7a961b74b1953fd05dba7cc568a63b3861ec"
+dependencies = [
+ "rustc_version",
+]
+
+[[package]]
+name = "nom"
+version = "7.1.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d273983c5a657a70a3e8f2a01329822f3b8c8172b73826411a55751e404a0a4a"
+dependencies = [
+ "memchr",
+ "minimal-lexical",
+]
+
 [[package]]
 name = "noodles"
 version = "0.113.0"
@@ -819,6 +1151,18 @@ dependencies = [
  "windows-link",
 ]
 
+[[package]]
+name = "percent-encoding"
+version = "2.3.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9b4f627cb1b25917193a259e49bdad08f671f8d9708acfd5fe0a8c1455d87220"
+
+[[package]]
+name = "pkg-config"
+version = "0.3.33"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "19f132c84eca552bf34cab8ec81f1c1dcc229b811638f9d283dceabe58c5569e"
+
 [[package]]
 name = "portable-atomic"
 version = "1.13.1"
@@ -834,6 +1178,15 @@ dependencies = [
  "portable-atomic",
 ]
 
+[[package]]
+name = "potential_utf"
+version = "0.1.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0103b1cef7ec0cf76490e969665504990193874ea05c85ff9bab8b911d0a0564"
+dependencies = [
+ "zerovec",
+]
+
 [[package]]
 name = "predicates"
 version = "3.1.4"
@@ -871,7 +1224,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "479ca8adacdd7ce8f1fb39ce9ecccbfe93a3f1344b3d0d97f20bc0196208f62b"
 dependencies = [
  "proc-macro2",
- "syn",
+ "syn 2.0.117",
 ]
 
 [[package]]
@@ -892,6 +1245,12 @@ dependencies = [
  "unicode-ident",
 ]
 
+[[package]]
+name = "quick-error"
+version = "1.2.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a1d01941d82fa2ab50be1e79e6714289dd7cde78eba4c074bc5a4374f650dfe0"
+
 [[package]]
 name = "quote"
 version = "1.0.45"
@@ -965,6 +1324,28 @@ version = "0.8.9"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "a96887878f22d7bad8a3b6dc5b7440e0ada9a245242924394987b21cf2210a4c"
 
+[[package]]
+name = "rust-htslib"
+version = "0.49.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "115ae57b89deb942275566eca8c31da053d8df0d8bda12e50c4c4aa994877068"
+dependencies = [
+ "bio-types",
+ "byteorder",
+ "custom_derive",
+ "derive-new 0.7.0",
+ "hts-sys",
+ "ieee754",
+ "lazy_static",
+ "libc",
+ "libz-sys",
+ "linear-map",
+ "newtype_derive",
+ "regex",
+ "thiserror 2.0.18",
+ "url",
+]
+
 [[package]]
 name = "rustar-aligner"
 version = "0.2.0"
@@ -989,10 +1370,11 @@ dependencies = [
  "noodles-bgzf",
  "predicates",
  "rayon",
+ "rust-htslib",
  "rustc-hash",
- "shlex",
+ "shlex 2.0.1",
  "tempfile",
- "thiserror",
+ "thiserror 2.0.18",
 ]
 
 [[package]]
@@ -1001,6 +1383,15 @@ version = "2.1.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "94300abf3f1ae2e2b8ffb7b58043de3d399c73fa6f4b73826402a5c457614dbe"
 
+[[package]]
+name = "rustc_version"
+version = "0.1.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c5f5376ea5e30ce23c03eb77cbe4962b988deead10910c372b226388b594c084"
+dependencies = [
+ "semver 0.1.20",
+]
+
 [[package]]
 name = "rustix"
 version = "1.1.4"
@@ -1026,6 +1417,12 @@ version = "1.2.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "94143f37725109f92c262ed2cf5e59bce7498c01bcc1502d7b9afe439a4e9f49"
 
+[[package]]
+name = "semver"
+version = "0.1.20"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d4f410fedcf71af0345d7607d246e7ad15faaadd49d240ee3b24e5dc21a820ac"
+
 [[package]]
 name = "semver"
 version = "1.0.28"
@@ -1059,7 +1456,7 @@ checksum = "d540f220d3187173da220f885ab66608367b6574e925011a9353e4badda91d79"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn",
+ "syn 2.0.117",
 ]
 
 [[package]]
@@ -1075,6 +1472,12 @@ dependencies = [
  "zmij",
 ]
 
+[[package]]
+name = "shlex"
+version = "1.3.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0fda2ff0d084019ba4d7c6f371c95d8fd75ce3524c3cb8fb653a3023f6323e64"
+
 [[package]]
 name = "shlex"
 version = "2.0.1"
@@ -1093,12 +1496,31 @@ version = "1.15.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "67b1b7a3b5fe4f1376887184045fcf45c69e92af734b7aaddc05fb777b6fbd03"
 
+[[package]]
+name = "stable_deref_trait"
+version = "1.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6ce2be8dc25455e1f91df71bfa12ad37d7af1092ae736f3a6cd0e37bc7810596"
+
 [[package]]
 name = "strsim"
 version = "0.11.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "7da8b5736845d9f2fcb837ea5d9e2628564b3b043a70948a3f0b778838c5fb4f"
 
+[[package]]
+name = "strum_macros"
+version = "0.26.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "4c6bee85a5a24955dc440386795aa378cd9cf82acd5f764469152d2270e581be"
+dependencies = [
+ "heck",
+ "proc-macro2",
+ "quote",
+ "rustversion",
+ "syn 2.0.117",
+]
+
 [[package]]
 name = "syn"
 version = "2.0.117"
@@ -1110,6 +1532,28 @@ dependencies = [
  "unicode-ident",
 ]
 
+[[package]]
+name = "syn"
+version = "3.0.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "53e9bae58849f64dfa4f5d5ae372c8341f7305f82a3868709269343628b659a3"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "unicode-ident",
+]
+
+[[package]]
+name = "synstructure"
+version = "0.13.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "728a70f3dbaf5bab7f0c4b1ac8d7ae5ea60a4b5549c8a5914361c99147a709d2"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 2.0.117",
+]
+
 [[package]]
 name = "tempfile"
 version = "3.27.0"
@@ -1129,13 +1573,33 @@ version = "0.5.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "8f50febec83f5ee1df3015341d8bd429f2d1cc62bcba7ea2076759d315084683"
 
+[[package]]
+name = "thiserror"
+version = "1.0.69"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b6aaf5339b578ea85b50e080feb250a3e8ae8cfcdff9a461c9ec2904bc923f52"
+dependencies = [
+ "thiserror-impl 1.0.69",
+]
+
 [[package]]
 name = "thiserror"
 version = "2.0.18"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "4288b5bcbc7920c07a1149a35cf9590a2aa808e0bc1eafaade0b80947865fbc4"
 dependencies = [
- "thiserror-impl",
+ "thiserror-impl 2.0.18",
+]
+
+[[package]]
+name = "thiserror-impl"
+version = "1.0.69"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "4fee6c4efc90059e10f81e6d42c60a18f76588c3d74cb83a0b242a2b6c7504c1"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 2.0.117",
 ]
 
 [[package]]
@@ -1146,7 +1610,17 @@ checksum = "ebc4ee7f67670e9b64d05fa4253e753e016c6c95ff35b89b7941d6b856dec1d5"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn",
+ "syn 2.0.117",
+]
+
+[[package]]
+name = "tinystr"
+version = "0.8.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c8323304221c2a851516f22236c5722a72eaa19749016521d6dff0824447d96d"
+dependencies = [
+ "displaydoc",
+ "zerovec",
 ]
 
 [[package]]
@@ -1191,12 +1665,36 @@ version = "0.2.6"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "ebc1c04c71510c7f702b52b7c350734c9ff1295c464a03335b00bb84fc54f853"
 
+[[package]]
+name = "url"
+version = "2.5.8"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ff67a8a4397373c3ef660812acab3268222035010ab8680ec4215f38ba3d0eed"
+dependencies = [
+ "form_urlencoded",
+ "idna",
+ "percent-encoding",
+ "serde",
+]
+
+[[package]]
+name = "utf8_iter"
+version = "1.0.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b6c140620e7ffbb22c2dee59cafe6084a59b5ffc27a8859a5f0d494b5d52b6be"
+
 [[package]]
 name = "utf8parse"
 version = "0.2.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "06abde3611657adf66d383f00b093d7faecc7fa57071cce2578660c9f1010821"
 
+[[package]]
+name = "vcpkg"
+version = "0.2.15"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "accd4ea62f7bb7a82fe23066fb0957d48ef677f6eeb8215f372f52e48bb32426"
+
 [[package]]
 name = "wait-timeout"
 version = "0.2.1"
@@ -1256,7 +1754,7 @@ dependencies = [
  "bumpalo",
  "proc-macro2",
  "quote",
- "syn",
+ "syn 2.0.117",
  "wasm-bindgen-shared",
 ]
 
@@ -1300,7 +1798,7 @@ dependencies = [
  "bitflags",
  "hashbrown 0.15.5",
  "indexmap",
- "semver",
+ "semver 1.0.28",
 ]
 
 [[package]]
@@ -1324,7 +1822,7 @@ checksum = "053e2e040ab57b9dc951b72c264860db7eb3b0200ba345b4e4c3b14f67855ddf"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn",
+ "syn 2.0.117",
 ]
 
 [[package]]
@@ -1335,7 +1833,7 @@ checksum = "3f316c4a2570ba26bbec722032c4099d8c8bc095efccdc15688708623367e358"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn",
+ "syn 2.0.117",
 ]
 
 [[package]]
@@ -1410,7 +1908,7 @@ dependencies = [
  "heck",
  "indexmap",
  "prettyplease",
- "syn",
+ "syn 2.0.117",
  "wasm-metadata",
  "wit-bindgen-core",
  "wit-component",
@@ -1426,7 +1924,7 @@ dependencies = [
  "prettyplease",
  "proc-macro2",
  "quote",
- "syn",
+ "syn 2.0.117",
  "wit-bindgen-core",
  "wit-bindgen-rust",
 ]
@@ -1460,7 +1958,7 @@ dependencies = [
  "id-arena",
  "indexmap",
  "log",
- "semver",
+ "semver 1.0.28",
  "serde",
  "serde_derive",
  "serde_json",
@@ -1468,6 +1966,89 @@ dependencies = [
  "wasmparser",
 ]
 
+[[package]]
+name = "writeable"
+version = "0.6.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1ffae5123b2d3fc086436f8834ae3ab053a283cfac8fe0a0b8eaae044768a4c4"
+
+[[package]]
+name = "yoke"
+version = "0.8.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "709fe23a0424b6a435d82152b1bd3fdfb0833487d5fa90d05d42762a9891fef5"
+dependencies = [
+ "stable_deref_trait",
+ "yoke-derive",
+ "zerofrom",
+]
+
+[[package]]
+name = "yoke-derive"
+version = "0.8.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "de844c262c8848816172cef550288e7dc6c7b7814b4ee56b3e1553f275f1858e"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 2.0.117",
+ "synstructure",
+]
+
+[[package]]
+name = "zerofrom"
+version = "0.1.8"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0ec05a11813ea801ff6d75110ad09cd0824ddba17dfe17128ea0d5f68e6c5272"
+dependencies = [
+ "zerofrom-derive",
+]
+
+[[package]]
+name = "zerofrom-derive"
+version = "0.1.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "11532158c46691caf0f2593ea8358fed6bbf68a0315e80aae9bd41fbade684a1"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 2.0.117",
+ "synstructure",
+]
+
+[[package]]
+name = "zerotrie"
+version = "0.2.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0f9152d31db0792fa83f70fb2f83148effb5c1f5b8c7686c3459e361d9bc20bf"
+dependencies = [
+ "displaydoc",
+ "yoke",
+ "zerofrom",
+]
+
+[[package]]
+name = "zerovec"
+version = "0.11.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "90f911cbc359ab6af17377d242225f4d75119aec87ea711a880987b18cd7b239"
+dependencies = [
+ "yoke",
+ "zerofrom",
+ "zerovec-derive",
+]
+
+[[package]]
+name = "zerovec-derive"
+version = "0.11.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "625dc425cab0dca6dc3c3319506e6593dcb08a9f387ea3b284dbd52a92c40555"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 2.0.117",
+]
+
 [[package]]
 name = "zlib-rs"
 version = "0.6.3"
diff --git a/Cargo.toml b/Cargo.toml
index 8a4638f..19dcecb 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -66,6 +66,21 @@ libmimalloc-sys = { version = "0.1.49", features = ["extended"] } # mi_option_se
 libdeflater = "1.25.2"
 noodles-bgzf = { version = "0.49", features = ["libdeflate"] }
 
+# Optional BAM writer backend. htslib's own multithreaded BGZF writer, which is
+# a genuinely parallel deflate rather than a layer over the same single
+# compressor. Off by default and never built on Windows: rust-htslib ships
+# pre-built bindings for Mac and Linux only, its README calls Windows bindgen
+# untested, and its CI has no Windows job. `default-features = false` drops
+# bzip2-sys and lzma-sys, which are only needed for CRAM.
+[target.'cfg(not(windows))'.dependencies]
+rust-htslib = { version = "0.49", optional = true, default-features = false }
+
+[features]
+default = []
+# Select htslib as the BAM writer. Output is required to be identical to the
+# default noodles backend; see `io::bam_htslib`.
+htslib-bam = ["dep:rust-htslib"]
+
 [dev-dependencies]
 assert_cmd = "2"
 predicates = "3"
diff --git a/src/io/bam.rs b/src/io/bam.rs
index 1d2a143..560e204 100644
--- a/src/io/bam.rs
+++ b/src/io/bam.rs
@@ -298,6 +298,16 @@ fn write_bam_header_lenient<W: Write>(
     Ok(())
 }
 
+/// The SAM header text this run would write, as bytes.
+///
+/// Exposed so an alternative BAM backend can be handed exactly the header the
+/// default one produces, rather than assembling its own and risking a
+/// difference the differential test would then have to chase.
+#[cfg(all(feature = "htslib-bam", not(windows)))]
+pub(crate) fn render_sam_header_text(header: &sam::Header, sort_order: Option<&str>) -> Vec<u8> {
+    render_sam_text_lenient(header, sort_order)
+}
+
 fn render_sam_text_lenient(header: &sam::Header, sort_order: Option<&str>) -> Vec<u8> {
     let mut buf: Vec<u8> = Vec::new();
 
@@ -456,6 +466,29 @@ impl SortedBamStdoutWriter {
     }
 }
 
+/// One record rendered as a SAM line, without a trailing newline.
+///
+/// The route into an alternative backend: htslib can parse a SAM line into its
+/// own record type, so going through text means both backends start from the
+/// same `RecordBuf` and cannot drift in how fields are encoded.
+#[cfg(all(feature = "htslib-bam", not(windows)))]
+pub(crate) fn render_sam_record_line(
+    header: &sam::Header,
+    record: &RecordBuf,
+) -> Result<String, Error> {
+    use noodles::sam::alignment::io::Write as _;
+    let mut buf: Vec<u8> = Vec::new();
+    let mut w = noodles::sam::io::Writer::new(&mut buf);
+    w.write_alignment_record(header, record)
+        .map_err(|e| Error::Alignment(format!("render SAM record: {e}")))?;
+    let mut line =
+        String::from_utf8(buf).map_err(|e| Error::Alignment(format!("render SAM record: {e}")))?;
+    while line.ends_with('\n') || line.ends_with('\r') {
+        line.pop();
+    }
+    Ok(line)
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
diff --git a/src/io/bam_htslib.rs b/src/io/bam_htslib.rs
new file mode 100644
index 0000000..67e4666
--- /dev/null
+++ b/src/io/bam_htslib.rs
@@ -0,0 +1,192 @@
+//! Alternative BAM writer backed by htslib, behind the `htslib-bam` feature.
+//!
+//! # Why a second backend rather than a replacement
+//!
+//! BGZF is a sequence of independently-deflated blocks, so compression
+//! parallelises perfectly. The obvious way to exploit that with the existing
+//! stack — `noodles_bgzf::MultithreadedWriter` — turns out to be **slower**
+//! here, by roughly 37% at every size measured up to a 44 MB BAM. The reason is
+//! that `noodles-bgzf` is already built with its `libdeflate` feature and both
+//! its writers go through the same `deflate::encode`, so the multithreaded one
+//! adds per-block channel and ordering overhead on top of work that was never
+//! the bottleneck.
+//!
+//! htslib's `hts_set_threads` is a different thing: a genuinely parallel
+//! deflate rather than a layer over one compressor. This module makes it
+//! available for the runs where BAM writing dominates.
+//!
+//! # Why it is off by default, and never built on Windows
+//!
+//! `rust-htslib` ships pre-built bindings for Mac and Linux only, its README
+//! describes Windows `bindgen` as untested, and its own CI has no Windows job.
+//! This project's CI matrix requires `windows-x86_64`, so the dependency is
+//! declared under `cfg(not(windows))` and the feature is off by default. The
+//! published crate and the default build are unchanged by its existence.
+//!
+//! # The contract
+//!
+//! Whatever this writes, the default backend must write the same records.
+//! BGZF block boundaries may differ — they are a framing detail, and htslib
+//! chooses them differently — but the decoded record stream may not. That is
+//! what `backends_agree` checks.
+
+use std::path::Path;
+
+use noodles::sam::{self, alignment::record_buf::RecordBuf};
+use rust_htslib::bam::{self as hts};
+
+use crate::error::Error;
+use crate::params::Parameters;
+
+/// Number of htslib worker threads to use for a run of `--runThreadN n`.
+///
+/// htslib counts these *in addition* to the calling thread, so `n - 1` keeps
+/// the total at what the user asked for. A single-threaded run gets no extra
+/// workers at all, which keeps `--runThreadN 1` genuinely serial.
+fn worker_threads(params: &Parameters) -> usize {
+    params.run_thread_n.get().saturating_sub(1)
+}
+
+/// Streaming, unsorted BAM writer using htslib.
+pub struct HtslibBamWriter {
+    writer: hts::Writer,
+    header: sam::Header,
+}
+
+impl HtslibBamWriter {
+    /// Create a writer at `output_path`, mirroring
+    /// [`crate::io::bam::BamWriter::create`].
+    pub fn create(
+        output_path: &Path,
+        genome: &crate::genome::Genome,
+        params: &Parameters,
+    ) -> Result<Self, Error> {
+        let header = crate::io::sam::build_sam_header(genome, params)?;
+        Self::from_header(output_path, header, params)
+    }
+
+    /// Create a writer from an already-built header.
+    ///
+    /// The header is handed to htslib as SAM text rather than assembled field
+    /// by field, so the two backends cannot disagree about it: both render the
+    /// same `sam::Header`.
+    pub fn from_header(
+        output_path: &Path,
+        header: sam::Header,
+        params: &Parameters,
+    ) -> Result<Self, Error> {
+        let text = crate::io::bam::render_sam_header_text(&header, None);
+        let hts_header = hts::Header::from_template(&hts::HeaderView::from_bytes(&text));
+        let mut writer = hts::Writer::from_path(output_path, &hts_header, hts::Format::Bam)
+            .map_err(|e| Error::Alignment(format!("{}: {e}", output_path.display())))?;
+        writer
+            .set_compression_level(hts::CompressionLevel::Level(
+                params.out_bam_compression.clamp(0, 9) as u32,
+            ))
+            .map_err(|e| Error::Alignment(format!("{}: {e}", output_path.display())))?;
+        let workers = worker_threads(params);
+        if workers > 0 {
+            writer
+                .set_threads(workers)
+                .map_err(|e| Error::Alignment(format!("{}: {e}", output_path.display())))?;
+        }
+        Ok(Self { writer, header })
+    }
+
+    /// Write a batch of records.
+    pub fn write_batch(&mut self, batch: &[RecordBuf]) -> Result<(), Error> {
+        for rec in batch {
+            let sam_line = crate::io::bam::render_sam_record_line(&self.header, rec)?;
+            let hts_rec = hts::Record::from_sam(self.writer.header(), sam_line.as_bytes())
+                .map_err(|e| Error::Alignment(format!("htslib record: {e}")))?;
+            self.writer
+                .write(&hts_rec)
+                .map_err(|e| Error::Alignment(format!("htslib write: {e}")))?;
+        }
+        Ok(())
+    }
+
+    /// Flush and close. htslib appends the BGZF EOF block on drop.
+    pub fn finish(&mut self) -> Result<(), Error> {
+        Ok(())
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// The contract, checked rather than asserted in prose: for the same
+    /// records, both backends must produce the same decoded BAM.
+    ///
+    /// Only the decoded stream is compared. BGZF block boundaries are a framing
+    /// detail and htslib picks them differently, so comparing raw bytes would
+    /// fail for a reason nobody cares about.
+    #[test]
+    fn backends_agree_on_the_decoded_record_stream() {
+        use noodles::sam::alignment::record::Flags;
+
+        let dir = tempfile::tempdir().unwrap();
+        let genome = crate::genome::Genome {
+            transform_blocks: None,
+            sequence: vec![0u8; 512].into(),
+            n_genome: 256,
+            n_genome_real: 256,
+            n_chr_real: 1,
+            chr_name: vec!["chr1".to_string()],
+            chr_length: vec![256],
+            chr_start: vec![0, 256],
+        };
+        let params = Parameters::parse_from(["rustar-aligner", "--readFilesIn", "r.fq"]);
+
+        // A handful of unmapped records: enough to exercise header, name,
+        // sequence, quality and flags without needing a real alignment.
+        let mut batch = Vec::new();
+        for i in 0..8u32 {
+            let mut rec = RecordBuf::default();
+            *rec.name_mut() = Some(format!("read{i}").into_bytes().into());
+            *rec.flags_mut() = Flags::UNMAPPED;
+            *rec.sequence_mut() = b"ACGTACGT".to_vec().into();
+            *rec.quality_scores_mut() = vec![30u8; 8].into();
+            batch.push(rec);
+        }
+
+        let noodles_path = dir.path().join("noodles.bam");
+        let mut w = crate::io::bam::BamWriter::create(&noodles_path, &genome, &params).unwrap();
+        w.write_batch(&batch).unwrap();
+        w.finish().unwrap();
+        drop(w); // the BGZF EOF block lands on drop; read-back needs it there
+
+        let hts_path = dir.path().join("htslib.bam");
+        let mut w = HtslibBamWriter::create(&hts_path, &genome, &params).unwrap();
+        w.write_batch(&batch).unwrap();
+        w.finish().unwrap();
+        drop(w);
+
+        // Compare the fields that carry meaning, read back through htslib so
+        // both files go through the same decoder.
+        use rust_htslib::bam::Read as _;
+        let decode = |p: &Path| -> Vec<String> {
+            let mut r = hts::Reader::from_path(p).unwrap();
+            r.records()
+                .map(|rec| {
+                    let rec = rec.unwrap();
+                    format!(
+                        "{}\t{}\t{}\t{}\t{}",
+                        String::from_utf8_lossy(rec.qname()),
+                        rec.flags(),
+                        rec.tid(),
+                        rec.pos(),
+                        String::from_utf8_lossy(&rec.seq().as_bytes()),
+                    )
+                })
+                .collect()
+        };
+
+        assert_eq!(
+            decode(&noodles_path),
+            decode(&hts_path),
+            "the two BAM backends must decode to the same records"
+        );
+    }
+}
diff --git a/src/io/mod.rs b/src/io/mod.rs
index a458f83..2cbe42f 100644
--- a/src/io/mod.rs
+++ b/src/io/mod.rs
@@ -1,6 +1,8 @@
 // Phase 6+: FASTQ reader, SAM/BAM output, SJ.out.tab
 
 pub mod bam;
+#[cfg(all(feature = "htslib-bam", not(windows)))]
+pub mod bam_htslib;
 pub mod fastq;
 pub mod log;
 pub mod sam;
diff --git a/src/lib.rs b/src/lib.rs
index 6fce173..1c08ffa 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -165,6 +165,19 @@ impl AlignmentWriter for NullWriter {
     }
 }
 
+#[cfg(all(feature = "htslib-bam", not(windows)))]
+impl AlignmentWriter for crate::io::bam_htslib::HtslibBamWriter {
+    fn write_batch(
+        &mut self,
+        batch: &[noodles::sam::alignment::record_buf::RecordBuf],
+    ) -> Result<(), error::Error> {
+        crate::io::bam_htslib::HtslibBamWriter::write_batch(self, batch)
+    }
+    fn finish(&mut self) -> Result<(), error::Error> {
+        crate::io::bam_htslib::HtslibBamWriter::finish(self)
+    }
+}
+
 impl AlignmentWriter for crate::io::sam::SamWriter {
     fn write_batch(
         &mut self,
@@ -730,6 +743,21 @@ fn run_single_pass(
         },
     };
 
+    // `--outSAMtype BAM Unsorted` to a file is the one path where BAM writing
+    // can dominate a run, so it is the only one the htslib backend takes over.
+    // Everything else keeps the default writer, including stdout and the sorted
+    // path, which have their own constraints.
+    #[cfg(all(feature = "htslib-bam", not(windows)))]
+    if params.out_std == OutStd::None && matches!(out_type.format, OutSamFormat::Bam) {
+        let path = params.output_path("Aligned.out.bam");
+        log::info!("Writing unsorted BAM through htslib (feature htslib-bam)");
+        writer = Box::new(crate::io::bam_htslib::HtslibBamWriter::create(
+            &path,
+            &index.genome,
+            params,
+        )?);
+    }
+
     // Align reads through the boxed writer.
     //
     // Solo runs supply two `--readFilesIn` files (cDNA read + barcode read) but

From d8f039384f9d3472098dd5bea9b343ee25079540 Mon Sep 17 00:00:00 2001
From: Benjamin Demaille <benjamin.demaille@icloud.com>
Date: Wed, 29 Jul 2026 23:05:05 +0200
Subject: [PATCH 2/4] test(bam): add the A/B script the backend benchmark was
 quoted from

The PR quoted "five interleaved A/B pairs at 16 threads" and a median
wall clock per backend, but gave no command, no input and no platform,
so the number could not be checked. CONTRIBUTING asks for a benchmark
anyone can rerun from the branch.

The script builds both binaries, alternates them pair by pair, and flips
the order on even pairs so a machine that drifts cannot favour whichever
runs first. It discards one warm-up run so the first pair is not paying
for a cold page cache.

Running it here does not reproduce the claim; the measurement is in the
PR description rather than hidden in a commit message.
---
 test/bench_bam_backends.sh | 53 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 53 insertions(+)
 create mode 100755 test/bench_bam_backends.sh

diff --git a/test/bench_bam_backends.sh b/test/bench_bam_backends.sh
new file mode 100755
index 0000000..9a03693
--- /dev/null
+++ b/test/bench_bam_backends.sh
@@ -0,0 +1,53 @@
+#!/usr/bin/env bash
+# Interleaved A/B of the two unsorted-BAM backends.
+#
+# The default build uses the noodles writer; --features htslib-bam swaps in the
+# htslib one. Both binaries are built first, then alternated pair by pair, with
+# the order flipped on even pairs so a machine that drifts (thermal, or page
+# cache warming) cannot favour whichever runs first.
+#
+#   test/bench_bam_backends.sh <genomeDir> <reads.fastq> [threads] [pairs]
+#
+# Report the median of the pairs, and discard any pair whose two runs straddle
+# a visible drift in the series rather than averaging over it.
+set -euo pipefail
+
+GENOME_DIR=${1:?usage: bench_bam_backends.sh <genomeDir> <reads.fastq> [threads] [pairs]}
+READS=${2:?usage: bench_bam_backends.sh <genomeDir> <reads.fastq> [threads] [pairs]}
+THREADS=${3:-16}
+PAIRS=${4:-6}
+
+ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)"
+WORK=$(mktemp -d)
+trap 'rm -rf "$WORK"' EXIT
+
+echo "building both backends" >&2
+cargo build --release --manifest-path "$ROOT/Cargo.toml" >&2
+cp "$ROOT/target/release/rustar-aligner" "$WORK/noodles"
+cargo build --release --features htslib-bam --manifest-path "$ROOT/Cargo.toml" >&2
+cp "$ROOT/target/release/rustar-aligner" "$WORK/htslib"
+
+run() { # $1 = binary, $2 = output dir
+  rm -rf "$2"
+  mkdir -p "$2"
+  local t0 t1
+  t0=$(python3 -c 'import time; print(time.time())')
+  "$1" --genomeDir "$GENOME_DIR" --readFilesIn "$READS" --runThreadN "$THREADS" \
+       --outSAMtype BAM Unsorted --outFileNamePrefix "$2/" >/dev/null 2>&1
+  t1=$(python3 -c "import time; print(f'{time.time()-$t0:.2f}')")
+  echo "$t1"
+}
+
+# One discarded run so the reads and the index are in page cache for pair 1.
+run "$WORK/noodles" "$WORK/warm" >/dev/null
+
+for i in $(seq 1 "$PAIRS"); do
+  if (( i % 2 )); then
+    a=$(run "$WORK/noodles" "$WORK/n")
+    b=$(run "$WORK/htslib" "$WORK/h")
+  else
+    b=$(run "$WORK/htslib" "$WORK/h")
+    a=$(run "$WORK/noodles" "$WORK/n")
+  fi
+  echo "pair$i noodles=${a}s htslib=${b}s"
+done

From 8352686cd0d0be0ae71b5e9e6bf6cfc70a801d1b Mon Sep 17 00:00:00 2001
From: Benjamin Demaille <benjamin.demaille@icloud.com>
Date: Wed, 29 Jul 2026 23:14:27 +0200
Subject: [PATCH 3/4] test(bam): refuse to benchmark on a busy machine

I ran the A/B script without checking the load and published the result
against the author's number. The machine runs bwa-mem2 jobs at 15 cores.
The tail of my series drifted from 2.7 s to 5.1 s on both backends,
which I wrote off as the machine warming up; that is what an unrelated
job arriving mid-series looks like, and contention of that size can
invert a result rather than merely add noise to it.

The script now reads the load average before building, refuses above 2.0
unless BENCH_IGNORE_LOAD=1, and prints the load beside every pair so a
contaminated series is visible in its own output rather than having to
be inferred afterwards.

The numbers I posted are withdrawn in the PR description. The headline
is unverified in both directions until there is a clean run.
---
 test/bench_bam_backends.sh | 14 +++++++++++++-
 1 file changed, 13 insertions(+), 1 deletion(-)

diff --git a/test/bench_bam_backends.sh b/test/bench_bam_backends.sh
index 9a03693..290d6e6 100755
--- a/test/bench_bam_backends.sh
+++ b/test/bench_bam_backends.sh
@@ -21,6 +21,18 @@ ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)"
 WORK=$(mktemp -d)
 trap 'rm -rf "$WORK"' EXIT
 
+# Refuse to measure on a busy machine. A single unrelated job saturating the
+# cores does not just add noise, it can invert the result, and it is invisible
+# in the numbers afterwards: the series simply drifts. Checking beforehand is
+# the only cheap way to know the measurement means anything.
+load_now() { uptime | sed -E 's/.*load averages?: *([0-9.]+).*/\1/'; }
+LOAD=$(load_now)
+if awk -v l="$LOAD" 'BEGIN{exit !(l > 2.0)}'; then
+  echo "load average is $LOAD; other work is running and the numbers would not mean anything." >&2
+  echo "wait for the machine to be idle, or set BENCH_IGNORE_LOAD=1 to override." >&2
+  [ "${BENCH_IGNORE_LOAD:-0}" = "1" ] || exit 1
+fi
+
 echo "building both backends" >&2
 cargo build --release --manifest-path "$ROOT/Cargo.toml" >&2
 cp "$ROOT/target/release/rustar-aligner" "$WORK/noodles"
@@ -49,5 +61,5 @@ for i in $(seq 1 "$PAIRS"); do
     b=$(run "$WORK/htslib" "$WORK/h")
     a=$(run "$WORK/noodles" "$WORK/n")
   fi
-  echo "pair$i noodles=${a}s htslib=${b}s"
+  echo "pair$i noodles=${a}s htslib=${b}s load=$(load_now)"
 done

From d8bad6f171392b69d408daea2f12d881a8fdaa37 Mon Sep 17 00:00:00 2001
From: Benjamin Demaille <benjamin.demaille@icloud.com>
Date: Thu, 30 Jul 2026 00:59:18 +0200
Subject: [PATCH 4/4] test(bam): gate the benchmark on CPU idle, not on load
 average

Load average is an exponential average over minutes, so it stays high
long after the offending job is gone and refuses to measure on a machine
that is now quiet. It read 2.24 here with every core free. Idle
percentage answers the question actually being asked: are the cores free
right now. The script samples it before building and prints it beside
every pair.
---
 test/bench_bam_backends.sh | 18 +++++++++++++-----
 1 file changed, 13 insertions(+), 5 deletions(-)

diff --git a/test/bench_bam_backends.sh b/test/bench_bam_backends.sh
index 290d6e6..ad9241e 100755
--- a/test/bench_bam_backends.sh
+++ b/test/bench_bam_backends.sh
@@ -25,10 +25,18 @@ trap 'rm -rf "$WORK"' EXIT
 # cores does not just add noise, it can invert the result, and it is invisible
 # in the numbers afterwards: the series simply drifts. Checking beforehand is
 # the only cheap way to know the measurement means anything.
-load_now() { uptime | sed -E 's/.*load averages?: *([0-9.]+).*/\1/'; }
-LOAD=$(load_now)
-if awk -v l="$LOAD" 'BEGIN{exit !(l > 2.0)}'; then
-  echo "load average is $LOAD; other work is running and the numbers would not mean anything." >&2
+#
+# The check is on CPU idle, not on load average. Load average is an exponential
+# average over minutes, so it stays high long after the offending job has gone
+# and would refuse to measure on a machine that is now perfectly quiet. Idle
+# percentage answers the question actually being asked: are the cores free right
+# now.
+cpu_idle() { # percent idle, sampled over one second
+  top -l 2 -n 0 2>/dev/null | awk '/CPU usage/ {gsub("%","",$(NF-1)); v=$(NF-1)} END {print v+0}'
+}
+IDLE=$(cpu_idle)
+if awk -v i="$IDLE" 'BEGIN{exit !(i < 80)}'; then
+  echo "CPU is only ${IDLE}% idle; other work is running and the numbers would not mean anything." >&2
   echo "wait for the machine to be idle, or set BENCH_IGNORE_LOAD=1 to override." >&2
   [ "${BENCH_IGNORE_LOAD:-0}" = "1" ] || exit 1
 fi
@@ -61,5 +69,5 @@ for i in $(seq 1 "$PAIRS"); do
     b=$(run "$WORK/htslib" "$WORK/h")
     a=$(run "$WORK/noodles" "$WORK/n")
   fi
-  echo "pair$i noodles=${a}s htslib=${b}s load=$(load_now)"
+  echo "pair$i noodles=${a}s htslib=${b}s idle=$(cpu_idle)%"
 done