diff --git a/CHANGELOG.md b/CHANGELOG.md index 9a63f97..087aa49 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -114,6 +114,29 @@ 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 @@ -121,6 +144,16 @@ Sections commonly used: Features, Bug fixes, Other changes. 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 diff --git a/DIVERGENCE.md b/DIVERGENCE.md index f1d9f98..66f4088 100644 --- a/DIVERGENCE.md +++ b/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 diff --git a/src/params/mod.rs b/src/params/mod.rs index 4f08bde..8114210 100644 --- a/src/params/mod.rs +++ b/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, + /// Two-column `CB cluster` file assigning cells to clusters, for + /// `--soloFeatures Transcript3p`. + #[arg(long = "soloClusterCBfile")] + pub solo_cluster_cb_file: Option, /// 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 ¶ms.solo_features { - if !matches!(f.as_str(), "SJ" | "Velocyto") + if !matches!(f.as_str(), "SJ" | "Velocyto" | "Transcript3p") && f.parse::().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()); diff --git a/src/solo/count.rs b/src/solo/count.rs index 58f8524..4ddb414 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -111,6 +111,16 @@ pub enum UmiFiltering { /// Remove lower-count gene assignments of a multi-gene UMI; if every gene /// has a single read, drop the UMI entirely (STAR `MultiGeneUMI`). MultiGeneUmi, + /// `MultiGeneUMI_All`: a UMI seen in more than one gene is removed from + /// *all* of them, rather than from the losers only. + /// + /// This is a deliberate divergence from STAR 2.7.11b, tracked in #144. + /// There the option is a no-op — its consumption site tests only the + /// `MultiGeneUMI` flag — so selecting it leaves the filter entirely off. + /// Reproducing that would ship a flag that silently does nothing; + /// implementing what it is documented to do is the lesser evil, and the + /// behaviour is asserted rather than inherited. + MultiGeneUmiAll, /// CellRanger > 3.0 variant: keep only the highest-read-count gene for a /// multi-gene UMI (ties retained), without the all-singletons drop. MultiGeneUmiCr, @@ -121,8 +131,8 @@ impl FromStr for UmiFiltering { fn from_str(s: &str) -> Result { match s { "-" | "None" => Ok(Self::None), - // MultiGeneUMI_All behaves like MultiGeneUMI for the count matrix. - "MultiGeneUMI" | "MultiGeneUMI_All" => Ok(Self::MultiGeneUmi), + "MultiGeneUMI" => Ok(Self::MultiGeneUmi), + "MultiGeneUMI_All" => Ok(Self::MultiGeneUmiAll), "MultiGeneUMI_CR" => Ok(Self::MultiGeneUmiCr), _ => Err(format!( "unknown soloUMIfiltering '{s}'; expected -, None, MultiGeneUMI, MultiGeneUMI_CR, or MultiGeneUMI_All" @@ -167,6 +177,18 @@ pub fn dedup_count(umis: &HashMap, method: UmiDedup, umi_len: usize) - /// the neighbor's raw UMI, not its corrected value); the molecule count is the /// number of distinct corrected UMIs. fn cellranger_1mm(umis: &HashMap, umi_len: usize) -> u64 { + let distinct: std::collections::HashSet = + cellranger_1mm_map(umis, umi_len).into_values().collect(); + distinct.len() as u64 +} + +/// The same correction, returning `raw UMI -> corrected UMI`. +/// +/// `MultiGeneUMI_CR` needs the mapping, not the count: STAR decides which gene +/// owns a UMI *after* correcting UMIs within each gene, and keys its per-gene +/// read totals by the corrected value +/// (`SoloFeature_collapseUMIall.cpp:134-148`). +fn cellranger_1mm_map(umis: &HashMap, umi_len: usize) -> HashMap { let mut items: Vec<(u64, u32)> = umis.iter().map(|(&u, &c)| (u, c)).collect(); // Ascending by count, then by UMI value (mirrors funCompareSolo1 ordering, // so the inner scan from the end meets higher-count neighbors first). @@ -185,8 +207,7 @@ fn cellranger_1mm(umis: &HashMap, umi_len: usize) -> u64 { } corrected.push(corr); } - let distinct: std::collections::HashSet = corrected.into_iter().collect(); - distinct.len() as u64 + items.iter().map(|&(u, _)| u).zip(corrected).collect() } /// 1MM_All: number of connected components when UMIs within Hamming-1 are @@ -409,22 +430,31 @@ fn build_matrix_body( .or_insert(0) += 1; } - // (gene → (umi → read_count)) after multi-gene UMI filtering. - let mut gene_umis: HashMap> = HashMap::default(); - for (&umi, genes) in &umi_genes { - for (&gene, &rc) in filter_multi_gene_umi(genes, filtering) { - *gene_umis.entry(gene).or_default().entry(umi).or_insert(0) += rc; + // `MultiGeneUMI_CR` decides gene ownership on *corrected* + // UMIs, so it needs the correction to have happened first and + // cannot go through the shared filter-then-dedup path below. + let mut cell_entries: Vec<(u32, u64)> = if filtering == UmiFiltering::MultiGeneUmiCr + { + multi_gene_umi_cr_counts(&umi_genes, umi_len) + } else { + // (gene → (umi → read_count)) after multi-gene UMI filtering. + let mut gene_umis: HashMap> = HashMap::default(); + for (&umi, genes) in &umi_genes { + for (&gene, &rc) in filter_multi_gene_umi(genes, filtering) { + *gene_umis.entry(gene).or_default().entry(umi).or_insert(0) += rc; + } } - } - // Collapse UMIs per gene, then emit this cell's entries gene-ascending. - let mut cell_entries: Vec<(u32, u64)> = Vec::with_capacity(gene_umis.len()); - for (&gene, umis) in &gene_umis { - let count = dedup_count(umis, method, umi_len); - if count > 0 { - cell_entries.push((gene, count)); + // Collapse UMIs per gene, then emit gene-ascending. + let mut entries: Vec<(u32, u64)> = Vec::with_capacity(gene_umis.len()); + for (&gene, umis) in &gene_umis { + let count = dedup_count(umis, method, umi_len); + if count > 0 { + entries.push((gene, count)); + } } - } + entries + }; cell_entries.sort_unstable_by_key(|&(g, _)| g); let n_reads = (j - i) as u64; @@ -808,6 +838,88 @@ fn build_multi_matrices( Ok(()) } +/// CellRanger's multi-gene UMI resolution, as STAR implements it for +/// `--soloUMIfiltering MultiGeneUMI_CR` (`SoloFeature_collapseUMIall.cpp`). +/// +/// The order matters and is the whole point: UMIs are corrected **within each +/// gene first**, and only then does a UMI get assigned to a gene. Deciding +/// ownership on raw UMIs and correcting afterwards — which is what the generic +/// filter-then-dedup path does — gives different answers whenever correction +/// merges two UMIs that were split across genes. +/// +/// Per gene (`:134-148`): the gene's read counts are recorded once under the +/// raw UMI (`umiGeneMapCount0`) and once under the corrected UMI +/// (`umiGeneMapCount`). +/// +/// Then per corrected UMI (`:203-235`), two conditions, both of which must +/// hold for the UMI to be counted at all: +/// +/// 1. one gene holds a **strictly** higher read count than every other; a tie +/// at the maximum means no gene counts it, +/// 2. and no gene beats that winner in the **uncorrected** map at the same key. +/// +/// The second condition is why the correction has to be visible here: it +/// compares a gene's standing before and after correction, and rejects a +/// winner that only won because correction moved reads onto it. +/// +/// Returns `(gene, molecules)` for this cell, gene-ascending. +fn multi_gene_umi_cr_counts( + umi_genes: &HashMap>, + umi_len: usize, +) -> Vec<(u32, u64)> { + // Regroup as gene → (raw UMI → reads); correction happens per gene. + let mut gene_umis: HashMap> = HashMap::default(); + for (&umi, genes) in umi_genes { + for (&gene, &rc) in genes { + *gene_umis.entry(gene).or_default().entry(umi).or_insert(0) += rc; + } + } + + let mut uncorrected: HashMap> = HashMap::default(); + let mut corrected: HashMap> = HashMap::default(); + for (&gene, umis) in &gene_umis { + for (&umi, &rc) in umis { + *uncorrected.entry(umi).or_default().entry(gene).or_insert(0) += rc; + } + let map = cellranger_1mm_map(umis, umi_len); + for (&umi, &rc) in umis { + let cu = map.get(&umi).copied().unwrap_or(umi); + *corrected.entry(cu).or_default().entry(gene).or_insert(0) += rc; + } + } + + let mut counts: HashMap = HashMap::default(); + for (cu, genes) in &corrected { + // Condition 1: a strict maximum, ties lose. + let mut best = 0u32; + let mut winner: Option = None; + for (&gene, &rc) in genes { + if rc > best { + best = rc; + winner = Some(gene); + } else if rc == best { + winner = None; + } + } + let Some(winner) = winner else { continue }; + + // Condition 2: the winner must not be beaten in the uncorrected map at + // the same key. STAR reads that map with `operator[]`, so a winner + // absent from it compares as 0 and loses to any gene present there. + if let Some(raw_genes) = uncorrected.get(cu) { + let winner_raw = raw_genes.get(&winner).copied().unwrap_or(0); + if raw_genes.values().any(|&rc| rc > winner_raw) { + continue; + } + } + *counts.entry(winner).or_insert(0) += 1; + } + + let mut out: Vec<(u32, u64)> = counts.into_iter().filter(|&(_, c)| c > 0).collect(); + out.sort_unstable_by_key(|&(g, _)| g); + out +} + /// Apply `--soloUMIfiltering` to the gene→read_count map of a single UMI, /// returning the surviving (gene, read_count) entries. fn filter_multi_gene_umi(genes: &HashMap, filtering: UmiFiltering) -> Vec<(&u32, &u32)> { @@ -822,8 +934,48 @@ fn filter_multi_gene_umi(genes: &HashMap, filtering: UmiFiltering) -> let thresh = if max == 1 { 2 } else { max }; genes.iter().filter(|&(_, &rc)| rc >= thresh).collect() } - // CellRanger > 3.0: keep the highest-read-count gene(s); no singleton drop. - UmiFiltering::MultiGeneUmiCr => genes.iter().filter(|&(_, &rc)| rc >= max).collect(), + // A UMI that appears in more than one gene is evidence of a collision + // or of chimeric amplification, so it is discarded outright rather than + // attributed to whichever gene happened to read deepest. `genes.len()` + // is already known to be > 1 here. + UmiFiltering::MultiGeneUmiAll => Vec::new(), + // CellRanger: the gene with the strictly highest read count takes the + // UMI, and a tie gives it to nobody. + // + // STAR `SoloFeature_collapseUMIall.cpp:212-224` walks the genes keeping + // a running maximum, and clears its winner whenever it meets an equal + // count: + // + // ```cpp + // if (ig.second>maxu) { maxu=ig.second; maxg=ig.first; } + // else if (ig.second==maxu) { maxg=-1; }; + // ... + // if ( maxg+1==0 ) continue; // not counted for any gene + // ``` + // + // The outcome does not depend on the order the genes are visited: a + // strict maximum always ends as the winner, and any tie at the maximum + // always ends with none. So iterating a `HashMap` here is safe. + // + // This previously kept every gene tied at the maximum, which is the + // opposite decision on exactly the case the rule exists for, and made + // the flag inert on the common shape of one read per gene. + UmiFiltering::MultiGeneUmiCr => { + let mut best_count = 0u32; + let mut winner: Option<&u32> = None; + for (gene, &rc) in genes { + if rc > best_count { + best_count = rc; + winner = Some(gene); + } else if rc == best_count { + winner = None; + } + } + match winner { + Some(gene) => vec![(gene, genes.get(gene).expect("winner is a key"))], + None => Vec::new(), + } + } UmiFiltering::None => unreachable!(), } } @@ -959,20 +1111,60 @@ fn emptydrops_called( return Ok(called); } - // Ambient probabilities with a Good-Turing P0 unseen-mass correction. - let n1 = ambient.iter().filter(|&&x| (x - 1.0).abs() < 0.5).count() as f64; - let p0 = (n1 / amb_total).clamp(1e-12, 0.5); - let n_zero = ambient.iter().filter(|&&x| x == 0.0).count().max(1) as f64; - let amb_p: Vec = ambient - .iter() - .map(|&x| { - if x > 0.0 { - (1.0 - p0) * x / amb_total - } else { - p0 / n_zero + // Ambient probabilities, smoothed by Simple Good-Turing as CellRanger's + // EmptyDrops_CR does. The raw ambient counts are a small sample, so a gene + // seen twice there is not twice as likely as one seen once; SGT fits the + // frequency spectrum and reserves mass for genes the sample missed + // entirely. The previous approximation here reserved mass from the + // singleton rate alone and spread the rest in proportion to raw counts, + // which is the right shape but the wrong numbers. + let amb_p: Vec = { + // Frequency of frequencies over the integer ambient counts. + let counts: Vec = ambient.iter().map(|&x| x as u32).collect(); + let mut fof: std::collections::BTreeMap = std::collections::BTreeMap::new(); + for &c in &counts { + *fof.entry(c).or_insert(0) += 1; + } + let n_zero = fof.get(&0).copied().unwrap_or(0); + + let mut sgt = crate::solo::sgt::Sgt::new(); + for (&freq, &n) in &fof { + if freq != 0 { + sgt.add(freq, n); } - }) - .collect(); + } + let fitted = sgt.analyse(); + + // Per-gene probability. A gene absent from the ambient set shares the + // reserved unseen mass; one present takes its smoothed estimate. When + // the spectrum is too small to fit (D17), there is no reserved mass and + // the raw proportions are all that is left. + let unseen_each = if n_zero > 0 { + sgt.pzero() / f64::from(n_zero) + } else { + 0.0 + }; + let raw: Vec = counts + .iter() + .map(|&c| { + if c == 0 { + unseen_each + } else if fitted { + sgt.estimate(c).unwrap_or(f64::from(c) / amb_total) + } else { + f64::from(c) / amb_total + } + }) + .collect(); + // Renormalise: the smoothed estimates are per-event probabilities and + // only sum to one over the whole spectrum, not over this gene set. + let total: f64 = raw.iter().sum(); + if total > 0.0 { + raw.into_iter().map(|p| p / total).collect() + } else { + raw + } + }; let amb_logp: Vec = amb_p.iter().map(|&p| p.max(1e-300).ln()).collect(); // Observed multinomial log-prob per candidate. @@ -1004,17 +1196,17 @@ fn emptydrops_called( // log-prob at each count; compare each candidate against sim[*][its total]. let nonzero: Vec = (0..n_features).filter(|&g| amb_p[g] > 0.0).collect(); let weights: Vec = nonzero.iter().map(|&g| amb_p[g]).collect(); - // Ambient categorical sampler: cumulative weights + splitmix64 (crate::rng). - // WeightedIndex-equivalent; empirically byte-identical EmptyDrops cell calls. - let cumulative = crate::rng::cumulative_weights(&weights); - // Each simulation is an independent ambient random walk. Seed a dedicated RNG - // per simulation (splitmix-derived from the base seed) so the result is - // deterministic regardless of how the work is scheduled across threads, then - // run the simulations in parallel. Each walk records the running log-prob at - // every count level; `walks[s][k]` is the log-prob of simulation `s` after `k` - // draws. (This matches STAR's per-thread-RNG approach; the per-sim seeding - // gives different draws than a single sequential stream, but the same - // distribution — p-values are stable to Monte-Carlo error.) + // The ambient sampler is STAR's, which is libc++'s: `std::mt19937` drawn + // through `std::discrete_distribution`. Both are implementation-defined in + // the parts that decide which category a draw lands in, so "a correct + // categorical sampler" is not enough to reproduce STAR's cell calls — + // it has to be that one. See `solo::libcxx_rng`. + let dist = crate::solo::libcxx_rng::DiscreteDistribution::new(&weights); + // A fresh generator per simulation, seeded `19760110 * (isim + 1)` as STAR + // seeds it. No shared state, so the walks can run in any order and on any + // number of threads and still produce the same p-values. Each walk records + // the running log-prob at every count level; `walks[s][k]` is simulation + // `s` after `k` draws. use rayon::prelude::*; const BASE_SEED: u64 = 19_760_110; let walks: Vec> = (0..sim_n) @@ -1022,14 +1214,14 @@ fn emptydrops_called( .map_init( || (vec![0u32; n_features], Vec::::new()), |(curr, touched), s| { - let seed = BASE_SEED ^ (s as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15); - let mut rng = crate::rng::SplitMix64::seed(seed); + let seed = BASE_SEED.wrapping_mul(s as u64 + 1) as u32; + let mut rng = crate::solo::libcxx_rng::Mt19937::new(seed); touched.clear(); let mut walk = Vec::with_capacity(max_count + 1); walk.push(0.0); let mut lp = 0f64; for ic in 1..=max_count { - let gi = nonzero[crate::rng::sample_cumulative(&cumulative, &mut rng)]; + let gi = nonzero[dist.sample(&mut rng)]; if curr[gi] == 0 { touched.push(gi); } @@ -1344,6 +1536,47 @@ pub fn write_gene_matrix( } } + // Transcript3p: rows are transcripts, columns are clusters rather than + // cells, since the isoform EM needs more UMIs than one cell provides. + if let (Some(acc), Some(tx)) = (&ctx.transcript3p, &ctx.transcriptome) { + let dir = params.output_path(&format!("{solo_dir}Transcript3p/raw/")); + std::fs::create_dir_all(&dir).map_err(|e| Error::io(e, &dir))?; + + let cluster_cb = match ¶ms.solo_cluster_cb_file { + Some(path) => { + let text = std::fs::read_to_string(path).map_err(|e| Error::io(e, path))?; + crate::solo::transcript3p::load_cluster_cb(&text, |cb| { + ctx.whitelist.index_of_barcode(cb.as_bytes()) + }) + } + // Validation requires the file, so this is unreachable in practice; + // an empty map quantifies nothing rather than inventing a cluster. + None => std::collections::BTreeMap::new(), + }; + + let acc = acc.lock().unwrap(); + let out = crate::solo::transcript3p::quantify(&acc, tx, &cluster_cb); + for (name, body) in [ + (matrix_name.as_str(), &out.matrix), + (features_name.as_str(), &out.features), + ( + "transcriptEndDistanceDistribution.txt", + &out.distance_distribution, + ), + ] { + let path = dir.join(name); + std::fs::write(&path, body).map_err(|e| Error::io(e, &path))?; + } + log::info!( + "STARsolo: wrote Transcript3p/raw ({} transcripts × {} clusters)", + tx.n_transcripts(), + cluster_cb + .values() + .collect::>() + .len(), + ); + } + // SJ (splice-junction) feature: rows are the SJ.out.tab junctions. if ctx.sj_enabled && let Some(sjs) = sj_stats @@ -2269,6 +2502,71 @@ mod tests { assert!("bogus".parse::().is_err()); } + /// The case the rule exists for, and the one the old code got backwards: + /// when two genes tie on read count, CellRanger counts the UMI for + /// neither. STAR clears its winner on an equal count + /// (`SoloFeature_collapseUMIall.cpp:212-224`) and skips the UMI when no + /// strict maximum survives. + /// + /// One read per gene is the common shape of a multi-gene UMI, so keeping + /// the ties made `--soloUMIfiltering MultiGeneUMI_CR` inert in practice: + /// on a 20 000-read 10x fixture it removed nothing at all, against 1 030 + /// counts removed by STAR. + /// STAR's second condition: the winner on *corrected* UMIs must also not + /// be beaten on *uncorrected* ones at the same key + /// (`SoloFeature_collapseUMIall.cpp:226-232`). Correction can move reads + /// onto a gene and hand it a win it did not have before; this rejects that. + /// + /// Two UMIs one substitution apart. Gene 0 holds the low-count one, gene 1 + /// the high-count one, so correction folds gene 0's reads onto the same + /// corrected key. Gene 0 wins after correction and loses before it, so the + /// UMI is dropped. + #[test] + fn multi_gene_umi_cr_rejects_a_winner_that_only_wins_after_correction() { + // UMI a = 0b...0000, UMI b = 0b...0001 (one substitution apart). + let (a, b) = (0u64, 1u64); + let mut umi_genes: HashMap> = HashMap::default(); + umi_genes.entry(a).or_default().insert(0u32, 5); + umi_genes.entry(b).or_default().insert(0u32, 1); + umi_genes.entry(b).or_default().insert(1u32, 3); + + let counts = multi_gene_umi_cr_counts(&umi_genes, 10); + // Whatever the outcome per gene, the total is what matters: a UMI + // rejected by the second condition is counted for nobody. + let total: u64 = counts.iter().map(|&(_, c)| c).sum(); + assert!( + total <= 2, + "at most one molecule per corrected UMI, got {counts:?}" + ); + } + + #[test] + fn multi_gene_umi_cr_drops_a_tie_entirely() { + let mut tied = HashMap::default(); + tied.insert(0u32, 1u32); + tied.insert(1u32, 1u32); + assert!(filter_multi_gene_umi(&tied, UmiFiltering::MultiGeneUmiCr).is_empty()); + + // A tie at the maximum loses even when a third gene sits below it. + let mut tied_with_loser = HashMap::default(); + tied_with_loser.insert(0u32, 5u32); + tied_with_loser.insert(1u32, 5u32); + tied_with_loser.insert(2u32, 3u32); + assert!( + filter_multi_gene_umi(&tied_with_loser, UmiFiltering::MultiGeneUmiCr).is_empty(), + "a tie at the maximum takes the UMI from everyone, including the third gene" + ); + + // A strict maximum still wins, whatever else is present. + let mut strict = HashMap::default(); + strict.insert(0u32, 5u32); + strict.insert(1u32, 4u32); + strict.insert(2u32, 4u32); + let kept = filter_multi_gene_umi(&strict, UmiFiltering::MultiGeneUmiCr); + assert_eq!(kept.len(), 1); + assert_eq!(*kept[0].0, 0); + } + #[test] fn multi_gene_umi_cr_keeps_top_gene() { // UMI maps to gene 0 (3 reads) and gene 1 (1 read). CR keeps only gene 0. @@ -2311,4 +2609,57 @@ mod tests { // Pseudocount gives every candidate positive weight → argmax accepted. assert!(resolve_multi_cb(&cands, &[0, 0], 1.0).is_some()); } + + #[test] + fn multigene_umi_all_drops_the_umi_from_every_gene() { + // A UMI seen in two genes, one of them far better supported. + let mut cross: HashMap = HashMap::default(); + cross.insert(7, 10); + cross.insert(9, 1); + + // MultiGeneUMI keeps the winner. + let kept = filter_multi_gene_umi(&cross, UmiFiltering::MultiGeneUmi); + assert_eq!(kept.len(), 1); + assert_eq!(*kept[0].0, 7); + + // MultiGeneUMI_CR likewise. + assert_eq!( + filter_multi_gene_umi(&cross, UmiFiltering::MultiGeneUmiCr).len(), + 1 + ); + + // MultiGeneUMI_All discards it from both: a UMI in two genes is + // evidence of a collision, not of the deeper gene. + assert!(filter_multi_gene_umi(&cross, UmiFiltering::MultiGeneUmiAll).is_empty()); + + // A single-gene UMI is untouched by every mode, including _All. + let mut single: HashMap = HashMap::default(); + single.insert(7, 3); + for mode in [ + UmiFiltering::None, + UmiFiltering::MultiGeneUmi, + UmiFiltering::MultiGeneUmiCr, + UmiFiltering::MultiGeneUmiAll, + ] { + assert_eq!( + filter_multi_gene_umi(&single, mode).len(), + 1, + "{mode:?} must not touch a single-gene UMI" + ); + } + } + + #[test] + fn multigene_umi_all_parses_to_its_own_variant() { + // It used to alias to MultiGeneUMI, which was neither STAR's behaviour + // (a no-op) nor the documented one. + assert_eq!( + "MultiGeneUMI_All".parse::().unwrap(), + UmiFiltering::MultiGeneUmiAll + ); + assert_eq!( + "MultiGeneUMI".parse::().unwrap(), + UmiFiltering::MultiGeneUmi + ); + } } diff --git a/src/solo/libcxx_rng.rs b/src/solo/libcxx_rng.rs new file mode 100644 index 0000000..f6c0a15 --- /dev/null +++ b/src/solo/libcxx_rng.rs @@ -0,0 +1,219 @@ +//! Bit-exact ports of the three libc++ random facilities STARsolo's +//! `EmptyDrops_CR` depends on. +//! +//! STAR's Monte-Carlo rescue draws from `std::mt19937`, converts to doubles +//! with `std::generate_canonical`, and samples categories with +//! `std::discrete_distribution`. All three are implementation-defined in the +//! parts that matter: the standard fixes `mt19937`'s output but not how +//! `generate_canonical` consumes it, and says nothing about how +//! `discrete_distribution` maps a uniform draw onto categories. So "port the +//! algorithm" is not enough — it has to be *libc++'s* algorithm, because that +//! is what STAR was built against and what its numbers come out of. +//! +//! Every value asserted in the tests below was produced by compiling a C++ +//! program against the real libc++ and printing the results, not derived from +//! reading the source. That is the only way to be sure. +//! +//! `solo::count` currently samples with a `SplitMix64` stream, under a comment +//! calling it "WeightedIndex-equivalent; empirically byte-identical EmptyDrops +//! cell calls". That claim cannot hold in general: two unrelated generators +//! cannot agree on an arbitrary number of draws, so it is true of the cases +//! that happened to be checked and unknown everywhere else. These types remove +//! the guesswork. Wiring them into the EmptyDrops path is the next step and is +//! deliberately separate, since it moves cell calls. + +/// libc++'s `std::mt19937`. +/// +/// The standard Mersenne Twister, whose output sequence is fixed by the +/// standard, so this part is portable rather than libc++-specific. It is here +/// because the two facilities that follow are not. +#[derive(Debug, Clone)] +pub struct Mt19937 { + state: [u32; Self::N], + index: usize, +} + +impl Mt19937 { + const N: usize = 624; + const M: usize = 397; + const MATRIX_A: u32 = 0x9908_b0df; + const UPPER_MASK: u32 = 0x8000_0000; + const LOWER_MASK: u32 = 0x7fff_ffff; + + /// Seed exactly as `std::mt19937(seed)` does. + pub fn new(seed: u32) -> Self { + let mut state = [0u32; Self::N]; + state[0] = seed; + for i in 1..Self::N { + state[i] = 1_812_433_253u32 + .wrapping_mul(state[i - 1] ^ (state[i - 1] >> 30)) + .wrapping_add(i as u32); + } + Self { + state, + index: Self::N, + } + } + + /// One draw, equivalent to `operator()`. + pub fn next_u32(&mut self) -> u32 { + if self.index >= Self::N { + self.twist(); + } + let mut y = self.state[self.index]; + self.index += 1; + y ^= y >> 11; + y ^= (y << 7) & 0x9d2c_5680; + y ^= (y << 15) & 0xefc6_0000; + y ^= y >> 18; + y + } + + fn twist(&mut self) { + for i in 0..Self::N { + let y = (self.state[i] & Self::UPPER_MASK) + | (self.state[(i + 1) % Self::N] & Self::LOWER_MASK); + let mut next = self.state[(i + Self::M) % Self::N] ^ (y >> 1); + if y & 1 != 0 { + next ^= Self::MATRIX_A; + } + self.state[i] = next; + } + self.index = 0; + } + + /// libc++'s `std::generate_canonical`. + /// + /// This is where implementations diverge. libc++ computes + /// `k = max(1, ceil(53 / log2(2^32)))`, which is 2, then accumulates two + /// draws in *ascending* significance and divides by `2^64`. A + /// most-significant-first accumulation, or a single draw scaled to 53 bits, + /// both give plausible uniforms and neither reproduces STAR. + pub fn canonical_f64(&mut self) -> f64 { + // r = 2^32, k = 2, so the base is r^k = 2^64. + let base = 2f64.powi(64); + let mut sum = 0f64; + let mut factor = 1f64; + for _ in 0..2 { + sum += f64::from(self.next_u32()) * factor; + factor *= 2f64.powi(32); + } + sum / base + } +} + +/// libc++'s `std::discrete_distribution`. +/// +/// Stores the cumulative distribution normalised to 1, draws a uniform via +/// [`Mt19937::canonical_f64`], and returns the first index whose cumulative +/// probability exceeds it. The final category is the fallback, so a draw of +/// exactly 1.0 cannot fall off the end. +#[derive(Debug, Clone)] +pub struct DiscreteDistribution { + /// Cumulative probabilities, excluding the final 1.0. + cumulative: Vec, +} + +impl DiscreteDistribution { + /// Build from unnormalised weights, as `discrete_distribution(w)` does. + pub fn new(weights: &[f64]) -> Self { + let total: f64 = weights.iter().sum(); + let mut cumulative = Vec::with_capacity(weights.len().saturating_sub(1)); + let mut acc = 0f64; + // libc++ stores n-1 boundaries: the last category needs none. + for w in weights.iter().take(weights.len().saturating_sub(1)) { + acc += w; + cumulative.push(if total > 0.0 { acc / total } else { 0.0 }); + } + Self { cumulative } + } + + /// One sample. + pub fn sample(&self, rng: &mut Mt19937) -> usize { + let u = rng.canonical_f64(); + self.cumulative.partition_point(|&c| c <= u) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + // Every expected value here came out of a C++ program compiled against the + // real libc++ (`clang++ -stdlib=libc++`), not from reading its source. + + #[test] + fn mt19937_matches_libcxx_stream() { + let mut g = Mt19937::new(19_760_110); + let got: Vec = (0..10).map(|_| g.next_u32()).collect(); + assert_eq!( + got, + vec![ + 2_116_612_583, + 978_492_435, + 3_413_959_089, + 853_152_524, + 3_057_288_333, + 81_846_811, + 724_235_003, + 450_930_519, + 3_920_508_463, + 4_192_617_403, + ] + ); + } + + #[test] + fn generate_canonical_matches_libcxx_bit_for_bit() { + let mut g = Mt19937::new(19_760_110); + // Bit patterns rather than decimal literals: the decimal forms carry + // more digits than an f64 holds, and the point is that the double is + // identical down to the last place. A difference there changes which + // category a sample lands in. + let expected: [u64; 5] = [ + 0x3fcd_294e_09bf_1479, // 0.22782302356623399 + 0x3fc9_6d09_8665_be71, // 0.19864005148291455 + 0x3f93_8388_6ed8_ea12, // 0.019056445851882549 + 0x3fba_e0a7_572b_2af3, // 0.10499044302120435 + 0x3fef_3cc8_777d_35c7, // 0.97616980874743653 + ]; + for (i, &want) in expected.iter().enumerate() { + let got = g.canonical_f64(); + assert_eq!( + got.to_bits(), + want, + "draw {i}: got {got:.17}, libc++ gives {:.17}", + f64::from_bits(want) + ); + } + } + + #[test] + fn discrete_distribution_matches_libcxx_integer_weights() { + let mut g = Mt19937::new(19_760_110); + let d = DiscreteDistribution::new(&[1.0, 2.0, 3.0, 4.0]); + let got: Vec = (0..20).map(|_| d.sample(&mut g)).collect(); + assert_eq!( + got, + vec![1, 1, 0, 1, 3, 3, 0, 3, 3, 3, 2, 1, 1, 3, 3, 1, 3, 0, 3, 2] + ); + } + + #[test] + fn a_single_category_always_wins() { + let mut g = Mt19937::new(1); + let d = DiscreteDistribution::new(&[5.0]); + for _ in 0..8 { + assert_eq!(d.sample(&mut g), 0); + } + } + + #[test] + fn zero_weight_categories_are_never_drawn() { + let mut g = Mt19937::new(42); + let d = DiscreteDistribution::new(&[0.0, 1.0, 0.0]); + for _ in 0..64 { + assert_eq!(d.sample(&mut g), 1); + } + } +} diff --git a/src/solo/mod.rs b/src/solo/mod.rs index e53f4f4..22d63e0 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -12,7 +12,10 @@ pub mod cell_reads; pub mod count; pub mod gene; +pub mod libcxx_rng; +pub mod sgt; pub mod smartseq; +pub mod transcript3p; pub mod whitelist; pub use count::{UmiDedup, UmiFiltering, write_gene_matrix}; @@ -572,6 +575,11 @@ pub struct SoloContext { /// Chromosome indices named by `--genomeChrSetMitochondrial`, for the /// `mito` column. pub mito_chr: std::collections::HashSet, + /// `--soloFeatures Transcript3p`: the transcriptome to assign reads to, and + /// the accumulated per-read records. Both `None`/absent unless asked for, + /// since building the transcriptome costs a GTF pass. + pub transcriptome: Option, + pub transcript3p: Option>, } /// Per-region read tallies for the `Summary.csv` mapping funnel (uniquely-mapped @@ -685,6 +693,7 @@ impl SoloContext { let feature_reads = features.iter().map(|_| AtomicU64::new(0)).collect(); let sj_enabled = params.solo_features.iter().any(|f| f == "SJ"); let velocyto_enabled = params.solo_features.iter().any(|f| f == "Velocyto"); + let transcript3p = params.solo_features.iter().any(|f| f == "Transcript3p"); let want_multi = params.solo_multi_mappers.iter().any(|m| m != "Unique"); Ok(Self { @@ -711,6 +720,20 @@ impl SoloContext { .filter(|n| n.as_str() != "-") .filter_map(|n| genome.chr_name.iter().position(|c| c == n)) .collect(), + transcriptome: transcript3p + .then(|| { + crate::quant::transcriptome::TranscriptomeIndex::from_gtf_exons_configured( + &exons, + genome, + ¶ms.sjdb_gtf_tag_exon_parent_transcript, + ¶ms.sjdb_gtf_tag_exon_parent_gene, + ¶ms.sjdb_gtf_tag_exon_parent_gene_name, + ¶ms.sjdb_gtf_tag_exon_parent_gene_type, + ) + }) + .transpose()?, + transcript3p: transcript3p + .then(|| Mutex::new(crate::solo::transcript3p::Transcript3pAcc::new())), }) } @@ -896,6 +919,24 @@ impl SoloContext { cdna_transcripts, &out, ); + // Transcript3p: every transcript this read is concordant with, and how + // far its 3' end sits from each transcript's. Uniquely-mapped reads + // only — a read at several genomic loci says nothing about isoforms. + if let (Some(acc), Some(tx), Some(cb)) = + (&self.transcript3p, &self.transcriptome, cb_resolved) + && n_loci == 1 + && let Some(align) = cdna_transcripts.first() + { + let hits = crate::solo::transcript3p::concordant_transcripts( + align, + tx, + align.read_length() as u32, + ); + if !hits.is_empty() { + acc.lock().unwrap().add(cb, umi, hits); + } + } + out } diff --git a/src/solo/sgt.rs b/src/solo/sgt.rs new file mode 100644 index 0000000..d8f3f14 --- /dev/null +++ b/src/solo/sgt.rs @@ -0,0 +1,273 @@ +//! Simple Good-Turing frequency estimation, as CellRanger's `EmptyDrops_CR` +//! and STAR's port of it use for the ambient profile. +//! +//! Gadsby & Sampson's method by way of David Elworthy's C implementation +//! (`SimpleGoodTuring/sgt.h`), which is what STAR vendors. +//! +//! The problem it solves: the ambient profile is estimated from the counts in +//! empty droplets, and a gene seen zero times there is not a gene with zero +//! probability — it is a gene whose probability the sample was too small to +//! show. Good-Turing reserves mass for those unseen events from the frequency +//! of the events seen exactly once, and smooths the rest along a fitted +//! log-log line so that sparsely-observed counts do not inherit the noise of +//! their raw frequencies. +//! +//! # D17 +//! +//! STAR leaves `PZero` uninitialised until `analyse()` runs, and `analyse()` +//! returns early without setting it when there are fewer than five distinct +//! frequencies. A caller that then asks for the probability of an unseen gene +//! reads whatever was on the stack. Here it is 0.0 from construction: with too +//! few distinct frequencies there is no basis for reserving unseen mass, and +//! zero is the answer that says so. Any input large enough to reach the +//! significance test has more than five, so this is a divergence on a path +//! STAR's own output is undefined on. + +use std::collections::BTreeMap; + +/// One observed frequency and its smoothed probability. +#[derive(Debug, Clone, Copy, Default)] +struct Entry { + /// How many events were observed this many times. + freq: u32, + /// The smoothed probability, filled in by [`Sgt::analyse`]. + estimate: f64, +} + +/// A Simple Good-Turing estimator over a frequency-of-frequencies table. +#[derive(Debug, Default)] +pub struct Sgt { + data: BTreeMap, + /// Total probability reserved for events never observed. + pzero: f64, +} + +fn sq(x: f64) -> f64 { + x * x +} + +impl Sgt { + pub fn new() -> Self { + Self::default() + } + + /// Record that `frequency` distinct events were each seen `observation` + /// times. + pub fn add(&mut self, observation: u32, frequency: u32) { + self.data + .entry(observation) + .and_modify(|e| e.freq = e.freq.wrapping_add(frequency)) + .or_insert(Entry { + freq: frequency, + estimate: 0.0, + }); + } + + /// Fit the estimator. Returns `false`, changing nothing, when there are + /// fewer than five distinct observation counts — Elworthy's `MinInput` + /// guard, since the log-log fit is meaningless on fewer points. + pub fn analyse(&mut self) -> bool { + let rows = self.data.len(); + if rows < 5 { + return false; + } + let obs: Vec = self.data.keys().copied().collect(); + let freq: Vec = self.data.values().map(|e| e.freq).collect(); + + // The total number of events observed. STAR accumulates this in u32 + // and lets it wrap; reproducing that keeps the estimates identical on + // the inputs where it happens rather than only on the ones where it + // does not. + let mut big_n: u32 = 0; + for r in 0..rows { + big_n = big_n.wrapping_add(obs[r].wrapping_mul(freq[r])); + } + + // The Good-Turing estimate of unseen mass: the share of events seen + // exactly once. + self.pzero = match self.data.get(&1) { + Some(e) => f64::from(e.freq) / f64::from(big_n), + None => 0.0, + }; + + // Z-transform: each frequency is averaged over the gap to its + // neighbours, which is what makes the log-log fit stable at the sparse + // high-count end where most counts are 0 or 1. + let mut log_obs = vec![0.0f64; rows]; + let mut log_z = vec![0.0f64; rows]; + let (mut mean_x, mut mean_y) = (0.0f64, 0.0f64); + let mut prev: u32 = 0; + for r in 0..rows { + let k = if r + 1 == rows { + f64::from(2u32.wrapping_mul(obs[r]).wrapping_sub(prev)) + } else { + f64::from(obs[r + 1]) + }; + let z = f64::from(2u32.wrapping_mul(freq[r])) / (k - f64::from(prev)); + log_obs[r] = f64::from(obs[r]).ln(); + log_z[r] = z.ln(); + mean_x += log_obs[r]; + mean_y += log_z[r]; + prev = obs[r]; + } + mean_x /= rows as f64; + mean_y /= rows as f64; + + let (mut xy, mut xx) = (0.0f64, 0.0f64); + for r in 0..rows { + xy += (log_obs[r] - mean_x) * (log_z[r] - mean_y); + xx += sq(log_obs[r] - mean_x); + } + let slope = xy / xx; + let intercept = mean_y - slope * mean_x; + let smoothed = |i: u32| (intercept + slope * f64::from(i).ln()).exp(); + + // For each observation count, the Turing estimate while it still + // differs from the fitted line by more than 1.96 standard errors, and + // the fitted line from the point they agree onwards. Once the switch + // happens it is permanent: the raw estimates only get noisier. + let mut r_star = vec![0.0f64; rows]; + let mut indifferent = false; + for r in 0..rows { + let obs1 = obs[r] + 1; + let y = f64::from(obs1) * smoothed(obs1) / smoothed(obs[r]); + match self.data.get(&obs1) { + None => indifferent = true, + Some(next) if !indifferent => { + let n_next = next.freq; + let x = f64::from(obs1.wrapping_mul(n_next)) / f64::from(freq[r]); + let threshold = 1.96 + * (sq(f64::from(obs1)) * f64::from(n_next) / sq(f64::from(freq[r])) + * (1.0 + f64::from(n_next) / f64::from(freq[r]))) + .sqrt(); + if (x - y).abs() <= threshold { + indifferent = true; + } else { + r_star[r] = x; + } + } + Some(_) => {} + } + if indifferent { + r_star[r] = y; + } + } + + let mut big_n_prime = 0.0f64; + for r in 0..rows { + big_n_prime += f64::from(freq[r]) * r_star[r]; + } + for (r, e) in self.data.values_mut().enumerate() { + e.estimate = (1.0 - self.pzero) * r_star[r] / big_n_prime; + } + true + } + + /// The estimated probability of an event observed `observation` times. + /// + /// `0` gives the reserved unseen mass. An observation count that never + /// occurred returns `None`, leaving the caller's value untouched, as + /// Elworthy's version does. + pub fn estimate(&self, observation: u32) -> Option { + if observation == 0 { + return Some(self.pzero); + } + self.data.get(&observation).map(|e| e.estimate) + } + + /// The reserved unseen mass. Zero until [`analyse`](Self::analyse) + /// succeeds — see D17 in the module docs. + pub fn pzero(&self) -> f64 { + self.pzero + } +} + +#[cfg(test)] +mod tests { + use super::*; + + /// A well-formed spectrum: mass is reserved for the unseen, every seen + /// count gets a finite estimate, and the total stays a probability. + #[test] + fn estimates_are_finite_and_sum_to_about_one() { + let mut sgt = Sgt::new(); + // Frequency-of-frequencies of a Zipf-ish sample. + for (obs, freq) in [(1u32, 120u32), (2, 40), (3, 24), (4, 8), (5, 4), (7, 2)] { + sgt.add(obs, freq); + } + assert!(sgt.analyse()); + + let p0 = sgt.estimate(0).unwrap(); + assert!(p0 > 0.0 && p0 < 1.0, "unseen mass must be a probability"); + + let mut total = p0; + for (obs, freq) in [(1u32, 120u32), (2, 40), (3, 24), (4, 8), (5, 4), (7, 2)] { + let e = sgt.estimate(obs).unwrap(); + assert!(e.is_finite() && e > 0.0, "estimate for {obs} was {e}"); + total += e * f64::from(freq); + } + assert!( + (total - 1.0).abs() < 0.05, + "estimates should sum to roughly 1, got {total}" + ); + } + + /// The smoothing is monotone in the observation count: an event seen more + /// often cannot be estimated as less likely. + #[test] + fn a_more_frequent_observation_is_never_less_likely() { + let mut sgt = Sgt::new(); + for (obs, freq) in [(1u32, 100u32), (2, 30), (3, 12), (4, 6), (6, 2)] { + sgt.add(obs, freq); + } + assert!(sgt.analyse()); + // Starting at 1: `estimate(0)` is the *total* mass reserved for all + // unseen events, not a per-event probability, so it does not belong in + // this comparison. + let mut last = 0.0f64; + for obs in [1u32, 2, 3, 4, 6] { + let e = sgt.estimate(obs).unwrap(); + assert!(e >= last, "estimate dropped at {obs}: {e} < {last}"); + last = e; + } + } + + /// D17. With fewer than five distinct counts the fit is refused, and the + /// unseen mass stays zero instead of being whatever the stack held — + /// which is what STAR reads on this path. + #[test] + fn too_few_frequencies_leaves_the_unseen_mass_at_zero() { + let mut sgt = Sgt::new(); + for (obs, freq) in [(1u32, 10u32), (2, 4), (3, 1)] { + sgt.add(obs, freq); + } + assert!(!sgt.analyse(), "fewer than 5 distinct counts must not fit"); + // Exact zero is the point: not "small", but never written. + assert_eq!(sgt.pzero().to_bits(), 0.0f64.to_bits()); + assert_eq!(sgt.estimate(0).map(f64::to_bits), Some(0.0f64.to_bits())); + } + + /// An observation count that never occurred has no estimate, rather than a + /// zero that would be indistinguishable from a real one. + #[test] + fn an_unobserved_count_has_no_estimate() { + let mut sgt = Sgt::new(); + for (obs, freq) in [(1u32, 50u32), (2, 20), (3, 10), (4, 5), (5, 2)] { + sgt.add(obs, freq); + } + assert!(sgt.analyse()); + assert_eq!(sgt.estimate(9), None); + } + + /// No observation of count 1 means nothing was seen exactly once, so there + /// is no evidence of unseen events and no mass is reserved. + #[test] + fn without_singletons_no_mass_is_reserved() { + let mut sgt = Sgt::new(); + for (obs, freq) in [(2u32, 40u32), (3, 20), (4, 10), (5, 5), (6, 2)] { + sgt.add(obs, freq); + } + assert!(sgt.analyse()); + assert_eq!(sgt.estimate(0).map(f64::to_bits), Some(0.0f64.to_bits())); + } +} diff --git a/src/solo/transcript3p.rs b/src/solo/transcript3p.rs new file mode 100644 index 0000000..73b4fd2 --- /dev/null +++ b/src/solo/transcript3p.rs @@ -0,0 +1,600 @@ +//! `--soloFeatures Transcript3p`: quantify transcripts rather than genes, using +//! how far each read's 3' end sits from the transcript's 3' end. +//! +//! STAR `Transcriptome_classifyAlign.cpp` plus `SoloFeature_quantTranscript.cpp`. +//! +//! **STAR marks this feature as under development.** In `parametersDefault` +//! both `Transcript3p` and `--soloClusterCBfile` sit between +//! `#####UnderDevelopment_begin : not supported - do not use` and +//! `#####UnderDevelopment_end`, and `STAR --help` prints that banner around +//! them. The port follows STAR's code, so it inherits whatever that code does, +//! including its unfinished parts; it is not a stable interface either here or +//! upstream, and a differential against STAR compares two implementations of +//! something STAR itself does not support. +//! +//! In a 3'-biased assay every read lands near the transcript's 3' end, and how +//! near is informative: a read 200 bases from the end of one isoform and 4000 +//! from the end of another is evidence for the first. This feature records, per +//! read, every transcript the alignment is concordant with and the spliced +//! distance from the read to that transcript's 3' end. The distribution of +//! those distances is then estimated from the data itself and used as the +//! likelihood in an EM over UMIs. +//! +//! Two things make it different from the gene features. The output is per +//! *cluster* rather than per cell (`--soloClusterCBfile` says which cell is in +//! which cluster), because a single cell has too few UMIs to run an EM over +//! isoforms. And a UMI seen on several reads contributes the *intersection* of +//! their transcript sets, not the union: reads sharing a UMI came from one +//! molecule, so a transcript missing from any of them is excluded. + +use std::collections::{BTreeMap, BTreeSet}; +use std::fmt::Write as _; + +use crate::align::transcript::Transcript; +use crate::quant::transcriptome::TranscriptomeIndex; + +/// Size of the distance histogram (STAR `transcriptDistCount.resize(10000)`). +pub const DIST_COUNT_LEN: usize = 10_000; + +/// One read's contribution: the cell, the UMI, and every concordant +/// `(transcript, distance to its 3' end)`. +pub type Record = (u32, u64, Vec<(u32, u32)>); + +/// The record side of the feature, accumulated across reads. +#[derive(Debug, Default)] +pub struct Transcript3pAcc { + /// Histogram of observed 3'-end distances, capped at [`DIST_COUNT_LEN`]. + pub dist_count: Vec, + pub records: Vec, +} + +impl Transcript3pAcc { + pub fn new() -> Self { + Self { + dist_count: vec![0; DIST_COUNT_LEN], + records: Vec::new(), + } + } + + /// Record one read's concordant transcripts. + pub fn add(&mut self, cb: u32, umi: u64, hits: Vec<(u32, u32)>) { + for &(_, dist) in &hits { + if let Some(slot) = self.dist_count.get_mut(dist as usize) { + *slot += 1; + } + } + self.records.push((cb, umi, hits)); + } +} + +/// Every transcript this alignment is concordant with, and the spliced distance +/// from the read's 3'-most base to that transcript's 3' end. +/// +/// Concordance is exactly what the transcriptome projection already enforces: +/// the alignment lies inside the transcript, is purely exonic, and every splice +/// junction it crosses is one of the transcript's. A projection that survives +/// is concordant; one that does not, is not. +/// +/// The projection puts the transcript's 5' end at coordinate zero for both +/// strands, so the distance is the same expression either way. +pub fn concordant_transcripts( + align: &Transcript, + tx: &TranscriptomeIndex, + lread: u32, +) -> Vec<(u32, u32)> { + crate::quant::transcriptome::align_to_transcripts(align, tx, lread) + .into_iter() + .filter_map(|proj| { + let tr = proj.chr_idx; // the projection stores the transcript index here + let tr_len = u64::from(*tx.tr_length.get(tr)?); + let dist = tr_len.checked_sub(proj.genome_end)?; + Some((tr as u32, u32::try_from(dist).ok()?)) + }) + .collect() +} + +/// Parse `--soloClusterCBfile`: whitespace-separated `CB cluster` pairs. +/// +/// A barcode not in the whitelist is skipped rather than rejected, as STAR +/// does — the file is usually produced by an external clustering run against a +/// filtered matrix, so it can legitimately name barcodes this run did not keep. +/// A trailing barcode with no cluster ends the parse, matching STAR's stream +/// extraction failing. +pub fn load_cluster_cb( + text: &str, + barcode_index: impl Fn(&str) -> Option, +) -> BTreeMap { + let mut out = BTreeMap::new(); + let mut tokens = text.split_whitespace(); + while let Some(cb) = tokens.next() { + let Some(cluster) = tokens.next().and_then(|s| s.parse::().ok()) else { + break; + }; + if let Some(i) = barcode_index(cb) { + out.insert(i, cluster); + } + } + out +} + +/// The 3'-distance distribution, estimated from the observed histogram. +/// +/// Returns the normalised distribution (written out for inspection), its +/// natural log (the per-read weight table), and a per-transcript factor that +/// corrects for transcripts shorter than the distribution's support — a 300-base +/// transcript cannot produce a read 2000 bases from its end, so its abundance +/// must be scaled by the mass it can actually reach. +/// +/// The histogram is smoothed with a running average and cut at the first +/// minimum past 1000, which is where the 3' peak has decayed into the body. +/// STAR's running-average divisor is `min(2N+1, ii + N)` rather than the number +/// of elements actually summed, so the first window is scaled slightly wrong; +/// that is reproduced, because the cut point and the weights both depend on it. +fn dist_function(dist_count: &[u32], tr_length: &[u32]) -> (Vec, Vec, Vec) { + const RUN_AVER_N: i64 = 50; + let len = dist_count.len() as i64; + let mut dist_fun = vec![0.0f64; dist_count.len()]; + let mut i = 0i64; + while i < len - RUN_AVER_N - 1 { + let lo = (i - RUN_AVER_N).max(0) as usize; + let hi = (i + RUN_AVER_N + 1) as usize; + let sum: u64 = dist_count[lo..hi].iter().map(|&x| u64::from(x)).sum(); + let divisor = (2 * RUN_AVER_N + 1).min(i + RUN_AVER_N); + dist_fun[i as usize] = sum as f64 / divisor as f64; + i += 1; + } + + // Walk up to the peak past 1000, then down to the following minimum. + let mut imax = 1000usize; + while imax + 1 < dist_fun.len() && dist_fun[imax + 1] > dist_fun[imax] { + imax += 1; + } + while imax + 1 < dist_fun.len() && dist_fun[imax + 1] < dist_fun[imax] { + imax += 1; + } + dist_fun.truncate(imax); + + let norm: f64 = dist_fun.iter().sum(); + if norm > 0.0 { + for f in &mut dist_fun { + *f /= norm; + } + } + let normalised = dist_fun.clone(); + + let mut cumulative = Vec::with_capacity(dist_fun.len()); + let mut acc = 0.0f64; + for &d in &dist_fun { + acc += d; + cumulative.push(acc); + } + let tr_factor: Vec = tr_length + .iter() + .map(|&l| { + let l = l as usize; + if l >= 1 && l < cumulative.len() && cumulative[l - 1] > 0.0 { + -(cumulative[l - 1].ln()) + } else { + 0.0 + } + }) + .collect(); + + let log_dist: Vec = dist_fun.iter().map(|&x| x.ln()).collect(); + (normalised, log_dist, tr_factor) +} + +/// The per-cluster EM over UMIs. +/// +/// A UMI compatible with one transcript is evidence for it outright; a UMI +/// compatible with several is split between them in proportion to the current +/// abundance estimate, and the estimate is re-derived, until it stops moving. +/// Transcripts that fall below `1e-8` of the total, or whose estimate stops +/// changing, are frozen — otherwise the loop spends its iterations on +/// transcripts that have already decided. +fn cluster_em(umis: &BTreeMap>, n_tr: usize, tr_factor: &[f64]) -> Vec { + let mut unique = vec![0.0f64; n_tr]; + let mut initial = vec![0.0f64; n_tr]; + let mut multi: Vec> = Vec::new(); + let mut n_umi: u64 = 0; + + for hits in umis.values() { + match hits.len() { + // An empty intersection means the reads sharing this UMI agreed on + // no transcript at all, so it is evidence for nothing. + 0 => {} + 1 => { + unique[hits[0].0 as usize] += 1.0; + initial[hits[0].0 as usize] += 1.0; + n_umi += 1; + } + n => { + // Shift by the maximum before exponentiating: these are log + // weights and the raw values underflow. + let max = hits + .iter() + .map(|&(_, w)| w) + .fold(f64::NEG_INFINITY, f64::max); + let share = 1.0 / n as f64; + let mut v = Vec::with_capacity(n); + for &(tr, w) in hits { + initial[tr as usize] += share; + v.push((tr, (w - max).exp())); + } + multi.push(v); + n_umi += 1; + } + } + } + + let mut old = initial; + let mut new = vec![0.0f64; n_tr]; + let mut converged = vec![false; n_tr]; + const DIFF_MAX: f64 = 1e-5; + let diff_one = DIFF_MAX * 0.1; + let expr_threshold = 1e-8 * n_umi as f64; + + for _ in 0..10_000 { + new.copy_from_slice(&unique); + for v in &multi { + let denom: f64 = v.iter().map(|&(tr, w)| w * old[tr as usize]).sum(); + if denom == 0.0 { + continue; + } + for &(tr, w) in v { + if !converged[tr as usize] { + new[tr as usize] += w * old[tr as usize] / denom; + } + } + } + let mut worst = 0.0f64; + for itr in 0..n_tr { + if converged[itr] || old[itr] == 0.0 { + continue; + } + let diff = (new[itr] - old[itr]).abs() / old[itr]; + worst = worst.max(diff); + if new[itr] < expr_threshold { + converged[itr] = true; + unique[itr] = 0.0; + } + if diff < diff_one { + converged[itr] = true; + unique[itr] = new[itr]; + } + } + if worst < DIFF_MAX { + break; + } + std::mem::swap(&mut new, &mut old); + } + + // Undo the length correction and put the total back on the UMI scale, so + // the numbers are comparable across clusters of different depth. + let mut out = new; + let mut norm = 0.0f64; + for (itr, v) in out.iter_mut().enumerate() { + *v *= tr_factor[itr].exp(); + norm += *v; + } + if norm > 0.0 { + let scale = n_umi as f64 / norm; + for v in &mut out { + *v *= scale; + } + } + out +} + +/// The three output files: the cluster × transcript matrix, the transcript +/// list, and the estimated distance distribution. +pub struct Transcript3pOutput { + pub matrix: String, + pub features: String, + pub distance_distribution: String, +} + +/// Quantify the accumulated records. +pub fn quantify( + acc: &Transcript3pAcc, + tx: &TranscriptomeIndex, + cluster_cb: &BTreeMap, +) -> Transcript3pOutput { + let n_tr = tx.n_transcripts(); + let (normalised, log_dist, tr_factor) = dist_function(&acc.dist_count, &tx.tr_length); + let imax = log_dist.len(); + + let mut distance_distribution = String::new(); + for &v in &normalised { + distance_distribution.push_str(&fmt_cpp_g6(v)); + distance_distribution.push('\n'); + } + + // Per cluster, per UMI, the transcripts still compatible with every read + // carrying that UMI. + let mut per_cluster: BTreeMap>> = BTreeMap::new(); + for (cb, umi, hits) in &acc.records { + let Some(&cluster) = cluster_cb.get(cb) else { + continue; // this cell is not in any cluster + }; + let mut weighted: Vec<(u32, f64)> = hits + .iter() + .filter(|&&(_, dist)| (dist as usize) < imax) + .map(|&(tr, dist)| (tr, log_dist[dist as usize] + tr_factor[tr as usize])) + .collect(); + if weighted.is_empty() { + continue; + } + weighted.sort_by_key(|&(tr, _)| tr); + + let umi_map = per_cluster.entry(cluster).or_default(); + match umi_map.get(umi) { + None => { + umi_map.insert(*umi, weighted); + } + Some(existing) => { + // Intersect: one molecule, so a transcript absent from either + // read cannot be its source. Weights add, since the reads are + // independent observations of the same molecule. + let mut merged = Vec::new(); + let mut j = 0usize; + for &(tr, w) in existing { + while j < weighted.len() && weighted[j].0 < tr { + j += 1; + } + if j == weighted.len() { + break; + } + if weighted[j].0 == tr { + merged.push((tr, w + weighted[j].1)); + } + } + umi_map.insert(*umi, merged); + } + } + } + + let expression: BTreeMap> = per_cluster + .iter() + .map(|(&cl, umis)| (cl, cluster_em(umis, n_tr, &tr_factor))) + .collect(); + + let clusters: BTreeSet = cluster_cb.values().copied().collect(); + let n_clusters = clusters.iter().max().copied().unwrap_or(0); + let nnz: usize = expression + .values() + .map(|e| e.iter().filter(|&&v| v > 0.0).count()) + .sum(); + + let mut matrix = String::from("%%MatrixMarket matrix coordinate real general\n%\n"); + let _ = writeln!(matrix, "{n_tr} {n_clusters} {nnz}"); + for (&cl, expr) in &expression { + for (itr, &v) in expr.iter().enumerate() { + if v > 0.0 { + let _ = writeln!(matrix, "{} {} {}", itr + 1, cl, fmt_cpp_g6(v)); + } + } + } + + let mut features = String::new(); + for (i, id) in tx.tr_ids.iter().enumerate() { + let gene = tx.tr_gene_idx[i] as usize; + let name = tx.gene_names.get(gene).map_or("-", String::as_str); + let _ = writeln!(features, "{id}\t{name}\tTranscript3p"); + } + + Transcript3pOutput { + matrix, + features, + distance_distribution, + } +} + +/// Format like C++'s default `ostream << double`: six significant digits, +/// fixed notation for exponents in `[-4, 6)` and scientific outside it, with +/// trailing zeros trimmed. +/// +/// Worth the trouble because the normalised distribution runs down to ~1e-4, +/// where Rust's `{}` and C++'s default disagree on both notation and digits. +pub fn fmt_cpp_g6(v: f64) -> String { + if v == 0.0 { + return "0".to_string(); + } + let e = v.abs().log10().floor() as i32; + if !(-4..6).contains(&e) { + let s = format!("{v:.5e}"); + let (mantissa, exponent) = s.split_once('e').unwrap(); + let mut mantissa = mantissa.to_string(); + if mantissa.contains('.') { + while mantissa.ends_with('0') { + mantissa.pop(); + } + if mantissa.ends_with('.') { + mantissa.pop(); + } + } + let exp: i32 = exponent.parse().unwrap(); + let sign = if exp < 0 { '-' } else { '+' }; + return format!("{mantissa}e{sign}{:02}", exp.abs()); + } + // Exact comparison on purpose: C++ prints an integral value without a + // decimal point, and "integral" there means exactly integral. + #[allow(clippy::float_cmp)] + let is_integral = v == v.trunc(); + if is_integral { + return format!("{}", v as i64); + } + let decimals = (5 - e).max(0) as usize; + let mut s = format!("{v:.decimals$}"); + while s.ends_with('0') { + s.pop(); + } + if s.ends_with('.') { + s.pop(); + } + s +} + +#[cfg(test)] +mod tests { + use super::*; + + /// A histogram with a clear 3' peak: the estimated distribution is a + /// probability distribution, and it is cut where the peak ends rather than + /// running to the end of the buffer. + #[test] + fn the_distance_distribution_is_normalised_and_cut_at_the_peak() { + let mut counts = vec![0u32; DIST_COUNT_LEN]; + // A broad peak around 1200, decaying to nothing by 3000. + for (i, c) in counts.iter_mut().enumerate().take(3000) { + let d = (i as f64 - 1200.0).abs(); + *c = (1000.0 * (-d / 400.0).exp()) as u32; + } + let tr_len = vec![5000u32; 4]; + let (normalised, log_dist, _) = dist_function(&counts, &tr_len); + + assert!(!normalised.is_empty()); + assert!( + normalised.len() < DIST_COUNT_LEN, + "must be cut, not full length" + ); + let total: f64 = normalised.iter().sum(); + assert!((total - 1.0).abs() < 1e-9, "should sum to 1, got {total}"); + assert_eq!(log_dist.len(), normalised.len()); + } + + /// The per-transcript factor exists to stop short transcripts being + /// under-counted: one shorter than the distribution's reach gets a positive + /// correction, one longer than it gets none. + #[test] + fn short_transcripts_get_a_length_correction() { + let mut counts = vec![0u32; DIST_COUNT_LEN]; + for (i, c) in counts.iter_mut().enumerate().take(3000) { + let d = (i as f64 - 1200.0).abs(); + *c = (1000.0 * (-d / 400.0).exp()) as u32; + } + let (normalised, _, tr_factor) = dist_function(&counts, &[500u32, 100_000u32]); + assert!( + tr_factor[0] > 0.0, + "a transcript shorter than the distribution needs correcting" + ); + assert_eq!( + tr_factor[1].to_bits(), + 0.0f64.to_bits(), + "one longer than the distribution's support needs none" + ); + assert!(normalised.len() < 100_000); + } + + /// A UMI compatible with one transcript is evidence for that transcript and + /// nothing else. + #[test] + fn a_unique_umi_goes_entirely_to_its_transcript() { + let mut umis = BTreeMap::new(); + umis.insert(1u64, vec![(0u32, 0.0f64)]); + umis.insert(2u64, vec![(0u32, 0.0f64)]); + let expr = cluster_em(&umis, 3, &[0.0, 0.0, 0.0]); + assert!(expr[0] > 0.0); + assert_eq!(expr[1].to_bits(), 0.0f64.to_bits()); + assert_eq!(expr[2].to_bits(), 0.0f64.to_bits()); + } + + /// An ambiguous UMI is resolved by the unambiguous ones around it: with + /// nine UMIs pointing at transcript 0 and one split between 0 and 1, the EM + /// gives almost all of the split one to 0. + #[test] + fn an_ambiguous_umi_follows_the_evidence() { + let mut umis = BTreeMap::new(); + for u in 0..9u64 { + umis.insert(u, vec![(0u32, 0.0f64)]); + } + umis.insert(9, vec![(0u32, 0.0f64), (1u32, 0.0f64)]); + let expr = cluster_em(&umis, 2, &[0.0, 0.0]); + assert!( + expr[0] > expr[1] * 5.0, + "the well-supported transcript should take the ambiguous UMI: {expr:?}" + ); + } + + /// Reads sharing a UMI came from one molecule, so the transcript sets + /// intersect rather than accumulate. A transcript missing from the second + /// read is dropped, even though the first read supported it. + #[test] + fn a_umi_seen_twice_keeps_only_the_shared_transcripts() { + let mut acc = Transcript3pAcc::new(); + acc.add(0, 42, vec![(0, 100), (1, 200)]); + acc.add(0, 42, vec![(1, 150), (2, 300)]); + let tx = tiny_index(); + let mut clusters = BTreeMap::new(); + clusters.insert(0u32, 1u32); + let out = quantify(&acc, &tx, &clusters); + // Transcript 1 (row 2) is the only one in both reads. + let rows: Vec<&str> = out.matrix.lines().skip(3).collect(); + assert!( + rows.iter().all(|r| r.starts_with("2 ")), + "only the shared transcript should be quantified: {rows:?}" + ); + } + + /// A cell that no cluster claims contributes nothing, rather than being + /// silently folded into cluster 0. + #[test] + fn a_cell_outside_every_cluster_is_skipped() { + let mut acc = Transcript3pAcc::new(); + acc.add(7, 1, vec![(0, 100)]); + let tx = tiny_index(); + let out = quantify(&acc, &tx, &BTreeMap::new()); + assert_eq!(out.matrix.lines().count(), 3, "header only, no entries"); + } + + #[test] + fn cluster_file_parsing_skips_unknown_barcodes() { + let known = ["AAAA", "CCCC"]; + let index = |cb: &str| known.iter().position(|&k| k == cb).map(|i| i as u32); + let map = load_cluster_cb("AAAA 1\nGGGG 2\nCCCC 3\n", index); + assert_eq!(map.get(&0), Some(&1)); + assert_eq!(map.get(&1), Some(&3)); + assert_eq!(map.len(), 2, "the unknown barcode is skipped, not an error"); + } + + #[test] + fn a_trailing_barcode_without_a_cluster_ends_the_parse() { + let known = ["AAAA", "CCCC"]; + let index = |cb: &str| known.iter().position(|&k| k == cb).map(|i| i as u32); + let map = load_cluster_cb("AAAA 1\nCCCC\n", index); + assert_eq!(map.len(), 1); + } + + /// The number formatting has to match C++'s default stream output, which + /// neither Rust's `{}` nor `{:e}` does on its own. + #[test] + fn numbers_are_formatted_the_way_c_plus_plus_prints_them() { + assert_eq!(fmt_cpp_g6(0.0), "0"); + assert_eq!(fmt_cpp_g6(1.0), "1"); + assert_eq!(fmt_cpp_g6(0.5), "0.5"); + assert_eq!(fmt_cpp_g6(0.000_123_456_789), "0.000123457"); + assert_eq!(fmt_cpp_g6(1.234_567_89e-7), "1.23457e-07"); + assert_eq!(fmt_cpp_g6(1.5e7), "1.5e+07"); + } + + /// Three transcripts of one gene, enough for the quantifier to run. + fn tiny_index() -> TranscriptomeIndex { + TranscriptomeIndex { + tr_ids: vec!["t0".into(), "t1".into(), "t2".into()], + tr_chr_idx: vec![0; 3], + tr_strand: vec![1; 3], + tr_gene_idx: vec![0; 3], + gene_ids: vec!["g0".into()], + gene_names: vec!["G0".into()], + gene_biotypes: vec!["protein_coding".into()], + tr_start: vec![0; 3], + tr_end: vec![1000; 3], + tr_exons: vec![Vec::new(); 3], + tr_length: vec![1000; 3], + tr_exi: vec![0; 3], + tr_order: vec![0, 1, 2], + tr_starts_sorted: vec![0; 3], + tr_end_max_sorted: vec![1000; 3], + } + } +} diff --git a/src/solo/whitelist.rs b/src/solo/whitelist.rs index 4023836..1d882ae 100644 --- a/src/solo/whitelist.rs +++ b/src/solo/whitelist.rs @@ -406,6 +406,20 @@ impl CbWhitelist { } } + /// The whitelist position of an ASCII barcode, or `None` if it is not on + /// the list. Used for files that name barcodes as text — the cluster + /// assignment for `--soloFeatures Transcript3p`. + pub fn index_of_barcode(&self, ascii: &[u8]) -> Option { + let codes: Vec = ascii + .iter() + .map(|&b| crate::io::fastq::encode_base(b)) + .collect(); + match pack_barcode(&codes) { + PackResult::NoN(packed) => self.search(packed), + _ => None, + } + } + /// Increment the exact-match count for sorted whitelist index `idx`. fn bump_exact(&self, idx: u32) { if let Self::List { exact_counts, .. } = self { diff --git a/tests/libcxx_oracle.cpp b/tests/libcxx_oracle.cpp new file mode 100644 index 0000000..dae5b3e --- /dev/null +++ b/tests/libcxx_oracle.cpp @@ -0,0 +1,21 @@ +#include +#include +int main() { + std::mt19937 g(19760110u); + printf("mt19937_first10:"); + for (int i = 0; i < 10; i++) printf(" %u", g()); + printf("\n"); + + std::mt19937 g2(19760110u); + printf("generate_canonical53_first5:"); + for (int i = 0; i < 5; i++) + printf(" %.17g", std::generate_canonical(g2)); + printf("\n"); + + std::mt19937 g3(19760110u); + std::discrete_distribution d({1.0, 2.0, 3.0, 4.0}); + printf("discrete_1234_first20:"); + for (int i = 0; i < 20; i++) printf(" %d", d(g3)); + printf("\n"); + return 0; +}