solo: cbMinP posterior threshold, oneExact guard, adapter-anchored geometry (replaces #150)

CHANGELOG.md

@@ -13,6 +13,13 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 ### Features
 
+- **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
   verified against real STARsolo (#90).
@@ -99,6 +106,13 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 ### Bug fixes
 
+- 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
   lies within the gene's exons, rather than merely overlapping one. This

DIVERGENCE.md

@@ -30,6 +30,37 @@ This is the reason faithfulness is reported **tie-adjusted**. On the 10k yeast b
 
 **Source.** `src/rng.rs`, `src/align/read_align.rs` (`per_read_seed`, `shuffle_tied_prefix`), `src/params/mod.rs` (`MultimapperOrder`). STAR: `ReadAlign_multMapSelect.cpp`, `ReadAlignChunk` RNG seeding.
 
+### 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<pSolo.umiL;ii++)
+        homoPolymer[jj]=(homoPolymer[jj]<<2)+jj;
+};
+```
+
+but `CB_UMI_Complex` does not know its UMI length until it has extracted one, and sets it at read-processing time (`SoloReadBarcode_getCBandUMI.cpp:353-354`):
+
+```cpp
+if ( pSolo.umiL == 0 )
+    pSolo.umiL = umiSeq.size();
+```
+
+So when `umiL` is still 0 at construction, the loop body never runs and all four constants stay `0`. The check (`SoloReadBarcode_getCBandUMI.cpp:139`) then compares the packed UMI against `0` four times. Poly-A packs to `0`, because `A` is nucleotide 0, so poly-A is still caught; poly-C, poly-G and poly-T pack to non-zero values and match nothing.
+
+**What rustar-aligner does.** Rejects all four homopolymers, in every barcode geometry.
+
+**Why.** A UMI of one repeated base carries no information whichever base it is, so rejecting only poly-A is not a threshold choice to be matched but an accident of when a length becomes known. STAR's own `CB_UMI_Simple` path, where `umiL` is set from the parameters before construction, rejects all four; the complex path is the outlier. Matching the accident would mean counting poly-T UMIs as real molecules.
+
+**Impact.** Under `--soloType CB_UMI_Complex`, reads whose UMI is a non-A homopolymer are dropped here and counted by STAR. On real data these are rare and are sequencing artefacts. Note this is not a new behaviour in this PR: rustar-aligner has always rejected all four. What is new is that it is now deliberate, locked by a test, and recorded here rather than being an undocumented difference.
+
+**Source.** `src/solo/whitelist.rs` (`check_umi`, `umi_valid_rejects_every_homopolymer`). STAR: `SoloReadBarcode.cpp:16-21`, `SoloReadBarcode_getCBandUMI.cpp:139` and `:353-354`.
+
 ---
 
 ## 2. Cases where rustar-aligner outperforms STAR

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<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,
@@ -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<String>,
 
+    /// 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<String>,
@@ -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());
+    }
 }

src/solo/count.rs

@@ -258,11 +258,26 @@ fn directional(umis: &HashMap<u64, u32>, 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,
+            &params.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,
+                &params.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,
+            &params.solo_out_format_features_gene_field3,
             gzip,
         )?;
         write_barcodes(
@@ -1812,12 +1840,17 @@ fn write_cellranger_summary(
     Ok(())
 }
 
-/// `features.tsv`: `gene_id <TAB> gene_name <TAB> "Gene Expression"` (CellRanger
-/// v3 layout). We have no gene names, so the id is repeated.
+/// `features.tsv`: `gene_id <TAB> gene_name <TAB> <field3>` (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 <TAB> gene_name <TAB> "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,70 @@ 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);
+    }
+
+    /// 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);
     }
 }