Solo: refuse what the read loops cannot do, check barcode and UMI input

CHANGELOG.md

@@ -11,6 +11,22 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 ## [Unreleased]
 
+### Bug fixes
+
+- Solo runs refuse options their read loops do not implement (`--quantMode
+  GeneCounts`/`TranscriptomeSAM`, `--quantGeneCoverage`, `--quantGeneSplicing`,
+  `--outReadsUnmapped Fastx`, `--outFilterType BySJout`, `--chimSegmentMin`,
+  `--waspOutputMode`, and more than one `--soloUMIdedup`) instead of writing
+  empty or zero output with exit 0.
+- A whitelist whose barcodes are not as long as the barcode read gives is now
+  fatal, as in STAR. It used to give an empty matrix, or false matches.
+- `CB_UMI_Complex`: a missing or bad position is fatal, and the UMI width comes
+  from `--soloUMIposition` rather than `--soloUMIlen`.
+- `--soloUMIlen` above 32 is refused rather than corrupting UMIs.
+- `--runMode soloCellFiltering` no longer writes gzip bytes into `features.tsv`.
+- `--sjdbGTFtagExonParentGeneName` applies when the GTF is given at mapping.
+- `CellReads.stats` counts each feature's reads in that feature's file.
+
 ### Other changes
 
 - Third-party notices completed. `LICENSES/` gains the licences of code

DIVERGENCE.md

@@ -171,6 +171,16 @@ A rustar-aligner extension with no STAR equivalent, off by default (`No`). With
 
 **Source.** `src/solo/count.rs` (`write_cellranger_summary`), `src/solo/gene.rs` (`Region`).
 
+### 1.10 Solo options that are refused rather than implemented
+
+**What STAR does.** STARsolo runs accept `--quantMode GeneCounts TranscriptomeSAM`, `--outReadsUnmapped Fastx`, `--outFilterType BySJout`, `--chimSegmentMin`, `--waspOutputMode`, and several `--soloUMIdedup` methods (one matrix per method). These descriptions come from STAR's documentation; solo runs of the Homebrew binary could not be compared on macOS, where it stops at `geneInfo.tab`.
+
+**What rustar-aligner does.** With `--soloType` it exits with an error naming the option. The same goes for this tool's own `--quantGeneCoverage` and `--quantGeneSplicing`.
+
+**Why.** The solo read loops have no code for these. They were accepted and wrote empty or all-zero files with exit 0, or dropped all methods but the first. An error is better than a result that looks real. Each one can be lifted when the loop implements it.
+
+**Source.** `src/params/mod.rs` (STARsolo validation).
+
 ---
 
 ## 2. Cases where rustar-aligner outperforms STAR

src/lib.rs

@@ -1051,7 +1051,7 @@ fn run_single_pass(
         }
         if let Some(info) = sctx.tag_read_info() {
             let info = info.lock().unwrap();
-            let umi_len = params.solo_umi_len as usize;
+            let umi_len = sctx.umi_len;
             writer.finish_with(&mut |record| {
                 crate::io::sam::apply_solo_read_info(record, &info, &sctx.whitelist, umi_len);
             })?;

src/params/mod.rs

@@ -2311,8 +2311,57 @@ impl Parameters {
 
         // ── STARsolo validation ─────────────────────────────────────────
         if params.run_mode() == RunMode::AlignReads && params.solo_enabled() {
+            // The solo read loops have no code for these, so they would write
+            // empty or all-zero files and exit 0 (DIVERGENCE §1.10).
+            let unsupported = [
+                (params.quant_gene_counts(), "--quantMode GeneCounts"),
+                (
+                    params.quant_transcriptome_sam(),
+                    "--quantMode TranscriptomeSAM",
+                ),
+                (params.quant_gene_coverage(), "--quantGeneCoverage"),
+                (params.quant_gene_splicing(), "--quantGeneSplicing"),
+                (
+                    params.out_reads_unmapped == OutReadsUnmapped::Fastx,
+                    "--outReadsUnmapped Fastx",
+                ),
+                (
+                    params.out_filter_type == OutFilterType::BySJout,
+                    "--outFilterType BySJout",
+                ),
+                (params.chim_segment_min > 0, "--chimSegmentMin"),
+                (
+                    params.wasp_output_mode != WaspOutputMode::None,
+                    "--waspOutputMode",
+                ),
+            ];
+            if let Some((_, flag)) = unsupported.iter().find(|(on, _)| *on) {
+                return Err(command.error(
+                    ErrorKind::ArgumentConflict,
+                    format!("{flag} is not implemented with --soloType, so it cannot be used here"),
+                ));
+            }
+            // Only the first method reaches the matrix (DIVERGENCE §1.10).
+            if params.solo_umi_dedup.len() > 1 {
+                return Err(command.error(
+                    ErrorKind::InvalidValue,
+                    format!(
+                        "--soloUMIdedup takes one method here, got {}",
+                        params.solo_umi_dedup.len()
+                    ),
+                ));
+            }
             // CB_UMI_Complex needs one CB position + whitelist per segment.
             if params.solo_type == SoloType::CbUmiComplex {
+                if params.solo_umi_position.is_empty() {
+                    return Err(command.error(
+                        ErrorKind::MissingRequiredArgument,
+                        "--soloType CB_UMI_Complex requires --soloUMIposition",
+                    ));
+                }
+                if let Err(e) = crate::solo::check_positions(&params) {
+                    return Err(command.error(ErrorKind::InvalidValue, e.to_string()));
+                }
                 if params.solo_cb_position.is_empty() {
                     return Err(command.error(
                         ErrorKind::MissingRequiredArgument,
@@ -2537,6 +2586,20 @@ impl Parameters {
                     "--soloCBlen and --soloUMIlen must be > 0 for soloType CB_UMI_Simple",
                 ));
             }
+            // UMIs are packed in a u64 too. CB_UMI_Complex takes its width from
+            // --soloUMIposition, checked above.
+            if params.solo_type != SoloType::CbUmiComplex
+                && params.solo_umi_len as usize > crate::solo::whitelist::CB_LEN_MAX
+            {
+                return Err(command.error(
+                    ErrorKind::InvalidValue,
+                    format!(
+                        "--soloUMIlen {} exceeds the maximum of {}",
+                        params.solo_umi_len,
+                        crate::solo::whitelist::CB_LEN_MAX
+                    ),
+                ));
+            }
             // Cell barcode cannot exceed a u64 packing (32 bases).
             if params.solo_cb_len as usize > crate::solo::whitelist::CB_LEN_MAX {
                 return Err(command.error(

src/quant/mod.rs

@@ -69,7 +69,12 @@ impl GeneAnnotation {
     /// Build from GTF exon records using `gene_tag` as the gene grouping attribute.
     ///
     /// `gene_tag` is STAR's `sjdbGTFtagExonParentGene` (default `"gene_id"`).
-    pub fn from_gtf_exons_configured(exons: &[GtfRecord], genome: &Genome, gene_tag: &str) -> Self {
+    pub fn from_gtf_exons_configured(
+        exons: &[GtfRecord],
+        genome: &Genome,
+        gene_tag: &str,
+        gene_name_tags: &[String],
+    ) -> Self {
         let n_chrs = genome.n_chr_real;
         let mut gene_ids: Vec<String> = Vec::new();
         let mut gene_names: Vec<String> = Vec::new();
@@ -96,11 +101,9 @@ impl GeneAnnotation {
                 gene_is_reverse.push(is_rev);
                 // gene_name (symbol) for features.tsv col 2; STAR falls back to
                 // the gene_id when the GTF record omits gene_name.
-                let gene_name = exon
-                    .attributes
-                    .get("gene_name")
-                    .cloned()
-                    .unwrap_or_else(|| gene_id.clone());
+                let gene_name =
+                    crate::quant::transcriptome::find_attr_list(&exon.attributes, gene_name_tags)
+                        .unwrap_or_else(|| gene_id.clone());
                 gene_names.push(gene_name);
                 gene_ids.push(gene_id);
                 gene_span.push(None);
@@ -290,7 +293,7 @@ impl GeneAnnotation {
 
     /// Build from GTF exon records using default `"gene_id"` attribute (backward-compatible).
     pub fn from_gtf_exons(exons: &[GtfRecord], genome: &Genome) -> Self {
-        Self::from_gtf_exons_configured(exons, genome, "gene_id")
+        Self::from_gtf_exons_configured(exons, genome, "gene_id", &["gene_name".to_string()])
     }
 
     pub fn n_genes(&self) -> usize {
@@ -629,6 +632,7 @@ pub fn resolve_gene_annotation(
             &exons,
             genome,
             &params.sjdb_gtf_tag_exon_parent_gene,
+            &params.sjdb_gtf_tag_exon_parent_gene_name,
         ));
     }
     match transcriptome {

src/quant/transcriptome.rs

@@ -31,7 +31,7 @@ const MISSING_GENE_TYPE: &str = "MissingGeneType";
 /// STAR's `find_attr_list` (`GTF.cpp`): the last of `tags`, in list order, that has an entry in
 /// `attrs`. `--sjdbGTFtagExonParentGeneName`/`GeneType` accept several candidate keys; when more
 /// than one matches, the last one in the list wins.
-fn find_attr_list(attrs: &HashMap<String, String>, tags: &[String]) -> Option<String> {
+pub(crate) fn find_attr_list(attrs: &HashMap<String, String>, tags: &[String]) -> Option<String> {
     tags.iter()
         .filter_map(|t| attrs.get(t))
         .next_back()

src/solo/count.rs

@@ -1534,7 +1534,7 @@ pub fn write_gene_matrix(
     } else {
         0.0
     };
-    let umi_len = params.solo_umi_len as usize;
+    let umi_len = ctx.umi_len;
 
     let solo_dir = params
         .solo_out_file_names
@@ -1654,7 +1654,7 @@ pub fn write_gene_matrix(
 
         // `--soloCellReadStats CB`: the per-cell read summary, alongside the
         // raw matrix because its UMI and gene columns are the raw totals.
-        if let Some(cell_stats) = &ctx.cell_read_stats {
+        if let Some(cell_stats) = ctx.cell_read_stats.as_ref().map(|s| &s[fi]) {
             let umi_gene: std::collections::BTreeMap<u32, (u32, u32)> = mstats
                 .cells
                 .iter()
@@ -1875,7 +1875,7 @@ pub fn write_gene_matrix(
             sorted.len(),
             gzip,
         )?;
-        let umi_len = params.solo_umi_len as usize;
+        let umi_len = ctx.umi_len;
         // STAR's per-line filters (`SoloReadFeature_inputRecords.cpp:23-104`),
         // counted only for a junction with a row: STAR skips the others first.
         let has_row = |start: u64, end: u64| row.contains_key(&(start, end));
@@ -2799,7 +2799,13 @@ pub fn run_cell_filtering(params: &crate::params::Parameters) -> anyhow::Result<
     std::fs::write(&bc_path, bc).map_err(|e| Error::io(e, &bc_path))?;
 
     let feat_path = dir.join(format!("{prefix}features.tsv"));
-    std::fs::copy(&features_path, &feat_path).map_err(|e| Error::io(e, &feat_path))?;
+    // The filtered matrix is written plain, whatever the raw one was.
+    let mut feat_out = std::fs::File::create(&feat_path).map_err(|e| Error::io(e, &feat_path))?;
+    std::io::copy(
+        &mut crate::solo::whitelist::open_maybe_gzip(&features_path)?,
+        &mut feat_out,
+    )
+    .map_err(|e| Error::io(e, &feat_path))?;
 
     let mtx_path = dir.join(format!("{prefix}matrix.mtx"));
     let nnz = finalize_matrix(