From 4cc1489a9768f7883c199241971408f24e0933ba Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 09:55:10 +0200 Subject: [PATCH 1/4] feat(solo): --soloFeatures Transcript3p, with --soloClusterCBfile MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Quantifies transcripts rather than genes, from where each read's 3' end sits relative to each transcript's. In a 3'-biased assay that distance is what separates isoforms: a read 200 bases from the end of one and 4000 from the end of another is evidence for the first. The distribution of those distances is estimated from the run's own histogram, smoothed and cut where the 3' peak decays into the body, and used as the likelihood in an EM over UMIs. Concordance needed no new code. `align_to_transcripts` already refuses to project an alignment that leaves the transcript, touches an intron, or crosses a junction the transcript does not have — which is exactly STAR's `Concordant` (`Transcriptome_classifyAlign.cpp`). A projection that survives is concordant; one that does not, is not. The projection also puts the 5' end at coordinate zero for both strands, so the distance to the 3' end is one expression rather than two. Two behaviours worth stating because they are not the obvious ones: Output is per cluster, not per cell, and `--soloClusterCBfile` is required. A single cell does not have enough UMIs to resolve isoforms, so the EM would be fitting noise. Asking for the feature without a clustering is refused rather than run. A UMI seen on several reads contributes the *intersection* of their transcript sets. Those reads came from one molecule, so a transcript missing from any of them cannot be its source. Taking the union would let a single stray read resurrect an isoform every other read excluded. Two of STAR's quirks are reproduced rather than corrected, because the cut point and every weight depend on them: the running-average divisor is `min(2N+1, i + N)` rather than the number of elements actually summed, and the transcript length factor is taken from the cumulative distribution at `trLen - 1` (`SoloFeature_quantTranscript.cpp`). Numbers are formatted the way C++'s default stream prints them — six significant digits, fixed inside `[1e-4, 1e6)` and scientific outside — since the normalised distribution runs down to ~1e-4 where Rust's `{}` and C++'s default disagree on both notation and digit count. --- CHANGELOG.md | 24 ++ src/params/mod.rs | 21 +- src/solo/count.rs | 41 +++ src/solo/mod.rs | 40 +++ src/solo/transcript3p.rs | 598 +++++++++++++++++++++++++++++++++++++++ src/solo/whitelist.rs | 14 + 6 files changed, 735 insertions(+), 3 deletions(-) create mode 100644 src/solo/transcript3p.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index b6e4893..aed5608 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -97,6 +97,30 @@ Sections commonly used: Features, Bug fixes, Other changes. glibc's malloc and per-thread heaps that return whole segments to the OS when abandoned, so allocator cache size stays bounded. +- **`--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 - **STARsolo `Gene` assignment now requires exon concordance**, matching diff --git a/src/params/mod.rs b/src/params/mod.rs index 0536a85..8c1a415 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1113,6 +1113,11 @@ pub struct Parameters { #[arg(long = "soloCBmatchWLtype", default_value = "1MM_multi")] pub solo_cb_match_wl_type: String, + /// 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()])] pub solo_cell_filter: Vec, @@ -1615,15 +1620,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" ), )); } @@ -1710,6 +1715,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(), diff --git a/src/solo/count.rs b/src/solo/count.rs index 9af3edf..b72487b 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -1323,6 +1323,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 diff --git a/src/solo/mod.rs b/src/solo/mod.rs index 2253dee..ef49a6a 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -12,6 +12,7 @@ pub mod count; pub mod gene; pub mod smartseq; +pub mod transcript3p; pub mod whitelist; pub use count::{UmiDedup, UmiFiltering, write_gene_matrix}; @@ -564,6 +565,11 @@ pub struct SoloContext { /// `--soloMultiMappers` includes a non-`Unique` method → capture gene- /// ambiguous reads for distribution into `UniqueAndMult-*.mtx`. pub want_multi: bool, + /// `--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 @@ -677,6 +683,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 { @@ -695,6 +702,20 @@ impl SoloContext { velocyto_enabled, velocyto_records: Mutex::new(Vec::new()), want_multi, + 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())), }) } @@ -871,6 +892,25 @@ impl SoloContext { fo }) .collect(); + + // 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/transcript3p.rs b/src/solo/transcript3p.rs new file mode 100644 index 0000000..bd63349 --- /dev/null +++ b/src/solo/transcript3p.rs @@ -0,0 +1,598 @@ +//! `--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`. +//! +//! 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)); + } + + pub fn merge(&mut self, other: Self) { + for (a, b) in self.dist_count.iter_mut().zip(&other.dist_count) { + *a += b; + } + self.records.extend(other.records); + } +} + +/// 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 { From 9d1f7871e484ee1601e0c8425ac4fcb0db5137b3 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 11:31:18 +0200 Subject: [PATCH 2/4] docs(changelog): keep only this PR's entry CONTRIBUTING.md requires the description to match the code; the entries for the other themes split out of #152 belong to their own PRs. --- CHANGELOG.md | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index aed5608..b04d78c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -97,16 +97,6 @@ Sections commonly used: Features, Bug fixes, Other changes. glibc's malloc and per-thread heaps that return whole segments to the OS when abandoned, so allocator cache size stays bounded. -- **`--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 From 295968acf18d5a4b5b209a85b45be55bd287aed6 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 11:53:49 +0200 Subject: [PATCH 3/4] refactor(solo): drop Transcript3pAcc::merge, which nothing calls Records are accumulated under a mutex, so there are no partials to merge. It was dead from the moment it was written; CONTRIBUTING.md rules out shipping it. --- src/solo/transcript3p.rs | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/solo/transcript3p.rs b/src/solo/transcript3p.rs index bd63349..14e475a 100644 --- a/src/solo/transcript3p.rs +++ b/src/solo/transcript3p.rs @@ -56,13 +56,6 @@ impl Transcript3pAcc { } self.records.push((cb, umi, hits)); } - - pub fn merge(&mut self, other: Self) { - for (a, b) in self.dist_count.iter_mut().zip(&other.dist_count) { - *a += b; - } - self.records.extend(other.records); - } } /// Every transcript this alignment is concordant with, and the spliced distance From bda9a345111962626918ffb9ab3abb14152a4577 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 21:43:33 +0200 Subject: [PATCH 4/4] docs(solo): note that STAR marks Transcript3p under development parametersDefault puts both Transcript3p and --soloClusterCBfile between "#####UnderDevelopment_begin : not supported - do not use" and "#####UnderDevelopment_end", and STAR --help prints that banner around them. The module said none of this. It matters for how the port is read: it follows STAR's code, so it inherits the unfinished parts of that code, and a differential against STAR compares two implementations of something STAR does not support. A reviewer should be told that before deciding to take it. --- src/solo/transcript3p.rs | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/solo/transcript3p.rs b/src/solo/transcript3p.rs index 14e475a..73b4fd2 100644 --- a/src/solo/transcript3p.rs +++ b/src/solo/transcript3p.rs @@ -3,6 +3,15 @@ //! //! 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