diff --git a/CHANGELOG.md b/CHANGELOG.md index 14b4346..de5e978 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 66f4088..59c6eca 100644 --- a/DIVERGENCE.md +++ b/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. --- diff --git a/src/params/mod.rs b/src/params/mod.rs index 3b50731..b8ce61e 100644 --- a/src/params/mod.rs +++ b/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 = ::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(&[ diff --git a/src/solo/count.rs b/src/solo/count.rs index 4ddb414..6704a1b 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -506,6 +506,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 @@ -1392,20 +1413,45 @@ pub fn write_gene_matrix( &ctx.gene_ann.gene_names, gzip, )?; - write_barcodes( - &raw_dir.join(&barcodes_name), - &ctx.whitelist, - sorted.len(), - gzip, - )?; + // `--soloOutRawBarcodes Observed` narrows the raw matrix to the + // barcodes that actually carry a count, which is what CellRanger's + // `raw_feature_bc_matrix` holds. STARsolo's raw matrix has a column per + // whitelist barcode, so the default keeps that. + let observed: Option> = if params.solo_out_raw_barcodes == "Observed" { + Some(observed_barcodes(&body)?) + } else { + None + }; + let (raw_cols, raw_remap) = match &observed { + Some(cbs) => { + let map: HashMap = cbs + .iter() + .enumerate() + .map(|(col, &cb)| (cb, col as u32 + 1)) + .collect(); + (cbs.len(), Some(map)) + } + None => (sorted.len(), None), + }; + match &observed { + Some(cbs) => { + write_barcodes_subset(&raw_dir.join(&barcodes_name), &ctx.whitelist, cbs, gzip)?; + } + None => write_barcodes( + &raw_dir.join(&barcodes_name), + &ctx.whitelist, + sorted.len(), + gzip, + )?, + } finalize_matrix( &body, &raw_dir.join(&matrix_name), gzip, n_genes, - sorted.len(), + raw_cols, mstats.nnz, - None, + raw_remap.as_ref(), )?; log::info!( "STARsolo: wrote {}/raw matrix ({} genes × {} barcodes, {} entries){}",