solo: cbMinP posterior threshold, oneExact guard, adapter-anchored geometry

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,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);
     }
 }