From 10f7c405701f70469c64c3e54e029f9495f1e74b Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 28 Jul 2026 23:58:06 +0200 Subject: [PATCH 1/6] 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 9af3edf..55dba72 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -258,11 +258,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], @@ -280,10 +295,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) } // --------------------------------------------------------------------------- @@ -2091,12 +2116,29 @@ 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); } } From f98c903c975a8f0a6f4d50303b22de6f0a45f614 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 00:10:52 +0200 Subject: [PATCH 2/6] 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 0536a85..571948c 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1123,6 +1123,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 55dba72..cea1ad1 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -1223,6 +1223,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( @@ -1279,6 +1280,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)?; @@ -1397,6 +1399,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( @@ -1837,12 +1840,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". @@ -1851,7 +1859,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 4023836..f63aa3d 100644 --- a/src/solo/whitelist.rs +++ b/src/solo/whitelist.rs @@ -737,4 +737,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 fde11336424688b79084ab83938cbbff28aa8866 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 00:18:37 +0200 Subject: [PATCH 3/6] 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 571948c..50b08cd 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1054,6 +1054,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, @@ -1559,6 +1570,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. @@ -2644,4 +2708,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 3227a24621f685ce031e2cb59f427291531213bc Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 01:46:06 +0200 Subject: [PATCH 4/6] 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 50b08cd..bc8a913 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1099,6 +1099,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 2253dee..d149e9e 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -49,11 +49,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 { @@ -87,6 +160,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. @@ -98,6 +201,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 { @@ -117,6 +243,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) @@ -146,6 +275,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(); @@ -1097,4 +1259,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 20494db54478f9ae9eb3c06d4267037fc18c5782 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 09:20:15 +0200 Subject: [PATCH 5/6] 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. --- 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 d149e9e..0c4d959 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -104,27 +104,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> { @@ -1323,15 +1344,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 c8a5d7db4a3f4207eb3d57bd82b1a0a13c548ac9 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 09:20:15 +0200 Subject: [PATCH 6/6] 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. --- src/params/mod.rs | 101 ---------------------------------------------- 1 file changed, 101 deletions(-) diff --git a/src/params/mod.rs b/src/params/mod.rs index bc8a913..645421d 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1054,12 +1054,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")] @@ -1581,46 +1575,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. @@ -2720,61 +2674,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 = [