Batch PRs #173, #175, #156, #158, #152

CHANGELOG.md

@@ -114,13 +114,46 @@ Sections commonly used: Features, Bug fixes, Other changes.
   and a matrix produced elsewhere should be callable too. It streams
   the matrix into the same form the align path produces, so the filters
   are the identical code rather than a second implementation.
+- **`--soloCellFilter EmptyDrops_CR` now uses CellRanger's actual
+  statistics.** The ambient profile is smoothed with Simple Good-Turing,
+  as CellRanger and STAR do, instead of an approximation that reserved
+  unseen mass from the singleton rate and spread the remainder in
+  proportion to raw counts. The Monte-Carlo null is drawn with libc++'s
+  `std::mt19937` and `std::discrete_distribution`, seeded
+  `19760110 * (isim + 1)` per simulation as STAR seeds it, replacing a
+  SplitMix64 stream that could not agree with STAR's over an arbitrary
+  number of draws. Cell calls move as a result.
+
+- **`--soloFeatures Transcript3p`** quantifies transcripts rather than
+  genes, using how far each read's 3' end sits from each transcript's.
+  In a 3'-biased assay that distance discriminates between isoforms: a
+  read 200 bases from the end of one and 4000 from the end of another
+  is evidence for the first. The distance distribution is estimated
+  from the data, then used as the likelihood in an EM over UMIs. Output
+  is per *cluster* rather than per cell — `--soloClusterCBfile` (new,
+  and required for this feature) says which cell is in which cluster,
+  because one cell has too few UMIs to resolve isoforms. Reads sharing
+  a UMI contribute the intersection of their transcript sets, not the
+  union: they came from one molecule. Writes `matrix.mtx`,
+  `features.tsv` and `transcriptEndDistanceDistribution.txt` under
+  `Solo.out/Transcript3p/raw/`.
 
 ### Bug fixes
 
 - `--runThreadN 1` ran on every logical core instead of on one. The
   rayon pool was configured only above 1, and skipping it leaves rayon's
   default of one worker per core. Output is unchanged; the run now uses
   the thread count asked for.
+- `--soloUMIfiltering MultiGeneUMI_CR` kept every gene tied at the
+  highest read count; CellRanger gives a tied UMI to no gene at all.
+  Since one read per gene is the ordinary shape of a multi-gene UMI, the
+  flag removed nothing in practice. On a 20k-read 10x fixture the count
+  matrix moves from 16 465 to 15 414 against STAR's 15 423.
+- `--soloUMIfiltering MultiGeneUMI_All` was aliased to `MultiGeneUMI`,
+  which is neither STAR's behaviour nor the documented one: in STAR
+  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).
 
 - **STARsolo `Gene` assignment now requires exon concordance**, matching
   STARsolo: a read counts toward a gene only when every aligned block

DIVERGENCE.md

@@ -30,6 +30,40 @@ 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 `--soloUMIfiltering MultiGeneUMI_All` filters, rather than doing nothing
+
+**What STAR does.** Nothing, in effect — but not because the rule is unimplemented. `SoloFeature_collapseUMIall.cpp:79-88` implements exactly the documented behaviour, zeroing every gene for any UMI seen in more than one:
+
+```cpp
+if (pSolo.umiFiltering.MultiGeneUMI_All) {
+    for (auto &iu : umiGeneMapCount)
+        if (iu.second.size()>1)
+            for (auto &ig : iu.second) ig.second=0; //kill all genes for this UMI
+};
+```
+
+The site that acts on those zeroed counts, however, gates on a different flag (`:116`, `if (pSolo.umiFiltering.MultiGeneUMI && umiGeneMapCount[...]==0)`), and `MultiGeneUMI` and `MultiGeneUMI_All` are set in mutually exclusive branches (`ParametersSolo.cpp:457-462`). Selecting `MultiGeneUMI_All` therefore zeroes the counts and then never reads them, and the run reports unfiltered counts.
+
+**What rustar-aligner does.** The documented behaviour: a UMI seen in more than one gene is removed from all of them.
+
+**Why.** This is a one-line wiring bug in STAR, not a design decision: STAR's own code, immediately above, computes the documented result and then discards it. Matching the binary would mean shipping a flag that silently does nothing to anyone who read either the documentation or STAR's source, which is what #144 was raised about. So the divergence is from STAR's behaviour but *not* from its intent. Single-gene UMIs are untouched, which the tests check across every mode.
+
+**Impact.** Confined to `--soloUMIfiltering MultiGeneUMI_All`. The default (`-`) and the other filtering modes produce identical counts. Inverting the choice is a one-line change, since the test asserts the behaviour either way.
+
+**Source.** `src/solo/count.rs` (`UmiFiltering::MultiGeneUmiAll`, `filter_multi_gene_umi`), locked by `multigene_umi_all_drops_the_umi_from_every_gene` and `multigene_umi_all_parses_to_its_own_variant`. STAR: `SoloFeature_collapseUMIall.cpp`, `ParametersSolo.cpp`.
+
+### 1.3 `EmptyDrops_CR` Simple-Good-Turing with fewer than five distinct frequencies
+
+**What STAR does.** The ambient profile for `--soloCellFilter EmptyDrops_CR` is smoothed with Simple Good-Turing (Elworthy's `SimpleGoodTuring/sgt.h`). `analyse()` returns early, doing nothing, when the frequency spectrum has fewer than five distinct counts — Elworthy's `MinInput` guard. `PZero`, the mass reserved for genes unseen in the ambient droplets, is neither assigned in that case nor initialised at construction, so a caller that asks for it reads whatever the stack held.
+
+**What rustar-aligner does.** `PZero` is zero from construction.
+
+**Why.** There is nothing to reproduce: the value STAR reads is not a decision it made. With fewer than five distinct frequencies there is no basis for reserving unseen mass, and zero says so. Reproducing STAR would mean writing code whose correct behaviour is to emit an uninitialised value, and a test asserting it.
+
+**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`.
+
 ---
 
 ## 2. Cases where rustar-aligner outperforms STAR
@@ -75,7 +109,7 @@ On the 10k yeast PE benchmark, 4 reads differ in alignment score (AS) because ST
 
 **Impact.** Past libc++'s load factor the map rehashes, and the order then depends on the bucket count, which depends on how many distinct barcodes were seen; beyond that size the order diverges. The **values never do** — only which line they appear on. Reading the file by barcode rather than by position is unaffected either way.
 
-**Source.** `src/solo/cell_reads.rs`, locked by `rows_are_emitted_in_reverse_first_appearance_order` and `merging_partials_preserves_order_and_sums`. STAR: `SoloFeature_statsOutput.cpp`.
+**Source.** `src/solo/cell_reads.rs`, locked by `rows_are_emitted_in_reverse_first_appearance_order`. STAR: `SoloFeature_statsOutput.cpp`.
 
 ---
 
@@ -91,9 +125,7 @@ For `--quantMode TranscriptomeSAM`, rustar-aligner builds the per-transcript exo
 
 ### 4.2 In-tree RNG generator
 
-rustar-aligner uses an in-tree splitmix64 (`src/rng.rs`) rather than the `rand` crate, avoiding the `getrandom`/`zerocopy`/`ppv-lite86` dependency chain. This is the generator underlying §1.1; it is called out separately because it is a dependency/implementation choice independent of the tie-break policy.
-
----
+rustar-aligner uses an in-tree splitmix64 (`src/rng.rs`) rather than the `rand` crate, avoiding the `getrandom`/`zerocopy`/`ppv-lite86` dependency chain. This is the generator underlying §1.1; it is called out separately because it is a dependency/implementation choice independent of the tie-break policy. It is not the only in-tree generator: `--soloCellFilter EmptyDrops_CR` samples with a bit-exact libc++ `mt19937` (`src/solo/libcxx_rng.rs`) so its Monte-Carlo null matches STAR's — a convergence with STAR rather than a divergence from it.
 
 ## 5. Known residual single-read differences
 

src/params/mod.rs

@@ -1129,6 +1129,10 @@ pub struct Parameters {
     /// `CellReads.stats`. `-` (the default) names none.
     #[arg(long = "genomeChrSetMitochondrial", num_args = 1.., default_values_t = vec!["-".to_string()])]
     pub genome_chr_set_mitochondrial: Vec<String>,
+    /// Two-column `CB cluster` file assigning cells to clusters, for
+    /// `--soloFeatures Transcript3p`.
+    #[arg(long = "soloClusterCBfile")]
+    pub solo_cluster_cb_file: Option<PathBuf>,
 
     /// Cell-calling / matrix filtering: None, CellRanger2.2, EmptyDrops_CR, TopCells.
     #[arg(long = "soloCellFilter", num_args = 1.., default_values_t = vec!["CellRanger2.2".to_string(), "3000".to_string(), "0.99".to_string(), "10".to_string()])]
@@ -1663,15 +1667,15 @@ impl Parameters {
                     ));
                 }
             }
-            // Gene / GeneFull / SJ / Velocyto are implemented.
+            // Gene / GeneFull / SJ / Velocyto / Transcript3p are implemented.
             for f in &params.solo_features {
-                if !matches!(f.as_str(), "SJ" | "Velocyto")
+                if !matches!(f.as_str(), "SJ" | "Velocyto" | "Transcript3p")
                     && f.parse::<crate::solo::SoloFeature>().is_err()
                 {
                     return Err(command.error(
                         ErrorKind::InvalidValue,
                         format!(
-                            "unsupported --soloFeatures '{f}'; supported: Gene, GeneFull, SJ, Velocyto"
+                            "unsupported --soloFeatures '{f}'; supported: Gene, GeneFull, SJ, Velocyto, Transcript3p"
                         ),
                     ));
                 }
@@ -1758,6 +1762,23 @@ impl Parameters {
                     ));
                 }
             }
+            // STAR refuses `MultiGeneUMI_CR` unless the dedup is exactly
+            // `1MM_CR` — one value, that value (`ParametersSolo.cpp:463-468`).
+            // The rule exists because the filter decides ownership from the
+            // corrected-UMI map, which only the CellRanger dedup builds.
+            if params
+                .solo_umi_filtering
+                .iter()
+                .any(|f| f == "MultiGeneUMI_CR")
+                && (params.solo_umi_dedup.len() > 1
+                    || params.solo_umi_dedup.first().map(String::as_str) != Some("1MM_CR"))
+            {
+                return Err(command.error(
+                    ErrorKind::InvalidValue,
+                    "--soloUMIfiltering MultiGeneUMI_CR only works with --soloUMIdedup 1MM_CR\n\
+                     SOLUTION: rerun with --soloUMIfiltering MultiGeneUMI_CR --soloUMIdedup 1MM_CR",
+                ));
+            }
             // --soloCellReadStats: `CB` is the only value STAR defines.
             if !matches!(params.solo_cell_read_stats.as_str(), "CB" | "None") {
                 return Err(command.error(
@@ -1768,6 +1789,16 @@ impl Parameters {
                     ),
                 ));
             }
+            // Transcript3p quantifies per cluster, so it needs the clustering.
+            if params.solo_features.iter().any(|f| f == "Transcript3p")
+                && params.solo_cluster_cb_file.is_none()
+            {
+                return Err(command.error(
+                    ErrorKind::MissingRequiredArgument,
+                    "--soloFeatures Transcript3p requires --soloClusterCBfile: the EM runs \
+                     per cluster of cells, since one cell has too few UMIs to resolve isoforms",
+                ));
+            }
             // Validate --clipAdapterType.
             if !matches!(
                 params.clip_adapter_type.as_str(),
@@ -2629,6 +2660,64 @@ mod tests {
         assert!(AlignEndsType::from_str("Bogus").is_err());
     }
 
+    /// STAR refuses `MultiGeneUMI_CR` unless the dedup is exactly `1MM_CR`
+    /// (`ParametersSolo.cpp:463-468`): the filter decides ownership from the
+    /// corrected-UMI map, which only the CellRanger dedup builds. We accepted
+    /// the combination silently and counted with an uncorrected map.
+    #[test]
+    fn multi_gene_umi_cr_requires_the_cellranger_dedup() {
+        // The solo validation block only runs in solo mode, which is also the
+        // only mode where these flags mean anything.
+        let base = [
+            "--readFilesIn",
+            "cdna.fq",
+            "bc.fq",
+            "--soloType",
+            "CB_UMI_Simple",
+            "--sjdbGTFfile",
+            "genes.gtf",
+            "--soloCBwhitelist",
+            "wl.txt",
+            "--soloUMIfiltering",
+            "MultiGeneUMI_CR",
+        ];
+
+        // Paired with 1MM_CR: accepted.
+        let mut ok = base.to_vec();
+        ok.extend_from_slice(&["--soloUMIdedup", "1MM_CR"]);
+        assert!(try_parse(&ok).is_ok());
+
+        // Default dedup (1MM_All) and any other single value: refused.
+        assert!(try_parse(&base).is_err(), "default dedup should be refused");
+        let mut wrong = base.to_vec();
+        wrong.extend_from_slice(&["--soloUMIdedup", "Exact"]);
+        assert!(try_parse(&wrong).is_err());
+
+        // More than one dedup value is refused even when 1MM_CR is among them,
+        // matching STAR's `typesIn.size()>1` half of the condition.
+        let mut multi = base.to_vec();
+        multi.extend_from_slice(&["--soloUMIdedup", "1MM_CR", "Exact"]);
+        assert!(try_parse(&multi).is_err());
+
+        // The pairing rule applies only to MultiGeneUMI_CR.
+        assert!(
+            try_parse(&[
+                "--readFilesIn",
+                "cdna.fq",
+                "bc.fq",
+                "--soloType",
+                "CB_UMI_Simple",
+                "--sjdbGTFfile",
+                "genes.gtf",
+                "--soloCBwhitelist",
+                "wl.txt",
+                "--soloUMIfiltering",
+                "MultiGeneUMI"
+            ])
+            .is_ok()
+        );
+    }
+
     #[test]
     fn out_sam_order_accepts_star_values_rejects_others() {
         assert!(try_parse(&["--readFilesIn", "r.fq", "--outSAMorder", "Paired"]).is_ok());