feat(solo): CellRanger behaviour by default on 10x geometry (changes default output, stacked on #175)

CHANGELOG.md

@@ -44,6 +44,17 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 - Read names are cut at `--readNameSeparator` (default `/`), as STAR does. A
   read named `foo/1` was previously emitted as `foo/1` where STAR emits `foo`.
+- 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

DIVERGENCE.md

@@ -110,6 +110,61 @@ On the 10k yeast PE benchmark, 4 reads differ in alignment score (AS) because ST
 **Impact.** Past libc++'s load factor the map rehashes, and the order then depends on the bucket count, which depends on how many distinct barcodes were seen; beyond that size the order diverges. The **values never do** — only which line they appear on. Reading the file by barcode rather than by position is unaffected either way.
 
 **Source.** `src/solo/cell_reads.rs`, locked by `rows_are_emitted_in_reverse_first_appearance_order`. STAR: `SoloFeature_statsOutput.cpp`.
+### 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.
 
 ---
 

src/params/mod.rs

@@ -1290,6 +1290,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
@@ -1515,6 +1530,8 @@ impl Parameters {
         let matches = command.clone().get_matches_from(args.iter());
         let mut params = <Self as clap::FromArgMatches>::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()
@@ -2174,8 +2191,225 @@ 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;
+        }
+        // MultiGeneUMI_CR decides ownership from the corrected-UMI map that
+        // only the CellRanger dedup builds, and STAR refuses the pair
+        // otherwise. So when the user has picked a different dedup, this
+        // default stays out of the way rather than composing into a
+        // combination the validation then rejects.
+        if id == "solo_umi_filtering"
+            && given("solo_umi_dedup")
+            && params.solo_umi_dedup.first().map(String::as_str) != Some("1MM_CR")
+        {
+            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).
@@ -2893,6 +3127,15 @@ mod tests {
             "genes.gtf",
             "--soloCBwhitelist",
             "wl.txt",
+            // Deliberately not 10x geometry: on a 10x run this build defaults
+            // the dedup to 1MM_CR, which would satisfy the rule on its own and
+            // hide what this test is about (see the 10x case at the end).
+            "--soloCBlen",
+            "12",
+            "--soloUMIlen",
+            "8",
+            "--soloUMIstart",
+            "13",
             "--soloUMIfiltering",
             "MultiGeneUMI_CR",
         ];
@@ -2914,6 +3157,24 @@ mod tests {
         multi.extend_from_slice(&["--soloUMIdedup", "1MM_CR", "Exact"]);
         assert!(try_parse(&multi).is_err());
 
+        // On 10x geometry the CellRanger defaults supply 1MM_CR themselves, so
+        // the same flags are accepted rather than refused.
+        let tenx = [
+            "--readFilesIn",
+            "cdna.fq",
+            "bc.fq",
+            "--soloType",
+            "CB_UMI_Simple",
+            "--sjdbGTFfile",
+            "genes.gtf",
+            "--soloCBwhitelist",
+            "wl.txt",
+            "--soloUMIfiltering",
+            "MultiGeneUMI_CR",
+        ];
+        let p = try_parse(&tenx).expect("10x defaults supply the CellRanger dedup");
+        assert_eq!(p.solo_umi_dedup, vec!["1MM_CR".to_string()]);
+
         // The pairing rule applies only to MultiGeneUMI_CR.
         assert!(
             try_parse(&[