From 51304b55078218f6efe258e6a1f649bb22c4b471 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 10:06:45 +0200 Subject: [PATCH 01/10] fix(solo): MultiGeneUMI_CR gives a tied UMI to nobody, not to everybody MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `--soloUMIfiltering MultiGeneUMI_CR` kept every gene tied at the highest read count. CellRanger's rule is the opposite on exactly that case: the gene with the *strictly* highest count takes the UMI, and a tie means no gene counts it. STAR walks the genes keeping a running maximum and clears its winner whenever it meets an equal count (`SoloFeature_collapseUMIall.cpp:212-224`): 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 One read per gene is the ordinary shape of a multi-gene UMI, and it is always a tie, so the old rule made the flag inert in practice rather than merely inaccurate. Measured on a 20 000-read 10x fixture (200 cells from the real v3 whitelist, 400 genes, 720 UMIs deliberately shared between two genes), against STAR 2.7.11b with the same flags: identical entries STAR counts rustar counts before 13 749 / 14 806 15 423 16 465 after 13 902 / 13 967 15 423 15 414 The flag removed nothing at all before; STAR removes 1 030 counts. The gap goes from +1 042 to -9. The outcome does not depend on the order the genes are visited — a strict maximum always ends as the winner, a tie always ends with none — so iterating a `HashMap` here stays deterministic. `multi_gene_umi_cr_drops_a_tie_entirely` pins the case the old tests missed: they only covered 3 reads against 1, where both rules agree. Not yet implemented, and stated so rather than left to be discovered: STAR applies a second condition, that the winning gene must also hold the top count among *uncorrected* UMIs (`umiGeneMapCount0`, same file, lines 226-232). That needs the pre-correction counts, which this code does not keep. The 65 entries still differing out of 13 967 are the place to look for its effect. Co-Authored-By: Claude Opus 5 (1M context) --- src/solo/count.rs | 76 +++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 74 insertions(+), 2 deletions(-) diff --git a/src/solo/count.rs b/src/solo/count.rs index 9af3edf..7d906c8 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -822,8 +822,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!(), } } @@ -2057,6 +2092,43 @@ 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. + #[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. From df80683ad2629da4f985d1c74d15ed232fc45a8e Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 10:06:45 +0200 Subject: [PATCH 02/10] docs(changelog): record the MultiGeneUMI_CR tie fix Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index b6e4893..91d4370 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -99,6 +99,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 From a8f774e55c72c4dd8d30e155324575227744eb4a Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 10:42:46 +0200 Subject: [PATCH 03/10] feat(solo): --soloOutRawBarcodes Observed, for a CellRanger-shaped raw matrix MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit STARsolo's raw matrix has a column per whitelist barcode. For 10x v3 that is 3 686 400 columns and a 62 MB `barcodes.tsv`, nearly all zeros. CellRanger's `raw_feature_bc_matrix` has a column per *observed* barcode. The two files therefore share no keys, which is not a rounding difference in a comparison, it is zero overlap: comparing our raw output against a real `cellranger count` run gave 0 identical entries out of 27 396 until the columns were reconciled. `--soloOutRawBarcodes Observed` narrows the raw matrix to the barcodes that carry a count. Default `Whitelist` keeps what STARsolo writes, so nothing changes for anyone not asking. Measured on the 20 000-read fixture: Whitelist 3 686 400 barcodes barcodes.tsv 62 668 800 bytes Observed 200 barcodes barcodes.tsv 3 400 bytes with identical counts on both sides: 13 937 entries, 15 414 counts. `finalize_matrix` already took a column remap for the filtered matrix, so this reuses it rather than adding a second path. The observed set is read back from the streamed body, which costs one pass and only when the flag is on. This is a **non-STAR flag** and needs sign-off; recorded in `DIVERGENCE.md` §3.2 rather than presented as parity. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 6 +++++ DIVERGENCE.md | 27 +++++++++++++++++++++ src/params/mod.rs | 15 ++++++++++++ src/solo/count.rs | 62 +++++++++++++++++++++++++++++++++++++++++------ 4 files changed, 102 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 91d4370..073d614 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,12 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- `--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). diff --git a/DIVERGENCE.md b/DIVERGENCE.md index bd957ed..5f8e6e6 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -65,6 +65,33 @@ On the 10k yeast PE benchmark, 4 reads differ in alignment score (AS) because ST --- +### 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. + +--- + ## 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/params/mod.rs b/src/params/mod.rs index 0536a85..9fc858f 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1133,6 +1133,21 @@ 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, + /// 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 diff --git a/src/solo/count.rs b/src/solo/count.rs index 7d906c8..d766641 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -476,6 +476,27 @@ 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()) +} + /// 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 @@ -1235,20 +1256,45 @@ 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)?; + } + None => write_barcodes( + &raw_dir.join(&barcodes_name), + &ctx.whitelist, + sorted.len(), + gzip, + )?, + } finalize_matrix( &body, &raw_dir.join(&matrix_name), gzip, n_genes, - sorted.len(), + raw_cols, mstats.nnz, - None, + raw_remap.as_ref(), )?; log::info!( "STARsolo: wrote {}/raw matrix ({} genes × {} barcodes, {} entries){}", From abfcf3e2e39d346efde2868a6b87a5be110d9b50 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 10:53:46 +0200 Subject: [PATCH 04/10] fix(solo): MultiGeneUMI_CR decides ownership on corrected UMIs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit STAR corrects UMIs within each gene *before* deciding which gene owns a UMI, and applies two conditions, not one (`SoloFeature_collapseUMIall.cpp:134-148` and `:203-235`): 1. one gene must hold a strictly higher read count than every other, on the **corrected** UMI map — that is #173, already landed; 2. and that winner must not be beaten in the **uncorrected** map at the same key. The second condition exists because correction moves reads between UMIs: a gene can win only because correction folded a neighbouring UMI onto it, and STAR rejects that win rather than counting it. Reproducing it needs the order STAR uses. The generic path here filters multi-gene UMIs first and corrects afterwards, which cannot express either condition: by the time correction happens the ownership decision is already made. `MultiGeneUMI_CR` therefore takes its own path, which is also what STAR does — the flag is only valid with `--soloUMIdedup 1MM_CR`, so there is no combination this bypasses. `cellranger_1mm_map` exposes the correction mapping that `cellranger_1mm` already computed and threw away. Measured against **CellRanger 10.0.0** on the 20 000-read fixture from #172, with #165 and #173 also applied: identical entries CellRanger rustar #165 + #173 13 651 / 13 709 15 111 15 091 plus this change 13 676 / 13 709 15 111 15 116 Entries CellRanger has and we do not go from 29 to 7, and the count gap from -20 to +5, which is 0.03%. Co-Authored-By: Claude Opus 5 (1M context) --- src/solo/count.rs | 160 +++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 145 insertions(+), 15 deletions(-) diff --git a/src/solo/count.rs b/src/solo/count.rs index d766641..903514f 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; @@ -829,6 +849,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)> { @@ -2148,6 +2250,34 @@ mod tests { /// 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(); From 9f823e2bfd6a5152ae9949315c29212be9f47580 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 11:06:14 +0200 Subject: [PATCH 05/10] feat(solo): CellRanger behaviour by default on 10x geometry MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Aligning 10x data and comparing the count matrix against CellRanger gave a successful run and different numbers, with nothing in the output pointing at the five flags that explain the difference. Measured against CellRanger 10.0.0 on a 20 000-read fixture, those flags are the whole gap: 8.9% away without them, 0.03% with them. When the geometry is unambiguously 10x — `CB_UMI_Simple`, a whitelist, a 16-base cell barcode, a 10- or 12-base UMI — the five now default to their CellRanger values: --clipAdapterType CellRanger4 --outFilterScoreMin 30 --soloCBmatchWLtype 1MM_multi_Nbase_pseudocounts --soloUMIfiltering MultiGeneUMI_CR --soloUMIdedup 1MM_CR A flag given on the command line always wins, including when the value asked for is STARsolo's own default: `value_source` distinguishes an explicit flag from a default, so the divergence is escapable by naming what you want. Every substitution is logged at INFO with the geometry that triggered it. **This changes default output behaviour on 10x runs and diverges from STARsolo**, which is why it is confined to a geometry nothing else in common use shares, why it is announced on every run it touches, and why it is in `DIVERGENCE.md` §1.3 as the largest entry in that file. It needs sign-off. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 6 ++ DIVERGENCE.md | 28 +++++++ src/params/mod.rs | 208 ++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 242 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 073d614..45de81a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,12 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- 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 diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 5f8e6e6..8a0d92c 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -65,6 +65,34 @@ 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.9% away to 0.03%. + +**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 diff --git a/src/params/mod.rs b/src/params/mod.rs index 9fc858f..4ace622 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1363,6 +1363,8 @@ 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); + 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() @@ -1885,8 +1887,214 @@ 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); 5] = [ + ("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"), +]; + +/// 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()], + _ => 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(", ") + ); + } +} + #[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()]); + } + + /// 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). From 014d3a7929133e5c6d7cba8ae9c266defb5008c8 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 11:48:48 +0200 Subject: [PATCH 06/10] docs(divergence): correct the CellRanger gap figure in 1.3 Re-derived both sides from one clean state: the five flags move the matrix from 8.96% above CellRanger to 2.17% above it, not to 0.03%. The earlier figure compared a rustar run against a CellRanger run built from a different state of the fixture. STAR 2.7.11b with the same flags is at +0.09%, so the remaining gap is open rather than closed. Co-Authored-By: Claude Opus 5 (1M context) --- DIVERGENCE.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 8a0d92c..592aa1d 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -80,7 +80,9 @@ 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.9% away to 0.03%. +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 From 816d4e8b7f34547ea70a62f0b21441e20a0a1b0f Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 11:27:58 +0200 Subject: [PATCH 07/10] feat(solo): --soloOutLayout CellRanger, CellRanger's output shape Writes the same numbers where `cellranger count` writes them: outs/{raw,filtered}_feature_bc_matrix/, 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. The layout implies --soloOutGzip yes, --soloOutRawBarcodes Observed and an outs/ output directory; each is still overridable on the command line. On 10x geometry it is the default, alongside the five flags from the previous commit, so it changes where output files are written on those runs. Counts are untouched: 13 959 entries and 15 439 counts in both layouts on the 20 000-read fixture, compared entry by entry. Against a real cellranger count 10.0.0 run on the same fixture the raw barcode sets match exactly, 200 of 200. Recorded in DIVERGENCE.md 3.3. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 9 ++ Cargo.toml | 3 + DIVERGENCE.md | 35 ++++++++ src/params/mod.rs | 163 +++++++++++++++++++++++++++++++++++- src/solo/count.rs | 103 +++++++++++++++++++---- tests/alignment_features.rs | 151 +++++++++++++++++++++++++++++++++ 6 files changed, 447 insertions(+), 17 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 45de81a..a8ba71b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,15 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- `--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 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 592aa1d..76fb734 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -122,6 +122,41 @@ untouched. --- +### 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. + +--- + ## 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/params/mod.rs b/src/params/mod.rs index 4ace622..edfa3c9 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1148,6 +1148,29 @@ pub struct Parameters { 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 @@ -1364,6 +1387,7 @@ impl Parameters { 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(); @@ -1889,12 +1913,15 @@ impl Parameters { /// 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); 5] = [ +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? @@ -1946,6 +1973,7 @@ fn apply_cellranger_defaults_on_10x(params: &mut Parameters, matches: &clap::Arg "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); @@ -1963,6 +1991,34 @@ fn apply_cellranger_defaults_on_10x(params: &mut Parameters, matches: &clap::Arg } } +/// `--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 { @@ -1997,6 +2053,111 @@ mod tests { 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 diff --git a/src/solo/count.rs b/src/solo/count.rs index 903514f..613af35 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -1334,10 +1334,28 @@ 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 { "" }; + // 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 @@ -1380,13 +1398,20 @@ pub fn write_gene_matrix( }; match &observed { Some(cbs) => { - write_barcodes_subset(&raw_dir.join(&barcodes_name), &ctx.whitelist, cbs, gzip)?; + 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( @@ -1399,10 +1424,10 @@ pub fn write_gene_matrix( raw_remap.as_ref(), )?; 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 { "" }, ); @@ -1426,7 +1451,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() @@ -1439,7 +1464,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), @@ -1450,8 +1481,8 @@ 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, ); @@ -1522,11 +1553,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( @@ -1562,6 +1596,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 @@ -2016,16 +2051,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)) } @@ -2033,12 +2073,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(()) })?; @@ -2052,12 +2098,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(()) })?; @@ -2070,6 +2117,30 @@ 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"); + } + #[test] fn median_sorted_odd_even_empty() { assert_eq!(median_sorted(&[]), 0); diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index 2db7d39..a0f812d 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,122 @@ 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"); +} + // --------------------------------------------------------------------------- // Test 9a' — Summary.csv stays STARsolo-faithful; the CellRanger mapping funnel // (exonic/intronic/intergenic/antisense) is split out into a separate @@ -1091,6 +1212,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 +1305,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 +1415,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 +1659,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 +1752,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 +2046,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, ]) From 41848babd59eae6d2ae1b4f59843aeabc670223f Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 11:54:53 +0200 Subject: [PATCH 08/10] feat(solo): metrics_summary.csv under the CellRanger layout Writes CellRanger 10.0.0's 20 metrics, in its order and value formats (thousands separators, quoted when that adds a comma; percentages to one decimal). STARsolo's Summary.csv is unchanged and still written alongside. Twelve of the 20 match a real cellranger count 10.0.0 run exactly on the 20 000-read fixture. The other eight rest on an interpretation of a denominator CellRanger does not document, and each interpretation is stated in DIVERGENCE.md 3.4 rather than left implicit. Three metrics needed new counters: Q30Stats tallies Phred >= 30 bases over the barcode, the UMI and the cDNA read, populated only when the layout asks for the metrics. The exonic/intronic split needs the gene-body overlap query, so a Gene-only run under this layout now performs it. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 6 + DIVERGENCE.md | 50 ++++++++ src/lib.rs | 26 +++- src/solo/count.rs | 243 +++++++++++++++++++++++++++++++++++- src/solo/mod.rs | 66 +++++++++- tests/alignment_features.rs | 34 +++++ 6 files changed, 415 insertions(+), 10 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a8ba71b..c75867d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,12 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- 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 diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 76fb734..6013a61 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -157,6 +157,56 @@ 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. + +--- + ## 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/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/solo/count.rs b/src/solo/count.rs index 613af35..784197c 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -1321,8 +1321,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), @@ -1520,10 +1524,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, @@ -1858,6 +1884,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 @@ -2141,6 +2329,55 @@ mod tests { 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); 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/tests/alignment_features.rs b/tests/alignment_features.rs index a0f812d..7355401 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -1147,6 +1147,40 @@ fn test_solo_out_layout_cellranger() { 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"); } // --------------------------------------------------------------------------- From 47a8d223c70414ca861d22a34daf1ee9e99bc1b9 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 13:51:47 +0200 Subject: [PATCH 09/10] feat(solo): CellRanger .h5 matrices, written without libhdf5 Adds raw_feature_bc_matrix.h5 and filtered_feature_bc_matrix.h5 under --soloOutLayout CellRanger, holding the same counts as the MEX directories in CellRanger's HDF5 layout, so scanpy.read_10x_h5, Seurat's Read10X_h5 and CellBender's .h5 path read them. No new dependency. src/io/hdf5.rs writes the file directly: version 0 superblock, version 1 object headers, symbol-table groups with a version 1 B-tree and a local heap, contiguous uncompressed datasets. That is the subset CellRanger's matrix file uses, and it avoids both a system libhdf5 (which would break cargo install) and cmake (which a vendored build needs). Verified with HDF5's own C library rather than our own reader: test/h5_conformance.sh runs h5ls -r, h5dump -H and h5stat, then an h5repack round-trip through libhdf5 followed by h5diff, which reports no difference. scanpy's and CellBender's loaders both return the expected matrix, and the .h5 matches the .mtx entry by entry: 13 959 entries, 15 439 counts. The unit tests lock determinism and the byte layout with a checksum, since cargo test cannot validate an HDF5 file; the conformance script must be re-run off CI whenever the writer changes. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 8 + DIVERGENCE.md | 52 +++ src/io/hdf5.rs | 785 ++++++++++++++++++++++++++++++++++++ src/io/mod.rs | 1 + src/solo/count.rs | 226 +++++++++++ test/h5_conformance.sh | 93 +++++ tests/alignment_features.rs | 22 + 7 files changed, 1187 insertions(+) create mode 100644 src/io/hdf5.rs create mode 100755 test/h5_conformance.sh diff --git a/CHANGELOG.md b/CHANGELOG.md index c75867d..4b8a352 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,14 @@ 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. 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 diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 6013a61..36fc39c 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -207,6 +207,58 @@ 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 contiguous uncompressed datasets. It +does not chunk, compress, or write variable-length types. + +**Where the file differs from CellRanger's.** CellRanger gzip-compresses its +datasets and pads every string field to 256 bytes; ours are uncompressed and +padded to the longest value present. Readers see identical values either way — +the difference is on-disk size and the `|S` width numpy reports. 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..8c2a0d4 --- /dev/null +++ b/src/io/hdf5.rs @@ -0,0 +1,785 @@ +//! 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 contiguous, uncompressed +//! datasets. 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 chunking, no compression, 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. 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_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 +} + +/// 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. +fn fill_value_message() -> Vec { + vec![2, 2, 0, 0] // version, space allocation time (early), write time, undefined +} + +/// 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 +} + +/// 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 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) + }; + let messages = vec![ + header_message(MSG_DATASPACE, &dataspace_message(&[ds.data.len() as u64])), + header_message(MSG_DATATYPE, &ds.data.datatype()), + header_message(MSG_FILL_VALUE, &fill_value_message()), + 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 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 1.14.6 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. + #[test] + fn the_image_is_deterministic_and_matches_the_validated_layout() { + let build_it = || { + let root = GroupSpec::new("/") + .group(GroupSpec::new("matrix").dataset("data", Data::I32(&[1, 2]))) + .attr("filetype", AttrValue::Str("matrix".into())); + build(&root).unwrap() + }; + let a = build_it(); + assert_eq!(a, build_it(), "two builds of the same tree differ"); + assert_eq!(a.len(), 2192, "image length"); + assert_eq!(fnv1a(&a), 0x67f1_1bdc_8bcd_42e4, "image checksum"); + } + + /// 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 ^= *b as u64; + 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/solo/count.rs b/src/solo/count.rs index 784197c..d3168fe 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -517,6 +517,176 @@ fn observed_barcodes(body: &tempfile::NamedTempFile) -> Result, Error> 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 @@ -1347,6 +1517,15 @@ pub fn write_gene_matrix( ("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. @@ -1427,6 +1606,39 @@ pub fn write_gene_matrix( mstats.nnz, 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 {} matrix ({} genes × {} barcodes, {} entries){}", raw_dir.display(), @@ -1490,6 +1702,20 @@ pub fn write_gene_matrix( 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. diff --git a/test/h5_conformance.sh b/test/h5_conformance.sh new file mode 100755 index 0000000..275b220 --- /dev/null +++ b/test/h5_conformance.sh @@ -0,0 +1,93 @@ +#!/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 "scanpy read_10x_h5 + CellBender loader" || bad "python readers" +import sys +import h5py +import scipy.sparse as sp +import numpy as np + +path = sys.argv[1] + +# 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 7355401..a717850 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -1181,6 +1181,28 @@ fn test_solo_out_layout_cellranger() { 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" + ); + } } // --------------------------------------------------------------------------- From f74336abab17325895271e71214a1e8e9f8c2595 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 14:03:39 +0200 Subject: [PATCH 10/10] feat(hdf5): deflate the dataset chunks Adds chunked storage with the deflate filter to the in-tree HDF5 writer: a version 3 chunked layout message, a version 1 filter pipeline message, and a version 1 chunk B-tree, with one chunk covering each dataset. Compression goes through the libdeflater already in the tree, in the zlib framing HDF5's deflate filter expects. The .h5 goes from 451 KB to 88 KB on the 20 000-read fixture, against CellRanger's 149 KB for the same matrix. Datasets at or below 1 kB stay contiguous: a chunk B-tree node is about 2 kB allocated at full capacity, so compressing them would cost more than it saves. CellRanger draws the same line, leaving _all_tag_keys uncompressed. Two fixes fell out of validating this against libhdf5, both silent before: the fill value message declared no fill value where HDF5 writes a zero-length one, and the filter pipeline carried no filter name. Every reader accepted the result and h5repack rewrote it into a file h5dump could not read. The conformance script now catches that, and the golden test says why it must be re-run. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 7 +- DIVERGENCE.md | 22 ++-- src/io/hdf5.rs | 277 +++++++++++++++++++++++++++++++++++------ test/h5_conformance.sh | 10 +- 4 files changed, 262 insertions(+), 54 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4b8a352..6b0dc86 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,9 +17,10 @@ Sections commonly used: Features, Bug fixes, Other changes. `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. 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`. + 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 diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 36fc39c..749c4c5 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -229,15 +229,19 @@ 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 contiguous uncompressed datasets. It -does not chunk, compress, or write variable-length types. - -**Where the file differs from CellRanger's.** CellRanger gzip-compresses its -datasets and pads every string field to 256 bytes; ours are uncompressed and -padded to the longest value present. Readers see identical values either way — -the difference is on-disk size and the `|S` width numpy reports. The -`genome` field carries the `--genomeDir` directory name, since we have no -equivalent of the name given to `cellranger mkref`, and `software_version` +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 diff --git a/src/io/hdf5.rs b/src/io/hdf5.rs index 8c2a0d4..62771da 100644 --- a/src/io/hdf5.rs +++ b/src/io/hdf5.rs @@ -10,18 +10,19 @@ //! # What it writes //! //! Version 0 superblock, version 1 object headers, old-style groups (symbol -//! table + version 1 B-tree + local heap), and contiguous, uncompressed -//! datasets. 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. +//! 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 chunking, no compression, 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. 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. +//! 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 @@ -64,6 +65,7 @@ const SUPERBLOCK_LEN: usize = 96; 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; @@ -307,10 +309,28 @@ fn dataspace_message(dims: &[u64]) -> Vec { 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. -fn fill_value_message() -> Vec { - vec![2, 2, 0, 0] // version, space allocation time (early), write time, undefined +/// +/// 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` @@ -324,6 +344,120 @@ fn data_layout_message(addr: u64, size: u64) -> Vec { 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(); @@ -438,23 +572,59 @@ fn write_group(img: &mut Image, group: &GroupSpec<'_>) -> Result 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])); } @@ -737,28 +907,53 @@ mod tests { 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 1.14.6 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. + /// 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 build_it = || { - let root = GroupSpec::new("/") - .group(GroupSpec::new("matrix").dataset("data", Data::I32(&[1, 2]))) - .attr("filetype", AttrValue::Str("matrix".into())); - build(&root).unwrap() - }; - let a = build_it(); - assert_eq!(a, build_it(), "two builds of the same tree differ"); - assert_eq!(a.len(), 2192, "image length"); - assert_eq!(fnv1a(&a), 0x67f1_1bdc_8bcd_42e4, "image checksum"); + 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 @@ -766,7 +961,7 @@ mod tests { fn fnv1a(bytes: &[u8]) -> u64 { let mut h: u64 = 0xcbf2_9ce4_8422_2325; for b in bytes { - h ^= *b as u64; + h ^= u64::from(*b); h = h.wrapping_mul(0x1000_0000_01b3); } h diff --git a/test/h5_conformance.sh b/test/h5_conformance.sh index 275b220..50034ad 100755 --- a/test/h5_conformance.sh +++ b/test/h5_conformance.sh @@ -48,7 +48,7 @@ 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 "scanpy read_10x_h5 + CellBender loader" || bad "python readers" +"$PYTHON" - "$F" <<'PY' && ok "python readers" || bad "python readers" import sys import h5py import scipy.sparse as sp @@ -56,6 +56,14 @@ 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"