From 340d2570cedd2690f31bb2e539d40c5b575e1373 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 28 Jul 2026 23:58:06 +0200 Subject: [PATCH 1/9] fix(solo): apply STAR's cbMinP posterior threshold and the oneExact guard MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `resolve_multi_cb` computed the count+quality posterior correctly, computed the normalising total, and then ignored it: any candidate with positive weight won by argmax. Two of STAR's conditions were missing. `cbMinP`. STAR requires the winner to hold at least 97.5% of the total posterior mass. Without it a barcode is assigned whenever any candidate has positive weight, including when two candidates are nearly tied — which is precisely the case the threshold exists to reject. A read whose barcode cannot be singled out should be discarded, not guessed. `oneExact`. Unless pseudocounts are in use, the winner must have been observed as an exact match at least once. A candidate whose only support is a pseudocount is eligible, not evidence. **This changes default counting behaviour**, and the existing unit test had to change with it: it asserted that priors of 10 against 3 resolve to the higher one, which is a 77% posterior share and now correctly resolves to nothing. The test was encoding the missing threshold. It now covers both sides — a decisive 1000-against-3 prior is accepted, a tie under pseudocounts is not. Verification I could not do here, and would want before this merges: the solo differential against real STARsolo (`test/solo_diff_docker.sh`). Reads that previously landed on an ambiguous barcode now land nowhere, so raw matrix counts will move, and the direction should be confirmed against the oracle rather than argued from the source. Everything in-tree is green — 560 lib + 26 integration — but that is not the same thing. Co-Authored-By: Claude Opus 5 (1M context) --- src/solo/count.rs | 62 +++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 52 insertions(+), 10 deletions(-) diff --git a/src/solo/count.rs b/src/solo/count.rs index 4ddb414..4f89e2c 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -279,11 +279,26 @@ fn directional(umis: &HashMap, umi_len: usize, dir_count_add: i64) -> // Cell-barcode multi-match resolution (deferred 1MM_multi) // --------------------------------------------------------------------------- +/// STAR's posterior threshold for accepting a corrected cell barcode +/// (`SoloReadBarcode` `cbMinP`). The winner must hold at least this share of +/// the total posterior mass, otherwise the read is discarded rather than +/// assigned to a barcode the evidence does not actually single out. +const CB_MIN_P: f64 = 0.975; + /// Resolve a 1MM_multi cell barcode to a single whitelist index using the /// count+quality posterior: weight = `(exactCount[cand] + pseudocount) · 10^(−q/10)` /// where `q` is the mismatch-position Phred score. `pseudocount` is 1 for the -/// `*_pseudocounts` match types (CellRanger ≥ 3.0). Returns the argmax, or -/// `None` if no candidate has positive weight. +/// `*_pseudocounts` match types (CellRanger ≥ 3.0). +/// +/// Returns the argmax, subject to two conditions STAR imposes and this did not: +/// +/// - the winner's share of the total posterior must reach [`CB_MIN_P`]. Taking +/// the argmax unconditionally assigns a barcode whenever any candidate has +/// positive weight, including when two candidates are all but tied, which is +/// exactly the case the threshold exists to reject; +/// - unless pseudocounts are in use, the winner must have been observed as an +/// exact match at least once (`oneExact`). A candidate whose only support is +/// a pseudocount is not evidence. fn resolve_multi_cb( candidates: &[crate::solo::whitelist::CbCandidate], exact_counts: &[u64], @@ -301,10 +316,20 @@ fn resolve_multi_cb( _ => best = Some((c.wl_index, weight)), } } - match best { - Some((idx, w)) if total > 0.0 && w > 0.0 => Some(idx), - _ => None, + let (idx, w) = best?; + if total <= 0.0 || w <= 0.0 { + return None; + } + // Posterior threshold: the winner has to actually win. + if w / total < CB_MIN_P { + return None; + } + // `oneExact`: required by every match type except the pseudocount ones, + // where the pseudocount is deliberately standing in for absent evidence. + if pseudocount == 0.0 && exact_counts.get(idx as usize).copied().unwrap_or(0) == 0 { + return None; } + Some(idx) } // --------------------------------------------------------------------------- @@ -2601,13 +2626,30 @@ mod tests { mismatch_qual: b'I', }, ]; - // Same quality → higher exact-count prior wins. - assert_eq!(resolve_multi_cb(&cands, &[10, 3], 0.0), Some(0)); - assert_eq!(resolve_multi_cb(&cands, &[3, 10], 0.0), Some(1)); + // Same quality, so the exact-count priors decide. The winner must also + // clear STAR's posterior threshold: 10 against 3 is only a 77% share, + // which is exactly the ambiguous case cbMinP exists to reject. + assert_eq!(resolve_multi_cb(&cands, &[10, 3], 0.0), None); + assert_eq!(resolve_multi_cb(&cands, &[3, 10], 0.0), None); + + // A decisive prior does clear it (1000/1003 = 99.7%). + assert_eq!(resolve_multi_cb(&cands, &[1000, 3], 0.0), Some(0)); + assert_eq!(resolve_multi_cb(&cands, &[3, 1000], 0.0), Some(1)); + // No prior signal and no pseudocount → rejected. assert_eq!(resolve_multi_cb(&cands, &[0, 0], 0.0), None); - // Pseudocount gives every candidate positive weight → argmax accepted. - assert!(resolve_multi_cb(&cands, &[0, 0], 1.0).is_some()); + + // Pseudocounts alone leave the two candidates tied at 50%, far below + // the threshold, so they are still rejected: a pseudocount makes a + // candidate *eligible*, it does not make it *evidence*. + assert_eq!(resolve_multi_cb(&cands, &[0, 0], 1.0), None); + + // With pseudocounts and a decisive prior, accepted. + assert_eq!(resolve_multi_cb(&cands, &[1000, 0], 1.0), Some(0)); + + // oneExact: without pseudocounts a winner that was never seen exactly + // is refused even when it dominates the posterior. + assert_eq!(resolve_multi_cb(&cands[..1], &[0, 0], 0.0), None); } #[test] From 70cda91f28d3a8b17111248a563ad1d221b1a78f Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 00:10:52 +0200 Subject: [PATCH 2/9] feat(solo): --soloOutFormatFeaturesGeneField3, and lock the homopolymer-UMI rule MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `--soloOutFormatFeaturesGeneField3` sets the third column of `features.tsv`, STAR's default being `Gene Expression`. The sentinel `-` drops the column, giving the two-column form some downstream tools expect. Previously the string was hard-coded. Also adds `umi_valid_rejects_every_homopolymer`. No code change was needed — `is_homopolymer` already rejects a UMI of any single repeated base — but the rule deserves a test that says why it is deliberate. STAR has a quirk (`umiL = 0`) under which a non-A homopolymer UMI slips through its filter, so poly-C, poly-G and poly-T survive there and are counted. That is a defect rather than a rule: a UMI of one repeated base carries no information whichever base it happens to be. This diverges from 2.7.11b knowingly, and the test now records that instead of leaving it to be rediscovered as a difference in a differential run. Co-Authored-By: Claude Opus 5 (1M context) --- src/params/mod.rs | 8 ++++++++ src/solo/count.rs | 19 ++++++++++++++++--- src/solo/whitelist.rs | 27 +++++++++++++++++++++++++++ 3 files changed, 51 insertions(+), 3 deletions(-) diff --git a/src/params/mod.rs b/src/params/mod.rs index 3b50731..220d28d 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1280,6 +1280,14 @@ pub struct Parameters { #[arg(long = "soloMultiMappers", num_args = 1.., default_values_t = vec!["Unique".to_string()])] pub solo_multi_mappers: Vec, + /// Third column of `features.tsv`. STAR's default is `Gene Expression`; + /// the sentinel `-` suppresses the column entirely. + #[arg( + long = "soloOutFormatFeaturesGeneField3", + default_value = "Gene Expression" + )] + pub solo_out_format_features_gene_field3: String, + /// Output directory name for solo matrices (relative to `--outFileNamePrefix`). #[arg(long = "soloOutFileNames", num_args = 1.., default_values_t = vec!["Solo.out/".to_string(), "features.tsv".to_string(), "barcodes.tsv".to_string(), "matrix.mtx".to_string()])] pub solo_out_file_names: Vec, diff --git a/src/solo/count.rs b/src/solo/count.rs index 4f89e2c..3f53103 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -1415,6 +1415,7 @@ pub fn write_gene_matrix( &raw_dir.join(&features_name), &ctx.gene_ann.gene_ids, &ctx.gene_ann.gene_names, + ¶ms.solo_out_format_features_gene_field3, gzip, )?; write_barcodes( @@ -1492,6 +1493,7 @@ pub fn write_gene_matrix( &filt_dir.join(&features_name), &ctx.gene_ann.gene_ids, &ctx.gene_ann.gene_names, + ¶ms.solo_out_format_features_gene_field3, gzip, )?; write_barcodes_subset(&filt_dir.join(&barcodes_name), &ctx.whitelist, &cbs, gzip)?; @@ -1651,6 +1653,7 @@ pub fn write_gene_matrix( &velo_dir.join(&features_name), &ctx.gene_ann.gene_ids, &ctx.gene_ann.gene_names, + ¶ms.solo_out_format_features_gene_field3, gzip, )?; write_barcodes( @@ -2091,12 +2094,17 @@ fn write_cellranger_summary( Ok(()) } -/// `features.tsv`: `gene_id gene_name "Gene Expression"` (CellRanger -/// v3 layout). We have no gene names, so the id is repeated. +/// `features.tsv`: `gene_id gene_name ` (CellRanger v3 +/// layout). +/// +/// `field3` is `--soloOutFormatFeaturesGeneField3`, whose STAR default is +/// `Gene Expression`. The sentinel `-` drops the column, giving the two-column +/// form some downstream tools expect. fn write_features( path: &Path, gene_ids: &[String], gene_names: &[String], + field3: &str, gzip: bool, ) -> Result<(), Error> { // CellRanger v3 layout: gene_id gene_name "Gene Expression". @@ -2105,7 +2113,12 @@ fn write_features( // fallback is already baked into `gene_names`. write_file(path, gzip, |w| { for (id, name) in gene_ids.iter().zip(gene_names.iter()) { - writeln!(w, "{id}\t{name}\tGene Expression").map_err(|e| Error::io(e, path))?; + if field3 == "-" { + writeln!(w, "{id}\t{name}") + } else { + writeln!(w, "{id}\t{name}\t{field3}") + } + .map_err(|e| Error::io(e, path))?; } Ok(()) })?; diff --git a/src/solo/whitelist.rs b/src/solo/whitelist.rs index 1d882ae..15b2a8a 100644 --- a/src/solo/whitelist.rs +++ b/src/solo/whitelist.rs @@ -751,4 +751,31 @@ mod tests { assert!(n.mm1_multi_nbase && n.pseudocounts); assert!(CbMatchType::from_str("bogus").is_err()); } + + #[test] + fn umi_valid_rejects_every_homopolymer() { + // STAR has a quirk (`umiL = 0`) under which a non-A homopolymer UMI + // slips through its filter. That is a defect, not a rule: a UMI of one + // repeated base carries no information whichever base it is. All four + // are rejected here, deliberately diverging from 2.7.11b. + // Sequences are base codes: 0..=3 for ACGT, 4 for N. + for (code, base) in [(0u8, 'A'), (1, 'C'), (2, 'G'), (3, 'T')] { + let umi = vec![code; 10]; + assert!( + matches!(check_umi(&umi), UmiCheck::Homopolymer), + "poly-{base} should be rejected" + ); + } + + // A UMI that merely starts with a run is fine. + assert!(matches!( + check_umi(&[0, 0, 0, 0, 0, 0, 0, 0, 0, 1]), + UmiCheck::Ok(_) + )); + // N anywhere still fails on its own grounds. + assert!(matches!( + check_umi(&[0, 0, 0, 0, 4, 0, 0, 0, 0, 0]), + UmiCheck::NinUmi + )); + } } From 613d4d3dd6add039d28aeca6454483eee0ec5b1e Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 00:18:37 +0200 Subject: [PATCH 3/9] feat(solo): --soloChemistry geometry presets, and refuse --soloCBtype String MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `--soloChemistry` sets the CB/UMI geometry from a 10x chemistry name so it does not have to be spelled out four flags at a time: SC3Pv1 (14-base CB, 10-base UMI), SC3Pv2 (16/10), SC3Pv3, SC3Pv4 and SC5P (16/12), plus the CellRanger aliases. Unknown names are refused. The preset wins over an explicit `--soloCBstart` and friends, and warns when it overrode one. Half-applying a preset — taking its CB length but the flag's UMI start, say — would produce a geometry neither the user nor the preset asked for, and would fail silently in the counts rather than loudly at startup. `--soloCBtype` is declared and accepts `Sequence`. `String` is refused rather than accepted and ignored: barcodes here are 2-bit packed ACGT, and a non-ACGT barcode packed that way is nonsense rather than an approximation. Supporting it means a different whitelist representation throughout, which is its own change. Co-Authored-By: Claude Opus 5 (1M context) --- src/params/mod.rs | 151 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 151 insertions(+) diff --git a/src/params/mod.rs b/src/params/mod.rs index 220d28d..175ad07 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1196,6 +1196,17 @@ pub struct Parameters { #[arg(long = "soloCBwhitelist", num_args = 1.., default_values_t = vec!["None".to_string()])] pub solo_cb_whitelist: Vec, + /// 10x chemistry preset. Sets the CB/UMI geometry so it does not have to be + /// spelled out with `--soloCBstart` and friends. `-` (the default) leaves + /// the geometry to those flags. + #[arg(long = "soloChemistry", default_value = "-")] + pub solo_chemistry: String, + + /// How cell barcodes are represented: `Sequence` (2-bit packed ACGT, the + /// default). + #[arg(long = "soloCBtype", default_value = "Sequence")] + pub solo_cb_type: String, + /// 1-based start position of the cell barcode in the barcode read. #[arg(long = "soloCBstart", default_value_t = 1)] pub solo_cb_start: u32, @@ -1826,6 +1837,59 @@ impl Parameters { } } + // --soloChemistry presets. STAR applies these before validating the + // geometry, so a preset and an explicit --soloCBstart cannot disagree: + // the preset wins and says so, rather than half-applying. + if params.solo_chemistry != "-" { + let geometry = match params.solo_chemistry.as_str() { + // (cb_start, cb_len, umi_start, umi_len), 1-based as STAR takes them. + "CR_2" | "CR_3" | "CR_3.1" | "CR_4" | "SC3Pv1" => Some((1, 14, 15, 10)), + "SC3Pv2" => Some((1, 16, 17, 10)), + "SC3Pv3" | "SC3Pv4" | "SC5P" => Some((1, 16, 17, 12)), + _ => None, + }; + let Some((cb_start, cb_len, umi_start, umi_len)) = geometry else { + return Err(command.error( + ErrorKind::InvalidValue, + format!( + "unknown --soloChemistry '{}'; expected SC3Pv1, SC3Pv2, SC3Pv3, \ + SC3Pv4, SC5P, or -", + params.solo_chemistry + ), + )); + }; + for (flag, given) in [ + ("--soloCBstart", params.solo_cb_start), + ("--soloCBlen", params.solo_cb_len), + ("--soloUMIstart", params.solo_umi_start), + ("--soloUMIlen", params.solo_umi_len), + ] { + if given != 0 { + log::warn!( + "--soloChemistry {} overrides the geometry given by {flag}", + params.solo_chemistry + ); + } + } + params.solo_cb_start = cb_start; + params.solo_cb_len = cb_len; + params.solo_umi_start = umi_start; + params.solo_umi_len = umi_len; + } + + // --soloCBtype: only 2-bit packed ACGT barcodes are supported. `String` + // needs a different whitelist representation entirely, so refuse rather + // than silently pack a non-ACGT barcode into nonsense. + if params.solo_cb_type != "Sequence" { + return Err(command.error( + ErrorKind::InvalidValue, + format!( + "--soloCBtype {} is not supported; expected Sequence", + params.solo_cb_type + ), + )); + } + // ── STARsolo validation ───────────────────────────────────────── if params.run_mode() == RunMode::AlignReads && params.solo_enabled() { // CB_UMI_Complex needs one CB position + whitelist per segment. @@ -3006,4 +3070,91 @@ mod tests { assert_eq!(p.out_sam_strand_field, "None"); assert!(!p.out_sam_attributes.contains(SamAttributes::XS)); } + + #[test] + fn solo_chemistry_presets_set_the_geometry() { + let base = [ + "--readFilesIn", + "cdna.fq", + "bc.fq", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--sjdbGTFfile", + "genes.gtf", + ]; + let with = |extra: &[&str]| { + let mut a = base.to_vec(); + a.extend_from_slice(extra); + try_parse(&a) + }; + + // v2: 16-base CB, 10-base UMI. + let p = with(&["--soloChemistry", "SC3Pv2"]).unwrap(); + assert_eq!( + ( + p.solo_cb_start, + p.solo_cb_len, + p.solo_umi_start, + p.solo_umi_len + ), + (1, 16, 17, 10) + ); + + // v3 lengthened the UMI to 12. + let p = with(&["--soloChemistry", "SC3Pv3"]).unwrap(); + assert_eq!( + ( + p.solo_cb_start, + p.solo_cb_len, + p.solo_umi_start, + p.solo_umi_len + ), + (1, 16, 17, 12) + ); + + // v1 used a 14-base CB. + let p = with(&["--soloChemistry", "SC3Pv1"]).unwrap(); + assert_eq!((p.solo_cb_len, p.solo_umi_len), (14, 10)); + + // The preset wins over an explicit geometry rather than half-applying. + let p = with(&["--soloChemistry", "SC3Pv3", "--soloCBlen", "14"]).unwrap(); + assert_eq!(p.solo_cb_len, 16); + + // Unknown presets are refused. + assert!(with(&["--soloChemistry", "SC9Pv9"]).is_err()); + } + + #[test] + fn solo_cb_type_string_is_refused_rather_than_mispacked() { + let base = [ + "--readFilesIn", + "cdna.fq", + "bc.fq", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--sjdbGTFfile", + "genes.gtf", + "--soloCBstart", + "1", + "--soloCBlen", + "16", + "--soloUMIstart", + "17", + "--soloUMIlen", + "10", + ]; + let with = |extra: &[&str]| { + let mut a = base.to_vec(); + a.extend_from_slice(extra); + try_parse(&a) + }; + assert!(with(&["--soloCBtype", "Sequence"]).is_ok()); + // `String` needs a whitelist representation this does not have; packing + // a non-ACGT barcode into 2 bits per base would produce nonsense. + assert!(with(&["--soloCBtype", "String"]).is_err()); + } } From 9841e4ace74f8d73c23f990777a86b2a6f755a33 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 01:46:06 +0200 Subject: [PATCH 4/9] feat(solo): adapter-anchored CB_UMI_Complex geometry MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `--soloCBposition` and `--soloUMIposition` take four fields, `startAnchor_startDist_endAnchor_endDist`, where the anchor codes are 0 = read start, 1 = read end, 2 = adapter start, 3 = adapter end. Only code 0 was accepted; the rest were a parse error. Codes 2 and 3 are the reason this needs a new layout rather than a wider parser: the adapter moves from read to read, so the offsets cannot be resolved once at startup the way every other layout's are. `SoloBarcodeLayout::Anchored` carries the specs with their anchors and the adapter, and resolves per read. Adds `--soloAdapterSequence` and `--soloAdapterMismatchesNmax`. The adapter search takes the leftmost position within the mismatch budget, which is what STAR takes — the first acceptable match, not the best-scoring one. A read where the adapter is not found, or where an anchor resolves outside the read, yields no barcode. Guessing would be worse than declining: a barcode read at an unknown offset is not salvageable, and a wrong barcode is a read attributed to the wrong cell rather than a read lost. Tests cover the adapter moving (both fields follow it), one mismatch inside the budget, no adapter at all, an anchor running off the front of the read, and the leftmost-match rule. Co-Authored-By: Claude Opus 5 (1M context) --- src/params/mod.rs | 11 ++ src/solo/mod.rs | 265 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 276 insertions(+) diff --git a/src/params/mod.rs b/src/params/mod.rs index 175ad07..d5deb50 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1241,6 +1241,17 @@ pub struct Parameters { #[arg(long = "genomeLoad", default_value = "NoSharedMemory")] pub genome_load: String, + /// Adapter sequence that `--soloCBposition` / `--soloUMIposition` anchor + /// codes 2 and 3 are measured from. `-` (the default) means no adapter, so + /// only the read-start and read-end anchors are usable. + #[arg(long = "soloAdapterSequence", default_value = "-")] + pub solo_adapter_sequence: String, + + /// Mismatches tolerated when locating `--soloAdapterSequence` in the + /// barcode read. + #[arg(long = "soloAdapterMismatchesNmax", default_value_t = 1)] + pub solo_adapter_mismatches_nmax: usize, + /// `CB_UMI_Complex` cell-barcode segment positions, one per segment, as /// `startAnchor_startDist_endAnchor_endDist`. Only read-start anchoring /// (`anchor = 0`, fixed positions) is supported, e.g. `0_0_0_7 0_8_0_15`. diff --git a/src/solo/mod.rs b/src/solo/mod.rs index 22d63e0..3177278 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -53,11 +53,84 @@ pub enum SoloBarcodeLayout { cb_segments: Vec<(usize, usize)>, umi: (usize, usize), }, + /// Multi-segment CB whose positions are measured from an adapter rather + /// than from the read ends. + /// + /// The adapter moves from read to read, so unlike every other layout the + /// offsets cannot be resolved once at startup; each read is searched. A + /// read where the adapter is not found within the mismatch budget yields no + /// barcode at all, which is what STAR does — a barcode read at an unknown + /// offset is not salvageable. + Anchored { + cb_segments: Vec<(Anchored, Anchored)>, + umi: (Anchored, Anchored), + adapter: Vec, + max_mismatches: usize, + }, } /// Parse a `--soloCBposition`/`--soloUMIposition` spec /// (`startAnchor_startDist_endAnchor_endDist`) into a 0-based `(start, len)`. /// Only read-start anchoring (`anchor = 0`) is supported. +/// One end of a `--soloCBposition` / `--soloUMIposition` field: which anchor +/// it is measured from, and how far. +/// +/// STAR's anchor codes: 0 = read start, 1 = read end, 2 = adapter start, +/// 3 = adapter end. Codes 2 and 3 make the position depend on where the +/// adapter turns out to be in each read, which is why resolution happens per +/// read rather than once at startup. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct Anchored { + pub anchor: u8, + pub dist: i64, +} + +impl Anchored { + /// Absolute offset in a read of length `read_len`, given where the adapter + /// was found. `None` when the anchor needs an adapter that was not found, + /// or when the offset falls outside the read. + fn resolve(self, read_len: usize, adapter: Option<(usize, usize)>) -> Option { + let base = match self.anchor { + 0 => 0i64, + 1 => read_len as i64 - 1, + 2 => adapter?.0 as i64, + 3 => adapter?.1 as i64 - 1, + _ => return None, + }; + let pos = base + self.dist; + if pos < 0 || pos as usize >= read_len { + None + } else { + Some(pos as usize) + } + } +} + +/// Locate `adapter` in `seq`, allowing up to `max_mm` mismatches. +/// +/// Returns the half-open `(start, end)` of the first position that fits, which +/// is what STAR takes: the leftmost acceptable match, not the best-scoring one. +fn find_adapter(seq: &[u8], adapter: &[u8], max_mm: usize) -> Option<(usize, usize)> { + if adapter.is_empty() || adapter.len() > seq.len() { + return None; + } + for start in 0..=(seq.len() - adapter.len()) { + let mut mm = 0usize; + for (a, b) in seq[start..start + adapter.len()].iter().zip(adapter) { + if a != b { + mm += 1; + if mm > max_mm { + break; + } + } + } + if mm <= max_mm { + return Some((start, start + adapter.len())); + } + } + None +} + fn parse_position(spec: &str) -> Result<(usize, usize), Error> { let f: Vec<&str> = spec.split('_').collect(); if f.len() != 4 { @@ -91,6 +164,36 @@ fn invalid_pos(spec: &str, why: &str) -> Error { )) } +/// Parse a full `startAnchor_startDist_endAnchor_endDist` field, keeping the +/// anchors rather than collapsing them to read-start offsets. +fn parse_position_anchored(spec: &str) -> Result<(Anchored, Anchored), Error> { + let f: Vec<&str> = spec.split('_').collect(); + if f.len() != 4 { + return Err(invalid_pos( + spec, + "expected startAnchor_startDist_endAnchor_endDist", + )); + } + let num = |x: &str| x.parse::().ok(); + match (num(f[0]), num(f[1]), num(f[2]), num(f[3])) { + (Some(sa), Some(sd), Some(ea), Some(ed)) + if (0..=3).contains(&sa) && (0..=3).contains(&ea) => + { + Ok(( + Anchored { + anchor: sa as u8, + dist: sd, + }, + Anchored { + anchor: ea as u8, + dist: ed, + }, + )) + } + _ => Err(invalid_pos(spec, "anchor must be 0, 1, 2 or 3")), + } +} + impl SoloBarcodeLayout { /// Build the layout from CLI parameters. `CB_UMI_Complex` parses /// `--soloCBposition`/`--soloUMIposition`; otherwise fixed Simple geometry. @@ -102,6 +205,29 @@ impl SoloBarcodeLayout { .filter_map(|s| parse_position(s).ok()) .collect(); let umi = parse_position(¶ms.solo_umi_position).unwrap_or((0, 0)); + + // An adapter turns the position specs into per-read offsets, so it + // selects a different layout rather than being a tweak to this one. + if params.solo_adapter_sequence != "-" { + let adapter: Vec = params + .solo_adapter_sequence + .bytes() + .map(crate::io::fastq::encode_base) + .collect(); + let anchored_cb: Vec<_> = params + .solo_cb_position + .iter() + .filter_map(|s| parse_position_anchored(s).ok()) + .collect(); + if let Ok(anchored_umi) = parse_position_anchored(¶ms.solo_umi_position) { + return Self::Anchored { + cb_segments: anchored_cb, + umi: anchored_umi, + adapter, + max_mismatches: params.solo_adapter_mismatches_nmax, + }; + } + } return Self::Complex { cb_segments, umi }; } Self::Simple { @@ -121,6 +247,9 @@ impl SoloBarcodeLayout { umi_start, umi_len, } => (cb_start + cb_len).max(umi_start + umi_len), + // Not knowable without the read: the adapter's position decides. + // `extract` does its own bounds checking instead. + Self::Anchored { .. } => 0, Self::Complex { cb_segments, umi } => cb_segments .iter() .map(|&(s, l)| s + l) @@ -150,6 +279,39 @@ impl SoloBarcodeLayout { umi_seq: seq[*umi_start..umi_start + umi_len].to_vec(), umi_qual: slice_or_empty(qual, *umi_start, *umi_len), }), + Self::Anchored { + cb_segments, + umi, + adapter, + max_mismatches, + } => { + let found = find_adapter(seq, adapter, *max_mismatches); + let span = |(a, b): &(Anchored, Anchored)| -> Option<(usize, usize)> { + let s0 = a.resolve(seq.len(), found)?; + let e0 = b.resolve(seq.len(), found)?; + (e0 >= s0).then_some((s0, e0 - s0 + 1)) + }; + let mut cb_seq = Vec::new(); + let mut cb_qual = Vec::new(); + for segment in cb_segments { + let (s0, l) = span(segment)?; + if s0 + l > seq.len() { + return None; + } + cb_seq.extend_from_slice(&seq[s0..s0 + l]); + cb_qual.extend_from_slice(&slice_or_empty(qual, s0, l)); + } + let (us, ul) = span(umi)?; + if us + ul > seq.len() { + return None; + } + Some(CellBarcode { + cb_seq, + cb_qual, + umi_seq: seq[us..us + ul].to_vec(), + umi_qual: slice_or_empty(qual, us, ul), + }) + } Self::Complex { cb_segments, umi } => { let mut cb_seq = Vec::new(); let mut cb_qual = Vec::new(); @@ -1212,4 +1374,107 @@ mod tests { let mut reader = SoloReadReader::open(cdna.path(), bc.path(), v2_layout(), None).unwrap(); assert!(reader.read_batch(10).is_err()); } + + #[test] + fn adapter_anchoring_finds_the_adapter_and_positions_around_it() { + // Adapter ACGT sits at offset 4. CB is the 4 bases before it + // (anchor 2, the adapter start), UMI the 4 after (anchor 3, adapter + // end). Both move with the adapter, which is the whole point. + let adapter = vec![0u8, 1, 2, 3]; + let layout = SoloBarcodeLayout::Anchored { + cb_segments: vec![( + Anchored { + anchor: 2, + dist: -4, + }, + Anchored { + anchor: 2, + dist: -1, + }, + )], + umi: ( + Anchored { anchor: 3, dist: 1 }, + Anchored { anchor: 3, dist: 4 }, + ), + adapter: adapter.clone(), + max_mismatches: 1, + }; + + // GGGG ACGT TTTT → CB = GGGG, UMI = TTTT + let read = EncodedRead { + name: "r".into(), + sequence: vec![2, 2, 2, 2, 0, 1, 2, 3, 3, 3, 3, 3], + quality: vec![b'I'; 12], + }; + let bc = layout.extract(&read).expect("adapter present"); + assert_eq!(bc.cb_seq, vec![2, 2, 2, 2]); + assert_eq!(bc.umi_seq, vec![3, 3, 3, 3]); + + // Shift the adapter two bases right; both fields follow it. + let read = EncodedRead { + name: "r".into(), + sequence: vec![1, 1, 2, 2, 2, 2, 0, 1, 2, 3, 3, 3, 3, 3], + quality: vec![b'I'; 14], + }; + let bc = layout.extract(&read).expect("adapter present, shifted"); + assert_eq!(bc.cb_seq, vec![2, 2, 2, 2]); + assert_eq!(bc.umi_seq, vec![3, 3, 3, 3]); + + // One mismatch inside the adapter is within budget. + let read = EncodedRead { + name: "r".into(), + sequence: vec![2, 2, 2, 2, 0, 1, 2, 0, 3, 3, 3, 3], + quality: vec![b'I'; 12], + }; + assert!(layout.extract(&read).is_some()); + + // No adapter at all: no barcode. A barcode read at an unknown offset + // is not salvageable, so guessing would be worse than declining. + let read = EncodedRead { + name: "r".into(), + sequence: vec![2; 12], + quality: vec![b'I'; 12], + }; + assert!(layout.extract(&read).is_none()); + } + + #[test] + fn find_adapter_takes_the_leftmost_acceptable_match() { + // Two copies of ACGT; STAR takes the first that fits, not the best. + let seq = [0u8, 1, 2, 3, 9, 0, 1, 2, 3]; + assert_eq!(find_adapter(&seq, &[0, 1, 2, 3], 0), Some((0, 4))); + // Budget 0 and no exact match anywhere. + assert_eq!(find_adapter(&[2u8; 8], &[0, 1, 2, 3], 0), None); + // An adapter longer than the read cannot match. + assert_eq!(find_adapter(&[0u8, 1], &[0, 1, 2, 3], 4), None); + } + + #[test] + fn anchors_outside_the_read_yield_no_barcode() { + // Anchor 2 with a distance that runs off the front. + let layout = SoloBarcodeLayout::Anchored { + cb_segments: vec![( + Anchored { + anchor: 2, + dist: -99, + }, + Anchored { + anchor: 2, + dist: -1, + }, + )], + umi: ( + Anchored { anchor: 3, dist: 1 }, + Anchored { anchor: 3, dist: 2 }, + ), + adapter: vec![0, 1, 2, 3], + max_mismatches: 0, + }; + let read = EncodedRead { + name: "r".into(), + sequence: vec![0, 1, 2, 3, 3, 3, 3], + quality: vec![b'I'; 7], + }; + assert!(layout.extract(&read).is_none()); + } } From 3f776e6c9914dc0c68c53ad341b130d35b1e1bbf Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 09:20:15 +0200 Subject: [PATCH 5/9] fix(solo): the adapter search takes the best-scoring match, not the leftmost STAR's `localAlignHammingDist` (`SequenceFuns.cpp:341`) scores every offset and keeps the one with the fewest mismatches. This took the first offset that fit the mismatch budget, which is the opposite rule: an early approximate match beat a later exact one. Since `CB_UMI_Complex` measures its barcode positions from the adapter, choosing the wrong offset shifts the whole barcode, so this decides what the read is counted as, not merely where a match is reported. Reported by @Psy-Fer on #150, with the upstream source. The tests did not catch it because their fixtures had a single candidate or an exact one, where leftmost and best-scoring give the same answer. There is now a test where they differ. Two further details from the same function, both of which were missing: Ties go to the leftmost of the best, because STAR replaces its incumbent only on a strictly smaller distance. An `N` in the adapter never counts as a mismatch. It is a wildcard in the query, not a base that fails to match, and treating it as a mismatch spends budget that STAR does not spend. STAR also seeds its best distance with the adapter length, so an offset that mismatches everywhere can never win. That only shows when the budget is as large as the adapter, where it keeps the answer `None` instead of an arbitrary offset. Co-Authored-By: Claude Opus 5 (1M context) --- src/solo/mod.rs | 69 ++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 59 insertions(+), 10 deletions(-) diff --git a/src/solo/mod.rs b/src/solo/mod.rs index 3177278..94a7546 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -108,27 +108,48 @@ impl Anchored { /// Locate `adapter` in `seq`, allowing up to `max_mm` mismatches. /// -/// Returns the half-open `(start, end)` of the first position that fits, which -/// is what STAR takes: the leftmost acceptable match, not the best-scoring one. +/// STAR's `localAlignHammingDist` (`SequenceFuns.cpp:341`), called from +/// `SoloReadBarcode_getCBandUMI.cpp:340`. Three details decide the result and +/// none of them are optional: +/// +/// * It is **best-scoring, not leftmost**. Every offset is scored and the one +/// with the fewest mismatches wins. An early offset that merely fits the +/// budget loses to a later exact match, and since the barcode positions are +/// measured from the adapter, picking the wrong offset shifts the whole +/// barcode. +/// * Ties go to the **leftmost** of the best, because STAR replaces its best +/// only on a strictly smaller distance. +/// * An `N` in the adapter never counts as a mismatch. It is a wildcard in the +/// query, not a base that fails to match. +/// +/// STAR seeds its best distance with the adapter length, so an offset that +/// mismatches at every position can never win; that only matters when the +/// budget is as large as the adapter, where it keeps the answer `None` instead +/// of an arbitrary offset. fn find_adapter(seq: &[u8], adapter: &[u8], max_mm: usize) -> Option<(usize, usize)> { if adapter.is_empty() || adapter.len() > seq.len() { return None; } + const N_CODE: u8 = 4; + let mut best = adapter.len(); + let mut best_start = None; for start in 0..=(seq.len() - adapter.len()) { let mut mm = 0usize; for (a, b) in seq[start..start + adapter.len()].iter().zip(adapter) { - if a != b { + if *b != N_CODE && a != b { mm += 1; - if mm > max_mm { - break; + if mm >= best { + break; // cannot beat the incumbent } } } - if mm <= max_mm { - return Some((start, start + adapter.len())); + if mm < best { + best = mm; + best_start = Some(start); } } - None + let start = best_start?; + (best <= max_mm).then_some((start, start + adapter.len())) } fn parse_position(spec: &str) -> Result<(usize, usize), Error> { @@ -1438,15 +1459,43 @@ mod tests { assert!(layout.extract(&read).is_none()); } + /// The discriminating case: an early offset that fits the budget against a + /// later one that fits it better. STAR scores every offset and keeps the + /// best, so the later exact match wins. Taking the first acceptable one + /// instead would anchor the barcode four bases off. + #[test] + fn find_adapter_takes_the_best_scoring_match_not_the_leftmost() { + // ACGT at 0 with one mismatch, ACGT exactly at 5. + let seq = [0u8, 1, 2, 0, 9, 0, 1, 2, 3]; + assert_eq!(find_adapter(&seq, &[0, 1, 2, 3], 1), Some((5, 9))); + } + + /// Among equally good offsets STAR keeps the first, because it replaces its + /// best only on a strictly smaller distance. #[test] - fn find_adapter_takes_the_leftmost_acceptable_match() { - // Two copies of ACGT; STAR takes the first that fits, not the best. + fn find_adapter_breaks_ties_leftmost() { let seq = [0u8, 1, 2, 3, 9, 0, 1, 2, 3]; assert_eq!(find_adapter(&seq, &[0, 1, 2, 3], 0), Some((0, 4))); + } + + /// An `N` in the adapter is a wildcard in STAR's query, not a base that + /// fails to match, so it costs nothing against any read base. + #[test] + fn an_n_in_the_adapter_is_not_a_mismatch() { + let seq = [0u8, 1, 2, 3]; + // N at position 1 matches C without spending the budget. + assert_eq!(find_adapter(&seq, &[0, 4, 2, 3], 0), Some((0, 4))); + } + + #[test] + fn find_adapter_rejects_what_cannot_match() { // Budget 0 and no exact match anywhere. assert_eq!(find_adapter(&[2u8; 8], &[0, 1, 2, 3], 0), None); // An adapter longer than the read cannot match. assert_eq!(find_adapter(&[0u8, 1], &[0, 1, 2, 3], 4), None); + // Every position mismatches everywhere: STAR seeds its best distance + // with the adapter length, so nothing wins even at a budget that large. + assert_eq!(find_adapter(&[2u8; 8], &[0, 0, 0, 0], 4), None); } #[test] From a28018507c07e45bae121ebad30d8909bae08d47 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 09:20:15 +0200 Subject: [PATCH 6/9] revert(solo): drop --soloChemistry, which is not a STAR parameter `--soloChemistry` does not exist in STAR 2.7.11b; it was carried over from STAR-rs. Adding it here widens the CLI beyond STAR's, which is a positioning decision for this project rather than something to arrive inside a PR about the barcode posterior. Raised by @Psy-Fer on #150. The geometry it set is still reachable through the STAR flags it was a shorthand for: --soloCBstart, --soloCBlen, --soloUMIstart, --soloUMIlen. `--soloCBtype String` stays refused, as before: it is a real STAR parameter and refusing an unimplemented value is the honest answer. Co-Authored-By: Claude Opus 5 (1M context) --- src/params/mod.rs | 101 ------------------------------------- tests/parameter_surface.rs | 5 -- 2 files changed, 106 deletions(-) diff --git a/src/params/mod.rs b/src/params/mod.rs index d5deb50..c66c783 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1196,12 +1196,6 @@ pub struct Parameters { #[arg(long = "soloCBwhitelist", num_args = 1.., default_values_t = vec!["None".to_string()])] pub solo_cb_whitelist: Vec, - /// 10x chemistry preset. Sets the CB/UMI geometry so it does not have to be - /// spelled out with `--soloCBstart` and friends. `-` (the default) leaves - /// the geometry to those flags. - #[arg(long = "soloChemistry", default_value = "-")] - pub solo_chemistry: String, - /// How cell barcodes are represented: `Sequence` (2-bit packed ACGT, the /// default). #[arg(long = "soloCBtype", default_value = "Sequence")] @@ -1848,46 +1842,6 @@ impl Parameters { } } - // --soloChemistry presets. STAR applies these before validating the - // geometry, so a preset and an explicit --soloCBstart cannot disagree: - // the preset wins and says so, rather than half-applying. - if params.solo_chemistry != "-" { - let geometry = match params.solo_chemistry.as_str() { - // (cb_start, cb_len, umi_start, umi_len), 1-based as STAR takes them. - "CR_2" | "CR_3" | "CR_3.1" | "CR_4" | "SC3Pv1" => Some((1, 14, 15, 10)), - "SC3Pv2" => Some((1, 16, 17, 10)), - "SC3Pv3" | "SC3Pv4" | "SC5P" => Some((1, 16, 17, 12)), - _ => None, - }; - let Some((cb_start, cb_len, umi_start, umi_len)) = geometry else { - return Err(command.error( - ErrorKind::InvalidValue, - format!( - "unknown --soloChemistry '{}'; expected SC3Pv1, SC3Pv2, SC3Pv3, \ - SC3Pv4, SC5P, or -", - params.solo_chemistry - ), - )); - }; - for (flag, given) in [ - ("--soloCBstart", params.solo_cb_start), - ("--soloCBlen", params.solo_cb_len), - ("--soloUMIstart", params.solo_umi_start), - ("--soloUMIlen", params.solo_umi_len), - ] { - if given != 0 { - log::warn!( - "--soloChemistry {} overrides the geometry given by {flag}", - params.solo_chemistry - ); - } - } - params.solo_cb_start = cb_start; - params.solo_cb_len = cb_len; - params.solo_umi_start = umi_start; - params.solo_umi_len = umi_len; - } - // --soloCBtype: only 2-bit packed ACGT barcodes are supported. `String` // needs a different whitelist representation entirely, so refuse rather // than silently pack a non-ACGT barcode into nonsense. @@ -3082,61 +3036,6 @@ mod tests { assert!(!p.out_sam_attributes.contains(SamAttributes::XS)); } - #[test] - fn solo_chemistry_presets_set_the_geometry() { - let base = [ - "--readFilesIn", - "cdna.fq", - "bc.fq", - "--soloType", - "CB_UMI_Simple", - "--soloCBwhitelist", - "wl.txt", - "--sjdbGTFfile", - "genes.gtf", - ]; - let with = |extra: &[&str]| { - let mut a = base.to_vec(); - a.extend_from_slice(extra); - try_parse(&a) - }; - - // v2: 16-base CB, 10-base UMI. - let p = with(&["--soloChemistry", "SC3Pv2"]).unwrap(); - assert_eq!( - ( - p.solo_cb_start, - p.solo_cb_len, - p.solo_umi_start, - p.solo_umi_len - ), - (1, 16, 17, 10) - ); - - // v3 lengthened the UMI to 12. - let p = with(&["--soloChemistry", "SC3Pv3"]).unwrap(); - assert_eq!( - ( - p.solo_cb_start, - p.solo_cb_len, - p.solo_umi_start, - p.solo_umi_len - ), - (1, 16, 17, 12) - ); - - // v1 used a 14-base CB. - let p = with(&["--soloChemistry", "SC3Pv1"]).unwrap(); - assert_eq!((p.solo_cb_len, p.solo_umi_len), (14, 10)); - - // The preset wins over an explicit geometry rather than half-applying. - let p = with(&["--soloChemistry", "SC3Pv3", "--soloCBlen", "14"]).unwrap(); - assert_eq!(p.solo_cb_len, 16); - - // Unknown presets are refused. - assert!(with(&["--soloChemistry", "SC9Pv9"]).is_err()); - } - #[test] fn solo_cb_type_string_is_refused_rather_than_mispacked() { let base = [ diff --git a/tests/parameter_surface.rs b/tests/parameter_surface.rs index df3dd94..314eaca 100644 --- a/tests/parameter_surface.rs +++ b/tests/parameter_surface.rs @@ -115,11 +115,6 @@ const NOT_YET_ACCEPTED: &[&str] = &[ "genomeTransformOutput", "genomeType", "sjdbInsertSave", - // STARsolo barcode chemistry. - "soloAdapterMismatchesNmax", - "soloAdapterSequence", - "soloCBtype", - "soloOutFormatFeaturesGeneField3", ]; fn star_parameter_names() -> Vec { From 079455bfe7bd77da53a1a1da5bb440a81ece0d7d Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 30 Jul 2026 14:38:12 +0200 Subject: [PATCH 7/9] docs(divergence): record that all four homopolymer UMIs are rejected rustar-aligner has always rejected poly-A/C/G/T UMIs in every geometry. STAR rejects only poly-A under CB_UMI_Complex, because the four packed constants are built in the SoloReadBarcode constructor from pSolo.umiL while that path does not learn its UMI length until it has extracted one (SoloReadBarcode.cpp:16-21 against getCBandUMI.cpp:353-354), leaving all four constants zero. Poly-A packs to zero and is still caught; the other three are not. CONTRIBUTING requires a deliberate divergence to be recorded with the STAR source it departs from, and this one was not. Verified against the 2.7.11b source rather than asserted. Co-Authored-By: Claude Opus 5 (1M context) --- DIVERGENCE.md | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 66f4088..9dcac3a 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -63,6 +63,36 @@ The site that acts on those zeroed counts, however, gates on a different flag (` **Impact.** Degenerate inputs only — any dataset large enough to reach the significance test has far more than five distinct frequencies. On those inputs an uninitialised read can place arbitrary mass on unseen genes, which makes the multinomial log-probabilities meaningless; zero keeps them defined. **Source.** `src/solo/sgt.rs`, locked by `solo::sgt::tests::too_few_frequencies_leaves_the_unseen_mass_at_zero` (asserting the exact bit pattern, since the point is that nothing was written). STAR: `SoloFeature_emptyDrops_CR.cpp`, `SimpleGoodTuring/sgt.h`. +### 1.2 Homopolymer UMIs are rejected whichever base repeats + +**What STAR does.** STAR means to reject a UMI made of one repeated base, and under `CB_UMI_Simple` it does. Under `CB_UMI_Complex` it rejects only poly-A, letting poly-C, poly-G and poly-T through. + +The mechanism is an initialisation-order defect, not a rule. The four packed constants are built in the `SoloReadBarcode` constructor (`SoloReadBarcode.cpp:16-21`): + +```cpp +for (uint32 jj=0;jj<4;jj++) { + homoPolymer[jj]=0; + for (uint32 ii=0; ii Date: Thu, 30 Jul 2026 14:42:40 +0200 Subject: [PATCH 8/9] test(solo): name the cbMinP threshold and oneExact guard as their own tests Both rules were asserted inside resolve_multi_prefers_higher_prior, whose name describes only the ordering. The threshold decides between correcting a barcode and dropping the read, and oneExact is not implied by it (a sole candidate holds the whole posterior and is still refused), so each gets a test that says what it checks. Co-Authored-By: Claude Opus 5 (1M context) --- src/solo/count.rs | 41 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/src/solo/count.rs b/src/solo/count.rs index 3f53103..6652ed0 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -2717,4 +2717,45 @@ mod tests { UmiFiltering::MultiGeneUmi ); } + + /// Two candidates of equal quality, so the exact-count priors decide the + /// posterior, and STAR accepts the winner only once its share reaches + /// `CB_MIN_P` (0.975, `SoloReadBarcode_getCBandUMI.cpp`). Named apart from + /// the ordering test above because this is the threshold rather than the + /// ordering: it decides between correcting a barcode and dropping the read. + #[test] + fn cb_posterior_below_min_p_is_rejected() { + use crate::solo::whitelist::CbCandidate; + let cands = vec![ + CbCandidate { + wl_index: 0, + mismatch_pos: 1, + mismatch_qual: b'I', + }, + CbCandidate { + wl_index: 1, + mismatch_pos: 2, + mismatch_qual: b'I', + }, + ]; + // 10/13 = 0.769: a clear winner on ordering, still short of the bar. + assert_eq!(resolve_multi_cb(&cands, &[10, 3], 0.0), None); + // 1000/1003 = 0.997 clears it. + assert_eq!(resolve_multi_cb(&cands, &[1000, 3], 0.0), Some(0)); + } + + /// STAR's `oneExact`: a correction needs at least one candidate that was + /// seen exactly in the whitelist counts. A sole candidate with no exact + /// support holds the whole posterior and is still refused, so this guard is + /// not implied by the threshold. + #[test] + fn a_cb_never_seen_exactly_is_refused_without_pseudocounts() { + use crate::solo::whitelist::CbCandidate; + let one = vec![CbCandidate { + wl_index: 0, + mismatch_pos: 1, + mismatch_qual: b'I', + }]; + assert_eq!(resolve_multi_cb(&one, &[0, 0], 0.0), None); + } } From ae7a0a027fc02b1c5e8c0dc57b370909357a157d Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 30 Jul 2026 20:12:18 +0200 Subject: [PATCH 9/9] docs(changelog): record the solo barcode flags and the adapter fix Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 14b4346..c00b574 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -44,6 +44,12 @@ 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`. +- **STARsolo barcode resolution**: `--soloCBtype`, + `--soloAdapterSequence` / `--soloAdapterMismatchesNmax` for + adapter-anchored `CB_UMI_Complex` geometry, and + `--soloOutFormatFeaturesGeneField3`. Barcode correction now applies + STAR's `cbMinP` posterior threshold (0.975) and its `oneExact` + guard, so a low-confidence correction is dropped rather than accepted. - **STARsolo single-cell quantification (`--soloType`)** — the 10x Chromium / plate-based count-matrix pipeline, ported from STAR and @@ -186,6 +192,12 @@ Sections commonly used: Features, Bug fixes, Other changes. 2.7.11b the variant is a no-op. It now removes a UMI seen in two or more genes from **all** of them, the behaviour the option name describes. Recorded in `DIVERGENCE.md` (closes #144). +- The STARsolo adapter search took the leftmost match inside the + mismatch budget; STAR's `localAlignHammingDist` takes the + best-scoring one. Since `CB_UMI_Complex` measures barcode positions + from the adapter, the wrong offset shifted the whole barcode. Ties now + go to the leftmost of the best, and an `N` in the adapter is a + wildcard rather than a mismatch. - **STARsolo `Gene` assignment now requires exon concordance**, matching STARsolo: a read counts toward a gene only when every aligned block