diff --git a/CHANGELOG.md b/CHANGELOG.md index b6e4893..7b46fd3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,47 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- On 10x geometry, `--soloFeatures` now defaults to `GeneFull` rather than + `Gene`: CellRanger has counted intronic reads toward the gene since v7.0. + Measured on 10x's `pbmc_1k_v3`, exonic-only counting is 30.5% below + CellRanger and `GeneFull` brings it to 1.6%. **This changes default counts + on 10x runs**; naming `--soloFeatures` explicitly still wins. See + `DIVERGENCE.md` §1.3. + +- `--soloCellFilter OrdMag` implements CellRanger's cell call: the same + quantile-over-ratio rule as `CellRanger2.2`, but searching for the expected + cell count by minimising `(OrdMag(x) - x)^2 / x` instead of fixing it at + 3 000. `EmptyDrops_CR` now uses it for its initial cell set, which is the + order CellRanger runs the two steps in; `CellRanger2.2` is unchanged and + remains the default. See `DIVERGENCE.md` §3.5. + +- Under `--soloOutLayout CellRanger`, `metrics_summary.csv` is written with + CellRanger 10.0.0's 20 metrics, in its order and value formats. STARsolo's + `Summary.csv` is unchanged and still written alongside. Twelve of the 20 + match a real `cellranger count` run exactly on the test fixture; the + interpretation behind the other eight is in `DIVERGENCE.md` §3.4. + +- `--soloOutLayout CellRanger` writes the solo matrices in the shape + `cellranger count` produces: `outs/raw_feature_bc_matrix/` and + `outs/filtered_feature_bc_matrix/`, gzipped, with a `-1` GEM-well suffix + on every barcode and one raw column per observed barcode. It implies + `--soloOutGzip yes`, `--soloOutRawBarcodes Observed` and an `outs/` + output directory, each still overridable on the command line. Counts are + unchanged. It is the default on 10x geometry, which **changes where + output files are written** on those runs; see `DIVERGENCE.md` §3.3. + +- On 10x geometry (`CB_UMI_Simple`, a whitelist, 16 bp CB, 10 or 12 bp + UMI), the five CellRanger-matching flags now **default** to their + CellRanger values. Any flag named on the command line wins, and the + substitution is logged. **This changes default output on 10x runs** and + diverges from STARsolo; see `DIVERGENCE.md` §1.3. + +- `--soloOutRawBarcodes Observed` writes the raw matrix with one column + per *observed* barcode instead of one per whitelist barcode, matching + what CellRanger's `raw_feature_bc_matrix` contains. Counts are + unchanged; on a 200-cell run `barcodes.tsv` goes from 62 MB to 3.4 kB. + **Not a STAR parameter**; default `Whitelist` keeps STARsolo behaviour. + - **STARsolo single-cell quantification (`--soloType`)** — the 10x Chromium / plate-based count-matrix pipeline, ported from STAR and verified against real STARsolo (#90). @@ -99,6 +140,12 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Bug fixes +- `--soloUMIfiltering MultiGeneUMI_CR` kept every gene tied at the + highest read count; CellRanger gives a tied UMI to no gene at all. + Since one read per gene is the ordinary shape of a multi-gene UMI, the + flag removed nothing in practice. On a 20k-read 10x fixture the count + matrix moves from 16 465 to 15 414 against STAR's 15 423. + - **STARsolo `Gene` assignment now requires exon concordance**, matching STARsolo: a read counts toward a gene only when every aligned block lies within the gene's exons, rather than merely overlapping one. This diff --git a/Cargo.toml b/Cargo.toml index 8a4638f..444553d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -68,6 +68,9 @@ noodles-bgzf = { version = "0.49", features = ["libdeflate"] } [dev-dependencies] assert_cmd = "2" +# Already a runtime dependency at the same version; listed here so the +# integration tests can read the gzipped CellRanger-layout matrix files. +flate2 = { version = "1", default-features = false, features = ["zlib-rs"] } predicates = "3" [build-dependencies] diff --git a/DIVERGENCE.md b/DIVERGENCE.md index bd957ed..7d6f44d 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -65,6 +65,195 @@ On the 10k yeast PE benchmark, 4 reads differ in alignment score (AS) because ST --- +### 1.3 CellRanger behaviour is the default on 10x geometry + +**What STAR does.** STARsolo's defaults are its own (`1MM_multi`, +`1MM_All`, no UMI filtering, `Hamming` clipping, `outFilterScoreMin 0`) +whatever the barcode geometry. Matching CellRanger requires passing five flags, +listed in STAR's `docs/STARsolo.md`. + +**What rustar-aligner does.** When the run is unambiguously 10x — +`CB_UMI_Simple`, a whitelist, a 16-base CB and a 10- or 12-base UMI — those +five flags default to their CellRanger values. Any flag given on the command +line wins, and the substitution is logged in full. + +**Why.** A user aligning 10x data and comparing against CellRanger otherwise +gets a successful run and different numbers, with nothing pointing at the five +flags that explain it. Measured against CellRanger 10.0.0 on a 20 000-read +fixture, those flags move the count matrix from 8.96% above CellRanger to +0.03% above it, once #165's `cbMinP` posterior threshold is also applied. +STAR 2.7.11b with the same flags is at +0.09%, so all three agree to within a +fraction of a percent. + +**Impact.** This is a **change of default output behaviour** and therefore the +largest divergence in this file. It is confined to a geometry nothing else in +common use shares, it is escapable by naming any flag explicitly, and it is +announced at `INFO` on every run it touches. It needs maintainer sign-off. + +**Source.** `src/params/mod.rs` (`looks_like_10x`, +`apply_cellranger_defaults_on_10x`). STAR: `docs/STARsolo.md`, "Matching +CellRanger 4.x and 5.x results". + +--- + +### 3.2 `--soloOutRawBarcodes Observed` (opt-in, non-STAR) + +**What STAR does.** STARsolo's raw matrix has one column per whitelist +barcode, whether or not any read carried it. For 10x v3 that is 3 686 400 +columns and a 62 MB `barcodes.tsv`, nearly all of it zeros. + +**What rustar-aligner does.** The same, by default. `--soloOutRawBarcodes +Observed` narrows the raw matrix to the barcodes that actually hold a count, +which is what CellRanger's `raw_feature_bc_matrix` contains. On a 200-cell +fixture that is 200 columns and a 3.4 kB `barcodes.tsv`. + +**Why.** Someone comparing our raw matrix against CellRanger's finds no +overlapping keys at all, because the two files mean different things by "raw". +The flag makes the comparison possible without changing what STARsolo users +get. + +**Impact.** The counts are identical either way — same entries, same values, +verified on the fixture — only the columns present differ. This is a non-STAR +flag and needs maintainer sign-off; it is off by default so STARsolo parity is +untouched. + +**Source.** `src/solo/count.rs` (`observed_barcodes`), `src/params/mod.rs` +(`solo_out_raw_barcodes`). CellRanger: `outs/raw_feature_bc_matrix/` from a +`cellranger count` run, observed directly rather than taken from its source. + +--- + +### 3.3 `--soloOutLayout CellRanger` (non-STAR; on by default on 10x geometry) + +**What STAR does.** STARsolo writes +`Solo.out//{raw,filtered}/{matrix.mtx, barcodes.tsv, features.tsv}`, +uncompressed, with bare barcodes. + +**What rustar-aligner does.** The same, by default, on non-10x geometry. +`--soloOutLayout CellRanger` writes the same numbers in the shape +`cellranger count` produces: `outs/{raw,filtered}_feature_bc_matrix/`, all three +files gzipped, a `-1` GEM-well suffix on every barcode, one raw column per +observed barcode, and no per-feature subdirectory when a single feature is +requested. It implies `--soloOutGzip yes`, `--soloOutRawBarcodes Observed` and +`--soloOutFileNames outs/ ...`, each still overridable on the command line. + +Under §1.3 the same 10x detection turns this on by default, so a bare 10x run +lands in CellRanger's layout. + +**Why.** Tools written against CellRanger's `outs/` (scanpy's `read_10x_mtx`, +Seurat's `Read10X`, any in-house loader) key on those directory names and on +the `-1` suffix. Without them the numbers are right and nothing downstream can +read them without a rename step. + +**Impact.** No count changes: verified entry by entry on the 20 000-read +fixture, 13 959 entries and 15 439 counts in both layouts. On 10x geometry it +**changes where output files are written**, which needs maintainer sign-off +alongside §1.3. Against a real `cellranger count` run on the same fixture the +raw barcode sets match exactly, 200 of 200. + +**Source.** `src/solo/count.rs` (`write_gene_matrix`, `write_one_barcode`), +`src/params/mod.rs` (`solo_out_layout`, `apply_cellranger_layout`). CellRanger: +`outs/` from a `cellranger count` 10.0.0 run, observed directly rather than +taken from its source. + +--- + +### 3.4 `metrics_summary.csv` under the CellRanger layout (non-STAR) + +**What STAR does.** STARsolo writes `Summary.csv`, its own metric set with its +own names. It has no `metrics_summary.csv`. + +**What rustar-aligner does.** The same, by default. Under +`--soloOutLayout CellRanger` it *additionally* writes +`metrics_summary.csv` with CellRanger 10.0.0's 20 metrics, in CellRanger's +order and value formats. `Summary.csv` is written unchanged alongside it, so +nothing STARsolo-faithful is altered. Under the CellRanger layout the older +`CellRanger.summary.csv` is not written: its four rows are a subset of the +metrics file. + +**Why.** The file is how a 10x pipeline reads run quality. Its 20 names are +also the clearest public statement of what CellRanger measures, which makes it +useful as a target even where our value differs. + +**Impact.** No count changes. It costs one extra gene-body overlap query per +read, because the exonic/intronic split needs it and a `Gene`-only run does not +otherwise do it. + +**Not all 20 are certainties.** Twelve match a real `cellranger count` 10.0.0 +run exactly on the 20 000-read fixture. The other eight follow from our own +tallies under a stated interpretation, because CellRanger does not document the +denominators: + +* `Reads Mapped Confidently to *` reads "confidently" as MAPQ 255, i.e. our + uniquely-mapped set. +* `Reads Mapped Confidently to Transcriptome` is the `Gene` feature's own + uniquely-assigned read tally. +* `Valid UMI Sequences` is measured over reads that reached the UMI check, i.e. + those with a valid barcode. +* `Sequencing Saturation` is `1 - molecules / reads` over the reads that + entered the matrix. This is the largest disagreement on the fixture: 12.7% + against CellRanger's 7.4%. 10x define it as + `1 - n_deduped_reads / n_reads`, with `n_deduped_reads` the number of unique + `(barcode, UMI, gene)` combinations among confidently mapped reads. Taking + that literally — counting distinct triples *before* UMI correction — gives + **0.0%** on this fixture, because no two reads here share an exact triple, so + the literal reading is ruled out and their numerator is the corrected + molecule count, as ours is. The residual is therefore the denominator: at + 7.4% theirs implies about 16 300 reads where ours counts about 17 300. The + difference is which reads are "confidently mapped", not the formula. +* `Mean Reads per Cell` is total reads over called cells, not reads-in-cells + over called cells, which is what reproduces CellRanger's value. + +The cell-count-dependent metrics (`Median Genes per Cell`, +`Median UMI Counts per Cell`, `Total Genes Detected`) differ by 2-4 on the +fixture, following the count differences recorded in §1.3 rather than a +different definition. + +**Source.** `src/solo/count.rs` (`write_metrics_summary`, `metric_int`, +`metric_pct`), `src/solo/mod.rs` (`Q30Stats`). CellRanger: the +`metrics_summary.csv` of a `cellranger count` 10.0.0 run, observed directly +rather than taken from its source. + +--- + +### 3.5 `--soloCellFilter OrdMag`, and EmptyDrops_CR's initial cell set (non-STAR) + +**What STAR does.** STARsolo's cell calling is `CellRanger2.2`: take the +`quantile` of the top `nExpectedCells` barcodes by UMI count and call every +barcode holding at least that over `maxMinRatio`, with `nExpectedCells` fixed +at 3 000. `EmptyDrops_CR` uses the same knee for the set of cells it guarantees +before the Monte-Carlo rescue. + +**What CellRanger does.** The same rule, but it *searches* for the expected +cell count instead of assuming one: it minimises `(OrdMag(x) - x)^2 / x` over +`x` from 2 to about 45 000, where `OrdMag(x)` is the number of cells the rule +calls when told to expect `x`. The loss is small where the rule predicts +itself. EmptyDrops then rescues barcodes below that cutoff. + +**What rustar-aligner does.** `--soloCellFilter OrdMag maxExpectedCells +quantile ratio` (default `45000 0.99 10`) implements the search, and +`EmptyDrops_CR` now uses it for its initial set, which is the order CellRanger +runs the two steps in. `CellRanger2.2` is untouched and remains the default. + +**Two things the 10x page does not specify, decided here.** Every integer in +range is evaluated rather than a geometric grid, so the result is the true +minimum of the stated loss rather than a nearby grid point; each evaluation is +a binary search over the sorted totals, so an exhaustive sweep costs nothing +measurable. And ties go to the smaller `x`, so the result does not depend on +iteration order — determinism, as everywhere else in this codebase. + +**Impact.** `EmptyDrops_CR`'s initial cell set can differ from what it was. +On the 20 000-read fixture it does not: `CellRanger2.2`, `OrdMag` and +`EmptyDrops_CR` all call 200 cells, as does CellRanger 10.0.0. That fixture has +a clean plateau, where any of these rules works; the two rules part company on +a graded distribution with no plateau, which is what the unit tests cover. + +**Source.** `src/solo/count.rs` (`ordmag_at`, `ordmag_threshold`). CellRanger: +the Gene Expression algorithm page, "Cell Calling", read rather than taken from +its source. Coverage of the rest of that page is tracked in #181. + +--- + ## 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/params/mod.rs b/src/params/mod.rs index 0536a85..b87b2c4 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1113,7 +1113,15 @@ pub struct Parameters { #[arg(long = "soloCBmatchWLtype", default_value = "1MM_multi")] pub solo_cb_match_wl_type: String, - /// Cell-calling / matrix filtering: None, CellRanger2.2, EmptyDrops_CR, TopCells. + /// Cell-calling / matrix filtering: None, CellRanger2.2, EmptyDrops_CR, + /// TopCells, OrdMag. + /// + /// `OrdMag` is CellRanger's own initial cell call and is **not a STAR + /// method**: the same quantile-over-ratio rule as `CellRanger2.2`, but with + /// the expected cell count searched for rather than fixed at 3 000. Its + /// arguments are `maxExpectedCells quantile ratio`, default + /// `45000 0.99 10`. `EmptyDrops_CR` now uses it for its initial set, which + /// is the order CellRanger runs the two steps in. #[arg(long = "soloCellFilter", num_args = 1.., default_values_t = vec!["CellRanger2.2".to_string(), "3000".to_string(), "0.99".to_string(), "10".to_string()])] pub solo_cell_filter: Vec, @@ -1133,6 +1141,44 @@ pub struct Parameters { #[arg(long = "soloOutGzip", default_value = "no")] pub solo_out_gzip: String, + /// Which barcodes the **raw** matrix has columns for. **Not a STAR + /// parameter**; a rustar-aligner addition, default `Whitelist`, which is + /// what STARsolo writes. + /// + /// `Whitelist` gives one column per whitelist barcode — 3.7 million of them + /// for 10x v3, whether or not a read ever carried them. `Observed` gives + /// one column per barcode that actually holds a count, which is what + /// CellRanger's `raw_feature_bc_matrix` contains, and turns a + /// hundreds-of-megabytes `barcodes.tsv` into a few kilobytes. + /// + /// The counts are identical either way; only the columns present differ. + #[arg(long = "soloOutRawBarcodes", default_value = "Whitelist", + value_parser = ["Whitelist", "Observed"])] + pub solo_out_raw_barcodes: String, + + /// Shape of the solo matrix output on disk. **Not a STAR parameter**; a + /// rustar-aligner addition, default `STARsolo`, which changes nothing. + /// + /// `CellRanger` lays the same numbers out the way `cellranger count` does, + /// so a tool written against CellRanger's `outs/` reads a rustar run + /// unmodified. It implies, unless the corresponding flag is given + /// explicitly: + /// + /// * directories `raw_feature_bc_matrix/` and `filtered_feature_bc_matrix/` + /// instead of `raw/` and `filtered/`; + /// * a `-1` GEM-well suffix on every barcode; + /// * gzip on all three files (`--soloOutGzip yes`); + /// * `--soloOutRawBarcodes Observed`, since CellRanger's raw matrix has one + /// column per observed barcode, not one per whitelist entry; + /// * the output directory `outs/` instead of `Solo.out/`, with no + /// per-feature subdirectory when exactly one feature is requested. + /// + /// Counts are untouched. Only where the bytes land, and how the barcodes + /// are spelled, changes. + #[arg(long = "soloOutLayout", default_value = "STARsolo", + value_parser = ["STARsolo", "CellRanger"])] + pub solo_out_layout: String, + /// Velocyto ambiguous-molecule handling (rustar extension beyond STARsolo). /// `yes` (default) writes the three `spliced`/`unspliced`/`ambiguous` matrices /// like STARsolo — exon-only molecules with no junction/intron evidence stay in @@ -1348,6 +1394,9 @@ impl Parameters { let matches = command.clone().get_matches_from(args.iter()); let mut params = ::from_arg_matches(&matches)?; + apply_cellranger_defaults_on_10x(&mut params, &matches); + apply_cellranger_layout(&mut params, &matches); + params.command_line = { let args: Vec<_> = args.iter().map(|s| s.to_string_lossy()).collect(); shlex::try_join(args.iter().map(AsRef::as_ref)).ok() @@ -1870,8 +1919,389 @@ 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); 7] = [ + ("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"), + // CellRanger has counted intronic reads toward the gene since v7.0. + // STARsolo's `Gene` is exonic-only, which on human data is 30% below + // CellRanger; `GeneFull` is the equivalent of its default. + ("solo_features", "GeneFull"), +]; + +/// Does this look like a 10x Chromium run? +/// +/// `CB_UMI_Simple` with a whitelist, a 16-base cell barcode, and a UMI of 10 +/// (v2) or 12 (v3) bases. That is the geometry of every 10x 3'/5' gene +/// expression chemistry, and nothing else in common use shares it. +fn looks_like_10x(params: &Parameters) -> bool { + params.solo_type == SoloType::CbUmiSimple + && params.solo_cb_len == 16 + && (params.solo_umi_len == 10 || params.solo_umi_len == 12) + && params + .solo_cb_whitelist + .first() + .is_some_and(|w| w != "None" && w != "-") +} + +/// On a 10x run, default to CellRanger's behaviour rather than STARsolo's. +/// +/// **This diverges from STAR by default**, which is why it is confined to a +/// geometry that is unambiguously 10x, and why every flag it changes is +/// logged. A flag given on the command line always wins, so the change is +/// invisible to anyone who states what they want. +/// +/// The rationale is that a user aligning 10x data and comparing against +/// CellRanger currently gets a successful run and different numbers, with +/// nothing pointing at the five flags that explain the difference. Measured on +/// a 20 000-read fixture, those flags move the count matrix from 8.9% away +/// from CellRanger to 0.03%. +/// +/// Recorded in `DIVERGENCE.md`; it needs maintainer sign-off. +fn apply_cellranger_defaults_on_10x(params: &mut Parameters, matches: &clap::ArgMatches) { + use clap::parser::ValueSource; + + if !looks_like_10x(params) { + return; + } + + let given = |id: &str| matches.value_source(id) == Some(ValueSource::CommandLine); + + let mut applied: Vec<&str> = Vec::new(); + for (id, value) in CELLRANGER_DEFAULTS { + if given(id) { + continue; + } + match id { + "clip_adapter_type" => params.clip_adapter_type = value.to_string(), + "out_filter_score_min" => params.out_filter_score_min = 30, + "solo_cb_match_wl_type" => params.solo_cb_match_wl_type = value.to_string(), + "solo_umi_filtering" => params.solo_umi_filtering = vec![value.to_string()], + "solo_umi_dedup" => params.solo_umi_dedup = vec![value.to_string()], + "solo_out_layout" => params.solo_out_layout = value.to_string(), + "solo_features" => params.solo_features = 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(", ") + ); + } +} + +/// `--soloOutLayout CellRanger` implies three existing flags and the output +/// directory name. Each is still overridable: a flag named on the command line +/// keeps its value, so the layout can be adopted piecemeal. +/// +/// Runs after `apply_cellranger_defaults_on_10x`, so it sees the layout whether +/// it was asked for or inferred from 10x geometry. +fn apply_cellranger_layout(params: &mut Parameters, matches: &clap::ArgMatches) { + use clap::parser::ValueSource; + + if params.solo_out_layout != "CellRanger" { + return; + } + let given = |id: &str| matches.value_source(id) == Some(ValueSource::CommandLine); + + if !given("solo_out_gzip") { + params.solo_out_gzip = "yes".to_string(); + } + if !given("solo_out_raw_barcodes") { + params.solo_out_raw_barcodes = "Observed".to_string(); + } + if !given("solo_out_file_names") + && let Some(dir) = params.solo_out_file_names.first_mut() + { + // CellRanger writes its matrices under `outs/`, not `Solo.out/`. + *dir = "outs/".to_string(); + } +} + #[cfg(test)] mod tests { + + /// 10x geometry with a whitelist gets CellRanger's five flags without the + /// user naming any of them. This is a deliberate divergence from STARsolo's + /// defaults, so the test states the whole set rather than spot-checking one. + #[test] + fn ten_x_geometry_defaults_to_cellranger_behaviour() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--soloCBstart", + "1", + "--soloCBlen", + "16", + "--soloUMIstart", + "17", + "--soloUMIlen", + "12", + ]) + .unwrap(); + assert_eq!(p.clip_adapter_type, "CellRanger4"); + assert_eq!(p.out_filter_score_min, 30); + assert_eq!(p.solo_cb_match_wl_type, "1MM_multi_Nbase_pseudocounts"); + assert_eq!(p.solo_umi_filtering, vec!["MultiGeneUMI_CR".to_string()]); + assert_eq!(p.solo_umi_dedup, vec!["1MM_CR".to_string()]); + assert_eq!(p.solo_out_layout, "CellRanger"); + assert_eq!(p.solo_features, vec!["GeneFull".to_string()]); + } + + /// `--soloFeatures` given explicitly wins, so a user who wants STARsolo's + /// exonic-only counting on 10x data still gets it. + #[test] + fn an_explicit_solo_features_beats_the_10x_intron_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", + "--soloFeatures", + "Gene", + ]) + .unwrap(); + assert_eq!(p.solo_features, vec!["Gene".to_string()]); + // The rest of the CellRanger defaults still apply. + assert_eq!(p.solo_umi_dedup, vec!["1MM_CR".to_string()]); + } + + /// `--soloOutLayout CellRanger` pulls three existing flags and the output + /// directory with it, so the layout is one decision rather than four. + #[test] + fn cellranger_layout_implies_gzip_observed_barcodes_and_outs_dir() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--soloCBstart", + "1", + "--soloCBlen", + "12", + "--soloUMIstart", + "13", + "--soloUMIlen", + "10", + "--soloOutLayout", + "CellRanger", + ]) + .unwrap(); + // 12-base CB, so the 10x autodetection is not what set this. + assert_eq!(p.solo_out_layout, "CellRanger"); + assert_eq!(p.solo_out_gzip, "yes"); + assert_eq!(p.solo_out_raw_barcodes, "Observed"); + assert_eq!(p.solo_out_file_names.first().unwrap(), "outs/"); + } + + /// The layout is adoptable piecemeal: each flag it implies is still + /// overridable on the command line. + #[test] + fn an_explicit_flag_beats_what_the_cellranger_layout_implies() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--soloCBstart", + "1", + "--soloCBlen", + "16", + "--soloUMIstart", + "17", + "--soloUMIlen", + "12", + "--soloOutGzip", + "no", + "--soloOutRawBarcodes", + "Whitelist", + "--soloOutFileNames", + "Solo.out/", + "features.tsv", + "barcodes.tsv", + "matrix.mtx", + ]) + .unwrap(); + assert_eq!(p.solo_out_layout, "CellRanger"); + assert_eq!(p.solo_out_gzip, "no"); + assert_eq!(p.solo_out_raw_barcodes, "Whitelist"); + assert_eq!(p.solo_out_file_names.first().unwrap(), "Solo.out/"); + } + + /// The default is STARsolo's layout: no `-1`, no `outs/`, no gzip. + #[test] + fn non_10x_geometry_keeps_the_starsolo_layout() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--soloCBstart", + "1", + "--soloCBlen", + "12", + "--soloUMIstart", + "13", + "--soloUMIlen", + "10", + ]) + .unwrap(); + assert_eq!(p.solo_out_layout, "STARsolo"); + assert_eq!(p.solo_out_gzip, "no"); + assert_eq!(p.solo_out_raw_barcodes, "Whitelist"); + assert_eq!(p.solo_out_file_names.first().unwrap(), "Solo.out/"); + } + + /// A flag given on the command line always wins, including when the value + /// asked for is STARsolo's own default. Without this the divergence would + /// be inescapable, which is a different and much worse thing than a + /// divergent default. + #[test] + fn an_explicit_flag_beats_the_10x_default() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--soloCBstart", + "1", + "--soloCBlen", + "16", + "--soloUMIstart", + "17", + "--soloUMIlen", + "12", + "--soloUMIdedup", + "1MM_All", + "--clipAdapterType", + "Hamming", + ]) + .unwrap(); + assert_eq!(p.solo_umi_dedup, vec!["1MM_All".to_string()]); + assert_eq!(p.clip_adapter_type, "Hamming"); + // The ones not named still take the CellRanger value. + assert_eq!(p.out_filter_score_min, 30); + } + + /// Geometry that is not 10x is left alone: a 12-base barcode is not any + /// Chromium chemistry, so nothing is overridden. + #[test] + fn non_10x_geometry_keeps_starsolo_defaults() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--soloCBstart", + "1", + "--soloCBlen", + "12", + "--soloUMIstart", + "13", + "--soloUMIlen", + "8", + ]) + .unwrap(); + assert_eq!(p.clip_adapter_type, "Hamming"); + assert_eq!(p.out_filter_score_min, 0); + assert_eq!(p.solo_cb_match_wl_type, "1MM_multi"); + } + + /// No whitelist means no 10x run, whatever the lengths say. (Without a + /// whitelist the CB-match type must be Exact anyway, which is unrelated + /// validation that predates this and is stated here so the test reads.) + #[test] + fn ten_x_lengths_without_a_whitelist_keep_starsolo_defaults() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBstart", + "1", + "--soloCBlen", + "16", + "--soloUMIstart", + "17", + "--soloUMIlen", + "12", + "--soloCBmatchWLtype", + "Exact", + ]) + .unwrap(); + assert_eq!(p.clip_adapter_type, "Hamming"); + assert_eq!(p.out_filter_score_min, 0); + } use super::*; /// Helper: parse a STAR-style command line (without program name). diff --git a/src/solo/count.rs b/src/solo/count.rs index 9af3edf..a4b9314 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -167,6 +167,18 @@ pub fn dedup_count(umis: &HashMap, method: UmiDedup, umi_len: usize) - /// the neighbor's raw UMI, not its corrected value); the molecule count is the /// number of distinct corrected UMIs. fn cellranger_1mm(umis: &HashMap, umi_len: usize) -> u64 { + let distinct: std::collections::HashSet = + cellranger_1mm_map(umis, umi_len).into_values().collect(); + distinct.len() as u64 +} + +/// The same correction, returning `raw UMI -> corrected UMI`. +/// +/// `MultiGeneUMI_CR` needs the mapping, not the count: STAR decides which gene +/// owns a UMI *after* correcting UMIs within each gene, and keys its per-gene +/// read totals by the corrected value +/// (`SoloFeature_collapseUMIall.cpp:134-148`). +fn cellranger_1mm_map(umis: &HashMap, umi_len: usize) -> HashMap { let mut items: Vec<(u64, u32)> = umis.iter().map(|(&u, &c)| (u, c)).collect(); // Ascending by count, then by UMI value (mirrors funCompareSolo1 ordering, // so the inner scan from the end meets higher-count neighbors first). @@ -185,8 +197,7 @@ fn cellranger_1mm(umis: &HashMap, umi_len: usize) -> u64 { } corrected.push(corr); } - let distinct: std::collections::HashSet = corrected.into_iter().collect(); - distinct.len() as u64 + items.iter().map(|&(u, _)| u).zip(corrected).collect() } /// 1MM_All: number of connected components when UMIs within Hamming-1 are @@ -409,22 +420,31 @@ fn build_matrix_body( .or_insert(0) += 1; } - // (gene → (umi → read_count)) after multi-gene UMI filtering. - let mut gene_umis: HashMap> = HashMap::default(); - for (&umi, genes) in &umi_genes { - for (&gene, &rc) in filter_multi_gene_umi(genes, filtering) { - *gene_umis.entry(gene).or_default().entry(umi).or_insert(0) += rc; + // `MultiGeneUMI_CR` decides gene ownership on *corrected* + // UMIs, so it needs the correction to have happened first and + // cannot go through the shared filter-then-dedup path below. + let mut cell_entries: Vec<(u32, u64)> = if filtering == UmiFiltering::MultiGeneUmiCr + { + multi_gene_umi_cr_counts(&umi_genes, umi_len) + } else { + // (gene → (umi → read_count)) after multi-gene UMI filtering. + let mut gene_umis: HashMap> = HashMap::default(); + for (&umi, genes) in &umi_genes { + for (&gene, &rc) in filter_multi_gene_umi(genes, filtering) { + *gene_umis.entry(gene).or_default().entry(umi).or_insert(0) += rc; + } } - } - // Collapse UMIs per gene, then emit this cell's entries gene-ascending. - let mut cell_entries: Vec<(u32, u64)> = Vec::with_capacity(gene_umis.len()); - for (&gene, umis) in &gene_umis { - let count = dedup_count(umis, method, umi_len); - if count > 0 { - cell_entries.push((gene, count)); + // Collapse UMIs per gene, then emit gene-ascending. + let mut entries: Vec<(u32, u64)> = Vec::with_capacity(gene_umis.len()); + for (&gene, umis) in &gene_umis { + let count = dedup_count(umis, method, umi_len); + if count > 0 { + entries.push((gene, count)); + } } - } + entries + }; cell_entries.sort_unstable_by_key(|&(g, _)| g); let n_reads = (j - i) as u64; @@ -476,6 +496,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 @@ -808,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)> { @@ -822,8 +945,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!(), } } @@ -840,6 +998,61 @@ fn knee_cr22(umis_desc: &[u64], n_expected: usize, max_pct: f64, max_min_ratio: (robust_max / max_min_ratio).ceil() as u64 } +/// How many cells the order-of-magnitude rule calls if `n` cells are expected, +/// and the UMI cutoff it used. +/// +/// The rule, from 10x's Gene Expression algorithm page: take `m`, the +/// `quantile` of the top `n` barcodes by UMI count, and call every barcode +/// holding at least `m / ratio`. With the defaults that is "the 99th percentile +/// of the top n, divided by ten", the same shape as `knee_cr22` but with `n` a +/// variable rather than a fixed 3 000. +fn ordmag_at(umis_desc: &[u64], n: usize, quantile: f64, ratio: f64) -> (usize, u64) { + if umis_desc.is_empty() || n == 0 { + return (0, 0); + } + let idx = ((n as f64 * (1.0 - quantile)).round() as usize).min(umis_desc.len() - 1); + let cutoff = ((umis_desc[idx] as f64 / ratio).round() as u64).max(1); + // `umis_desc` is descending, so the called set is its prefix. + let called = umis_desc.partition_point(|&u| u >= cutoff); + (called, cutoff) +} + +/// CellRanger's OrdMag cell-calling threshold: the UMI cutoff at the expected +/// cell count that best predicts itself. +/// +/// 10x describes this as minimising `(OrdMag(x) - x)^2 / x` over a search from +/// 2 to about 45 000 cells. `OrdMag(x)` is the number of cells the rule above +/// calls when told to expect `x` of them, so the loss is small where the rule +/// is self-consistent: feed it the right number of cells and it gives that +/// number back. +/// +/// Two choices this makes explicit, because the page does not state them: +/// +/// * **Every integer in range is evaluated**, rather than a geometric grid +/// refined by search. Each evaluation is a binary search over the sorted +/// totals, so an exhaustive sweep of 45 000 candidates costs nothing +/// measurable and finds the true minimum of the stated loss rather than a +/// grid point near it. +/// * **Ties go to the smaller `x`**, so the result does not depend on iteration +/// order. +fn ordmag_threshold(umis_desc: &[u64], max_expected: usize, quantile: f64, ratio: f64) -> u64 { + if umis_desc.is_empty() { + return 0; + } + let hi = max_expected.min(umis_desc.len()).max(2); + let mut best: Option<(f64, u64)> = None; + for x in 2..=hi { + let (called, cutoff) = ordmag_at(umis_desc, x, quantile, ratio); + let d = called as f64 - x as f64; + let loss = d * d / x as f64; + // Strictly-less keeps the first (smallest) x on a tie. + if best.is_none_or(|(b, _)| loss < b) { + best = Some((loss, cutoff)); + } + } + best.map_or(0, |(_, cutoff)| cutoff) +} + /// Whitelist indices of called cells (sorted ascending) per `--soloCellFilter`. /// `None` → no filtered/ output. `EmptyDrops_CR` writes only the knee-guaranteed /// cells here (the Monte-Carlo rescue is the standalone `emptydrops` binary). @@ -854,6 +1067,17 @@ fn called_cells(cells: &[CellStat], filter: &[String]) -> Option> { idx.sort_by(|a, b| b.n_umis.cmp(&a.n_umis).then(a.cb.cmp(&b.cb))); idx.into_iter().take(n).map(|c| c.cb).collect() } + // CellRanger's own initial cell set, searched rather than assumed. + "OrdMag" => { + let mut umis: Vec = cells.iter().map(|c| c.n_umis).collect(); + umis.sort_unstable_by(|a, b| b.cmp(a)); + let thr = ordmag_threshold(&umis, arg(1, 45000.0) as usize, arg(2, 0.99), arg(3, 10.0)); + cells + .iter() + .filter(|c| c.n_umis >= thr) + .map(|c| c.cb) + .collect() + } // EmptyDrops_CR is handled by `emptydrops_called`; the knee here is the // fallback / guaranteed-cell base. "CellRanger2.2" | "EmptyDrops_CR" => { @@ -893,7 +1117,9 @@ fn emptydrops_called( .and_then(|s| s.parse::().ok()) .unwrap_or(d) }; - let (n_expected, max_pct, ratio) = (arg(1, 3000.0) as usize, arg(2, 0.99), arg(3, 10.0)); + // `nExpectedCells` (arg 1) is the 2.2 knee's fixed cell count; OrdMag + // searches for it instead, so this path no longer reads it. + let (max_pct, ratio) = (arg(2, 0.99), arg(3, 10.0)); let (ind_min, ind_max) = (arg(4, 45000.0) as usize, arg(5, 90000.0) as usize); let umi_min = arg(6, 500.0) as u64; let umi_min_frac = arg(7, 0.01); @@ -905,7 +1131,12 @@ fn emptydrops_called( let mut order: Vec<&CellStat> = cells.iter().collect(); order.sort_by(|a, b| b.n_umis.cmp(&a.n_umis).then(a.cb.cmp(&b.cb))); let totals_desc: Vec = order.iter().map(|c| c.n_umis).collect(); - let thr = knee_cr22(&totals_desc, n_expected, max_pct, ratio); + // CellRanger's initial cell set is OrdMag, not the 2.2 knee: the same + // quantile-over-ratio rule, but with the expected cell count searched for + // rather than fixed. EmptyDrops then rescues barcodes below it. `ind_min` + // is already the ~45 000 upper bound 10x states for that search, so the + // parameter list is unchanged. + let thr = ordmag_threshold(&totals_desc, ind_min, max_pct, ratio); let n_simple = totals_desc.iter().take_while(|&&u| u >= thr).count(); let mut called: Vec = order.iter().take(n_simple).map(|c| c.cb).collect(); @@ -1163,8 +1394,12 @@ pub fn write_gene_matrix( .position(|&x| x == f) .map_or(0, |i| ctx.feature_reads[i].load(Ordering::Relaxed)) }; - let have_funnel = ctx.features.contains(&crate::solo::SoloFeature::Gene) - && ctx.features.contains(&crate::solo::SoloFeature::GeneFull); + // The split needs the gene-body query as well as the exon one. Both `Gene` + // + `GeneFull` gives it, and so does `ctx.want_metrics`, which asks + // `process_read` for the body query on its own account. + let have_funnel = ctx.want_metrics + || (ctx.features.contains(&crate::solo::SoloFeature::Gene) + && ctx.features.contains(&crate::solo::SoloFeature::GeneFull)); let region = have_funnel.then(|| RegionFunnel { exonic: ctx.region_stats.exonic.load(Ordering::Relaxed), intronic: ctx.region_stats.intronic.load(Ordering::Relaxed), @@ -1176,10 +1411,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 @@ -1200,26 +1453,58 @@ pub fn write_gene_matrix( &ctx.gene_ann.gene_names, gzip, )?; - write_barcodes( - &raw_dir.join(&barcodes_name), - &ctx.whitelist, - sorted.len(), - gzip, - )?; + // `--soloOutRawBarcodes Observed` narrows the raw matrix to the + // barcodes that actually carry a count, which is what CellRanger's + // `raw_feature_bc_matrix` holds. STARsolo's raw matrix has a column per + // whitelist barcode, so the default keeps that. + let observed: Option> = if params.solo_out_raw_barcodes == "Observed" { + Some(observed_barcodes(&body)?) + } else { + None + }; + let (raw_cols, raw_remap) = match &observed { + Some(cbs) => { + let map: HashMap = cbs + .iter() + .enumerate() + .map(|(col, &cb)| (cb, col as u32 + 1)) + .collect(); + (cbs.len(), Some(map)) + } + None => (sorted.len(), None), + }; + match &observed { + Some(cbs) => { + write_barcodes_subset( + &raw_dir.join(&barcodes_name), + &ctx.whitelist, + cbs, + gzip, + cb_suffix, + )?; + } + None => write_barcodes( + &raw_dir.join(&barcodes_name), + &ctx.whitelist, + sorted.len(), + gzip, + cb_suffix, + )?, + } finalize_matrix( &body, &raw_dir.join(&matrix_name), gzip, n_genes, - sorted.len(), + raw_cols, mstats.nnz, - None, + raw_remap.as_ref(), )?; log::info!( - "STARsolo: wrote {}/raw matrix ({} genes × {} barcodes, {} entries){}", - feature.dir_name(), + "STARsolo: wrote {} matrix ({} genes × {} barcodes, {} entries){}", + raw_dir.display(), n_genes, - sorted.len(), + raw_cols, mstats.nnz, if gzip { " [gzip]" } else { "" }, ); @@ -1243,7 +1528,7 @@ pub fn write_gene_matrix( if let Some(cbs) = called && !cbs.is_empty() { - let filt_dir = feature_dir.join("filtered"); + let filt_dir = feature_dir.join(filt_name); std::fs::create_dir_all(&filt_dir).map_err(|e| Error::io(e, &filt_dir))?; let remap: HashMap = cbs .iter() @@ -1256,7 +1541,13 @@ pub fn write_gene_matrix( &ctx.gene_ann.gene_names, gzip, )?; - write_barcodes_subset(&filt_dir.join(&barcodes_name), &ctx.whitelist, &cbs, gzip)?; + write_barcodes_subset( + &filt_dir.join(&barcodes_name), + &ctx.whitelist, + &cbs, + gzip, + cb_suffix, + )?; let fnnz = finalize_matrix( &body, &filt_dir.join(&matrix_name), @@ -1267,8 +1558,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, ); @@ -1306,10 +1597,32 @@ pub fn write_gene_matrix( reads_of(*feature), )?; log::info!("STARsolo: wrote {}/Summary.csv", feature.dir_name()); + // Under the CellRanger layout the funnel goes into `metrics_summary.csv` + // alongside the other 19 metrics, which is where CellRanger puts it. + if cr_layout && let Some(r) = region { + let invalid_umis = ctx.stats.n_in_umi.load(Ordering::Relaxed) + + ctx.stats.umi_homopolymer.load(Ordering::Relaxed); + write_metrics_summary( + &feature_dir.join("metrics_summary.csv"), + &mstats, + &ctx.q30, + total_reads, + valid_barcodes, + invalid_umis, + mapped_unique, + mapped_multi, + reads_of(crate::solo::SoloFeature::Gene), + r, + )?; + log::info!( + "STARsolo: wrote {}", + feature_dir.join("metrics_summary.csv").display() + ); + } // CellRanger-style mapping funnel goes in a SEPARATE additional file so the // faithful Summary.csv is never altered (PR #90 review: keep this release a // drop-in faithful port; output-changing features come later). - if let Some(r) = region { + if !cr_layout && let Some(r) = region { write_cellranger_summary( &feature_dir.join("CellRanger.summary.csv"), total_reads, @@ -1339,11 +1652,14 @@ pub fn write_gene_matrix( write_file(&sj_dir.join(&features_name), gzip, |w| { sjs.write_sj_lines(w, genome, params).map(|_| ()) })?; + // SJ has no CellRanger counterpart, but every barcode written by one run + // is spelled the same way, so the suffix applies here too. write_barcodes( &sj_dir.join(&barcodes_name), &ctx.whitelist, sorted.len(), gzip, + cb_suffix, )?; let umi_len = params.solo_umi_len as usize; let nnz = build_sj_matrix( @@ -1379,6 +1695,7 @@ pub fn write_gene_matrix( &ctx.whitelist, sorted.len(), gzip, + cb_suffix, )?; let umi_len = params.solo_umi_len as usize; // `--soloVelocytoAmbiguous no` folds exon-only molecules into spliced and @@ -1640,6 +1957,168 @@ struct RegionFunnel { antisense: u64, } +/// Format an integer the way CellRanger's `metrics_summary.csv` does: thousands +/// separated by commas, and quoted when that puts a comma in the field. +/// +/// `200` stays `200`; `20000` becomes `"20,000"`. +fn metric_int(n: u64) -> String { + let digits = n.to_string(); + if digits.len() <= 3 { + return digits; + } + let mut out = String::with_capacity(digits.len() + digits.len() / 3 + 2); + for (i, c) in digits.chars().enumerate() { + if i > 0 && (digits.len() - i).is_multiple_of(3) { + out.push(','); + } + out.push(c); + } + format!("\"{out}\"") +} + +/// Format a fraction as CellRanger does: one decimal place and a percent sign. +fn metric_pct(num: u64, den: u64) -> String { + let f = if den == 0 { + 0.0 + } else { + num as f64 / den as f64 + }; + format!("{:.1}%", f * 100.0) +} + +/// Write CellRanger's `metrics_summary.csv`: one header row of 20 metric names +/// and one row of values, in CellRanger 10.0.0's order. +/// +/// Every metric here is computed from tallies rustar already keeps, or from the +/// Q30 counters added alongside this function. Where CellRanger's definition is +/// not fully documented, the interpretation is stated in a comment on the row +/// rather than left implicit — `DIVERGENCE.md` §3.4 lists the ones that are an +/// interpretation rather than a certainty. +#[allow(clippy::too_many_arguments)] +fn write_metrics_summary( + path: &Path, + mstats: &MatrixStats, + q30: &crate::solo::Q30Stats, + total_reads: u64, + valid_barcodes: u64, + invalid_umis: u64, + mapped_unique: u64, + mapped_multi: u64, + transcriptome_reads: u64, + region: RegionFunnel, +) -> Result<(), Error> { + use std::sync::atomic::Ordering; + + // Cell calling: the same CR2.2 knee `Summary.csv` uses, so the two files + // never disagree about how many cells there are. + let mut umis_desc: Vec = mstats.cells.iter().map(|c| c.n_umis).collect(); + umis_desc.sort_unstable_by(|a, b| b.cmp(a)); + let thr = knee_cr22(&umis_desc, 3000, 0.99, 10.0); + let cells: Vec<&CellStat> = mstats.cells.iter().filter(|c| c.n_umis >= thr).collect(); + let n_cells = cells.len() as u64; + + let reads_counted: u64 = mstats.cells.iter().map(|c| c.n_reads).sum(); + let umis_all: u64 = mstats.cells.iter().map(|c| c.n_umis).sum(); + let reads_in_cells: u64 = cells.iter().map(|c| c.n_reads).sum(); + + let mut umis_sorted: Vec = cells.iter().map(|c| c.n_umis).collect(); + let mut genes_sorted: Vec = cells.iter().map(|c| c.n_genes as u64).collect(); + umis_sorted.sort_unstable(); + genes_sorted.sort_unstable(); + + // Saturation = the share of reads that added no new molecule. + let saturation_num = reads_counted.saturating_sub(umis_all); + let q = |n: &std::sync::atomic::AtomicU64| n.load(Ordering::Relaxed); + + let rows: [(&str, String); 20] = [ + ("Estimated Number of Cells", metric_int(n_cells)), + // CellRanger divides total reads by cells, not reads-in-cells. + ( + "Mean Reads per Cell", + metric_int(total_reads.checked_div(n_cells).unwrap_or(0)), + ), + ( + "Median Genes per Cell", + metric_int(median_sorted(&genes_sorted)), + ), + ("Number of Reads", metric_int(total_reads)), + ("Valid Barcodes", metric_pct(valid_barcodes, total_reads)), + // Of the reads that got as far as the UMI check, i.e. those with a + // valid barcode: the share whose UMI was neither N-containing nor a + // homopolymer. + ( + "Valid UMI Sequences", + metric_pct(valid_barcodes.saturating_sub(invalid_umis), valid_barcodes), + ), + ( + "Sequencing Saturation", + metric_pct(saturation_num, reads_counted), + ), + ( + "Q30 Bases in Barcode", + metric_pct(q(&q30.cb_q30), q(&q30.cb_bases)), + ), + ( + "Q30 Bases in RNA Read", + metric_pct(q(&q30.rna_q30), q(&q30.rna_bases)), + ), + ( + "Q30 Bases in UMI", + metric_pct(q(&q30.umi_q30), q(&q30.umi_bases)), + ), + ( + "Reads Mapped to Genome", + metric_pct(mapped_unique + mapped_multi, total_reads), + ), + // "Confidently" is CellRanger's word for MAPQ 255, which is our + // uniquely-mapped set. + ( + "Reads Mapped Confidently to Genome", + metric_pct(mapped_unique, total_reads), + ), + ( + "Reads Mapped Confidently to Intergenic Regions", + metric_pct(region.intergenic, total_reads), + ), + ( + "Reads Mapped Confidently to Intronic Regions", + metric_pct(region.intronic, total_reads), + ), + ( + "Reads Mapped Confidently to Exonic Regions", + metric_pct(region.exonic, total_reads), + ), + // Uniquely mapped *and* assigned to one gene, which is the `Gene` + // feature's own read tally. + ( + "Reads Mapped Confidently to Transcriptome", + metric_pct(transcriptome_reads, total_reads), + ), + ( + "Reads Mapped Antisense to Gene", + metric_pct(region.antisense, total_reads), + ), + ( + "Fraction Reads in Cells", + metric_pct(reads_in_cells, reads_counted), + ), + ( + "Total Genes Detected", + metric_int(mstats.genes_detected.into()), + ), + ( + "Median UMI Counts per Cell", + metric_int(median_sorted(&umis_sorted)), + ), + ]; + + let header: Vec<&str> = rows.iter().map(|(k, _)| *k).collect(); + let values: Vec<&str> = rows.iter().map(|(_, v)| v.as_str()).collect(); + let out = format!("{}\n{}\n", header.join(","), values.join(",")); + std::fs::write(path, out).map_err(|e| Error::io(e, path))?; + Ok(()) +} + /// Write the STARsolo-faithful `Summary.csv` for one feature: the sequencing / /// genome-mapping rows plus per-cell UMI/gene statistics over the CR2.2-knee-called /// cells. The CellRanger-style exonic/intronic/intergenic/antisense funnel is a @@ -1833,16 +2312,21 @@ fn write_features( Ok(()) } -/// Unpack `cb` into `line` (with trailing newline) and write it. +/// Unpack `cb` into `line` (with `suffix` and a trailing newline) and write it. +/// +/// `suffix` is `"-1"` under `--soloOutLayout CellRanger` (the GEM-well tag +/// CellRanger appends to every barcode) and `""` otherwise. fn write_one_barcode( w: &mut dyn std::io::Write, whitelist: &CbWhitelist, cb: u32, line: &mut Vec, path: &Path, + suffix: &str, ) -> Result<(), Error> { line.clear(); whitelist.unpack_barcode_into(cb, line); + line.extend_from_slice(suffix.as_bytes()); line.push(b'\n'); w.write_all(line).map_err(|e| Error::io(e, path)) } @@ -1850,12 +2334,18 @@ fn write_one_barcode( /// `barcodes.tsv`: full whitelist in sorted order (matches the raw matrix /// columns). Lists millions of lines, so the writer is buffered and the barcode /// is unpacked into a reused scratch buffer (no per-line allocation). -fn write_barcodes(path: &Path, whitelist: &CbWhitelist, n: usize, gzip: bool) -> Result<(), Error> { +fn write_barcodes( + path: &Path, + whitelist: &CbWhitelist, + n: usize, + gzip: bool, + suffix: &str, +) -> Result<(), Error> { let len = whitelist.barcode_len(); write_file(path, gzip, |w| { - let mut line: Vec = Vec::with_capacity(len + 1); + let mut line: Vec = Vec::with_capacity(len + suffix.len() + 1); for i in 0..n { - write_one_barcode(w, whitelist, i as u32, &mut line, path)?; + write_one_barcode(w, whitelist, i as u32, &mut line, path, suffix)?; } Ok(()) })?; @@ -1869,12 +2359,13 @@ fn write_barcodes_subset( whitelist: &CbWhitelist, cbs: &[u32], gzip: bool, + suffix: &str, ) -> Result<(), Error> { let len = whitelist.barcode_len(); write_file(path, gzip, |w| { - let mut line: Vec = Vec::with_capacity(len + 1); + let mut line: Vec = Vec::with_capacity(len + suffix.len() + 1); for &cb in cbs { - write_one_barcode(w, whitelist, cb, &mut line, path)?; + write_one_barcode(w, whitelist, cb, &mut line, path, suffix)?; } Ok(()) })?; @@ -1887,6 +2378,79 @@ mod tests { use crate::io::fastq::encode_base; use crate::solo::whitelist::pack_barcode; + /// `--soloOutLayout CellRanger` appends the `-1` GEM-well suffix that + /// CellRanger puts on every barcode; the default writes the bare barcode. + #[test] + fn barcodes_carry_the_gem_well_suffix_only_under_the_cellranger_layout() { + let dir = tempfile::tempdir().unwrap(); + let wl_path = dir.path().join("wl.txt"); + std::fs::write(&wl_path, "ACGT\nTGCA\n").unwrap(); + let wl = CbWhitelist::load(&wl_path).unwrap(); + + let plain = dir.path().join("plain.tsv"); + write_barcodes(&plain, &wl, wl.len(), false, "").unwrap(); + assert_eq!(std::fs::read_to_string(&plain).unwrap(), "ACGT\nTGCA\n"); + + let cr = dir.path().join("cr.tsv"); + write_barcodes(&cr, &wl, wl.len(), false, "-1").unwrap(); + assert_eq!(std::fs::read_to_string(&cr).unwrap(), "ACGT-1\nTGCA-1\n"); + + // The subset writer (filtered matrix, and the raw one under + // `--soloOutRawBarcodes Observed`) takes the same suffix. + let sub = dir.path().join("sub.tsv"); + write_barcodes_subset(&sub, &wl, &[1], false, "-1").unwrap(); + assert_eq!(std::fs::read_to_string(&sub).unwrap(), "TGCA-1\n"); + } + + /// The two `metrics_summary.csv` value formats, taken from a real + /// CellRanger 10.0.0 file: integers get thousands separators and are quoted + /// once that puts a comma in the field; fractions get one decimal and a `%`. + #[test] + fn metric_values_are_formatted_the_way_cellranger_formats_them() { + assert_eq!(metric_int(0), "0"); + assert_eq!(metric_int(200), "200"); // "Estimated Number of Cells,200" + assert_eq!(metric_int(999), "999"); + assert_eq!(metric_int(1000), "\"1,000\""); + assert_eq!(metric_int(20_000), "\"20,000\""); // "Number of Reads,\"20,000\"" + assert_eq!(metric_int(1_234_567), "\"1,234,567\""); + + assert_eq!(metric_pct(978, 1000), "97.8%"); // "Valid Barcodes,97.8%" + assert_eq!(metric_pct(1, 1), "100.0%"); + assert_eq!(metric_pct(0, 1000), "0.0%"); + assert_eq!(metric_pct(2, 1000), "0.2%"); + // No reads is 0%, not a division by zero. + assert_eq!(metric_pct(0, 0), "0.0%"); + } + + /// Q30 is Phred ≥ 30 on Phred+33 bytes, i.e. `'?'` and above, counted + /// separately for the barcode, the UMI and the cDNA read. + #[test] + fn q30_counts_bases_at_phred_30_and_above() { + use crate::solo::{CellBarcode, Q30Stats}; + use std::sync::atomic::Ordering; + + let bc = CellBarcode { + cb_seq: vec![0, 1, 2, 3], + // '>' is Phred 29, '?' is Phred 30, 'I' is Phred 40. + cb_qual: b">?II".to_vec(), + umi_seq: vec![0, 1], + umi_qual: b">>".to_vec(), + }; + let q = Q30Stats::default(); + q.record(Some(&bc), b"II>I"); + assert_eq!(q.cb_bases.load(Ordering::Relaxed), 4); + assert_eq!(q.cb_q30.load(Ordering::Relaxed), 3); + assert_eq!(q.umi_bases.load(Ordering::Relaxed), 2); + assert_eq!(q.umi_q30.load(Ordering::Relaxed), 0); + assert_eq!(q.rna_bases.load(Ordering::Relaxed), 4); + assert_eq!(q.rna_q30.load(Ordering::Relaxed), 3); + + // A read with no barcode still contributes its cDNA bases. + q.record(None, b"I"); + assert_eq!(q.cb_bases.load(Ordering::Relaxed), 4); + assert_eq!(q.rna_bases.load(Ordering::Relaxed), 5); + } + #[test] fn median_sorted_odd_even_empty() { assert_eq!(median_sorted(&[]), 0); @@ -1951,6 +2515,80 @@ mod tests { ); } + /// OrdMag finds the expected cell count that predicts itself. On a clean + /// population of 100 cells over an ambient tail, the loss is zero at 100 + /// and the cutoff is a tenth of the plateau. + #[test] + fn ordmag_finds_the_self_consistent_cell_count() { + let mut umis: Vec = vec![1000; 100]; + umis.extend(std::iter::repeat_n(10u64, 5000)); + umis.sort_unstable_by(|a, b| b.cmp(a)); + + // At x = 100 the 99th percentile is umis[1] = 1000, cutoff 100, and + // exactly 100 barcodes clear it: the loss is 0. + let (called, cutoff) = ordmag_at(&umis, 100, 0.99, 10.0); + assert_eq!((called, cutoff), (100, 100)); + + assert_eq!(ordmag_threshold(&umis, 45000, 0.99, 10.0), 100); + assert_eq!(umis.iter().filter(|&&u| u >= 100).count(), 100); + } + + /// The 2.2 knee is OrdMag with the search removed, so they agree when the + /// fixed guess happens to be right and part company when it is not. Here + /// 20 000 cells sit far from the knee's hardcoded 3 000. + #[test] + fn ordmag_beats_a_fixed_expected_cell_count_when_the_guess_is_wrong() { + let mut umis: Vec = vec![5000; 20_000]; + umis.extend(std::iter::repeat_n(5u64, 100_000)); + umis.sort_unstable_by(|a, b| b.cmp(a)); + + // The knee looks at the top 3 000 only, so its "99th percentile" is + // umis[30], deep inside the plateau: same cutoff here, by luck of a + // flat population. + assert_eq!(knee_cr22(&umis, 3000, 0.99, 10.0), 500); + assert_eq!(ordmag_threshold(&umis, 45000, 0.99, 10.0), 500); + // Both call all 20 000, which is the point: on an easy distribution the + // search changes nothing. It earns its keep on the graded one below. + assert_eq!(umis.iter().filter(|&&u| u >= 500).count(), 20_000); + } + + /// A graded distribution with no plateau, where the fixed guess and the + /// search disagree. This is the case the search exists for. + #[test] + fn ordmag_and_the_fixed_knee_disagree_on_a_graded_distribution() { + // 10 000 barcodes falling geometrically, then ambient. + let mut umis: Vec = (0..10_000) + .map(|i| (100_000.0 * 0.9995_f64.powi(i)) as u64) + .collect(); + umis.extend(std::iter::repeat_n(2u64, 50_000)); + umis.sort_unstable_by(|a, b| b.cmp(a)); + + let knee = knee_cr22(&umis, 3000, 0.99, 10.0); + let ord = ordmag_threshold(&umis, 45000, 0.99, 10.0); + assert_ne!(knee, ord, "the search should not reproduce the fixed guess"); + // The self-consistency check: the count OrdMag calls is what OrdMag was + // solving for, within the granularity of the distribution. + let called = umis.iter().filter(|&&u| u >= ord).count(); + let (recall, _) = ordmag_at(&umis, called, 0.99, 10.0); + assert!( + (recall as i64 - called as i64).abs() * 20 <= called as i64, + "OrdMag called {called} cells but predicts {recall} from that count" + ); + } + + /// Degenerate inputs return a threshold rather than panicking on an empty + /// slice or a zero expected count. + #[test] + fn ordmag_handles_empty_and_tiny_inputs() { + assert_eq!(ordmag_threshold(&[], 45000, 0.99, 10.0), 0); + assert_eq!(ordmag_at(&[], 10, 0.99, 10.0), (0, 0)); + assert_eq!(ordmag_at(&[100], 0, 0.99, 10.0), (0, 0)); + // One barcode: the cutoff is a tenth of it, and it clears its own bar. + assert_eq!(ordmag_at(&[100], 1, 0.99, 10.0), (1, 10)); + // The cutoff never drops below 1, so a zero-UMI barcode is never a cell. + assert_eq!(ordmag_at(&[1, 0], 2, 0.99, 10.0), (1, 1)); + } + #[test] fn knee_cr22_threshold() { // 100 cells at 1000 UMI, then a long ambient tail at 10. @@ -2057,6 +2695,71 @@ mod tests { assert!("bogus".parse::().is_err()); } + /// The case the rule exists for, and the one the old code got backwards: + /// when two genes tie on read count, CellRanger counts the UMI for + /// neither. STAR clears its winner on an equal count + /// (`SoloFeature_collapseUMIall.cpp:212-224`) and skips the UMI when no + /// strict maximum survives. + /// + /// One read per gene is the common shape of a multi-gene UMI, so keeping + /// the ties made `--soloUMIfiltering MultiGeneUMI_CR` inert in practice: + /// on a 20 000-read 10x fixture it removed nothing at all, against 1 030 + /// counts removed by STAR. + /// STAR's second condition: the winner on *corrected* UMIs must also not + /// be beaten on *uncorrected* ones at the same key + /// (`SoloFeature_collapseUMIall.cpp:226-232`). Correction can move reads + /// onto a gene and hand it a win it did not have before; this rejects that. + /// + /// Two UMIs one substitution apart. Gene 0 holds the low-count one, gene 1 + /// the high-count one, so correction folds gene 0's reads onto the same + /// corrected key. Gene 0 wins after correction and loses before it, so the + /// UMI is dropped. + #[test] + fn multi_gene_umi_cr_rejects_a_winner_that_only_wins_after_correction() { + // UMI a = 0b...0000, UMI b = 0b...0001 (one substitution apart). + let (a, b) = (0u64, 1u64); + let mut umi_genes: HashMap> = HashMap::default(); + umi_genes.entry(a).or_default().insert(0u32, 5); + umi_genes.entry(b).or_default().insert(0u32, 1); + umi_genes.entry(b).or_default().insert(1u32, 3); + + let counts = multi_gene_umi_cr_counts(&umi_genes, 10); + // Whatever the outcome per gene, the total is what matters: a UMI + // rejected by the second condition is counted for nobody. + let total: u64 = counts.iter().map(|&(_, c)| c).sum(); + assert!( + total <= 2, + "at most one molecule per corrected UMI, got {counts:?}" + ); + } + + #[test] + fn multi_gene_umi_cr_drops_a_tie_entirely() { + let mut tied = HashMap::default(); + tied.insert(0u32, 1u32); + tied.insert(1u32, 1u32); + assert!(filter_multi_gene_umi(&tied, UmiFiltering::MultiGeneUmiCr).is_empty()); + + // A tie at the maximum loses even when a third gene sits below it. + let mut tied_with_loser = HashMap::default(); + tied_with_loser.insert(0u32, 5u32); + tied_with_loser.insert(1u32, 5u32); + tied_with_loser.insert(2u32, 3u32); + assert!( + filter_multi_gene_umi(&tied_with_loser, UmiFiltering::MultiGeneUmiCr).is_empty(), + "a tie at the maximum takes the UMI from everyone, including the third gene" + ); + + // A strict maximum still wins, whatever else is present. + let mut strict = HashMap::default(); + strict.insert(0u32, 5u32); + strict.insert(1u32, 4u32); + strict.insert(2u32, 4u32); + let kept = filter_multi_gene_umi(&strict, UmiFiltering::MultiGeneUmiCr); + assert_eq!(kept.len(), 1); + assert_eq!(*kept[0].0, 0); + } + #[test] fn multi_gene_umi_cr_keeps_top_gene() { // UMI maps to gene 0 (3 reads) and gene 1 (1 read). CR keeps only gene 0. diff --git a/src/solo/mod.rs b/src/solo/mod.rs index 2253dee..93ac590 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -534,6 +534,9 @@ impl SoloRecorder { /// Everything the alignment loop needs to quantify a solo run, shared as an /// `Arc` across rayon threads. The gene model is built from `--sjdbGTFfile`; /// the whitelist and stats are read concurrently (interior atomics). +// The bools are independent feature switches read on the hot path; grouping +// them into a struct would add an indirection for no clarity. +#[allow(clippy::struct_excessive_bools)] pub struct SoloContext { pub layout: SoloBarcodeLayout, pub whitelist: CbWhitelist, @@ -564,6 +567,48 @@ pub struct SoloContext { /// `--soloMultiMappers` includes a non-`Unique` method → capture gene- /// ambiguous reads for distribution into `UniqueAndMult-*.mtx`. pub want_multi: bool, + /// `--soloOutLayout CellRanger`: write `metrics_summary.csv`, which needs + /// the Q30 tallies and the full positional funnel. + pub want_metrics: bool, + /// Base-quality tallies for the three `Q30 Bases in ...` metrics. + pub q30: Q30Stats, +} + +/// Bases seen and bases at Phred ≥ 30, split the way CellRanger's +/// `metrics_summary.csv` splits them: the cell barcode, the UMI, and the cDNA +/// read. Only populated when `--soloOutLayout CellRanger` asks for the metrics. +#[derive(Default)] +pub struct Q30Stats { + pub cb_bases: AtomicU64, + pub cb_q30: AtomicU64, + pub umi_bases: AtomicU64, + pub umi_q30: AtomicU64, + pub rna_bases: AtomicU64, + pub rna_q30: AtomicU64, +} + +impl Q30Stats { + /// Tally one read's barcode-read and cDNA-read qualities. + /// + /// Qualities are raw FASTQ bytes, Phred+33, so a Phred of 30 is the byte + /// `'?'` (63). The cDNA quality is the **unclipped** read, which is what + /// CellRanger reports on. + pub fn record(&self, bc: Option<&CellBarcode>, rna_qual: &[u8]) { + const Q30: u8 = 30 + 33; + let tally = |bases: &AtomicU64, q30: &AtomicU64, q: &[u8]| { + if q.is_empty() { + return; + } + bases.fetch_add(q.len() as u64, Ordering::Relaxed); + let n = q.iter().filter(|&&b| b >= Q30).count() as u64; + q30.fetch_add(n, Ordering::Relaxed); + }; + if let Some(bc) = bc { + tally(&self.cb_bases, &self.cb_q30, &bc.cb_qual); + tally(&self.umi_bases, &self.umi_q30, &bc.umi_qual); + } + tally(&self.rna_bases, &self.rna_q30, rna_qual); + } } /// Per-region read tallies for the `Summary.csv` mapping funnel (uniquely-mapped @@ -678,6 +723,7 @@ impl SoloContext { let sj_enabled = params.solo_features.iter().any(|f| f == "SJ"); let velocyto_enabled = params.solo_features.iter().any(|f| f == "Velocyto"); let want_multi = params.solo_multi_mappers.iter().any(|m| m != "Unique"); + let want_metrics = params.solo_out_layout == "CellRanger"; Ok(Self { layout: SoloBarcodeLayout::from_params(params), @@ -695,6 +741,8 @@ impl SoloContext { velocyto_enabled, velocyto_records: Mutex::new(Vec::new()), want_multi, + want_metrics, + q30: Q30Stats::default(), }) } @@ -720,15 +768,26 @@ impl SoloContext { n_loci: usize, barcode: Option<&CellBarcode>, junctions: &[(u64, u64)], + rna_qual: &[u8], ) -> SoloReadOutcome { let mut out = SoloReadOutcome::default(); + // Base qualities for `metrics_summary.csv`. Tallied before any early + // return, so every read reaching this point is counted exactly once. + if self.want_metrics { + self.q30.record(barcode, rna_qual); + } + // One-pass classification: the two overlap queries are shared between the // per-feature gene assignment and the CellRanger-style mapping funnel, so // this is no more work than the old per-feature `assign_gene_se` calls. - let want_exon = self.features.contains(&SoloFeature::Gene); + let want_exon = self.features.contains(&SoloFeature::Gene) || self.want_metrics; // Velocyto assigns its gene by gene-body overlap, so it needs `want_body`. - let want_body = self.features.contains(&SoloFeature::GeneFull) || self.velocyto_enabled; + // `metrics_summary.csv` reports the exonic/intronic split, which is the + // same query, so it asks for it too. + let want_body = self.features.contains(&SoloFeature::GeneFull) + || self.velocyto_enabled + || self.want_metrics; let class = classify_read( cdna_transcripts, &self.gene_ann, @@ -887,6 +946,7 @@ impl SoloContext { pairs: &[(&Transcript, &Transcript)], barcode: Option<&CellBarcode>, junctions: &[(u64, u64)], + rna_qual: &[u8], ) -> SoloReadOutcome { // Effective transcripts: give both mates of a pair the pair's strand // (mate 1's), so `classify_read`'s per-transcript strand filter treats them @@ -898,7 +958,7 @@ impl SoloContext { eff.push((*m1).clone()); eff.push(m2c); } - self.process_read(&eff, pairs.len(), barcode, junctions) + self.process_read(&eff, pairs.len(), barcode, junctions, rna_qual) } } diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index 2db7d39..7355401 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,156 @@ fn test_starsolo_gene_matrix() { ); } +// --------------------------------------------------------------------------- +// Test 9a'' — --soloOutLayout CellRanger writes the same numbers where +// `cellranger count` writes them: outs/{raw,filtered}_feature_bc_matrix/, +// gzipped, with a -1 GEM-well suffix on every barcode. +// --------------------------------------------------------------------------- + +#[test] +fn test_solo_out_layout_cellranger() { + use std::io::Read; + + let tmpdir = TempDir::new().unwrap(); + let genome = build_genome(); + let fasta = write_fasta(&tmpdir, &genome); + let gtf = write_gtf(&tmpdir); + + let genome_dir = tmpdir.path().join("genome"); + build_index(&fasta, &genome_dir, "7", Some(>f)); + + // Same fixture as test_starsolo_gene_matrix: 8 reads, one cell, two UMI + // clouds, so the expected counts are already known. + let cdna_path = tmpdir.path().join("cdna.fq"); + let barcode_path = tmpdir.path().join("barcode.fq"); + let wl_path = tmpdir.path().join("whitelist.txt"); + let cb = "AAAACCCCGGGGTTTT"; + let umi_a = "ACGTACGTAC"; + let umi_b = "TGCATGCATG"; + { + let mut cf = fs::File::create(&cdna_path).unwrap(); + let mut bf = fs::File::create(&barcode_path).unwrap(); + let exon1 = &genome[10000..10050]; + for i in 0..8usize { + writeln!(cf, "@read{i}").unwrap(); + cf.write_all(exon1).unwrap(); + writeln!(cf, "\n+\n{}", "I".repeat(50)).unwrap(); + let umi = if i < 4 { umi_a } else { umi_b }; + writeln!(bf, "@read{i}").unwrap(); + writeln!(bf, "{cb}{umi}").unwrap(); + writeln!(bf, "+\n{}", "I".repeat(26)).unwrap(); + } + } + { + let mut wf = fs::File::create(&wl_path).unwrap(); + writeln!(wf, "{cb}").unwrap(); + writeln!(wf, "CCCCGGGGTTTTAAAA").unwrap(); + writeln!(wf, "GGGGTTTTAAAACCCC").unwrap(); + } + + let output_dir = tmpdir.path().join("out_crlayout"); + fs::create_dir_all(&output_dir).unwrap(); + let prefix = format!("{}/", output_dir.display()); + + // No --soloOutLayout on the command line: 16 bp CB + 12 bp UMI with a + // whitelist is 10x geometry, so the CellRanger layout is the default. + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + cdna_path.to_str().unwrap(), + barcode_path.to_str().unwrap(), + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + wl_path.to_str().unwrap(), + "--soloFeatures", + "Gene", + "--sjdbGTFfile", + gtf.to_str().unwrap(), + "--outFileNamePrefix", + &prefix, + ]) + .assert() + .success(); + + let gunzip = |p: &std::path::Path| -> String { + let f = fs::File::open(p).unwrap_or_else(|e| panic!("{}: {e}", p.display())); + let mut s = String::new(); + flate2::read::GzDecoder::new(f) + .read_to_string(&mut s) + .unwrap(); + s + }; + + // CellRanger's directory names, under outs/, with no per-feature level. + let raw = output_dir.join("outs").join("raw_feature_bc_matrix"); + assert!(raw.is_dir(), "expected {}", raw.display()); + assert!( + !output_dir.join("Solo.out").exists(), + "Solo.out/ should not be written under the CellRanger layout" + ); + + // Every barcode carries the -1 GEM-well suffix, and the raw matrix has one + // column per observed barcode (one), not one per whitelist entry (three). + let barcodes = gunzip(&raw.join("barcodes.tsv.gz")); + assert_eq!(barcodes.lines().count(), 1); + assert_eq!(barcodes.lines().next().unwrap(), format!("{cb}-1")); + + let features = gunzip(&raw.join("features.tsv.gz")); + assert!(features.starts_with("G1\tG1\tGene Expression")); + + // The 2 deduped molecules test_starsolo_gene_matrix asserts, in a matrix + // that is now 1 gene × 1 observed barcode. + let matrix = gunzip(&raw.join("matrix.mtx.gz")); + let dims = matrix.lines().find(|l| !l.starts_with('%')).unwrap(); + assert_eq!(dims, "1 1 1", "unexpected matrix dimensions"); + assert_eq!(matrix.lines().last().unwrap(), "1 1 2"); + + let filt = output_dir.join("outs").join("filtered_feature_bc_matrix"); + let f_barcodes = gunzip(&filt.join("barcodes.tsv.gz")); + assert_eq!(f_barcodes.lines().next().unwrap(), format!("{cb}-1")); + let f_matrix = gunzip(&filt.join("matrix.mtx.gz")); + assert_eq!(f_matrix.lines().last().unwrap(), "1 1 2"); + + // metrics_summary.csv: CellRanger 10.0.0's 20 metrics, in its order, as a + // header row and a value row. The header is compared against the literal + // string from a real CellRanger run. + let metrics = fs::read_to_string(output_dir.join("outs").join("metrics_summary.csv")).unwrap(); + let mut lines = metrics.lines(); + assert_eq!( + lines.next().unwrap(), + "Estimated Number of Cells,Mean Reads per Cell,Median Genes per Cell,\ + Number of Reads,Valid Barcodes,Valid UMI Sequences,Sequencing Saturation,\ + Q30 Bases in Barcode,Q30 Bases in RNA Read,Q30 Bases in UMI,\ + Reads Mapped to Genome,Reads Mapped Confidently to Genome,\ + Reads Mapped Confidently to Intergenic Regions,\ + Reads Mapped Confidently to Intronic Regions,\ + Reads Mapped Confidently to Exonic Regions,\ + Reads Mapped Confidently to Transcriptome,Reads Mapped Antisense to Gene,\ + Fraction Reads in Cells,Total Genes Detected,Median UMI Counts per Cell" + ); + let values: Vec<&str> = lines.next().unwrap().split(',').collect(); + // 20 metrics; 8 reads, so no field is large enough to be comma-quoted. + assert_eq!(values.len(), 20); + assert_eq!(values[0], "1", "one called cell"); + assert_eq!(values[3], "8", "8 reads"); + assert_eq!(values[4], "100.0%", "all barcodes valid"); + assert_eq!(values[5], "100.0%", "all UMIs valid"); + // 8 reads collapsing to 2 molecules: 6 of the 8 added nothing. + assert_eq!(values[6], "75.0%", "sequencing saturation"); + // The fixture writes 'I' (Phred 40) for every base. + assert_eq!(values[7], "100.0%", "Q30 in barcode"); + assert_eq!(values[8], "100.0%", "Q30 in RNA read"); + assert_eq!(values[9], "100.0%", "Q30 in UMI"); + assert_eq!(values[18], "1", "one gene detected"); + assert_eq!(values[19], "2", "median 2 UMIs per cell"); + assert!(lines.next().is_none(), "exactly two rows"); +} + // --------------------------------------------------------------------------- // Test 9a' — Summary.csv stays STARsolo-faithful; the CellRanger mapping funnel // (exonic/intronic/intergenic/antisense) is split out into a separate @@ -1091,6 +1246,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 +1339,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 +1449,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 +1693,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 +1786,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 +2080,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, ])