solo: CellRanger .h5 matrices, written without libhdf5 (answers #177, stacked on #179)

CHANGELOG.md

@@ -13,6 +13,42 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 ### Features
 
+- Under `--soloOutLayout CellRanger`, `raw_feature_bc_matrix.h5` and
+  `filtered_feature_bc_matrix.h5` are written beside the MEX directories,
+  holding the same counts in CellRanger's HDF5 layout, so
+  `scanpy.read_10x_h5`, Seurat's `Read10X_h5` and CellBender's `.h5` path
+  read them. Datasets are deflate-compressed, through the `libdeflater`
+  already in the tree. Written by an in-tree HDF5 writer
+  (`src/io/hdf5.rs`): **no new dependency**, no libhdf5, no `cmake`. See
+  `DIVERGENCE.md` §3.5 and `test/h5_conformance.sh`.
+
+- Under `--soloOutLayout CellRanger`, `metrics_summary.csv` is written with
+  CellRanger 10.0.0's 20 metrics, in its order and value formats. STARsolo's
+  `Summary.csv` is unchanged and still written alongside. Twelve of the 20
+  match a real `cellranger count` run exactly on the test fixture; the
+  interpretation behind the other eight is in `DIVERGENCE.md` §3.4.
+
+- `--soloOutLayout CellRanger` writes the solo matrices in the shape
+  `cellranger count` produces: `outs/raw_feature_bc_matrix/` and
+  `outs/filtered_feature_bc_matrix/`, gzipped, with a `-1` GEM-well suffix
+  on every barcode and one raw column per observed barcode. It implies
+  `--soloOutGzip yes`, `--soloOutRawBarcodes Observed` and an `outs/`
+  output directory, each still overridable on the command line. Counts are
+  unchanged. It is the default on 10x geometry, which **changes where
+  output files are written** on those runs; see `DIVERGENCE.md` §3.3.
+
+- On 10x geometry (`CB_UMI_Simple`, a whitelist, 16 bp CB, 10 or 12 bp
+  UMI), the five CellRanger-matching flags now **default** to their
+  CellRanger values. Any flag named on the command line wins, and the
+  substitution is logged. **This changes default output on 10x runs** and
+  diverges from STARsolo; see `DIVERGENCE.md` §1.3.
+
+- `--soloOutRawBarcodes Observed` writes the raw matrix with one column
+  per *observed* barcode instead of one per whitelist barcode, matching
+  what CellRanger's `raw_feature_bc_matrix` contains. Counts are
+  unchanged; on a 200-cell run `barcodes.tsv` goes from 62 MB to 3.4 kB.
+  **Not a STAR parameter**; default `Whitelist` keeps STARsolo behaviour.
+
 - **STARsolo single-cell quantification (`--soloType`)** — the 10x
   Chromium / plate-based count-matrix pipeline, ported from STAR and
   verified against real STARsolo (#90).
@@ -99,6 +135,12 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 ### Bug fixes
 
+- `--soloUMIfiltering MultiGeneUMI_CR` kept every gene tied at the
+  highest read count; CellRanger gives a tied UMI to no gene at all.
+  Since one read per gene is the ordinary shape of a multi-gene UMI, the
+  flag removed nothing in practice. On a 20k-read 10x fixture the count
+  matrix moves from 16 465 to 15 414 against STAR's 15 423.
+
 - **STARsolo `Gene` assignment now requires exon concordance**, matching
   STARsolo: a read counts toward a gene only when every aligned block
   lies within the gene's exons, rather than merely overlapping one. This

Cargo.toml

@@ -68,6 +68,9 @@ noodles-bgzf = { version = "0.49", features = ["libdeflate"] }
 
 [dev-dependencies]
 assert_cmd = "2"
+# Already a runtime dependency at the same version; listed here so the
+# integration tests can read the gzipped CellRanger-layout matrix files.
+flate2 = { version = "1", default-features = false, features = ["zlib-rs"] }
 predicates = "3"
 
 [build-dependencies]

DIVERGENCE.md

@@ -65,6 +65,204 @@ On the 10k yeast PE benchmark, 4 reads differ in alignment score (AS) because ST
 
 ---
 
+### 1.3 CellRanger behaviour is the default on 10x geometry
+
+**What STAR does.** STARsolo's defaults are its own (`1MM_multi`,
+`1MM_All`, no UMI filtering, `Hamming` clipping, `outFilterScoreMin 0`)
+whatever the barcode geometry. Matching CellRanger requires passing five flags,
+listed in STAR's `docs/STARsolo.md`.
+
+**What rustar-aligner does.** When the run is unambiguously 10x —
+`CB_UMI_Simple`, a whitelist, a 16-base CB and a 10- or 12-base UMI — those
+five flags default to their CellRanger values. Any flag given on the command
+line wins, and the substitution is logged in full.
+
+**Why.** A user aligning 10x data and comparing against CellRanger otherwise
+gets a successful run and different numbers, with nothing pointing at the five
+flags that explain it. Measured against CellRanger 10.0.0 on a 20 000-read
+fixture, those flags move the count matrix from 8.96% above CellRanger to
+2.17% above it. The remaining 2.17% is an open divergence: STAR 2.7.11b with
+the same flags is at +0.09%, so this closes most of the gap and not all of it.
+
+**Impact.** This is a **change of default output behaviour** and therefore the
+largest divergence in this file. It is confined to a geometry nothing else in
+common use shares, it is escapable by naming any flag explicitly, and it is
+announced at `INFO` on every run it touches. It needs maintainer sign-off.
+
+**Source.** `src/params/mod.rs` (`looks_like_10x`,
+`apply_cellranger_defaults_on_10x`). STAR: `docs/STARsolo.md`, "Matching
+CellRanger 4.x and 5.x results".
+
+---
+
+### 3.2 `--soloOutRawBarcodes Observed` (opt-in, non-STAR)
+
+**What STAR does.** STARsolo's raw matrix has one column per whitelist
+barcode, whether or not any read carried it. For 10x v3 that is 3 686 400
+columns and a 62 MB `barcodes.tsv`, nearly all of it zeros.
+
+**What rustar-aligner does.** The same, by default. `--soloOutRawBarcodes
+Observed` narrows the raw matrix to the barcodes that actually hold a count,
+which is what CellRanger's `raw_feature_bc_matrix` contains. On a 200-cell
+fixture that is 200 columns and a 3.4 kB `barcodes.tsv`.
+
+**Why.** Someone comparing our raw matrix against CellRanger's finds no
+overlapping keys at all, because the two files mean different things by "raw".
+The flag makes the comparison possible without changing what STARsolo users
+get.
+
+**Impact.** The counts are identical either way — same entries, same values,
+verified on the fixture — only the columns present differ. This is a non-STAR
+flag and needs maintainer sign-off; it is off by default so STARsolo parity is
+untouched.
+
+**Source.** `src/solo/count.rs` (`observed_barcodes`), `src/params/mod.rs`
+(`solo_out_raw_barcodes`). CellRanger: `outs/raw_feature_bc_matrix/` from a
+`cellranger count` run, observed directly rather than taken from its source.
+
+---
+
+### 3.3 `--soloOutLayout CellRanger` (non-STAR; on by default on 10x geometry)
+
+**What STAR does.** STARsolo writes
+`Solo.out/<feature>/{raw,filtered}/{matrix.mtx, barcodes.tsv, features.tsv}`,
+uncompressed, with bare barcodes.
+
+**What rustar-aligner does.** The same, by default, on non-10x geometry.
+`--soloOutLayout CellRanger` writes the same numbers in the shape
+`cellranger count` produces: `outs/{raw,filtered}_feature_bc_matrix/`, all three
+files gzipped, a `-1` GEM-well suffix on every barcode, one raw column per
+observed barcode, and no per-feature subdirectory when a single feature is
+requested. It implies `--soloOutGzip yes`, `--soloOutRawBarcodes Observed` and
+`--soloOutFileNames outs/ ...`, each still overridable on the command line.
+
+Under §1.3 the same 10x detection turns this on by default, so a bare 10x run
+lands in CellRanger's layout.
+
+**Why.** Tools written against CellRanger's `outs/` (scanpy's `read_10x_mtx`,
+Seurat's `Read10X`, any in-house loader) key on those directory names and on
+the `-1` suffix. Without them the numbers are right and nothing downstream can
+read them without a rename step.
+
+**Impact.** No count changes: verified entry by entry on the 20 000-read
+fixture, 13 959 entries and 15 439 counts in both layouts. On 10x geometry it
+**changes where output files are written**, which needs maintainer sign-off
+alongside §1.3. Against a real `cellranger count` run on the same fixture the
+raw barcode sets match exactly, 200 of 200.
+
+**Source.** `src/solo/count.rs` (`write_gene_matrix`, `write_one_barcode`),
+`src/params/mod.rs` (`solo_out_layout`, `apply_cellranger_layout`). CellRanger:
+`outs/` from a `cellranger count` 10.0.0 run, observed directly rather than
+taken from its source.
+
+---
+
+### 3.4 `metrics_summary.csv` under the CellRanger layout (non-STAR)
+
+**What STAR does.** STARsolo writes `Summary.csv`, its own metric set with its
+own names. It has no `metrics_summary.csv`.
+
+**What rustar-aligner does.** The same, by default. Under
+`--soloOutLayout CellRanger` it *additionally* writes
+`metrics_summary.csv` with CellRanger 10.0.0's 20 metrics, in CellRanger's
+order and value formats. `Summary.csv` is written unchanged alongside it, so
+nothing STARsolo-faithful is altered. Under the CellRanger layout the older
+`CellRanger.summary.csv` is not written: its four rows are a subset of the
+metrics file.
+
+**Why.** The file is how a 10x pipeline reads run quality. Its 20 names are
+also the clearest public statement of what CellRanger measures, which makes it
+useful as a target even where our value differs.
+
+**Impact.** No count changes. It costs one extra gene-body overlap query per
+read, because the exonic/intronic split needs it and a `Gene`-only run does not
+otherwise do it.
+
+**Not all 20 are certainties.** Twelve match a real `cellranger count` 10.0.0
+run exactly on the 20 000-read fixture. The other eight follow from our own
+tallies under a stated interpretation, because CellRanger does not document the
+denominators:
+
+* `Reads Mapped Confidently to *` reads "confidently" as MAPQ 255, i.e. our
+  uniquely-mapped set.
+* `Reads Mapped Confidently to Transcriptome` is the `Gene` feature's own
+  uniquely-assigned read tally.
+* `Valid UMI Sequences` is measured over reads that reached the UMI check, i.e.
+  those with a valid barcode.
+* `Sequencing Saturation` is `1 - molecules / reads` over the reads that
+  entered the matrix. This is the largest disagreement on the fixture: 12.8%
+  against CellRanger's 7.4%.
+* `Mean Reads per Cell` is total reads over called cells, not reads-in-cells
+  over called cells, which is what reproduces CellRanger's value.
+
+The cell-count-dependent metrics (`Median Genes per Cell`,
+`Median UMI Counts per Cell`, `Total Genes Detected`) differ by 2-4 on the
+fixture, following the count differences recorded in §1.3 rather than a
+different definition.
+
+**Source.** `src/solo/count.rs` (`write_metrics_summary`, `metric_int`,
+`metric_pct`), `src/solo/mod.rs` (`Q30Stats`). CellRanger: the
+`metrics_summary.csv` of a `cellranger count` 10.0.0 run, observed directly
+rather than taken from its source.
+
+---
+
+### 3.5 `*_feature_bc_matrix.h5` written without libhdf5 (non-STAR)
+
+**What STAR does.** STARsolo writes MEX only. It has no HDF5 output.
+
+**What rustar-aligner does.** The same, by default. Under
+`--soloOutLayout CellRanger` it additionally writes
+`raw_feature_bc_matrix.h5` and `filtered_feature_bc_matrix.h5` beside the MEX
+directories, holding the same counts in CellRanger's HDF5 layout.
+
+**Why.** `scanpy.read_10x_h5`, Seurat's `Read10X_h5` and CellBender's `.h5`
+path all read that single file, and a good deal of 10x tooling reaches for it
+in preference to the directory.
+
+**How, and why that is the divergence.** The file is produced by
+`src/io/hdf5.rs`, a writer in this repository, not by libhdf5. A binding would
+mean either a system `libhdf5` — which would make `cargo install
+rustar-aligner` fail on a machine without it — or `cmake` and a multi-minute C
+build on all five CI targets, for one optional output file. The maintainer
+declined both, so the file is written directly.
+
+The writer emits the oldest and most widely-read variant of each structure:
+version 0 superblock, version 1 object headers, symbol-table groups with a
+version 1 B-tree and a local heap, and datasets stored either contiguously or
+as a single deflate-compressed chunk indexed by a version 1 chunk B-tree. It
+does not write variable-length types, and it never splits a dataset into more
+than one chunk.
+
+**Where the file differs from CellRanger's.** Both deflate their datasets at
+level 6. CellRanger pads every string field to 256 bytes and chunks at 65 536
+elements; ours are padded to the longest value present and stored as one chunk
+per dataset. Readers see identical values either way — the difference is
+on-disk size and the `|S<n>` width numpy reports. On the fixture ours is 88 216
+bytes against CellRanger's 149 472, the string padding accounting for most of
+it. The `genome` field carries the `--genomeDir` directory name, since we have
+no equivalent of the name given to `cellranger mkref`, and `software_version`
+says `rustar-aligner`, not a CellRanger version.
+
+**Impact.** No count changes: the `.h5` and the `.mtx` were compared entry by
+entry on the fixture and hold the same 13 959 entries and 15 439 counts.
+
+**How it is checked.** `cargo test` cannot validate an HDF5 file, so it locks
+determinism and the byte layout with a checksum, and the real check is
+`test/h5_conformance.sh`, which runs the HDF5 command-line tools and the
+readers. On the fixture: `h5ls -r`, `h5dump -H` and `h5stat` all succeed, an
+`h5repack` round-trip through libhdf5 followed by `h5diff` reports no
+difference, and scanpy's and CellBender's loaders both return the expected
+matrix. That script must be re-run, off CI, whenever `src/io/hdf5.rs` changes:
+only libhdf5 can say whether new bytes are still a valid file.
+
+**Source.** `src/io/hdf5.rs`, `src/solo/count.rs` (`write_matrix_h5`). Format:
+*HDF5 File Format Specification Version 3.0*. CellRanger: the `.h5` of a
+`cellranger count` 10.0.0 run, inspected with `h5py` rather than taken from its
+source.
+
+---
+
 ## 4. Implementation divergences (no intended output difference)
 
 These differ in *how* a result is produced, not *what* is produced. They are documented so a reviewer chasing a discrepancy knows the mechanism differs by design.