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 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, + /// 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, @@ -1230,6 +1235,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`. @@ -1280,6 +1296,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, @@ -1818,6 +1842,19 @@ impl Parameters { } } + // --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. @@ -2998,4 +3035,36 @@ mod tests { assert_eq!(p.out_sam_strand_field, "None"); assert!(!p.out_sam_attributes.contains(SamAttributes::XS)); } + + #[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()); + } } diff --git a/src/solo/count.rs b/src/solo/count.rs index 4ddb414..6652ed0 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) } // --------------------------------------------------------------------------- @@ -1390,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( @@ -1467,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)?; @@ -1626,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( @@ -2066,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". @@ -2080,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(()) })?; @@ -2601,13 +2639,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] @@ -2662,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); + } } diff --git a/src/solo/mod.rs b/src/solo/mod.rs index 22d63e0..94a7546 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -53,11 +53,105 @@ 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. +/// +/// 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 *b != N_CODE && a != b { + mm += 1; + if mm >= best { + break; // cannot beat the incumbent + } + } + } + if mm < best { + best = mm; + best_start = Some(start); + } + } + let start = best_start?; + (best <= max_mm).then_some((start, start + adapter.len())) +} + fn parse_position(spec: &str) -> Result<(usize, usize), Error> { let f: Vec<&str> = spec.split('_').collect(); if f.len() != 4 { @@ -91,6 +185,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 +226,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 +268,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 +300,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 +1395,135 @@ 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()); + } + + /// 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_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] + 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()); + } } 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 + )); + } } 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 {