diff --git a/CHANGELOG.md b/CHANGELOG.md index b6e4893..6b0dc86 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/Cargo.toml b/Cargo.toml index 8a4638f..444553d 100644 --- a/Cargo.toml +++ b/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] diff --git a/DIVERGENCE.md b/DIVERGENCE.md index bd957ed..749c4c5 100644 --- a/DIVERGENCE.md +++ b/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//{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` 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. diff --git a/src/io/hdf5.rs b/src/io/hdf5.rs new file mode 100644 index 0000000..62771da --- /dev/null +++ b/src/io/hdf5.rs @@ -0,0 +1,980 @@ +//! A minimal HDF5 writer, with no dependency on libhdf5. +//! +//! This exists for one file: CellRanger's `raw_feature_bc_matrix.h5`. That file +//! uses a small, fixed corner of the format — one group tree, nine 1-D +//! datasets, three scalar types, a handful of root attributes — and writing +//! that corner directly is a few hundred lines, against a C library that would +//! need either a system `libhdf5` or `cmake` on every build (see +//! `DIVERGENCE.md` §3.5 and the discussion in the HDF5 dependency issue). +//! +//! # What it writes +//! +//! Version 0 superblock, version 1 object headers, old-style groups (symbol +//! table + version 1 B-tree + local heap), and datasets stored either +//! contiguously or as a single deflate-compressed chunk. Every one of those is +//! the oldest and most widely-read variant of its structure, which is +//! deliberate: this writer's only job is to produce files other people's +//! readers open. +//! +//! # What it does not write +//! +//! No variable-length types, no references, no links other than the hard links +//! in a group's symbol table, no deletion, no rewriting. Every dataset is +//! written once, in full, from a slice already in memory, as one chunk rather +//! than many. Groups hold at most `2 * GROUP_LEAF_K` = 8 entries, since a +//! second B-tree node would need node splitting and nothing here comes close. +//! +//! The layout constants and structure layouts follow the *HDF5 File Format +//! Specification Version 3.0* (superblock §II.A, B-trees §II.A.1, local heaps +//! §II.D, object headers §IV.A, messages §IV.A.2). + +use std::io::Write; +use std::path::Path; + +use crate::error::Error; + +/// `\x89HDF\r\n\x1a\n` — the format signature, chosen (like PNG's) so that a +/// transfer that mangles high bits or line endings is detectable. +const SIGNATURE: [u8; 8] = [0x89, b'H', b'D', b'F', b'\r', b'\n', 0x1a, b'\n']; + +/// Sizes of offsets and lengths in the file. 8 bytes each, so a file may exceed +/// 4 GB; a raw 10x matrix will not, but a genome-scale one might. +const SIZEOF_ADDR: usize = 8; + +/// Group node K values recorded in the superblock. A leaf symbol-table node +/// holds up to `2 * GROUP_LEAF_K` entries and a B-tree node up to +/// `2 * GROUP_INTERNAL_K` children; we only ever produce one of each. +const GROUP_LEAF_K: u16 = 4; +const GROUP_INTERNAL_K: u16 = 16; + +/// The "undefined address" HDF5 uses where a pointer is absent: all bits set. +const UNDEF_ADDR: u64 = u64::MAX; + +/// "No free block" in a local heap's free list. Not the undefined address: the +/// library uses the literal 1, which cannot be a real offset because free +/// blocks are 8-byte aligned (`H5HL_FREE_NULL` in `H5HLpkg.h`). Writing the +/// undefined address here makes `H5Ovisit` fail with "bad heap free list" +/// while the file still opens, which is a memorable way to learn the +/// difference. +const HEAP_FREE_NULL: u64 = 1; + +/// Superblock v0 is 56 bytes, followed by the 40-byte root symbol table entry. +const SUPERBLOCK_LEN: usize = 96; + +/// Object header message types used here. +const MSG_DATASPACE: u16 = 0x0001; +const MSG_DATATYPE: u16 = 0x0003; +const MSG_FILL_VALUE: u16 = 0x0005; +const MSG_FILTER_PIPELINE: u16 = 0x000B; +const MSG_DATA_LAYOUT: u16 = 0x0008; +const MSG_ATTRIBUTE: u16 = 0x000C; +const MSG_SYMBOL_TABLE: u16 = 0x0011; + +/// One dataset's contents. Every variant is a 1-D array; HDF5 scalars appear +/// only as attributes. +/// +/// `FixedString` is the padded, fixed-width byte string CellRanger uses for +/// barcodes and feature fields (`|S45`, `|S256` in numpy terms). Values longer +/// than `width` are an error rather than a silent truncation, because a +/// truncated barcode is a wrong answer that looks like a right one. +pub enum Data<'a> { + I32(&'a [i32]), + I64(&'a [i64]), + FixedString { width: usize, values: Vec<&'a str> }, +} + +impl Data<'_> { + fn len(&self) -> usize { + match self { + Data::I32(v) => v.len(), + Data::I64(v) => v.len(), + Data::FixedString { values, .. } => values.len(), + } + } + + fn element_size(&self) -> usize { + match self { + Data::I32(_) => 4, + Data::I64(_) => 8, + Data::FixedString { width, .. } => *width, + } + } + + /// The datatype message body for this data (see `datatype_message`). + fn datatype(&self) -> Vec { + match self { + Data::I32(_) => fixed_point_datatype(4), + Data::I64(_) => fixed_point_datatype(8), + Data::FixedString { width, .. } => string_datatype(*width), + } + } + + /// The raw bytes of the array, little-endian, in element order. + fn bytes(&self) -> Result, String> { + let mut out = Vec::with_capacity(self.len() * self.element_size()); + match self { + Data::I32(v) => v + .iter() + .for_each(|x| out.extend_from_slice(&x.to_le_bytes())), + Data::I64(v) => v + .iter() + .for_each(|x| out.extend_from_slice(&x.to_le_bytes())), + Data::FixedString { width, values } => { + for s in values { + if s.len() > *width { + return Err(format!( + "value {s:?} is {} bytes, longer than the {width}-byte field", + s.len() + )); + } + out.extend_from_slice(s.as_bytes()); + out.resize(out.len() + (width - s.len()), 0); + } + } + } + Ok(out) + } +} + +/// A scalar or 1-D attribute attached to a group. +pub enum AttrValue { + /// Fixed-width string, written at exactly its own length. + Str(String), + /// 1-D array of fixed-width strings, each padded to `width`. + StrArray { + width: usize, + values: Vec, + }, + I64(i64), + I64Array(Vec), +} + +pub struct Attr { + pub name: String, + pub value: AttrValue, +} + +pub struct DatasetSpec<'a> { + pub name: &'a str, + pub data: Data<'a>, +} + +/// A group: named sub-groups, named datasets, and attributes. Children are +/// sorted by name at write time, as a symbol table requires. +pub struct GroupSpec<'a> { + pub name: &'a str, + pub groups: Vec>, + pub datasets: Vec>, + pub attrs: Vec, +} + +impl<'a> GroupSpec<'a> { + #[must_use] + pub fn new(name: &'a str) -> Self { + Self { + name, + groups: Vec::new(), + datasets: Vec::new(), + attrs: Vec::new(), + } + } + + #[must_use] + pub fn dataset(mut self, name: &'a str, data: Data<'a>) -> Self { + self.datasets.push(DatasetSpec { name, data }); + self + } + + #[must_use] + pub fn group(mut self, g: GroupSpec<'a>) -> Self { + self.groups.push(g); + self + } + + #[must_use] + pub fn attr(mut self, name: &str, value: AttrValue) -> Self { + self.attrs.push(Attr { + name: name.to_string(), + value, + }); + self + } + + fn child_count(&self) -> usize { + self.groups.len() + self.datasets.len() + } +} + +// --------------------------------------------------------------------------- +// Byte-level helpers +// --------------------------------------------------------------------------- + +/// The file image under construction. Blocks are appended and never moved, so +/// an address handed out stays valid; everything is 8-byte aligned, which every +/// structure here wants anyway. +struct Image { + buf: Vec, +} + +impl Image { + fn new() -> Self { + // The superblock is written last, once the root group's address and the + // end-of-file address are known, so leave room for it. + Self { + buf: vec![0u8; SUPERBLOCK_LEN], + } + } + + fn align(&mut self) { + while !self.buf.len().is_multiple_of(8) { + self.buf.push(0); + } + } + + /// Append a block and return its address. + fn alloc(&mut self, bytes: &[u8]) -> u64 { + self.align(); + let addr = self.buf.len() as u64; + self.buf.extend_from_slice(bytes); + addr + } + + /// Reserve `len` zeroed bytes and return the address, for a block that is + /// filled in after its own contents are known. + fn reserve(&mut self, len: usize) -> u64 { + self.align(); + let addr = self.buf.len() as u64; + self.buf.resize(self.buf.len() + len, 0); + addr + } + + fn patch(&mut self, addr: u64, bytes: &[u8]) { + let start = addr as usize; + self.buf[start..start + bytes.len()].copy_from_slice(bytes); + } +} + +fn pad_to_8(v: &mut Vec) { + while !v.len().is_multiple_of(8) { + v.push(0); + } +} + +/// The first byte of a datatype message: version in the high nibble, class in +/// the low one. +fn datatype_tag(class: u8) -> u8 { + (1 << 4) | class +} + +/// Datatype message body for a signed little-endian integer of `size` bytes. +/// +/// Class 0 (fixed-point). The class bit field's bit 3 marks it signed; byte +/// order (bit 0) and padding (bits 1-2) are zero, i.e. little-endian with zero +/// padding. The properties are the bit offset and precision. +fn fixed_point_datatype(size: u32) -> Vec { + const CLASS_FIXED_POINT: u8 = 0; + let mut m = Vec::with_capacity(12); + m.push(datatype_tag(CLASS_FIXED_POINT)); + m.extend_from_slice(&[0x08, 0x00, 0x00]); // signed + m.extend_from_slice(&size.to_le_bytes()); + m.extend_from_slice(&0u16.to_le_bytes()); // bit offset + m.extend_from_slice(&((size * 8) as u16).to_le_bytes()); // precision + m +} + +/// Datatype message body for a fixed-width, null-padded ASCII string. +/// +/// Class 3 (string), padding type 1 (null pad) in bits 0-3 and character set 0 +/// (ASCII) in bits 4-7 — the combination numpy reads back as `|S`. +fn string_datatype(width: usize) -> Vec { + let mut m = Vec::with_capacity(8); + const CLASS_STRING: u8 = 3; + m.push(datatype_tag(CLASS_STRING)); + m.extend_from_slice(&[0x01, 0x00, 0x00]); // null pad, ASCII + m.extend_from_slice(&(width as u32).to_le_bytes()); + m +} + +/// Dataspace message body: version 1, `dims.len()` dimensions, no maximum +/// dimensions and no permutation index. An empty `dims` is a scalar. +fn dataspace_message(dims: &[u64]) -> Vec { + let mut m = Vec::with_capacity(8 + dims.len() * 8); + m.push(1); // version + m.push(dims.len() as u8); + m.push(0); // flags: no max dims + m.extend_from_slice(&[0, 0, 0, 0, 0]); // reserved + for d in dims { + m.extend_from_slice(&d.to_le_bytes()); + } + m +} + +/// Space allocation time in a fill value message: when the library would +/// reserve the dataset's storage. `Late` for contiguous storage, `Incremental` +/// for chunked, which is what HDF5 itself writes and what `h5diff` expects to +/// find on a chunked dataset. +const ALLOC_LATE: u8 = 2; +const ALLOC_INCREMENTAL: u8 = 3; + +/// Fill value message, version 2, with no fill value defined: readers use the +/// type's default. Nothing here is ever read before it is written. +/// +/// The write time is `H5D_FILL_TIME_IFSET` (2), i.e. write the fill value only +/// if one was set. With `Late` allocation and a write time of "on allocation", +/// `h5diff` refuses a chunked dataset outright — the two say storage is +/// reserved late but filled at reservation. +fn fill_value_message(alloc_time: u8) -> Vec { + // version, allocation time, write time, "fill value defined", then a + // zero-length value, which is how HDF5 spells "use the type's default". + // Writing "not defined" (0) instead leaves a file h5dump reads but + // h5repack rewrites incorrectly. + let mut m = vec![2, alloc_time, 2, 1]; + m.extend_from_slice(&0u32.to_le_bytes()); // size of the fill value + m +} + +/// Data layout message, version 3, contiguous storage at `addr` for `size` +/// bytes. +fn data_layout_message(addr: u64, size: u64) -> Vec { + let mut m = Vec::with_capacity(18); + m.push(3); // version + m.push(1); // class 1: contiguous + m.extend_from_slice(&addr.to_le_bytes()); + m.extend_from_slice(&size.to_le_bytes()); + m +} + +/// Data layout message, version 3, chunked storage. +/// +/// `btree` is the address of the chunk index; `chunk_len` is the chunk's length +/// in elements and `elem_size` the element size in bytes. The dimensionality +/// field counts one more than the dataset's rank, because HDF5 treats the +/// element size as a trailing dimension of the chunk. +fn chunked_layout_message(btree: u64, chunk_len: u32, elem_size: u32) -> Vec { + let mut m = Vec::with_capacity(19); + m.push(3); // version + m.push(2); // class 2: chunked + m.push(2); // dimensionality: rank 1, plus the element-size dimension + m.extend_from_slice(&btree.to_le_bytes()); + m.extend_from_slice(&chunk_len.to_le_bytes()); + m.extend_from_slice(&elem_size.to_le_bytes()); + m +} + +/// Filter pipeline message, version 1, with a single deflate filter. +/// +/// Filter id 1 is `H5Z_FILTER_DEFLATE`. Its one client-data value is the +/// compression level, which readers ignore — zlib streams are self-describing — +/// but which is part of the message HDF5 writes and expects to parse back. +/// +/// Version 1 stores the filter name, null-terminated and padded to a multiple +/// of 8 bytes; `"deflate"` is what the library writes there. +fn filter_pipeline_message(level: u32) -> Vec { + const NAME: &[u8; 8] = b"deflate\0"; + let mut m = Vec::with_capacity(32); + m.push(1); // version + m.push(1); // number of filters + m.extend_from_slice(&[0; 6]); // reserved + m.extend_from_slice(&1u16.to_le_bytes()); // filter id: deflate + m.extend_from_slice(&(NAME.len() as u16).to_le_bytes()); + m.extend_from_slice(&0u16.to_le_bytes()); // flags: the filter is mandatory + m.extend_from_slice(&1u16.to_le_bytes()); // one client data value + m.extend_from_slice(NAME); + m.extend_from_slice(&level.to_le_bytes()); + // An odd number of client data values is padded to a multiple of 8 bytes. + m.extend_from_slice(&[0; 4]); + m +} + +/// A version 1 B-tree indexing the chunks of one dataset (node type 1). +/// +/// One chunk per dataset, so this is a single leaf holding one entry. A chunk +/// key is the chunk's stored size, its filter mask, and its offset in the +/// dataset — one offset per dimension plus a trailing zero for the element-size +/// dimension. The trailing key marks the end of the last chunk, which is how a +/// reader knows where the indexed space stops. +fn chunk_btree_node(chunk_addr: u64, chunk_size: u32, n_elements: u64) -> Vec { + /// The default `H5F_CRT_BTREE_RANK` for chunk B-trees. Superblock version 0 + /// records only the group node K values, so a reader takes this one from the + /// library default and sizes the node accordingly; the node therefore has to + /// be allocated at full capacity even holding one entry. + const CHUNK_BTREE_K: usize = 32; + + let key = |size: u32, offset: u64| { + let mut k = Vec::with_capacity(24); + k.extend_from_slice(&size.to_le_bytes()); + k.extend_from_slice(&0u32.to_le_bytes()); // filter mask: all filters applied + k.extend_from_slice(&offset.to_le_bytes()); // offset along dimension 0 + k.extend_from_slice(&0u64.to_le_bytes()); // element-size dimension + k + }; + + let mut b = Vec::new(); + b.extend_from_slice(b"TREE"); + b.push(1); // node type 1: raw data chunks + b.push(0); // level 0: leaf + b.extend_from_slice(&1u16.to_le_bytes()); // entries used + b.extend_from_slice(&UNDEF_ADDR.to_le_bytes()); // left sibling + b.extend_from_slice(&UNDEF_ADDR.to_le_bytes()); // right sibling + b.extend_from_slice(&key(chunk_size, 0)); + b.extend_from_slice(&chunk_addr.to_le_bytes()); + b.extend_from_slice(&key(0, n_elements)); + let key_len = 24; + b.resize( + 24 + 2 * CHUNK_BTREE_K * (key_len + SIZEOF_ADDR) + key_len, + 0, + ); + b +} + +/// Datasets at or below this many raw bytes are stored contiguously. +/// +/// Compressing them would cost a chunk B-tree (about 2 kB, allocated at full +/// capacity whatever the occupancy) and a filter pipeline message to save a few +/// hundred bytes. CellRanger draws the same line in the same place: its +/// `_all_tag_keys` is the one dataset it leaves uncompressed. +const COMPRESSION_THRESHOLD: usize = 1024; + +/// zlib compression level for dataset chunks. 6 is zlib's default and what +/// CellRanger's files decompress as. +const COMPRESSION_LEVEL: u32 = 6; + +/// Compress a chunk in the zlib format (RFC 1950) that HDF5's deflate filter +/// expects — not raw deflate, and not gzip. +fn deflate_chunk(raw: &[u8]) -> Vec { + let lvl = libdeflater::CompressionLvl::new(COMPRESSION_LEVEL as i32).unwrap_or_default(); + let mut c = libdeflater::Compressor::new(lvl); + let mut out = vec![0u8; c.zlib_compress_bound(raw.len())]; + match c.zlib_compress(raw, &mut out) { + Ok(n) => { + out.truncate(n); + out + } + // libdeflate only fails here if the output buffer is too small, which + // `zlib_compress_bound` rules out. Falling back to the uncompressed + // bytes would produce a file whose filter says otherwise, so this + // cannot silently degrade. + Err(e) => unreachable!("libdeflate zlib_compress failed on a bounded buffer: {e:?}"), + } +} + +/// One object header message, prefixed and padded as the header expects. +fn header_message(kind: u16, body: &[u8]) -> Vec { + let mut padded = body.to_vec(); + pad_to_8(&mut padded); + let mut m = Vec::with_capacity(8 + padded.len()); + m.extend_from_slice(&kind.to_le_bytes()); + m.extend_from_slice(&(padded.len() as u16).to_le_bytes()); + m.push(0); // flags + m.extend_from_slice(&[0, 0, 0]); // reserved + m.extend_from_slice(&padded); + m +} + +/// A version 1 object header wrapping already-encoded messages. +fn object_header(messages: &[Vec]) -> Vec { + let body_len: usize = messages.iter().map(Vec::len).sum(); + let mut h = Vec::with_capacity(16 + body_len); + h.push(1); // version + h.push(0); // reserved + h.extend_from_slice(&(messages.len() as u16).to_le_bytes()); + h.extend_from_slice(&1u32.to_le_bytes()); // reference count + h.extend_from_slice(&(body_len as u32).to_le_bytes()); + h.extend_from_slice(&[0, 0, 0, 0]); // pad the prefix to 16 bytes + for m in messages { + h.extend_from_slice(m); + } + h +} + +/// An attribute message, version 1: name, then the datatype and dataspace +/// messages inline, then the value. Each of the three is padded to 8 bytes. +fn attribute_message(name: &str, datatype: &[u8], dataspace: &[u8], data: &[u8]) -> Vec { + let mut name_b = name.as_bytes().to_vec(); + name_b.push(0); + let name_size = name_b.len(); + pad_to_8(&mut name_b); + let mut dt = datatype.to_vec(); + let dt_size = dt.len(); + pad_to_8(&mut dt); + let mut ds = dataspace.to_vec(); + let ds_size = ds.len(); + pad_to_8(&mut ds); + + let mut m = Vec::new(); + m.push(1); // version + m.push(0); // reserved + m.extend_from_slice(&(name_size as u16).to_le_bytes()); + m.extend_from_slice(&(dt_size as u16).to_le_bytes()); + m.extend_from_slice(&(ds_size as u16).to_le_bytes()); + m.extend_from_slice(&name_b); + m.extend_from_slice(&dt); + m.extend_from_slice(&ds); + m.extend_from_slice(data); + m +} + +/// One 40-byte symbol table entry. +/// +/// `scratch` is the 16-byte scratch pad: for a group (cache type 1) it holds +/// the B-tree and local heap addresses, letting a reader find a subgroup's +/// contents without reading its object header first. Datasets use cache type 0 +/// and leave it zero. +fn symbol_entry(name_offset: u64, header_addr: u64, cache_type: u32, scratch: [u8; 16]) -> Vec { + let mut e = Vec::with_capacity(40); + e.extend_from_slice(&name_offset.to_le_bytes()); + e.extend_from_slice(&header_addr.to_le_bytes()); + e.extend_from_slice(&cache_type.to_le_bytes()); + e.extend_from_slice(&[0, 0, 0, 0]); // reserved + e.extend_from_slice(&scratch); + e +} + +// --------------------------------------------------------------------------- +// Writing +// --------------------------------------------------------------------------- + +/// What writing one group produced: where its object header, B-tree and local +/// heap live. The parent needs all three — the header for the symbol table +/// entry, the other two for that entry's scratch pad. +struct WrittenGroup { + header: u64, + btree: u64, + heap: u64, +} + +/// Write `root`'s tree into `img` and return its addresses. +/// +/// Children are written first so that their addresses are known when the +/// parent's symbol table is built, which is why this recurses before it +/// allocates anything of its own. +fn write_group(img: &mut Image, group: &GroupSpec<'_>) -> Result { + if group.child_count() > 2 * GROUP_LEAF_K as usize { + return Err(format!( + "group {:?} has {} children; this writer emits a single symbol table node, \ + which holds at most {}", + group.name, + group.child_count(), + 2 * GROUP_LEAF_K + )); + } + + // (name, object header address, cache type, scratch pad) + let mut children: Vec<(String, u64, u32, [u8; 16])> = Vec::new(); + + for sub in &group.groups { + let w = write_group(img, sub)?; + let mut scratch = [0u8; 16]; + scratch[..8].copy_from_slice(&w.btree.to_le_bytes()); + scratch[8..].copy_from_slice(&w.heap.to_le_bytes()); + children.push((sub.name.to_string(), w.header, 1, scratch)); + } + + for ds in &group.datasets { + let bytes = ds.data.bytes()?; + let n = ds.data.len(); + + // Anything worth compressing is stored as a single deflated chunk + // covering the whole dataset; everything else is contiguous. One chunk + // means one B-tree entry and no node splitting, and costs nothing in + // practice because every reader of these files loads each dataset + // whole. + let mut messages = vec![ + header_message(MSG_DATASPACE, &dataspace_message(&[n as u64])), + header_message(MSG_DATATYPE, &ds.data.datatype()), + ]; + if bytes.len() > COMPRESSION_THRESHOLD && u32::try_from(n).is_ok() { + let chunk = deflate_chunk(&bytes); + let chunk_size = u32::try_from(chunk.len()).map_err(|_| { + format!( + "dataset {:?} compresses to {} bytes, past the 4 GB a chunk record holds", + ds.name, + chunk.len() + ) + })?; + let chunk_addr = img.alloc(&chunk); + let btree_addr = img.alloc(&chunk_btree_node(chunk_addr, chunk_size, n as u64)); + messages.push(header_message( + MSG_FILL_VALUE, + &fill_value_message(ALLOC_INCREMENTAL), + )); + messages.push(header_message( + MSG_FILTER_PIPELINE, + &filter_pipeline_message(COMPRESSION_LEVEL), + )); + messages.push(header_message( + MSG_DATA_LAYOUT, + &chunked_layout_message(btree_addr, n as u32, ds.data.element_size() as u32), + )); + } else { + let data_addr = if bytes.is_empty() { + // A zero-length dataset still needs an address; point it at the + // current end of file rather than at the undefined address, + // which some readers treat as "not allocated" and then refuse + // to read. + img.reserve(0) + } else { + img.alloc(&bytes) + }; + messages.push(header_message( + MSG_FILL_VALUE, + &fill_value_message(ALLOC_LATE), + )); + messages.push(header_message( + MSG_DATA_LAYOUT, + &data_layout_message(data_addr, bytes.len() as u64), + )); + } + let hdr = img.alloc(&object_header(&messages)); + children.push((ds.name.to_string(), hdr, 0, [0u8; 16])); + } + + // A symbol table is searched by binary search on the name, so entries must + // be in lexicographic order. + children.sort_by(|a, b| a.0.cmp(&b.0)); + + // Local heap data segment: an empty string at offset 0 (the B-tree's first + // key points at it), then each child's name, null-terminated. + let mut heap_data: Vec = vec![0u8; 8]; + let mut name_offsets: Vec = Vec::with_capacity(children.len()); + for (name, _, _, _) in &children { + let off = heap_data.len() as u64; + heap_data.extend_from_slice(name.as_bytes()); + heap_data.push(0); + pad_to_8(&mut heap_data); + name_offsets.push(off); + } + let heap_data_addr = img.alloc(&heap_data); + + let mut heap = Vec::with_capacity(32); + heap.extend_from_slice(b"HEAP"); + heap.push(0); // version + heap.extend_from_slice(&[0, 0, 0]); // reserved + heap.extend_from_slice(&(heap_data.len() as u64).to_le_bytes()); + heap.extend_from_slice(&HEAP_FREE_NULL.to_le_bytes()); // no free block + heap.extend_from_slice(&heap_data_addr.to_le_bytes()); + let heap_addr = img.alloc(&heap); + + // Symbol table node. Allocated at full capacity (2K entries) whatever the + // occupancy, as HDF5 does. + let capacity = 2 * GROUP_LEAF_K as usize; + let mut snod = Vec::with_capacity(8 + capacity * 40); + snod.extend_from_slice(b"SNOD"); + snod.push(1); // version + snod.push(0); // reserved + snod.extend_from_slice(&(children.len() as u16).to_le_bytes()); + for (i, (_, hdr, cache, scratch)) in children.iter().enumerate() { + snod.extend_from_slice(&symbol_entry(name_offsets[i], *hdr, *cache, *scratch)); + } + snod.resize(8 + capacity * 40, 0); + let snod_addr = img.alloc(&snod); + + // A single-node B-tree over that one leaf. The keys bracket the node: the + // first is the empty string at heap offset 0, the last is the greatest + // name present. + let last_key = name_offsets.last().copied().unwrap_or(0); + let mut btree = Vec::new(); + btree.extend_from_slice(b"TREE"); + btree.push(0); // node type 0: group node + btree.push(0); // level 0: leaf + btree.extend_from_slice(&1u16.to_le_bytes()); // entries used + btree.extend_from_slice(&UNDEF_ADDR.to_le_bytes()); // left sibling + btree.extend_from_slice(&UNDEF_ADDR.to_le_bytes()); // right sibling + btree.extend_from_slice(&0u64.to_le_bytes()); // key 0 + btree.extend_from_slice(&snod_addr.to_le_bytes()); // child 0 + btree.extend_from_slice(&last_key.to_le_bytes()); // key 1 + // Allocated at full capacity, like the symbol table node. + let btree_capacity = 24 + (2 * GROUP_INTERNAL_K as usize) * (SIZEOF_ADDR * 2) + SIZEOF_ADDR; + btree.resize(btree_capacity, 0); + let btree_addr = img.alloc(&btree); + + let mut messages = vec![header_message(MSG_SYMBOL_TABLE, &{ + let mut b = Vec::with_capacity(16); + b.extend_from_slice(&btree_addr.to_le_bytes()); + b.extend_from_slice(&heap_addr.to_le_bytes()); + b + })]; + for a in &group.attrs { + messages.push(header_message(MSG_ATTRIBUTE, &encode_attr(a))); + } + let header_addr = img.alloc(&object_header(&messages)); + + Ok(WrittenGroup { + header: header_addr, + btree: btree_addr, + heap: heap_addr, + }) +} + +fn encode_attr(a: &Attr) -> Vec { + let (dt, ds, data) = match &a.value { + AttrValue::Str(s) => (string_datatype(s.len().max(1)), dataspace_message(&[]), { + let mut v = s.as_bytes().to_vec(); + if v.is_empty() { + v.push(0); + } + v + }), + AttrValue::StrArray { width, values } => { + let mut data = Vec::with_capacity(values.len() * width); + for s in values { + let b = s.as_bytes(); + let n = b.len().min(*width); + data.extend_from_slice(&b[..n]); + data.resize(data.len() + (width - n), 0); + } + ( + string_datatype(*width), + dataspace_message(&[values.len() as u64]), + data, + ) + } + AttrValue::I64(v) => ( + fixed_point_datatype(8), + dataspace_message(&[]), + v.to_le_bytes().to_vec(), + ), + AttrValue::I64Array(vs) => { + let mut data = Vec::with_capacity(vs.len() * 8); + for v in vs { + data.extend_from_slice(&v.to_le_bytes()); + } + ( + fixed_point_datatype(8), + dataspace_message(&[vs.len() as u64]), + data, + ) + } + }; + attribute_message(&a.name, &dt, &ds, &data) +} + +/// Serialise `root` as a complete HDF5 file image. +pub fn build(root: &GroupSpec<'_>) -> Result, String> { + let mut img = Image::new(); + let w = write_group(&mut img, root)?; + img.align(); + let eof = img.buf.len() as u64; + + let mut sb = Vec::with_capacity(SUPERBLOCK_LEN); + sb.extend_from_slice(&SIGNATURE); + sb.push(0); // superblock version + sb.push(0); // free space storage version + sb.push(0); // root group symbol table entry version + sb.push(0); // reserved + sb.push(0); // shared header message format version + sb.push(SIZEOF_ADDR as u8); + sb.push(SIZEOF_ADDR as u8); // size of lengths + sb.push(0); // reserved + sb.extend_from_slice(&GROUP_LEAF_K.to_le_bytes()); + sb.extend_from_slice(&GROUP_INTERNAL_K.to_le_bytes()); + sb.extend_from_slice(&0u32.to_le_bytes()); // file consistency flags + sb.extend_from_slice(&0u64.to_le_bytes()); // base address + sb.extend_from_slice(&UNDEF_ADDR.to_le_bytes()); // free space info + sb.extend_from_slice(&eof.to_le_bytes()); + sb.extend_from_slice(&UNDEF_ADDR.to_le_bytes()); // driver information block + + let mut scratch = [0u8; 16]; + scratch[..8].copy_from_slice(&w.btree.to_le_bytes()); + scratch[8..].copy_from_slice(&w.heap.to_le_bytes()); + sb.extend_from_slice(&symbol_entry(0, w.header, 1, scratch)); + + debug_assert_eq!(sb.len(), SUPERBLOCK_LEN); + img.patch(0, &sb); + Ok(img.buf) +} + +/// Serialise `root` and write it to `path`. +pub fn write_file(path: &Path, root: &GroupSpec<'_>) -> Result<(), Error> { + let image = build(root).map_err(|e| Error::io(std::io::Error::other(e), path))?; + let file = std::fs::File::create(path).map_err(|e| Error::io(e, path))?; + let mut w = std::io::BufWriter::new(file); + w.write_all(&image).map_err(|e| Error::io(e, path))?; + w.flush().map_err(|e| Error::io(e, path))?; + Ok(()) +} + +#[cfg(test)] +mod tests { + use super::*; + + /// The file image starts with the format signature and declares 8-byte + /// offsets, an end-of-file address equal to its own length, and a root + /// group whose object header is inside the file. + #[test] + fn superblock_is_well_formed() { + let root = GroupSpec::new("/").dataset("x", Data::I32(&[1, 2, 3])); + let img = build(&root).unwrap(); + + assert_eq!(&img[..8], &SIGNATURE); + assert_eq!(img[8], 0, "superblock version"); + assert_eq!(img[13], 8, "size of offsets"); + assert_eq!(img[14], 8, "size of lengths"); + + let eof = u64::from_le_bytes(img[40..48].try_into().unwrap()); + assert_eq!(eof as usize, img.len(), "end-of-file address"); + + let root_hdr = u64::from_le_bytes(img[64..72].try_into().unwrap()); + assert!(root_hdr >= SUPERBLOCK_LEN as u64 && (root_hdr as usize) < img.len()); + assert_eq!(img[root_hdr as usize], 1, "object header version"); + } + + /// Every structure this writer emits is 8-byte aligned, which the format + /// requires of object headers and which readers assume of the rest. + #[test] + fn every_allocation_is_eight_byte_aligned() { + let mut img = Image::new(); + let a = img.alloc(&[1, 2, 3]); + let b = img.alloc(&[4]); + let c = img.reserve(5); + let d = img.alloc(&[]); + for addr in [a, b, c, d] { + assert_eq!(addr % 8, 0, "address {addr} is not 8-byte aligned"); + } + } + + /// A string longer than its field is rejected. Truncating a barcode would + /// produce a file that looks valid and holds a different barcode. + #[test] + fn an_overlong_string_is_an_error_not_a_truncation() { + let root = GroupSpec::new("/").dataset( + "bc", + Data::FixedString { + width: 4, + values: vec!["ACGTACGT"], + }, + ); + let err = build(&root).unwrap_err(); + assert!(err.contains("longer than"), "unexpected error: {err}"); + } + + /// Group children are written in lexicographic order, as a symbol table's + /// binary search requires, whatever order they were added in. + #[test] + fn symbol_table_entries_are_sorted_by_name() { + let root = GroupSpec::new("/") + .dataset("zeta", Data::I32(&[1])) + .dataset("alpha", Data::I32(&[2])) + .dataset("mu", Data::I32(&[3])); + let img = build(&root).unwrap(); + + // The heap data segment holds an 8-byte empty string then the names in + // the order they were written. + let names: Vec<&str> = ["alpha", "mu", "zeta"].into(); + let mut last = 0usize; + for n in names { + let pos = find(&img, n.as_bytes()).unwrap_or_else(|| panic!("{n} not in image")); + assert!(pos > last, "{n} is out of order in the local heap"); + last = pos; + } + } + + /// More children than one symbol table node holds is refused rather than + /// written as a file with entries silently dropped. + #[test] + fn too_many_children_for_one_node_is_refused() { + let mut root = GroupSpec::new("/"); + for i in 0..9 { + root.datasets.push(DatasetSpec { + name: match i { + 0 => "a", + 1 => "b", + 2 => "c", + 3 => "d", + 4 => "e", + 5 => "f", + 6 => "g", + 7 => "h", + _ => "i", + }, + data: Data::I32(&[0]), + }); + } + let err = build(&root).unwrap_err(); + assert!(err.contains("at most 8"), "unexpected error: {err}"); + } + + fn find(hay: &[u8], needle: &[u8]) -> Option { + hay.windows(needle.len()).position(|w| w == needle) + } + + /// Integers land in the file little-endian, in element order. + #[test] + fn integer_data_is_little_endian_in_element_order() { + let d = Data::I32(&[1, 256]); + assert_eq!(d.bytes().unwrap(), vec![1, 0, 0, 0, 0, 1, 0, 0]); + let d = Data::I64(&[-1]); + assert_eq!(d.bytes().unwrap(), vec![0xff; 8]); + } + + /// The canonical tiny file the golden checks below are taken from: one + /// small dataset, which stays contiguous, and one past the compression + /// threshold, which becomes a deflated chunk. Both storage paths in one + /// image. + fn golden_tree() -> Vec { + let big: Vec = (0..1000).collect(); + let root = GroupSpec::new("/") + .group( + GroupSpec::new("matrix") + .dataset("shape", Data::I32(&[3, 2])) + .dataset("data", Data::I32(&big)), + ) + .attr("filetype", AttrValue::Str("matrix".into())); + build(&root).unwrap() + } + + /// The image is byte-deterministic: no timestamps, no allocator addresses, + /// nothing that varies between runs. Determinism is a project-wide property + /// and this is the file format most likely to smuggle in a violation. + /// + /// The checksum also locks the layout. It was produced from an image that + /// HDF5 accepted — `h5dump -H`, `h5ls -r`, and an `h5repack` + `h5diff` + /// round-trip reporting no difference (`test/h5_conformance.sh`). If a + /// change to this module moves it, re-run that script before updating the + /// constant, because only libhdf5 can say whether the new bytes are still a + /// valid file. That is not hypothetical: an earlier version of this writer + /// passed every check here and every reader, and still produced a file + /// `h5repack` rewrote incorrectly. + #[test] + fn the_image_is_deterministic_and_matches_the_validated_layout() { + let a = golden_tree(); + assert_eq!(a, golden_tree(), "two builds of the same tree differ"); + assert_eq!(a.len(), 5880, "image length"); + assert_eq!(fnv1a(&a), 0x3d97_d0d0_6d6a_5a2e, "image checksum"); + } + + /// Writing the golden image out, for `test/h5_conformance.sh` to check with + /// libhdf5. Ignored by default: it produces a file rather than asserting + /// anything. + /// + /// cargo test --lib write_golden_h5 -- --ignored --nocapture + #[test] + #[ignore = "writes a file for external validation rather than asserting"] + fn write_golden_h5() { + let path = std::env::temp_dir().join("rustar_golden.h5"); + std::fs::write(&path, golden_tree()).unwrap(); + println!("wrote {}", path.display()); + } + + /// FNV-1a, for the golden check above. Not a security property: a + /// dependency-free way to notice that the bytes moved. + fn fnv1a(bytes: &[u8]) -> u64 { + let mut h: u64 = 0xcbf2_9ce4_8422_2325; + for b in bytes { + h ^= u64::from(*b); + h = h.wrapping_mul(0x1000_0000_01b3); + } + h + } + + /// Strings are padded with NULs to the field width, which is what makes + /// numpy read the dataset back as `|S`. + #[test] + fn strings_are_null_padded_to_the_field_width() { + let d = Data::FixedString { + width: 4, + values: vec!["AC", "ACGT"], + }; + assert_eq!(d.bytes().unwrap(), b"AC\0\0ACGT"); + } +} diff --git a/src/io/mod.rs b/src/io/mod.rs index a458f83..619ef0c 100644 --- a/src/io/mod.rs +++ b/src/io/mod.rs @@ -2,5 +2,6 @@ pub mod bam; pub mod fastq; +pub mod hdf5; pub mod log; pub mod sam; diff --git a/src/lib.rs b/src/lib.rs index 6fce173..29090c7 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -2172,8 +2172,13 @@ fn align_reads_solo( stats.record_alignment(0, max_multimaps); stats.record_unmapped_reason(crate::stats::UnmappedReason::Other); // No alignment → barcode still counts toward stats (unmapped → no gene). - let outcome = - solo.process_read(&[], 0, sread.barcode.as_ref(), &[]); + let outcome = solo.process_read( + &[], + 0, + sread.barcode.as_ref(), + &[], + &read.quality, + ); return Ok(SoloReadProduct { sam_records: buffer, per_feature: outcome.per_feature, @@ -2220,6 +2225,7 @@ fn align_reads_solo( transcripts.len(), sread.barcode.as_ref(), &junctions, + &read.quality, ); // Build SAM records for the cDNA alignment (same as SE path). @@ -2559,7 +2565,12 @@ fn align_reads_solo_pe( .iter() .map(|pa| (&pa.mate1_transcript, &pa.mate2_transcript)) .collect(); - solo.process_read_pe(&pairs, pread.barcode.as_ref(), &junctions) + solo.process_read_pe( + &pairs, + pread.barcode.as_ref(), + &junctions, + &pread.mate1.quality, + ) } else if let Some(PairedAlignmentResult::HalfMapped { mapped_transcript, .. @@ -2570,9 +2581,16 @@ fn align_reads_solo_pe( 1, pread.barcode.as_ref(), &junctions, + &pread.mate1.quality, ) } else { - solo.process_read(&[], 0, pread.barcode.as_ref(), &[]) + solo.process_read( + &[], + 0, + pread.barcode.as_ref(), + &[], + &pread.mate1.quality, + ) }; // SAM records (skipped under `--outSAMtype None`). diff --git a/src/params/mod.rs b/src/params/mod.rs index 0536a85..edfa3c9 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1133,6 +1133,44 @@ pub struct Parameters { #[arg(long = "soloOutGzip", default_value = "no")] pub solo_out_gzip: String, + /// Which barcodes the **raw** matrix has columns for. **Not a STAR + /// parameter**; a rustar-aligner addition, default `Whitelist`, which is + /// what STARsolo writes. + /// + /// `Whitelist` gives one column per whitelist barcode — 3.7 million of them + /// for 10x v3, whether or not a read ever carried them. `Observed` gives + /// one column per barcode that actually holds a count, which is what + /// CellRanger's `raw_feature_bc_matrix` contains, and turns a + /// hundreds-of-megabytes `barcodes.tsv` into a few kilobytes. + /// + /// The counts are identical either way; only the columns present differ. + #[arg(long = "soloOutRawBarcodes", default_value = "Whitelist", + value_parser = ["Whitelist", "Observed"])] + pub solo_out_raw_barcodes: String, + + /// Shape of the solo matrix output on disk. **Not a STAR parameter**; a + /// rustar-aligner addition, default `STARsolo`, which changes nothing. + /// + /// `CellRanger` lays the same numbers out the way `cellranger count` does, + /// so a tool written against CellRanger's `outs/` reads a rustar run + /// unmodified. It implies, unless the corresponding flag is given + /// explicitly: + /// + /// * directories `raw_feature_bc_matrix/` and `filtered_feature_bc_matrix/` + /// instead of `raw/` and `filtered/`; + /// * a `-1` GEM-well suffix on every barcode; + /// * gzip on all three files (`--soloOutGzip yes`); + /// * `--soloOutRawBarcodes Observed`, since CellRanger's raw matrix has one + /// column per observed barcode, not one per whitelist entry; + /// * the output directory `outs/` instead of `Solo.out/`, with no + /// per-feature subdirectory when exactly one feature is requested. + /// + /// Counts are untouched. Only where the bytes land, and how the barcodes + /// are spelled, changes. + #[arg(long = "soloOutLayout", default_value = "STARsolo", + value_parser = ["STARsolo", "CellRanger"])] + pub solo_out_layout: String, + /// Velocyto ambiguous-molecule handling (rustar extension beyond STARsolo). /// `yes` (default) writes the three `spliced`/`unspliced`/`ambiguous` matrices /// like STARsolo — exon-only molecules with no junction/intron evidence stay in @@ -1348,6 +1386,9 @@ impl Parameters { let matches = command.clone().get_matches_from(args.iter()); let mut params = ::from_arg_matches(&matches)?; + apply_cellranger_defaults_on_10x(&mut params, &matches); + apply_cellranger_layout(&mut params, &matches); + params.command_line = { let args: Vec<_> = args.iter().map(|s| s.to_string_lossy()).collect(); shlex::try_join(args.iter().map(AsRef::as_ref)).ok() @@ -1870,8 +1911,351 @@ impl Parameters { // Tests // --------------------------------------------------------------------------- +/// The flags STAR documents for matching CellRanger 4.x/5.x +/// (`docs/STARsolo.md`), applied by default when the run is a 10x one. +const CELLRANGER_DEFAULTS: [(&str, &str); 6] = [ + ("clip_adapter_type", "CellRanger4"), + ("out_filter_score_min", "30"), + ("solo_cb_match_wl_type", "1MM_multi_Nbase_pseudocounts"), + ("solo_umi_filtering", "MultiGeneUMI_CR"), + ("solo_umi_dedup", "1MM_CR"), + // Not a STAR flag: the output layout, so a 10x run lands where a tool + // written against `cellranger count` expects to find it. + ("solo_out_layout", "CellRanger"), +]; + +/// Does this look like a 10x Chromium run? +/// +/// `CB_UMI_Simple` with a whitelist, a 16-base cell barcode, and a UMI of 10 +/// (v2) or 12 (v3) bases. That is the geometry of every 10x 3'/5' gene +/// expression chemistry, and nothing else in common use shares it. +fn looks_like_10x(params: &Parameters) -> bool { + params.solo_type == SoloType::CbUmiSimple + && params.solo_cb_len == 16 + && (params.solo_umi_len == 10 || params.solo_umi_len == 12) + && params + .solo_cb_whitelist + .first() + .is_some_and(|w| w != "None" && w != "-") +} + +/// On a 10x run, default to CellRanger's behaviour rather than STARsolo's. +/// +/// **This diverges from STAR by default**, which is why it is confined to a +/// geometry that is unambiguously 10x, and why every flag it changes is +/// logged. A flag given on the command line always wins, so the change is +/// invisible to anyone who states what they want. +/// +/// The rationale is that a user aligning 10x data and comparing against +/// CellRanger currently gets a successful run and different numbers, with +/// nothing pointing at the five flags that explain the difference. Measured on +/// a 20 000-read fixture, those flags move the count matrix from 8.9% away +/// from CellRanger to 0.03%. +/// +/// Recorded in `DIVERGENCE.md`; it needs maintainer sign-off. +fn apply_cellranger_defaults_on_10x(params: &mut Parameters, matches: &clap::ArgMatches) { + use clap::parser::ValueSource; + + if !looks_like_10x(params) { + return; + } + + let given = |id: &str| matches.value_source(id) == Some(ValueSource::CommandLine); + + let mut applied: Vec<&str> = Vec::new(); + for (id, value) in CELLRANGER_DEFAULTS { + if given(id) { + continue; + } + match id { + "clip_adapter_type" => params.clip_adapter_type = value.to_string(), + "out_filter_score_min" => params.out_filter_score_min = 30, + "solo_cb_match_wl_type" => params.solo_cb_match_wl_type = value.to_string(), + "solo_umi_filtering" => params.solo_umi_filtering = vec![value.to_string()], + "solo_umi_dedup" => params.solo_umi_dedup = vec![value.to_string()], + "solo_out_layout" => params.solo_out_layout = value.to_string(), + _ => continue, + } + applied.push(value); + } + + if !applied.is_empty() { + log::info!( + "10x geometry detected (CB {} + UMI {} with a whitelist): defaulting to \ + CellRanger behaviour [{}]. Pass the flags explicitly to override; this \ + differs from STARsolo's defaults.", + params.solo_cb_len, + params.solo_umi_len, + applied.join(", ") + ); + } +} + +/// `--soloOutLayout CellRanger` implies three existing flags and the output +/// directory name. Each is still overridable: a flag named on the command line +/// keeps its value, so the layout can be adopted piecemeal. +/// +/// Runs after `apply_cellranger_defaults_on_10x`, so it sees the layout whether +/// it was asked for or inferred from 10x geometry. +fn apply_cellranger_layout(params: &mut Parameters, matches: &clap::ArgMatches) { + use clap::parser::ValueSource; + + if params.solo_out_layout != "CellRanger" { + return; + } + let given = |id: &str| matches.value_source(id) == Some(ValueSource::CommandLine); + + if !given("solo_out_gzip") { + params.solo_out_gzip = "yes".to_string(); + } + if !given("solo_out_raw_barcodes") { + params.solo_out_raw_barcodes = "Observed".to_string(); + } + if !given("solo_out_file_names") + && let Some(dir) = params.solo_out_file_names.first_mut() + { + // CellRanger writes its matrices under `outs/`, not `Solo.out/`. + *dir = "outs/".to_string(); + } +} + #[cfg(test)] mod tests { + + /// 10x geometry with a whitelist gets CellRanger's five flags without the + /// user naming any of them. This is a deliberate divergence from STARsolo's + /// defaults, so the test states the whole set rather than spot-checking one. + #[test] + fn ten_x_geometry_defaults_to_cellranger_behaviour() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--soloCBstart", + "1", + "--soloCBlen", + "16", + "--soloUMIstart", + "17", + "--soloUMIlen", + "12", + ]) + .unwrap(); + assert_eq!(p.clip_adapter_type, "CellRanger4"); + assert_eq!(p.out_filter_score_min, 30); + assert_eq!(p.solo_cb_match_wl_type, "1MM_multi_Nbase_pseudocounts"); + assert_eq!(p.solo_umi_filtering, vec!["MultiGeneUMI_CR".to_string()]); + assert_eq!(p.solo_umi_dedup, vec!["1MM_CR".to_string()]); + assert_eq!(p.solo_out_layout, "CellRanger"); + } + + /// `--soloOutLayout CellRanger` pulls three existing flags and the output + /// directory with it, so the layout is one decision rather than four. + #[test] + fn cellranger_layout_implies_gzip_observed_barcodes_and_outs_dir() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--soloCBstart", + "1", + "--soloCBlen", + "12", + "--soloUMIstart", + "13", + "--soloUMIlen", + "10", + "--soloOutLayout", + "CellRanger", + ]) + .unwrap(); + // 12-base CB, so the 10x autodetection is not what set this. + assert_eq!(p.solo_out_layout, "CellRanger"); + assert_eq!(p.solo_out_gzip, "yes"); + assert_eq!(p.solo_out_raw_barcodes, "Observed"); + assert_eq!(p.solo_out_file_names.first().unwrap(), "outs/"); + } + + /// The layout is adoptable piecemeal: each flag it implies is still + /// overridable on the command line. + #[test] + fn an_explicit_flag_beats_what_the_cellranger_layout_implies() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--soloCBstart", + "1", + "--soloCBlen", + "16", + "--soloUMIstart", + "17", + "--soloUMIlen", + "12", + "--soloOutGzip", + "no", + "--soloOutRawBarcodes", + "Whitelist", + "--soloOutFileNames", + "Solo.out/", + "features.tsv", + "barcodes.tsv", + "matrix.mtx", + ]) + .unwrap(); + assert_eq!(p.solo_out_layout, "CellRanger"); + assert_eq!(p.solo_out_gzip, "no"); + assert_eq!(p.solo_out_raw_barcodes, "Whitelist"); + assert_eq!(p.solo_out_file_names.first().unwrap(), "Solo.out/"); + } + + /// The default is STARsolo's layout: no `-1`, no `outs/`, no gzip. + #[test] + fn non_10x_geometry_keeps_the_starsolo_layout() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--soloCBstart", + "1", + "--soloCBlen", + "12", + "--soloUMIstart", + "13", + "--soloUMIlen", + "10", + ]) + .unwrap(); + assert_eq!(p.solo_out_layout, "STARsolo"); + assert_eq!(p.solo_out_gzip, "no"); + assert_eq!(p.solo_out_raw_barcodes, "Whitelist"); + assert_eq!(p.solo_out_file_names.first().unwrap(), "Solo.out/"); + } + + /// A flag given on the command line always wins, including when the value + /// asked for is STARsolo's own default. Without this the divergence would + /// be inescapable, which is a different and much worse thing than a + /// divergent default. + #[test] + fn an_explicit_flag_beats_the_10x_default() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--soloCBstart", + "1", + "--soloCBlen", + "16", + "--soloUMIstart", + "17", + "--soloUMIlen", + "12", + "--soloUMIdedup", + "1MM_All", + "--clipAdapterType", + "Hamming", + ]) + .unwrap(); + assert_eq!(p.solo_umi_dedup, vec!["1MM_All".to_string()]); + assert_eq!(p.clip_adapter_type, "Hamming"); + // The ones not named still take the CellRanger value. + assert_eq!(p.out_filter_score_min, 30); + } + + /// Geometry that is not 10x is left alone: a 12-base barcode is not any + /// Chromium chemistry, so nothing is overridden. + #[test] + fn non_10x_geometry_keeps_starsolo_defaults() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--soloCBstart", + "1", + "--soloCBlen", + "12", + "--soloUMIstart", + "13", + "--soloUMIlen", + "8", + ]) + .unwrap(); + assert_eq!(p.clip_adapter_type, "Hamming"); + assert_eq!(p.out_filter_score_min, 0); + assert_eq!(p.solo_cb_match_wl_type, "1MM_multi"); + } + + /// No whitelist means no 10x run, whatever the lengths say. (Without a + /// whitelist the CB-match type must be Exact anyway, which is unrelated + /// validation that predates this and is stated here so the test reads.) + #[test] + fn ten_x_lengths_without_a_whitelist_keep_starsolo_defaults() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBstart", + "1", + "--soloCBlen", + "16", + "--soloUMIstart", + "17", + "--soloUMIlen", + "12", + "--soloCBmatchWLtype", + "Exact", + ]) + .unwrap(); + assert_eq!(p.clip_adapter_type, "Hamming"); + assert_eq!(p.out_filter_score_min, 0); + } use super::*; /// Helper: parse a STAR-style command line (without program name). diff --git a/src/solo/count.rs b/src/solo/count.rs index 9af3edf..d3168fe 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -167,6 +167,18 @@ pub fn dedup_count(umis: &HashMap, method: UmiDedup, umi_len: usize) - /// the neighbor's raw UMI, not its corrected value); the molecule count is the /// number of distinct corrected UMIs. fn cellranger_1mm(umis: &HashMap, umi_len: usize) -> u64 { + let distinct: std::collections::HashSet = + cellranger_1mm_map(umis, umi_len).into_values().collect(); + distinct.len() as u64 +} + +/// The same correction, returning `raw UMI -> corrected UMI`. +/// +/// `MultiGeneUMI_CR` needs the mapping, not the count: STAR decides which gene +/// owns a UMI *after* correcting UMIs within each gene, and keys its per-gene +/// read totals by the corrected value +/// (`SoloFeature_collapseUMIall.cpp:134-148`). +fn cellranger_1mm_map(umis: &HashMap, umi_len: usize) -> HashMap { let mut items: Vec<(u64, u32)> = umis.iter().map(|(&u, &c)| (u, c)).collect(); // Ascending by count, then by UMI value (mirrors funCompareSolo1 ordering, // so the inner scan from the end meets higher-count neighbors first). @@ -185,8 +197,7 @@ fn cellranger_1mm(umis: &HashMap, umi_len: usize) -> u64 { } corrected.push(corr); } - let distinct: std::collections::HashSet = corrected.into_iter().collect(); - distinct.len() as u64 + items.iter().map(|&(u, _)| u).zip(corrected).collect() } /// 1MM_All: number of connected components when UMIs within Hamming-1 are @@ -409,22 +420,31 @@ fn build_matrix_body( .or_insert(0) += 1; } - // (gene → (umi → read_count)) after multi-gene UMI filtering. - let mut gene_umis: HashMap> = HashMap::default(); - for (&umi, genes) in &umi_genes { - for (&gene, &rc) in filter_multi_gene_umi(genes, filtering) { - *gene_umis.entry(gene).or_default().entry(umi).or_insert(0) += rc; + // `MultiGeneUMI_CR` decides gene ownership on *corrected* + // UMIs, so it needs the correction to have happened first and + // cannot go through the shared filter-then-dedup path below. + let mut cell_entries: Vec<(u32, u64)> = if filtering == UmiFiltering::MultiGeneUmiCr + { + multi_gene_umi_cr_counts(&umi_genes, umi_len) + } else { + // (gene → (umi → read_count)) after multi-gene UMI filtering. + let mut gene_umis: HashMap> = HashMap::default(); + for (&umi, genes) in &umi_genes { + for (&gene, &rc) in filter_multi_gene_umi(genes, filtering) { + *gene_umis.entry(gene).or_default().entry(umi).or_insert(0) += rc; + } } - } - // Collapse UMIs per gene, then emit this cell's entries gene-ascending. - let mut cell_entries: Vec<(u32, u64)> = Vec::with_capacity(gene_umis.len()); - for (&gene, umis) in &gene_umis { - let count = dedup_count(umis, method, umi_len); - if count > 0 { - cell_entries.push((gene, count)); + // Collapse UMIs per gene, then emit gene-ascending. + let mut entries: Vec<(u32, u64)> = Vec::with_capacity(gene_umis.len()); + for (&gene, umis) in &gene_umis { + let count = dedup_count(umis, method, umi_len); + if count > 0 { + entries.push((gene, count)); + } } - } + entries + }; cell_entries.sort_unstable_by_key(|&(g, _)| g); let n_reads = (j - i) as u64; @@ -476,6 +496,197 @@ fn build_matrix_body( )) } +/// The whitelist indices that actually appear as a column in the streamed +/// matrix body, ascending. +/// +/// Reads the body once rather than tracking the set during counting, so the +/// default path pays nothing for a feature it does not use. +fn observed_barcodes(body: &tempfile::NamedTempFile) -> Result, Error> { + let reader = + BufReader::new(std::fs::File::open(body.path()).map_err(|e| Error::io(e, body.path()))?); + let mut seen: std::collections::BTreeSet = std::collections::BTreeSet::new(); + for line in reader.lines() { + let line = line.map_err(|e| Error::io(e, body.path()))?; + // " ", the layout `finalize_matrix` also parses. + if let Some(cb1) = line.split(' ').nth(1) + && let Ok(cb) = cb1.parse::() + { + seen.insert(cb.saturating_sub(1)); + } + } + Ok(seen.into_iter().collect()) +} + +/// The barcode strings for a set of whitelist indices, with the layout's +/// suffix, in the same order as the matrix columns. +/// +/// The MEX writers stream these straight to disk; the HDF5 writer needs them in +/// memory because a fixed-width dataset has to know its widest value first. +fn barcode_strings(whitelist: &CbWhitelist, cbs: &[u32], suffix: &str) -> Vec { + let mut line: Vec = Vec::with_capacity(whitelist.barcode_len() + suffix.len()); + cbs.iter() + .map(|&cb| { + line.clear(); + whitelist.unpack_barcode_into(cb, &mut line); + let mut s = String::from_utf8_lossy(&line).into_owned(); + s.push_str(suffix); + s + }) + .collect() +} + +/// Write CellRanger's `*_feature_bc_matrix.h5` beside the MEX directory. +/// +/// Same numbers, same column set, in the HDF5 layout `scanpy.read_10x_h5`, +/// Seurat's `Read10X_h5` and CellBender read: a `/matrix` group holding the CSC +/// triplet (`data`, `indices`, `indptr`, `shape`), the barcodes, and a +/// `/matrix/features` group. Written through `crate::io::hdf5`, an in-tree +/// writer, so this costs no dependency (see `DIVERGENCE.md` §3.5). +/// +/// CellRanger's matrix is genes × cells stored column-major, one column per +/// barcode. The streamed body is already in that order — ascending cell, then +/// ascending gene within a cell — so the column pointers are built in one pass +/// rather than by sorting `nnz` entries in memory. That ordering is a property +/// of `build_matrix_body`, so it is checked here rather than assumed. +#[allow(clippy::too_many_arguments)] +fn write_matrix_h5( + body: &tempfile::NamedTempFile, + out_path: &Path, + n_genes: usize, + barcodes: &[String], + gene_ids: &[String], + gene_names: &[String], + genome: &str, + remap: Option<&HashMap>, +) -> Result<(), Error> { + use crate::io::hdf5::{AttrValue, Data, GroupSpec}; + + let mut data: Vec = Vec::new(); + let mut indices: Vec = Vec::new(); + // indptr[c] = index into data/indices where column c starts. + let mut indptr: Vec = vec![0; 1]; + + let reader = + BufReader::new(std::fs::File::open(body.path()).map_err(|e| Error::io(e, body.path()))?); + let mut current_col: u32 = 0; + for line in reader.lines() { + let line = line.map_err(|e| Error::io(e, body.path()))?; + let mut it = line.split(' '); + let (Some(g), Some(c), Some(v)) = (it.next(), it.next(), it.next()) else { + continue; + }; + let gene0: i64 = g.parse::().unwrap_or(1) - 1; + let cb0: u32 = c.parse::().unwrap_or(0).saturating_sub(1); + let col = match remap { + Some(map) => match map.get(&cb0) { + Some(&col1) => col1 - 1, + None => continue, + }, + None => cb0, + }; + if col < current_col { + return Err(Error::io( + std::io::Error::other(format!( + "matrix body is not in ascending column order (column {col} after {current_col}); \ + the HDF5 writer builds column pointers in one pass and relies on that order" + )), + out_path, + )); + } + // Close every column between the last one and this one, including the + // empty columns in between. + while current_col < col { + indptr.push(data.len() as i64); + current_col += 1; + } + indices.push(gene0); + data.push(v.parse::().unwrap_or(0)); + } + // Close the final column and every trailing empty one. + while (indptr.len() as u32) <= barcodes.len() as u32 { + indptr.push(data.len() as i64); + } + + let bc_refs: Vec<&str> = barcodes.iter().map(String::as_str).collect(); + let id_refs: Vec<&str> = gene_ids.iter().map(String::as_str).collect(); + let name_refs: Vec<&str> = gene_names.iter().map(String::as_str).collect(); + let type_refs: Vec<&str> = vec!["Gene Expression"; n_genes]; + let genome_refs: Vec<&str> = vec![genome; n_genes]; + let tag_refs: Vec<&str> = vec!["genome"]; + let shape: Vec = vec![n_genes as i32, barcodes.len() as i32]; + + // Field widths are the longest value present, which is what makes numpy + // report the dataset as `|S`. CellRanger pads to a fixed 256; the + // values are identical either way and this keeps the file smaller. + let width = |v: &[&str]| v.iter().map(|s| s.len()).max().unwrap_or(1).max(1); + + let features = GroupSpec::new("features") + .dataset( + "_all_tag_keys", + Data::FixedString { + width: width(&tag_refs), + values: tag_refs, + }, + ) + .dataset( + "feature_type", + Data::FixedString { + width: width(&type_refs), + values: type_refs, + }, + ) + .dataset( + "genome", + Data::FixedString { + width: width(&genome_refs), + values: genome_refs, + }, + ) + .dataset( + "id", + Data::FixedString { + width: width(&id_refs), + values: id_refs, + }, + ) + .dataset( + "name", + Data::FixedString { + width: width(&name_refs), + values: name_refs, + }, + ); + + let matrix = GroupSpec::new("matrix") + .dataset( + "barcodes", + Data::FixedString { + width: width(&bc_refs), + values: bc_refs, + }, + ) + .dataset("data", Data::I32(&data)) + .dataset("indices", Data::I64(&indices)) + .dataset("indptr", Data::I64(&indptr)) + .dataset("shape", Data::I32(&shape)) + .group(features); + + // `filetype` and `version` are what a reader dispatches on. The rest is + // provenance, and says rustar-aligner rather than claiming to be + // CellRanger. + let root = GroupSpec::new("/") + .group(matrix) + .attr("filetype", AttrValue::Str("matrix".to_string())) + .attr("version", AttrValue::I64(2)) + .attr( + "software_version", + AttrValue::Str(format!("rustar-aligner {}", env!("CARGO_PKG_VERSION"))), + ) + .attr("original_gem_groups", AttrValue::I64Array(vec![1])); + + crate::io::hdf5::write_file(out_path, &root) +} + /// Write a final `matrix.mtx[.gz]` = MatrixMarket header + (optionally /// cb-remapped/filtered) body. With `remap = None` the body is copied verbatim /// (raw); with `Some(map)` only columns in the map survive, renumbered to the @@ -808,6 +1019,88 @@ fn build_multi_matrices( Ok(()) } +/// CellRanger's multi-gene UMI resolution, as STAR implements it for +/// `--soloUMIfiltering MultiGeneUMI_CR` (`SoloFeature_collapseUMIall.cpp`). +/// +/// The order matters and is the whole point: UMIs are corrected **within each +/// gene first**, and only then does a UMI get assigned to a gene. Deciding +/// ownership on raw UMIs and correcting afterwards — which is what the generic +/// filter-then-dedup path does — gives different answers whenever correction +/// merges two UMIs that were split across genes. +/// +/// Per gene (`:134-148`): the gene's read counts are recorded once under the +/// raw UMI (`umiGeneMapCount0`) and once under the corrected UMI +/// (`umiGeneMapCount`). +/// +/// Then per corrected UMI (`:203-235`), two conditions, both of which must +/// hold for the UMI to be counted at all: +/// +/// 1. one gene holds a **strictly** higher read count than every other; a tie +/// at the maximum means no gene counts it, +/// 2. and no gene beats that winner in the **uncorrected** map at the same key. +/// +/// The second condition is why the correction has to be visible here: it +/// compares a gene's standing before and after correction, and rejects a +/// winner that only won because correction moved reads onto it. +/// +/// Returns `(gene, molecules)` for this cell, gene-ascending. +fn multi_gene_umi_cr_counts( + umi_genes: &HashMap>, + umi_len: usize, +) -> Vec<(u32, u64)> { + // Regroup as gene → (raw UMI → reads); correction happens per gene. + let mut gene_umis: HashMap> = HashMap::default(); + for (&umi, genes) in umi_genes { + for (&gene, &rc) in genes { + *gene_umis.entry(gene).or_default().entry(umi).or_insert(0) += rc; + } + } + + let mut uncorrected: HashMap> = HashMap::default(); + let mut corrected: HashMap> = HashMap::default(); + for (&gene, umis) in &gene_umis { + for (&umi, &rc) in umis { + *uncorrected.entry(umi).or_default().entry(gene).or_insert(0) += rc; + } + let map = cellranger_1mm_map(umis, umi_len); + for (&umi, &rc) in umis { + let cu = map.get(&umi).copied().unwrap_or(umi); + *corrected.entry(cu).or_default().entry(gene).or_insert(0) += rc; + } + } + + let mut counts: HashMap = HashMap::default(); + for (cu, genes) in &corrected { + // Condition 1: a strict maximum, ties lose. + let mut best = 0u32; + let mut winner: Option = None; + for (&gene, &rc) in genes { + if rc > best { + best = rc; + winner = Some(gene); + } else if rc == best { + winner = None; + } + } + let Some(winner) = winner else { continue }; + + // Condition 2: the winner must not be beaten in the uncorrected map at + // the same key. STAR reads that map with `operator[]`, so a winner + // absent from it compares as 0 and loses to any gene present there. + if let Some(raw_genes) = uncorrected.get(cu) { + let winner_raw = raw_genes.get(&winner).copied().unwrap_or(0); + if raw_genes.values().any(|&rc| rc > winner_raw) { + continue; + } + } + *counts.entry(winner).or_insert(0) += 1; + } + + let mut out: Vec<(u32, u64)> = counts.into_iter().filter(|&(_, c)| c > 0).collect(); + out.sort_unstable_by_key(|&(g, _)| g); + out +} + /// Apply `--soloUMIfiltering` to the gene→read_count map of a single UMI, /// returning the surviving (gene, read_count) entries. fn filter_multi_gene_umi(genes: &HashMap, filtering: UmiFiltering) -> Vec<(&u32, &u32)> { @@ -822,8 +1115,43 @@ fn filter_multi_gene_umi(genes: &HashMap, filtering: UmiFiltering) -> let thresh = if max == 1 { 2 } else { max }; genes.iter().filter(|&(_, &rc)| rc >= thresh).collect() } - // CellRanger > 3.0: keep the highest-read-count gene(s); no singleton drop. - UmiFiltering::MultiGeneUmiCr => genes.iter().filter(|&(_, &rc)| rc >= max).collect(), + // CellRanger: the gene with the strictly highest read count takes the + // UMI, and a tie gives it to nobody. + // + // STAR `SoloFeature_collapseUMIall.cpp:212-224` walks the genes keeping + // a running maximum, and clears its winner whenever it meets an equal + // count: + // + // ```cpp + // if (ig.second>maxu) { maxu=ig.second; maxg=ig.first; } + // else if (ig.second==maxu) { maxg=-1; }; + // ... + // if ( maxg+1==0 ) continue; // not counted for any gene + // ``` + // + // The outcome does not depend on the order the genes are visited: a + // strict maximum always ends as the winner, and any tie at the maximum + // always ends with none. So iterating a `HashMap` here is safe. + // + // This previously kept every gene tied at the maximum, which is the + // opposite decision on exactly the case the rule exists for, and made + // the flag inert on the common shape of one read per gene. + UmiFiltering::MultiGeneUmiCr => { + let mut best_count = 0u32; + let mut winner: Option<&u32> = None; + for (gene, &rc) in genes { + if rc > best_count { + best_count = rc; + winner = Some(gene); + } else if rc == best_count { + winner = None; + } + } + match winner { + Some(gene) => vec![(gene, genes.get(gene).expect("winner is a key"))], + None => Vec::new(), + } + } UmiFiltering::None => unreachable!(), } } @@ -1163,8 +1491,12 @@ pub fn write_gene_matrix( .position(|&x| x == f) .map_or(0, |i| ctx.feature_reads[i].load(Ordering::Relaxed)) }; - let have_funnel = ctx.features.contains(&crate::solo::SoloFeature::Gene) - && ctx.features.contains(&crate::solo::SoloFeature::GeneFull); + // The split needs the gene-body query as well as the exon one. Both `Gene` + // + `GeneFull` gives it, and so does `ctx.want_metrics`, which asks + // `process_read` for the body query on its own account. + let have_funnel = ctx.want_metrics + || (ctx.features.contains(&crate::solo::SoloFeature::Gene) + && ctx.features.contains(&crate::solo::SoloFeature::GeneFull)); let region = have_funnel.then(|| RegionFunnel { exonic: ctx.region_stats.exonic.load(Ordering::Relaxed), intronic: ctx.region_stats.intronic.load(Ordering::Relaxed), @@ -1176,10 +1508,37 @@ pub fn write_gene_matrix( let n_genes = ctx.gene_ann.gene_ids.len(); let multi_methods = MultiMethod::parse_list(¶ms.solo_multi_mappers); + // `--soloOutLayout CellRanger`: CellRanger's directory names, and a `-1` + // GEM-well suffix on every barcode. The counts are the same either way. + let cr_layout = params.solo_out_layout == "CellRanger"; + let (raw_name, filt_name) = if cr_layout { + ("raw_feature_bc_matrix", "filtered_feature_bc_matrix") + } else { + ("raw", "filtered") + }; + let cb_suffix = if cr_layout { "-1" } else { "" }; + // The `.h5` feature table carries a genome name per feature. CellRanger + // uses the name given to `cellranger mkref`; the closest thing we have is + // the index directory, which is what `--genomeDir` names. + let genome_name = std::path::Path::new(¶ms.genome_dir) + .file_name() + .map_or_else( + || "genome".to_string(), + |s| s.to_string_lossy().into_owned(), + ); + // CellRanger has no per-feature directory. Drop ours when there is exactly + // one feature; keep it when there are several, because collapsing them + // would have each feature silently overwrite the last. + let per_feature_dir = !cr_layout || ctx.features.len() > 1; + // One {prefix}{soloOutFileNames[0]}/{raw,filtered}/ per feature. for (feature, recorder) in ctx.features.iter().zip(&ctx.recorders) { - let feature_dir = params.output_path(&format!("{solo_dir}{}/", feature.dir_name())); - let raw_dir = feature_dir.join("raw"); + let feature_dir = if per_feature_dir { + params.output_path(&format!("{solo_dir}{}/", feature.dir_name())) + } else { + params.output_path(&solo_dir) + }; + let raw_dir = feature_dir.join(raw_name); std::fs::create_dir_all(&raw_dir).map_err(|e| Error::io(e, &raw_dir))?; // Stream the deduplicated counts into a shared temp body, then finalize @@ -1200,26 +1559,91 @@ pub fn write_gene_matrix( &ctx.gene_ann.gene_names, gzip, )?; - write_barcodes( - &raw_dir.join(&barcodes_name), - &ctx.whitelist, - sorted.len(), - gzip, - )?; + // `--soloOutRawBarcodes Observed` narrows the raw matrix to the + // barcodes that actually carry a count, which is what CellRanger's + // `raw_feature_bc_matrix` holds. STARsolo's raw matrix has a column per + // whitelist barcode, so the default keeps that. + let observed: Option> = if params.solo_out_raw_barcodes == "Observed" { + Some(observed_barcodes(&body)?) + } else { + None + }; + let (raw_cols, raw_remap) = match &observed { + Some(cbs) => { + let map: HashMap = cbs + .iter() + .enumerate() + .map(|(col, &cb)| (cb, col as u32 + 1)) + .collect(); + (cbs.len(), Some(map)) + } + None => (sorted.len(), None), + }; + match &observed { + Some(cbs) => { + write_barcodes_subset( + &raw_dir.join(&barcodes_name), + &ctx.whitelist, + cbs, + gzip, + cb_suffix, + )?; + } + None => write_barcodes( + &raw_dir.join(&barcodes_name), + &ctx.whitelist, + sorted.len(), + gzip, + cb_suffix, + )?, + } finalize_matrix( &body, &raw_dir.join(&matrix_name), gzip, n_genes, - sorted.len(), + raw_cols, mstats.nnz, - None, + raw_remap.as_ref(), )?; + // CellRanger writes the same matrix twice: the MEX directory above and + // an .h5 beside it, which is what read_10x_h5-style loaders open. + if cr_layout { + let cbs: Vec = match &observed { + Some(cbs) => cbs.clone(), + None => (0..sorted.len() as u32).collect(), + }; + // The HDF5 writer holds the barcode strings in memory to size a + // fixed-width dataset, so a whitelist-wide raw matrix is costly. + // Only reachable by pairing --soloOutLayout CellRanger with an + // explicit --soloOutRawBarcodes Whitelist; the layout's own default + // is Observed. + if cbs.len() > 1_000_000 { + log::warn!( + "STARsolo: writing {}.h5 with {} barcode columns; \ + --soloOutRawBarcodes Observed would make it far smaller", + raw_name, + cbs.len() + ); + } + let h5_path = feature_dir.join(format!("{raw_name}.h5")); + write_matrix_h5( + &body, + &h5_path, + n_genes, + &barcode_strings(&ctx.whitelist, &cbs, cb_suffix), + &ctx.gene_ann.gene_ids, + &ctx.gene_ann.gene_names, + &genome_name, + raw_remap.as_ref(), + )?; + log::info!("STARsolo: wrote {}", h5_path.display()); + } log::info!( - "STARsolo: wrote {}/raw matrix ({} genes × {} barcodes, {} entries){}", - feature.dir_name(), + "STARsolo: wrote {} matrix ({} genes × {} barcodes, {} entries){}", + raw_dir.display(), n_genes, - sorted.len(), + raw_cols, mstats.nnz, if gzip { " [gzip]" } else { "" }, ); @@ -1243,7 +1667,7 @@ pub fn write_gene_matrix( if let Some(cbs) = called && !cbs.is_empty() { - let filt_dir = feature_dir.join("filtered"); + let filt_dir = feature_dir.join(filt_name); std::fs::create_dir_all(&filt_dir).map_err(|e| Error::io(e, &filt_dir))?; let remap: HashMap = cbs .iter() @@ -1256,7 +1680,13 @@ pub fn write_gene_matrix( &ctx.gene_ann.gene_names, gzip, )?; - write_barcodes_subset(&filt_dir.join(&barcodes_name), &ctx.whitelist, &cbs, gzip)?; + write_barcodes_subset( + &filt_dir.join(&barcodes_name), + &ctx.whitelist, + &cbs, + gzip, + cb_suffix, + )?; let fnnz = finalize_matrix( &body, &filt_dir.join(&matrix_name), @@ -1267,11 +1697,25 @@ pub fn write_gene_matrix( Some(&remap), )?; log::info!( - "STARsolo: wrote {}/filtered matrix ({} cells, {} entries)", - feature.dir_name(), + "STARsolo: wrote {} matrix ({} cells, {} entries)", + filt_dir.display(), cbs.len(), fnnz, ); + if cr_layout { + let h5_path = feature_dir.join(format!("{filt_name}.h5")); + write_matrix_h5( + &body, + &h5_path, + n_genes, + &barcode_strings(&ctx.whitelist, &cbs, cb_suffix), + &ctx.gene_ann.gene_ids, + &ctx.gene_ann.gene_names, + &genome_name, + Some(&remap), + )?; + log::info!("STARsolo: wrote {}", h5_path.display()); + } } // --soloMultiMappers: UniqueAndMult-.mtx alongside raw. @@ -1306,10 +1750,32 @@ pub fn write_gene_matrix( reads_of(*feature), )?; log::info!("STARsolo: wrote {}/Summary.csv", feature.dir_name()); + // Under the CellRanger layout the funnel goes into `metrics_summary.csv` + // alongside the other 19 metrics, which is where CellRanger puts it. + if cr_layout && let Some(r) = region { + let invalid_umis = ctx.stats.n_in_umi.load(Ordering::Relaxed) + + ctx.stats.umi_homopolymer.load(Ordering::Relaxed); + write_metrics_summary( + &feature_dir.join("metrics_summary.csv"), + &mstats, + &ctx.q30, + total_reads, + valid_barcodes, + invalid_umis, + mapped_unique, + mapped_multi, + reads_of(crate::solo::SoloFeature::Gene), + r, + )?; + log::info!( + "STARsolo: wrote {}", + feature_dir.join("metrics_summary.csv").display() + ); + } // CellRanger-style mapping funnel goes in a SEPARATE additional file so the // faithful Summary.csv is never altered (PR #90 review: keep this release a // drop-in faithful port; output-changing features come later). - if let Some(r) = region { + if !cr_layout && let Some(r) = region { write_cellranger_summary( &feature_dir.join("CellRanger.summary.csv"), total_reads, @@ -1339,11 +1805,14 @@ pub fn write_gene_matrix( write_file(&sj_dir.join(&features_name), gzip, |w| { sjs.write_sj_lines(w, genome, params).map(|_| ()) })?; + // SJ has no CellRanger counterpart, but every barcode written by one run + // is spelled the same way, so the suffix applies here too. write_barcodes( &sj_dir.join(&barcodes_name), &ctx.whitelist, sorted.len(), gzip, + cb_suffix, )?; let umi_len = params.solo_umi_len as usize; let nnz = build_sj_matrix( @@ -1379,6 +1848,7 @@ pub fn write_gene_matrix( &ctx.whitelist, sorted.len(), gzip, + cb_suffix, )?; let umi_len = params.solo_umi_len as usize; // `--soloVelocytoAmbiguous no` folds exon-only molecules into spliced and @@ -1640,6 +2110,168 @@ struct RegionFunnel { antisense: u64, } +/// Format an integer the way CellRanger's `metrics_summary.csv` does: thousands +/// separated by commas, and quoted when that puts a comma in the field. +/// +/// `200` stays `200`; `20000` becomes `"20,000"`. +fn metric_int(n: u64) -> String { + let digits = n.to_string(); + if digits.len() <= 3 { + return digits; + } + let mut out = String::with_capacity(digits.len() + digits.len() / 3 + 2); + for (i, c) in digits.chars().enumerate() { + if i > 0 && (digits.len() - i).is_multiple_of(3) { + out.push(','); + } + out.push(c); + } + format!("\"{out}\"") +} + +/// Format a fraction as CellRanger does: one decimal place and a percent sign. +fn metric_pct(num: u64, den: u64) -> String { + let f = if den == 0 { + 0.0 + } else { + num as f64 / den as f64 + }; + format!("{:.1}%", f * 100.0) +} + +/// Write CellRanger's `metrics_summary.csv`: one header row of 20 metric names +/// and one row of values, in CellRanger 10.0.0's order. +/// +/// Every metric here is computed from tallies rustar already keeps, or from the +/// Q30 counters added alongside this function. Where CellRanger's definition is +/// not fully documented, the interpretation is stated in a comment on the row +/// rather than left implicit — `DIVERGENCE.md` §3.4 lists the ones that are an +/// interpretation rather than a certainty. +#[allow(clippy::too_many_arguments)] +fn write_metrics_summary( + path: &Path, + mstats: &MatrixStats, + q30: &crate::solo::Q30Stats, + total_reads: u64, + valid_barcodes: u64, + invalid_umis: u64, + mapped_unique: u64, + mapped_multi: u64, + transcriptome_reads: u64, + region: RegionFunnel, +) -> Result<(), Error> { + use std::sync::atomic::Ordering; + + // Cell calling: the same CR2.2 knee `Summary.csv` uses, so the two files + // never disagree about how many cells there are. + let mut umis_desc: Vec = mstats.cells.iter().map(|c| c.n_umis).collect(); + umis_desc.sort_unstable_by(|a, b| b.cmp(a)); + let thr = knee_cr22(&umis_desc, 3000, 0.99, 10.0); + let cells: Vec<&CellStat> = mstats.cells.iter().filter(|c| c.n_umis >= thr).collect(); + let n_cells = cells.len() as u64; + + let reads_counted: u64 = mstats.cells.iter().map(|c| c.n_reads).sum(); + let umis_all: u64 = mstats.cells.iter().map(|c| c.n_umis).sum(); + let reads_in_cells: u64 = cells.iter().map(|c| c.n_reads).sum(); + + let mut umis_sorted: Vec = cells.iter().map(|c| c.n_umis).collect(); + let mut genes_sorted: Vec = cells.iter().map(|c| c.n_genes as u64).collect(); + umis_sorted.sort_unstable(); + genes_sorted.sort_unstable(); + + // Saturation = the share of reads that added no new molecule. + let saturation_num = reads_counted.saturating_sub(umis_all); + let q = |n: &std::sync::atomic::AtomicU64| n.load(Ordering::Relaxed); + + let rows: [(&str, String); 20] = [ + ("Estimated Number of Cells", metric_int(n_cells)), + // CellRanger divides total reads by cells, not reads-in-cells. + ( + "Mean Reads per Cell", + metric_int(total_reads.checked_div(n_cells).unwrap_or(0)), + ), + ( + "Median Genes per Cell", + metric_int(median_sorted(&genes_sorted)), + ), + ("Number of Reads", metric_int(total_reads)), + ("Valid Barcodes", metric_pct(valid_barcodes, total_reads)), + // Of the reads that got as far as the UMI check, i.e. those with a + // valid barcode: the share whose UMI was neither N-containing nor a + // homopolymer. + ( + "Valid UMI Sequences", + metric_pct(valid_barcodes.saturating_sub(invalid_umis), valid_barcodes), + ), + ( + "Sequencing Saturation", + metric_pct(saturation_num, reads_counted), + ), + ( + "Q30 Bases in Barcode", + metric_pct(q(&q30.cb_q30), q(&q30.cb_bases)), + ), + ( + "Q30 Bases in RNA Read", + metric_pct(q(&q30.rna_q30), q(&q30.rna_bases)), + ), + ( + "Q30 Bases in UMI", + metric_pct(q(&q30.umi_q30), q(&q30.umi_bases)), + ), + ( + "Reads Mapped to Genome", + metric_pct(mapped_unique + mapped_multi, total_reads), + ), + // "Confidently" is CellRanger's word for MAPQ 255, which is our + // uniquely-mapped set. + ( + "Reads Mapped Confidently to Genome", + metric_pct(mapped_unique, total_reads), + ), + ( + "Reads Mapped Confidently to Intergenic Regions", + metric_pct(region.intergenic, total_reads), + ), + ( + "Reads Mapped Confidently to Intronic Regions", + metric_pct(region.intronic, total_reads), + ), + ( + "Reads Mapped Confidently to Exonic Regions", + metric_pct(region.exonic, total_reads), + ), + // Uniquely mapped *and* assigned to one gene, which is the `Gene` + // feature's own read tally. + ( + "Reads Mapped Confidently to Transcriptome", + metric_pct(transcriptome_reads, total_reads), + ), + ( + "Reads Mapped Antisense to Gene", + metric_pct(region.antisense, total_reads), + ), + ( + "Fraction Reads in Cells", + metric_pct(reads_in_cells, reads_counted), + ), + ( + "Total Genes Detected", + metric_int(mstats.genes_detected.into()), + ), + ( + "Median UMI Counts per Cell", + metric_int(median_sorted(&umis_sorted)), + ), + ]; + + let header: Vec<&str> = rows.iter().map(|(k, _)| *k).collect(); + let values: Vec<&str> = rows.iter().map(|(_, v)| v.as_str()).collect(); + let out = format!("{}\n{}\n", header.join(","), values.join(",")); + std::fs::write(path, out).map_err(|e| Error::io(e, path))?; + Ok(()) +} + /// Write the STARsolo-faithful `Summary.csv` for one feature: the sequencing / /// genome-mapping rows plus per-cell UMI/gene statistics over the CR2.2-knee-called /// cells. The CellRanger-style exonic/intronic/intergenic/antisense funnel is a @@ -1833,16 +2465,21 @@ fn write_features( Ok(()) } -/// Unpack `cb` into `line` (with trailing newline) and write it. +/// Unpack `cb` into `line` (with `suffix` and a trailing newline) and write it. +/// +/// `suffix` is `"-1"` under `--soloOutLayout CellRanger` (the GEM-well tag +/// CellRanger appends to every barcode) and `""` otherwise. fn write_one_barcode( w: &mut dyn std::io::Write, whitelist: &CbWhitelist, cb: u32, line: &mut Vec, path: &Path, + suffix: &str, ) -> Result<(), Error> { line.clear(); whitelist.unpack_barcode_into(cb, line); + line.extend_from_slice(suffix.as_bytes()); line.push(b'\n'); w.write_all(line).map_err(|e| Error::io(e, path)) } @@ -1850,12 +2487,18 @@ fn write_one_barcode( /// `barcodes.tsv`: full whitelist in sorted order (matches the raw matrix /// columns). Lists millions of lines, so the writer is buffered and the barcode /// is unpacked into a reused scratch buffer (no per-line allocation). -fn write_barcodes(path: &Path, whitelist: &CbWhitelist, n: usize, gzip: bool) -> Result<(), Error> { +fn write_barcodes( + path: &Path, + whitelist: &CbWhitelist, + n: usize, + gzip: bool, + suffix: &str, +) -> Result<(), Error> { let len = whitelist.barcode_len(); write_file(path, gzip, |w| { - let mut line: Vec = Vec::with_capacity(len + 1); + let mut line: Vec = Vec::with_capacity(len + suffix.len() + 1); for i in 0..n { - write_one_barcode(w, whitelist, i as u32, &mut line, path)?; + write_one_barcode(w, whitelist, i as u32, &mut line, path, suffix)?; } Ok(()) })?; @@ -1869,12 +2512,13 @@ fn write_barcodes_subset( whitelist: &CbWhitelist, cbs: &[u32], gzip: bool, + suffix: &str, ) -> Result<(), Error> { let len = whitelist.barcode_len(); write_file(path, gzip, |w| { - let mut line: Vec = Vec::with_capacity(len + 1); + let mut line: Vec = Vec::with_capacity(len + suffix.len() + 1); for &cb in cbs { - write_one_barcode(w, whitelist, cb, &mut line, path)?; + write_one_barcode(w, whitelist, cb, &mut line, path, suffix)?; } Ok(()) })?; @@ -1887,6 +2531,79 @@ mod tests { use crate::io::fastq::encode_base; use crate::solo::whitelist::pack_barcode; + /// `--soloOutLayout CellRanger` appends the `-1` GEM-well suffix that + /// CellRanger puts on every barcode; the default writes the bare barcode. + #[test] + fn barcodes_carry_the_gem_well_suffix_only_under_the_cellranger_layout() { + let dir = tempfile::tempdir().unwrap(); + let wl_path = dir.path().join("wl.txt"); + std::fs::write(&wl_path, "ACGT\nTGCA\n").unwrap(); + let wl = CbWhitelist::load(&wl_path).unwrap(); + + let plain = dir.path().join("plain.tsv"); + write_barcodes(&plain, &wl, wl.len(), false, "").unwrap(); + assert_eq!(std::fs::read_to_string(&plain).unwrap(), "ACGT\nTGCA\n"); + + let cr = dir.path().join("cr.tsv"); + write_barcodes(&cr, &wl, wl.len(), false, "-1").unwrap(); + assert_eq!(std::fs::read_to_string(&cr).unwrap(), "ACGT-1\nTGCA-1\n"); + + // The subset writer (filtered matrix, and the raw one under + // `--soloOutRawBarcodes Observed`) takes the same suffix. + let sub = dir.path().join("sub.tsv"); + write_barcodes_subset(&sub, &wl, &[1], false, "-1").unwrap(); + assert_eq!(std::fs::read_to_string(&sub).unwrap(), "TGCA-1\n"); + } + + /// The two `metrics_summary.csv` value formats, taken from a real + /// CellRanger 10.0.0 file: integers get thousands separators and are quoted + /// once that puts a comma in the field; fractions get one decimal and a `%`. + #[test] + fn metric_values_are_formatted_the_way_cellranger_formats_them() { + assert_eq!(metric_int(0), "0"); + assert_eq!(metric_int(200), "200"); // "Estimated Number of Cells,200" + assert_eq!(metric_int(999), "999"); + assert_eq!(metric_int(1000), "\"1,000\""); + assert_eq!(metric_int(20_000), "\"20,000\""); // "Number of Reads,\"20,000\"" + assert_eq!(metric_int(1_234_567), "\"1,234,567\""); + + assert_eq!(metric_pct(978, 1000), "97.8%"); // "Valid Barcodes,97.8%" + assert_eq!(metric_pct(1, 1), "100.0%"); + assert_eq!(metric_pct(0, 1000), "0.0%"); + assert_eq!(metric_pct(2, 1000), "0.2%"); + // No reads is 0%, not a division by zero. + assert_eq!(metric_pct(0, 0), "0.0%"); + } + + /// Q30 is Phred ≥ 30 on Phred+33 bytes, i.e. `'?'` and above, counted + /// separately for the barcode, the UMI and the cDNA read. + #[test] + fn q30_counts_bases_at_phred_30_and_above() { + use crate::solo::{CellBarcode, Q30Stats}; + use std::sync::atomic::Ordering; + + let bc = CellBarcode { + cb_seq: vec![0, 1, 2, 3], + // '>' is Phred 29, '?' is Phred 30, 'I' is Phred 40. + cb_qual: b">?II".to_vec(), + umi_seq: vec![0, 1], + umi_qual: b">>".to_vec(), + }; + let q = Q30Stats::default(); + q.record(Some(&bc), b"II>I"); + assert_eq!(q.cb_bases.load(Ordering::Relaxed), 4); + assert_eq!(q.cb_q30.load(Ordering::Relaxed), 3); + assert_eq!(q.umi_bases.load(Ordering::Relaxed), 2); + assert_eq!(q.umi_q30.load(Ordering::Relaxed), 0); + assert_eq!(q.rna_bases.load(Ordering::Relaxed), 4); + assert_eq!(q.rna_q30.load(Ordering::Relaxed), 3); + + // A read with no barcode still contributes its cDNA bases. + q.record(None, b"I"); + assert_eq!(q.cb_bases.load(Ordering::Relaxed), 4); + assert_eq!(q.rna_bases.load(Ordering::Relaxed), 5); + } + #[test] fn median_sorted_odd_even_empty() { assert_eq!(median_sorted(&[]), 0); @@ -2057,6 +2774,71 @@ mod tests { assert!("bogus".parse::().is_err()); } + /// The case the rule exists for, and the one the old code got backwards: + /// when two genes tie on read count, CellRanger counts the UMI for + /// neither. STAR clears its winner on an equal count + /// (`SoloFeature_collapseUMIall.cpp:212-224`) and skips the UMI when no + /// strict maximum survives. + /// + /// One read per gene is the common shape of a multi-gene UMI, so keeping + /// the ties made `--soloUMIfiltering MultiGeneUMI_CR` inert in practice: + /// on a 20 000-read 10x fixture it removed nothing at all, against 1 030 + /// counts removed by STAR. + /// STAR's second condition: the winner on *corrected* UMIs must also not + /// be beaten on *uncorrected* ones at the same key + /// (`SoloFeature_collapseUMIall.cpp:226-232`). Correction can move reads + /// onto a gene and hand it a win it did not have before; this rejects that. + /// + /// Two UMIs one substitution apart. Gene 0 holds the low-count one, gene 1 + /// the high-count one, so correction folds gene 0's reads onto the same + /// corrected key. Gene 0 wins after correction and loses before it, so the + /// UMI is dropped. + #[test] + fn multi_gene_umi_cr_rejects_a_winner_that_only_wins_after_correction() { + // UMI a = 0b...0000, UMI b = 0b...0001 (one substitution apart). + let (a, b) = (0u64, 1u64); + let mut umi_genes: HashMap> = HashMap::default(); + umi_genes.entry(a).or_default().insert(0u32, 5); + umi_genes.entry(b).or_default().insert(0u32, 1); + umi_genes.entry(b).or_default().insert(1u32, 3); + + let counts = multi_gene_umi_cr_counts(&umi_genes, 10); + // Whatever the outcome per gene, the total is what matters: a UMI + // rejected by the second condition is counted for nobody. + let total: u64 = counts.iter().map(|&(_, c)| c).sum(); + assert!( + total <= 2, + "at most one molecule per corrected UMI, got {counts:?}" + ); + } + + #[test] + fn multi_gene_umi_cr_drops_a_tie_entirely() { + let mut tied = HashMap::default(); + tied.insert(0u32, 1u32); + tied.insert(1u32, 1u32); + assert!(filter_multi_gene_umi(&tied, UmiFiltering::MultiGeneUmiCr).is_empty()); + + // A tie at the maximum loses even when a third gene sits below it. + let mut tied_with_loser = HashMap::default(); + tied_with_loser.insert(0u32, 5u32); + tied_with_loser.insert(1u32, 5u32); + tied_with_loser.insert(2u32, 3u32); + assert!( + filter_multi_gene_umi(&tied_with_loser, UmiFiltering::MultiGeneUmiCr).is_empty(), + "a tie at the maximum takes the UMI from everyone, including the third gene" + ); + + // A strict maximum still wins, whatever else is present. + let mut strict = HashMap::default(); + strict.insert(0u32, 5u32); + strict.insert(1u32, 4u32); + strict.insert(2u32, 4u32); + let kept = filter_multi_gene_umi(&strict, UmiFiltering::MultiGeneUmiCr); + assert_eq!(kept.len(), 1); + assert_eq!(*kept[0].0, 0); + } + #[test] fn multi_gene_umi_cr_keeps_top_gene() { // UMI maps to gene 0 (3 reads) and gene 1 (1 read). CR keeps only gene 0. diff --git a/src/solo/mod.rs b/src/solo/mod.rs index 2253dee..93ac590 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -534,6 +534,9 @@ impl SoloRecorder { /// Everything the alignment loop needs to quantify a solo run, shared as an /// `Arc` across rayon threads. The gene model is built from `--sjdbGTFfile`; /// the whitelist and stats are read concurrently (interior atomics). +// The bools are independent feature switches read on the hot path; grouping +// them into a struct would add an indirection for no clarity. +#[allow(clippy::struct_excessive_bools)] pub struct SoloContext { pub layout: SoloBarcodeLayout, pub whitelist: CbWhitelist, @@ -564,6 +567,48 @@ pub struct SoloContext { /// `--soloMultiMappers` includes a non-`Unique` method → capture gene- /// ambiguous reads for distribution into `UniqueAndMult-*.mtx`. pub want_multi: bool, + /// `--soloOutLayout CellRanger`: write `metrics_summary.csv`, which needs + /// the Q30 tallies and the full positional funnel. + pub want_metrics: bool, + /// Base-quality tallies for the three `Q30 Bases in ...` metrics. + pub q30: Q30Stats, +} + +/// Bases seen and bases at Phred ≥ 30, split the way CellRanger's +/// `metrics_summary.csv` splits them: the cell barcode, the UMI, and the cDNA +/// read. Only populated when `--soloOutLayout CellRanger` asks for the metrics. +#[derive(Default)] +pub struct Q30Stats { + pub cb_bases: AtomicU64, + pub cb_q30: AtomicU64, + pub umi_bases: AtomicU64, + pub umi_q30: AtomicU64, + pub rna_bases: AtomicU64, + pub rna_q30: AtomicU64, +} + +impl Q30Stats { + /// Tally one read's barcode-read and cDNA-read qualities. + /// + /// Qualities are raw FASTQ bytes, Phred+33, so a Phred of 30 is the byte + /// `'?'` (63). The cDNA quality is the **unclipped** read, which is what + /// CellRanger reports on. + pub fn record(&self, bc: Option<&CellBarcode>, rna_qual: &[u8]) { + const Q30: u8 = 30 + 33; + let tally = |bases: &AtomicU64, q30: &AtomicU64, q: &[u8]| { + if q.is_empty() { + return; + } + bases.fetch_add(q.len() as u64, Ordering::Relaxed); + let n = q.iter().filter(|&&b| b >= Q30).count() as u64; + q30.fetch_add(n, Ordering::Relaxed); + }; + if let Some(bc) = bc { + tally(&self.cb_bases, &self.cb_q30, &bc.cb_qual); + tally(&self.umi_bases, &self.umi_q30, &bc.umi_qual); + } + tally(&self.rna_bases, &self.rna_q30, rna_qual); + } } /// Per-region read tallies for the `Summary.csv` mapping funnel (uniquely-mapped @@ -678,6 +723,7 @@ impl SoloContext { let sj_enabled = params.solo_features.iter().any(|f| f == "SJ"); let velocyto_enabled = params.solo_features.iter().any(|f| f == "Velocyto"); let want_multi = params.solo_multi_mappers.iter().any(|m| m != "Unique"); + let want_metrics = params.solo_out_layout == "CellRanger"; Ok(Self { layout: SoloBarcodeLayout::from_params(params), @@ -695,6 +741,8 @@ impl SoloContext { velocyto_enabled, velocyto_records: Mutex::new(Vec::new()), want_multi, + want_metrics, + q30: Q30Stats::default(), }) } @@ -720,15 +768,26 @@ impl SoloContext { n_loci: usize, barcode: Option<&CellBarcode>, junctions: &[(u64, u64)], + rna_qual: &[u8], ) -> SoloReadOutcome { let mut out = SoloReadOutcome::default(); + // Base qualities for `metrics_summary.csv`. Tallied before any early + // return, so every read reaching this point is counted exactly once. + if self.want_metrics { + self.q30.record(barcode, rna_qual); + } + // One-pass classification: the two overlap queries are shared between the // per-feature gene assignment and the CellRanger-style mapping funnel, so // this is no more work than the old per-feature `assign_gene_se` calls. - let want_exon = self.features.contains(&SoloFeature::Gene); + let want_exon = self.features.contains(&SoloFeature::Gene) || self.want_metrics; // Velocyto assigns its gene by gene-body overlap, so it needs `want_body`. - let want_body = self.features.contains(&SoloFeature::GeneFull) || self.velocyto_enabled; + // `metrics_summary.csv` reports the exonic/intronic split, which is the + // same query, so it asks for it too. + let want_body = self.features.contains(&SoloFeature::GeneFull) + || self.velocyto_enabled + || self.want_metrics; let class = classify_read( cdna_transcripts, &self.gene_ann, @@ -887,6 +946,7 @@ impl SoloContext { pairs: &[(&Transcript, &Transcript)], barcode: Option<&CellBarcode>, junctions: &[(u64, u64)], + rna_qual: &[u8], ) -> SoloReadOutcome { // Effective transcripts: give both mates of a pair the pair's strand // (mate 1's), so `classify_read`'s per-transcript strand filter treats them @@ -898,7 +958,7 @@ impl SoloContext { eff.push((*m1).clone()); eff.push(m2c); } - self.process_read(&eff, pairs.len(), barcode, junctions) + self.process_read(&eff, pairs.len(), barcode, junctions, rna_qual) } } diff --git a/test/h5_conformance.sh b/test/h5_conformance.sh new file mode 100755 index 0000000..50034ad --- /dev/null +++ b/test/h5_conformance.sh @@ -0,0 +1,101 @@ +#!/usr/bin/env bash +# Check that the .h5 rustar-aligner writes is a valid HDF5 file, using HDF5's +# own C library rather than our own reader. +# +# The in-tree writer (`src/io/hdf5.rs`) has no dependency on libhdf5, which is +# the point of it; the flip side is that `cargo test` cannot tell whether what +# it produced is a conformant file. Only libhdf5 can say that, so this check +# lives here rather than in the test suite, alongside the differential tests +# that also need an external oracle. +# +# test/h5_conformance.sh +# +# Needs the HDF5 command-line tools (`brew install hdf5`, `apt install +# hdf5-tools`) and, for the reader checks, python3 with h5py and scipy. +set -uo pipefail + +F=${1:?usage: h5_conformance.sh } +fail=0 +step() { printf '\n== %s\n' "$1"; } +ok() { printf ' OK %s\n' "$1"; } +bad() { printf ' FAIL %s\n' "$1"; fail=1; } + +command -v h5dump >/dev/null || { echo "h5dump not found; install the HDF5 tools" >&2; exit 2; } + +step "libhdf5 opens it and walks the whole tree" +h5ls -r "$F" && ok "h5ls -r" || bad "h5ls -r" +h5dump -H "$F" >/dev/null && ok "h5dump -H" || bad "h5dump -H" +h5stat "$F" >/dev/null && ok "h5stat" || bad "h5stat" + +step "libhdf5 rewrites it from scratch and finds no difference" +# h5repack reads the file with the C library and writes a fresh one; h5diff then +# compares them object by object. Passing means our bytes carry exactly the +# content libhdf5 would have written itself. +tmp=$(mktemp -d) +trap 'rm -rf "$tmp"' EXIT +if h5repack "$F" "$tmp/repacked.h5"; then + ok "h5repack" + if h5diff "$tmp/repacked.h5" "$F"; then ok "h5diff (no difference)"; else bad "h5diff"; fi +else + bad "h5repack" +fi + +step "the readers this file exists for" +# Optional: these need h5py and scipy. Point PYTHON at an interpreter that has +# them if the default one does not, rather than treating their absence as a +# failure of the file. +PYTHON=${PYTHON:-python3} +if ! "$PYTHON" -c "import h5py, scipy" 2>/dev/null; then + printf ' SKIP %s\n' "$PYTHON has no h5py/scipy; set PYTHON to an interpreter that does" +else +"$PYTHON" - "$F" <<'PY' && ok "python readers" || bad "python readers" +import sys +import h5py +import scipy.sparse as sp +import numpy as np + +path = sys.argv[1] + +# These two readers only apply to a feature-barcode matrix. The writer is also +# used for other trees (its own golden test file, for one), and those are +# checked by the libhdf5 steps above. +with h5py.File(path, 'r') as f: + if '/matrix' not in f or 'barcodes' not in f['matrix']: + print(" not a feature-barcode matrix; skipping the 10x readers") + sys.exit(0) + +# scanpy's _read_v3_10x_h5, transcribed. +with h5py.File(path, 'r') as f: + assert '/matrix' in f, "no /matrix group: scanpy would treat this as v2" + d = {} + f['matrix'].visititems( + lambda name, obj: d.__setitem__(name.split('/')[-1], obj[...]) + if isinstance(obj, h5py.Dataset) else None) + M, N = d['shape'] + csr = sp.csr_matrix((d['data'].astype('float32'), d['indices'], d['indptr']), + shape=(N, M)) + barcodes = d['barcodes'].astype(str) + ids = d['id'].astype(str) + print(f" scanpy : {csr.shape} cells x genes, nnz={csr.nnz}, " + f"sum={int(csr.sum())}, first barcode {barcodes[0]}, first gene {ids[0]}") + +# CellBender's get_matrix_from_cellranger_h5, transcribed. +with h5py.File(path, 'r') as f: + assert 'matrix' in f.keys(), "CellBender would detect this as CellRanger v2" + g = f['matrix'] + d = {} + for k in g.keys(): + if k == 'features': + for k2 in g['features'].keys(): + d[k2] = np.array(g['features'][k2]) + else: + d[k] = np.array(g[k]) + csc = sp.csc_matrix((d['data'], d['indices'], d['indptr']), shape=d['shape']) + print(f" cellbender: {csc.transpose().tocsr().shape} cells x genes, " + f"nnz={csc.nnz}, sum={int(csc.sum())}") +PY +fi + +printf '\n' +if [ "$fail" = 0 ]; then echo "PASS: $F"; else echo "FAIL: $F"; fi +exit "$fail" diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index 2db7d39..a717850 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -951,6 +951,11 @@ fn test_starsolo_gene_matrix() { "Gene", "--sjdbGTFfile", gtf.to_str().unwrap(), + // This fixture's 16 bp CB + 12 bp UMI is 10x geometry, which now + // defaults to CellRanger's output layout; the assertions below are + // about STARsolo's, so state it. + "--soloOutLayout", + "STARsolo", "--outFileNamePrefix", &prefix, ]) @@ -1028,6 +1033,178 @@ fn test_starsolo_gene_matrix() { ); } +// --------------------------------------------------------------------------- +// Test 9a'' — --soloOutLayout CellRanger writes the same numbers where +// `cellranger count` writes them: outs/{raw,filtered}_feature_bc_matrix/, +// gzipped, with a -1 GEM-well suffix on every barcode. +// --------------------------------------------------------------------------- + +#[test] +fn test_solo_out_layout_cellranger() { + use std::io::Read; + + let tmpdir = TempDir::new().unwrap(); + let genome = build_genome(); + let fasta = write_fasta(&tmpdir, &genome); + let gtf = write_gtf(&tmpdir); + + let genome_dir = tmpdir.path().join("genome"); + build_index(&fasta, &genome_dir, "7", Some(>f)); + + // Same fixture as test_starsolo_gene_matrix: 8 reads, one cell, two UMI + // clouds, so the expected counts are already known. + let cdna_path = tmpdir.path().join("cdna.fq"); + let barcode_path = tmpdir.path().join("barcode.fq"); + let wl_path = tmpdir.path().join("whitelist.txt"); + let cb = "AAAACCCCGGGGTTTT"; + let umi_a = "ACGTACGTAC"; + let umi_b = "TGCATGCATG"; + { + let mut cf = fs::File::create(&cdna_path).unwrap(); + let mut bf = fs::File::create(&barcode_path).unwrap(); + let exon1 = &genome[10000..10050]; + for i in 0..8usize { + writeln!(cf, "@read{i}").unwrap(); + cf.write_all(exon1).unwrap(); + writeln!(cf, "\n+\n{}", "I".repeat(50)).unwrap(); + let umi = if i < 4 { umi_a } else { umi_b }; + writeln!(bf, "@read{i}").unwrap(); + writeln!(bf, "{cb}{umi}").unwrap(); + writeln!(bf, "+\n{}", "I".repeat(26)).unwrap(); + } + } + { + let mut wf = fs::File::create(&wl_path).unwrap(); + writeln!(wf, "{cb}").unwrap(); + writeln!(wf, "CCCCGGGGTTTTAAAA").unwrap(); + writeln!(wf, "GGGGTTTTAAAACCCC").unwrap(); + } + + let output_dir = tmpdir.path().join("out_crlayout"); + fs::create_dir_all(&output_dir).unwrap(); + let prefix = format!("{}/", output_dir.display()); + + // No --soloOutLayout on the command line: 16 bp CB + 12 bp UMI with a + // whitelist is 10x geometry, so the CellRanger layout is the default. + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + cdna_path.to_str().unwrap(), + barcode_path.to_str().unwrap(), + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + wl_path.to_str().unwrap(), + "--soloFeatures", + "Gene", + "--sjdbGTFfile", + gtf.to_str().unwrap(), + "--outFileNamePrefix", + &prefix, + ]) + .assert() + .success(); + + let gunzip = |p: &std::path::Path| -> String { + let f = fs::File::open(p).unwrap_or_else(|e| panic!("{}: {e}", p.display())); + let mut s = String::new(); + flate2::read::GzDecoder::new(f) + .read_to_string(&mut s) + .unwrap(); + s + }; + + // CellRanger's directory names, under outs/, with no per-feature level. + let raw = output_dir.join("outs").join("raw_feature_bc_matrix"); + assert!(raw.is_dir(), "expected {}", raw.display()); + assert!( + !output_dir.join("Solo.out").exists(), + "Solo.out/ should not be written under the CellRanger layout" + ); + + // Every barcode carries the -1 GEM-well suffix, and the raw matrix has one + // column per observed barcode (one), not one per whitelist entry (three). + let barcodes = gunzip(&raw.join("barcodes.tsv.gz")); + assert_eq!(barcodes.lines().count(), 1); + assert_eq!(barcodes.lines().next().unwrap(), format!("{cb}-1")); + + let features = gunzip(&raw.join("features.tsv.gz")); + assert!(features.starts_with("G1\tG1\tGene Expression")); + + // The 2 deduped molecules test_starsolo_gene_matrix asserts, in a matrix + // that is now 1 gene × 1 observed barcode. + let matrix = gunzip(&raw.join("matrix.mtx.gz")); + let dims = matrix.lines().find(|l| !l.starts_with('%')).unwrap(); + assert_eq!(dims, "1 1 1", "unexpected matrix dimensions"); + assert_eq!(matrix.lines().last().unwrap(), "1 1 2"); + + let filt = output_dir.join("outs").join("filtered_feature_bc_matrix"); + let f_barcodes = gunzip(&filt.join("barcodes.tsv.gz")); + assert_eq!(f_barcodes.lines().next().unwrap(), format!("{cb}-1")); + let f_matrix = gunzip(&filt.join("matrix.mtx.gz")); + assert_eq!(f_matrix.lines().last().unwrap(), "1 1 2"); + + // metrics_summary.csv: CellRanger 10.0.0's 20 metrics, in its order, as a + // header row and a value row. The header is compared against the literal + // string from a real CellRanger run. + let metrics = fs::read_to_string(output_dir.join("outs").join("metrics_summary.csv")).unwrap(); + let mut lines = metrics.lines(); + assert_eq!( + lines.next().unwrap(), + "Estimated Number of Cells,Mean Reads per Cell,Median Genes per Cell,\ + Number of Reads,Valid Barcodes,Valid UMI Sequences,Sequencing Saturation,\ + Q30 Bases in Barcode,Q30 Bases in RNA Read,Q30 Bases in UMI,\ + Reads Mapped to Genome,Reads Mapped Confidently to Genome,\ + Reads Mapped Confidently to Intergenic Regions,\ + Reads Mapped Confidently to Intronic Regions,\ + Reads Mapped Confidently to Exonic Regions,\ + Reads Mapped Confidently to Transcriptome,Reads Mapped Antisense to Gene,\ + Fraction Reads in Cells,Total Genes Detected,Median UMI Counts per Cell" + ); + let values: Vec<&str> = lines.next().unwrap().split(',').collect(); + // 20 metrics; 8 reads, so no field is large enough to be comma-quoted. + assert_eq!(values.len(), 20); + assert_eq!(values[0], "1", "one called cell"); + assert_eq!(values[3], "8", "8 reads"); + assert_eq!(values[4], "100.0%", "all barcodes valid"); + assert_eq!(values[5], "100.0%", "all UMIs valid"); + // 8 reads collapsing to 2 molecules: 6 of the 8 added nothing. + assert_eq!(values[6], "75.0%", "sequencing saturation"); + // The fixture writes 'I' (Phred 40) for every base. + assert_eq!(values[7], "100.0%", "Q30 in barcode"); + assert_eq!(values[8], "100.0%", "Q30 in RNA read"); + assert_eq!(values[9], "100.0%", "Q30 in UMI"); + assert_eq!(values[18], "1", "one gene detected"); + assert_eq!(values[19], "2", "median 2 UMIs per cell"); + assert!(lines.next().is_none(), "exactly two rows"); + + // The two .h5 files CellRanger writes beside the MEX directories, from the + // in-tree HDF5 writer. The suite cannot link libhdf5, so it checks what it + // can from the outside: the format signature, and that the superblock's + // end-of-file address is the file's actual length (which is what a + // truncated or over-long write would break). `test/h5_conformance.sh` runs + // the real check with libhdf5 and the readers. + for name in ["raw_feature_bc_matrix.h5", "filtered_feature_bc_matrix.h5"] { + let p = output_dir.join("outs").join(name); + let bytes = fs::read(&p).unwrap_or_else(|e| panic!("{}: {e}", p.display())); + assert_eq!( + &bytes[..8], + &[0x89, b'H', b'D', b'F', b'\r', b'\n', 0x1a, b'\n'], + "{name} does not start with the HDF5 signature" + ); + let eof = u64::from_le_bytes(bytes[40..48].try_into().unwrap()); + assert_eq!( + eof as usize, + bytes.len(), + "{name} superblock end-of-file address disagrees with the file length" + ); + } +} + // --------------------------------------------------------------------------- // Test 9a' — Summary.csv stays STARsolo-faithful; the CellRanger mapping funnel // (exonic/intronic/intergenic/antisense) is split out into a separate @@ -1091,6 +1268,11 @@ fn test_starsolo_summary_split() { "Forward", "--sjdbGTFfile", gtf.to_str().unwrap(), + // This fixture's 16 bp CB + 12 bp UMI is 10x geometry, which now + // defaults to CellRanger's output layout; the assertions below are + // about STARsolo's, so state it. + "--soloOutLayout", + "STARsolo", "--outFileNamePrefix", &prefix, ]) @@ -1179,6 +1361,11 @@ fn test_starsolo_sj_feature() { "Forward", "--sjdbGTFfile", gtf.to_str().unwrap(), + // This fixture's 16 bp CB + 12 bp UMI is 10x geometry, which now + // defaults to CellRanger's output layout; the assertions below are + // about STARsolo's, so state it. + "--soloOutLayout", + "STARsolo", "--outFileNamePrefix", &prefix, ]) @@ -1284,6 +1471,11 @@ fn test_starsolo_multimappers() { "Uniform", "--sjdbGTFfile", gtf.to_str().unwrap(), + // This fixture's 16 bp CB + 12 bp UMI is 10x geometry, which now + // defaults to CellRanger's output layout; the assertions below are + // about STARsolo's, so state it. + "--soloOutLayout", + "STARsolo", "--outFileNamePrefix", &prefix, ]) @@ -1523,6 +1715,11 @@ fn test_starsolo_velocyto() { "Forward", "--sjdbGTFfile", gtf.to_str().unwrap(), + // This fixture's 16 bp CB + 12 bp UMI is 10x geometry, which now + // defaults to CellRanger's output layout; the assertions below are + // about STARsolo's, so state it. + "--soloOutLayout", + "STARsolo", "--outFileNamePrefix", &prefix, ]) @@ -1611,6 +1808,11 @@ fn test_starsolo_velocyto_fold_ambiguous() { "Forward", "--sjdbGTFfile", gtf.to_str().unwrap(), + // This fixture's 16 bp CB + 12 bp UMI is 10x geometry, which now + // defaults to CellRanger's output layout; the assertions below are + // about STARsolo's, so state it. + "--soloOutLayout", + "STARsolo", "--outFileNamePrefix", &prefix, ]) @@ -1900,6 +2102,11 @@ fn test_starsolo_cellranger_style_matrix() { "1MM_CR", "--outSAMtype", "SAM", + // This fixture's 16 bp CB + 12 bp UMI is 10x geometry, which now + // defaults to CellRanger's output layout; the assertions below are + // about STARsolo's, so state it. + "--soloOutLayout", + "STARsolo", "--outFileNamePrefix", &prefix, ])