diff --git a/src/params/mod.rs b/src/params/mod.rs index 0536a85..645421d 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1054,6 +1054,11 @@ pub struct Parameters { #[arg(long = "soloCBwhitelist", num_args = 1.., default_values_t = vec!["None".to_string()])] pub solo_cb_whitelist: Vec, + /// 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, @@ -1088,6 +1093,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`. @@ -1123,6 +1139,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, @@ -1551,6 +1575,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. @@ -2636,4 +2673,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 9af3edf..cea1ad1 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) } // --------------------------------------------------------------------------- @@ -1198,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( @@ -1254,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)?; @@ -1372,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( @@ -1812,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". @@ -1826,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(()) })?; @@ -2091,12 +2129,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); } } diff --git a/src/solo/mod.rs b/src/solo/mod.rs index 2253dee..0c4d959 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -49,11 +49,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 { @@ -87,6 +181,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 +222,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 +264,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 +296,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 +1280,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 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 + )); + } }