diff --git a/CHANGELOG.md b/CHANGELOG.md index b6e4893..5d7fcd2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -97,6 +97,17 @@ 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. +- **`--soloCellReadStats CB`** writes `Solo.out//CellReads.stats`: + one row per cell barcode with fifteen counters describing what + happened to its reads — barcode match quality, unique or multi + genomic mapping, feature assignment, exonic/intronic and their + antisense counterparts, mitochondrial, and whether the read reached + the matrix — plus the per-cell UMI and gene totals. Reads whose + barcode never resolved are summed into a `CBnotInPasslist` row rather + than dropped, so the columns account for the whole input. + `--genomeChrSetMitochondrial` names the chromosomes behind the `mito` + column. + ### Bug fixes - **STARsolo `Gene` assignment now requires exon concordance**, matching diff --git a/DIVERGENCE.md b/DIVERGENCE.md index bd957ed..f1d9f98 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -65,6 +65,20 @@ On the 10k yeast PE benchmark, 4 reads differ in alignment score (AS) because ST --- +### 3.2 `CellReads.stats` row order + +**What STAR does.** `--soloCellReadStats CB` emits its rows by iterating a libc++ `std::unordered_map`, so the order is a hash-table walk rather than a sort. At the map sizes this produces, libc++ chains new entries at the head of their bucket and walks buckets in order, which comes out as the reverse of each barcode's first appearance in read order. + +**What rustar-aligner does.** Emits that same reverse-first-appearance order, including across threads: the per-read accumulator merges in read order, so a threaded run writes the same file as a serial one. + +**Why.** Reproducing the order where it is reproducible costs nothing and keeps a byte-comparison against STAR usable on the sizes where it can work at all. + +**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`. + +--- + ## 4. Implementation divergences (no intended output difference) These differ in *how* a result is produced, not *what* is produced. They are documented so a reviewer chasing a discrepancy knows the mechanism differs by design. diff --git a/src/params/mod.rs b/src/params/mod.rs index 0536a85..1b7bf48 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1113,6 +1113,17 @@ pub struct Parameters { #[arg(long = "soloCBmatchWLtype", default_value = "1MM_multi")] pub solo_cb_match_wl_type: String, + /// `CB`: write `Solo.out//CellReads.stats`, a per-cell-barcode + /// summary of what happened to the reads carrying it. `None` (the default) + /// writes nothing. + #[arg(long = "soloCellReadStats", default_value = "None")] + pub solo_cell_read_stats: String, + + /// Chromosome names treated as mitochondrial, for the `mito` column of + /// `CellReads.stats`. `-` (the default) names none. + #[arg(long = "genomeChrSetMitochondrial", num_args = 1.., default_values_t = vec!["-".to_string()])] + pub genome_chr_set_mitochondrial: Vec, + /// 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, @@ -1710,6 +1721,16 @@ impl Parameters { )); } } + // --soloCellReadStats: `CB` is the only value STAR defines. + if !matches!(params.solo_cell_read_stats.as_str(), "CB" | "None") { + return Err(command.error( + ErrorKind::InvalidValue, + format!( + "unknown --soloCellReadStats '{}'; expected CB or None", + params.solo_cell_read_stats + ), + )); + } // Validate --clipAdapterType. if !matches!( params.clip_adapter_type.as_str(), diff --git a/src/solo/cell_reads.rs b/src/solo/cell_reads.rs new file mode 100644 index 0000000..5d140a5 --- /dev/null +++ b/src/solo/cell_reads.rs @@ -0,0 +1,288 @@ +//! `--soloCellReadStats CB`: the per-cell-barcode read summary STARsolo writes +//! as `Solo.out//CellReads.stats`. +//! +//! One row per cell barcode, fifteen counters describing what happened to the +//! reads carrying it — how the barcode matched, whether the read mapped +//! uniquely, whether it landed on a feature, where in the gene, and whether it +//! reached the matrix — plus the per-cell UMI and gene totals. The reads whose +//! barcode never resolved are not dropped; they are summed into a single +//! `CBnotInPasslist` row, so the columns account for every read rather than +//! only the ones that succeeded. +//! +//! # D24: row order +//! +//! STAR iterates a libc++ `std::unordered_map` to emit these rows, so the order +//! is a hash-table walk, not a sort. For the small maps this produces, libc++ +//! chains new entries at the head of their bucket and walks buckets in order, +//! which comes out as the reverse of first appearance in read order. That is +//! what this reproduces. +//! +//! It is not reproducible in general: past the load factor libc++ rehashes, and +//! the order after a rehash depends on the bucket count, which depends on how +//! many distinct barcodes were seen. At that size the order diverges. The +//! **values never do** — only which line they appear on. A consumer that reads +//! this file by barcode rather than by position is unaffected either way, and +//! sorting by barcode is the only stable thing to do with it. +//! +//! Reads are folded in under a mutex held by `SoloContext`, so the order is the +//! order reads were processed in regardless of thread count. + +use std::collections::{BTreeMap, HashSet}; +use std::fmt::Write as _; + +/// What happened to one read, as the fourteen optional flags STAR tracks. +/// +/// `cbMatch` is not here: every read that reaches the accumulator matched a +/// barcode well enough to be attributed somewhere, so it is always set. +/// +/// Fourteen bools rather than a bitfield because they are written once per read +/// and read once per fold, and the names are the column names of the file. +#[allow(clippy::struct_excessive_bools)] +#[derive(Debug, Clone, Copy, Default)] +pub struct CellReadFlag { + /// The barcode was an exact whitelist hit. + pub cb_perfect: bool, + /// Corrected via a single one-mismatch neighbour. + pub cb_mm_unique: bool, + /// Corrected via several one-mismatch neighbours, resolved by the posterior. + pub cb_mm_multiple: bool, + /// The read mapped to exactly one genomic locus. + pub genome_u: bool, + /// The read mapped to several genomic loci. + pub genome_m: bool, + /// It fell on exactly one feature (gene). + pub feature_u: bool, + /// It fell on several features. + pub feature_m: bool, + /// Exonic, on the annotated strand. + pub exonic: bool, + /// Intronic, on the annotated strand. + pub intronic: bool, + /// Exonic, antisense to the annotation. + pub exonic_as: bool, + /// Intronic, antisense to the annotation. + pub intronic_as: bool, + /// On a chromosome named by `--genomeChrSetMitochondrial`. + pub mito: bool, + /// Counted into the unique-gene matrix. + pub counted_u: bool, + /// Counted through the multi-gene distribution. + pub counted_m: bool, +} + +/// The accumulator behind `CellReads.stats`. +#[derive(Debug, Default, Clone)] +pub struct CellReadStats { + /// Whitelist index to its fifteen counters. + cells: BTreeMap, + /// The single bucket for reads whose barcode did not resolve. + no_cb: [u64; 15], + /// Whitelist indices in first-appearance order; emitted reversed (D24). + order: Vec, + seen: HashSet, +} + +impl CellReadStats { + pub fn new() -> Self { + Self::default() + } + + fn fold(v: &mut [u64; 15], f: &CellReadFlag) { + v[0] += 1; // cbMatch: every read that got this far + for (i, set) in [ + f.cb_perfect, + f.cb_mm_unique, + f.cb_mm_multiple, + f.genome_u, + f.genome_m, + f.feature_u, + f.feature_m, + f.exonic, + f.intronic, + f.exonic_as, + f.intronic_as, + f.mito, + f.counted_u, + f.counted_m, + ] + .into_iter() + .enumerate() + { + if set { + v[i + 1] += 1; + } + } + } + + /// Record a read that resolved to whitelist cell `cb`. + pub fn add_cell(&mut self, cb: u32, flag: &CellReadFlag) { + if self.seen.insert(cb) { + self.order.push(cb); + } + Self::fold(self.cells.entry(cb).or_insert([0; 15]), flag); + } + + /// Record a read whose barcode did not resolve, or whose UMI was rejected. + pub fn add_no_cb(&mut self, flag: &CellReadFlag) { + Self::fold(&mut self.no_cb, flag); + } + + /// Render the file. `umi_gene` gives each cell its final + /// `(nUMIunique, nGenesUnique)`; a cell absent from it prints zeros. + /// `barcode_of` renders a whitelist index as its barcode string. + pub fn render( + &self, + barcode_of: impl Fn(u32) -> String, + umi_gene: &BTreeMap, + ) -> String { + let mut s = String::from( + "CB\tcbMatch\tcbPerfect\tcbMMunique\tcbMMmultiple\tgenomeU\tgenomeM\tfeatureU\t\ + featureM\texonic\tintronic\texonicAS\tintronicAS\tmito\tcountedU\tcountedM\t\ + nUMIunique\tnGenesUnique\tnUMImulti\tnGenesMulti\n", + ); + s.push_str("CBnotInPasslist"); + for v in &self.no_cb { + s.push('\t'); + s.push_str(&v.to_string()); + } + s.push_str("\t0\t0\t0\t0\n"); + // Reverse first-appearance order — see D24 in the module docs. + for &cb in self.order.iter().rev() { + s.push_str(&barcode_of(cb)); + for v in &self.cells[&cb] { + s.push('\t'); + s.push_str(&v.to_string()); + } + let (n_umi, n_gene) = umi_gene.get(&cb).copied().unwrap_or((0, 0)); + // nUMImulti and nGenesMulti stay zero: multi-gene UMIs are not + // collapsed into per-cell totals here. + let _ = writeln!(s, "\t{n_umi}\t{n_gene}\t0\t0"); + } + s + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn flag_perfect_counted() -> CellReadFlag { + CellReadFlag { + cb_perfect: true, + genome_u: true, + feature_u: true, + exonic: true, + counted_u: true, + ..Default::default() + } + } + + fn header_and_rows(text: &str) -> (Vec<&str>, Vec<&str>) { + let mut lines = text.lines(); + let header: Vec<&str> = lines.next().unwrap().split('\t').collect(); + (header, lines.collect()) + } + + /// Every column has a name, and every row has a value under each of them. + #[test] + fn each_row_fills_the_header() { + let mut st = CellReadStats::new(); + st.add_cell(0, &flag_perfect_counted()); + st.add_no_cb(&CellReadFlag::default()); + let out = st.render(|i| format!("CB{i}"), &BTreeMap::new()); + let (header, rows) = header_and_rows(&out); + assert_eq!(header.len(), 20); + assert_eq!(rows.len(), 2, "the passlist-miss row plus one cell"); + for row in rows { + assert_eq!(row.split('\t').count(), header.len()); + } + } + + /// A read is counted once under `cbMatch` and once under each flag it sets, + /// so the flag columns can exceed neither `cbMatch` nor each other's logic. + #[test] + fn flags_accumulate_per_read() { + let mut st = CellReadStats::new(); + for _ in 0..3 { + st.add_cell(7, &flag_perfect_counted()); + } + st.add_cell( + 7, + &CellReadFlag { + cb_mm_unique: true, + genome_m: true, + ..Default::default() + }, + ); + let out = st.render(|i| format!("CB{i}"), &BTreeMap::new()); + let row: Vec<&str> = out.lines().nth(2).unwrap().split('\t').collect(); + assert_eq!(row[0], "CB7"); + assert_eq!(row[1], "4", "cbMatch counts every read"); + assert_eq!(row[2], "3", "cbPerfect"); + assert_eq!(row[3], "1", "cbMMunique"); + assert_eq!(row[5], "3", "genomeU"); + assert_eq!(row[6], "1", "genomeM"); + assert_eq!(row[14], "3", "countedU"); + } + + /// Reads whose barcode never resolved are summed rather than dropped, so + /// the file accounts for the whole input. + #[test] + fn unresolved_reads_land_in_the_passlist_miss_row() { + let mut st = CellReadStats::new(); + st.add_cell(0, &flag_perfect_counted()); + for _ in 0..5 { + st.add_no_cb(&CellReadFlag { + genome_u: true, + ..Default::default() + }); + } + let out = st.render(|i| format!("CB{i}"), &BTreeMap::new()); + let row: Vec<&str> = out.lines().nth(1).unwrap().split('\t').collect(); + assert_eq!(row[0], "CBnotInPasslist"); + assert_eq!(row[1], "5"); + assert_eq!(row[5], "5", "genomeU"); + assert_eq!(&row[16..], ["0", "0", "0", "0"], "no UMI columns for it"); + } + + /// D24: rows come out in reverse first-appearance order, which is what + /// STAR's libc++ hash-map walk produces at these sizes. + #[test] + fn rows_are_emitted_in_reverse_first_appearance_order() { + let mut st = CellReadStats::new(); + for cb in [4u32, 1, 9] { + st.add_cell(cb, &flag_perfect_counted()); + } + st.add_cell(1, &flag_perfect_counted()); // seen again: order unchanged + let out = st.render(|i| format!("CB{i}"), &BTreeMap::new()); + let cbs: Vec<&str> = out + .lines() + .skip(2) + .map(|l| l.split('\t').next().unwrap()) + .collect(); + assert_eq!(cbs, ["CB9", "CB1", "CB4"]); + } + + /// The UMI and gene totals come from the final matrix, not from the read + /// counters, so a cell missing from that map prints zeros rather than + /// inheriting a read count. + #[test] + fn umi_and_gene_totals_come_from_the_matrix() { + let mut st = CellReadStats::new(); + st.add_cell(2, &flag_perfect_counted()); + st.add_cell(3, &flag_perfect_counted()); + let mut umi_gene = BTreeMap::new(); + umi_gene.insert(2u32, (17u32, 5u32)); + let out = st.render(|i| format!("CB{i}"), &umi_gene); + let rows: Vec> = out + .lines() + .skip(2) + .map(|l| l.split('\t').collect()) + .collect(); + // Reverse order: CB3 first, then CB2. + assert_eq!(rows[0][0], "CB3"); + assert_eq!(&rows[0][16..18], ["0", "0"]); + assert_eq!(rows[1][0], "CB2"); + assert_eq!(&rows[1][16..18], ["17", "5"]); + } +} diff --git a/src/solo/count.rs b/src/solo/count.rs index 9af3edf..cd45309 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -1224,6 +1224,27 @@ pub fn write_gene_matrix( if gzip { " [gzip]" } else { "" }, ); + // `--soloCellReadStats CB`: the per-cell read summary, alongside the + // raw matrix because its UMI and gene columns are the raw totals. + if let Some(cell_stats) = &ctx.cell_read_stats { + let umi_gene: std::collections::BTreeMap = mstats + .cells + .iter() + .map(|c| (c.cb, (c.n_umis as u32, c.n_genes))) + .collect(); + let path = feature_dir.join("CellReads.stats"); + let text = cell_stats.lock().unwrap().render( + |cb| { + ctx.whitelist + .barcode_string(cb) + .unwrap_or_else(|| cb.to_string()) + }, + &umi_gene, + ); + std::fs::write(&path, text).map_err(|e| Error::io(e, &path))?; + log::info!("STARsolo: wrote {}/CellReads.stats", feature.dir_name()); + } + // Filtered (cell-called) matrix per --soloCellFilter. EmptyDrops_CR runs // the Monte-Carlo rescue (needs the per-cell profiles in the body). let called = if params diff --git a/src/solo/mod.rs b/src/solo/mod.rs index 2253dee..e53f4f4 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -9,6 +9,7 @@ //! The barcode read is the SECOND `--readFilesIn` file (STAR convention: //! `--readFilesIn cDNA_read barcode_read`). It is never aligned — only parsed. +pub mod cell_reads; pub mod count; pub mod gene; pub mod smartseq; @@ -564,6 +565,13 @@ pub struct SoloContext { /// `--soloMultiMappers` includes a non-`Unique` method → capture gene- /// ambiguous reads for distribution into `UniqueAndMult-*.mtx`. pub want_multi: bool, + /// `--soloCellReadStats CB`: the per-cell read summary, or `None` when the + /// flag is off. Behind a mutex like the other per-read collections; the + /// lock is only taken when the flag asks for the file. + pub cell_read_stats: Option>, + /// Chromosome indices named by `--genomeChrSetMitochondrial`, for the + /// `mito` column. + pub mito_chr: std::collections::HashSet, } /// Per-region read tallies for the `Summary.csv` mapping funnel (uniquely-mapped @@ -695,6 +703,14 @@ impl SoloContext { velocyto_enabled, velocyto_records: Mutex::new(Vec::new()), want_multi, + cell_read_stats: (params.solo_cell_read_stats == "CB") + .then(|| Mutex::new(crate::solo::cell_reads::CellReadStats::new())), + mito_chr: params + .genome_chr_set_mitochondrial + .iter() + .filter(|n| n.as_str() != "-") + .filter_map(|n| genome.chr_name.iter().position(|c| c == n)) + .collect(), }) } @@ -871,9 +887,67 @@ impl SoloContext { fo }) .collect(); + + self.record_cell_read( + cb_resolved, + &cb_match, + n_loci, + &class, + cdna_transcripts, + &out, + ); out } + /// Fold one read into `CellReads.stats`, when `--soloCellReadStats CB` + /// asked for it. + /// + /// Called at the end of read processing so the counted flags reflect what + /// the read actually produced, rather than what it looked eligible for. + #[allow(clippy::too_many_arguments)] + fn record_cell_read( + &self, + cb_resolved: Option, + cb_match: &CbMatch, + n_loci: usize, + class: &crate::solo::gene::ReadClass, + transcripts: &[Transcript], + out: &SoloReadOutcome, + ) { + let Some(stats) = &self.cell_read_stats else { + return; + }; + let counted_u = out.per_feature.iter().any(|f| f.record.is_some()); + let counted_m = out.per_feature.iter().any(|f| f.multi_gene.is_some()); + let feature_u = counted_u || out.per_feature.iter().any(|f| f.multi.is_some()); + let flag = crate::solo::cell_reads::CellReadFlag { + cb_perfect: matches!(cb_match, CbMatch::Exact(_)), + cb_mm_unique: matches!(cb_match, CbMatch::Corrected(_)), + cb_mm_multiple: matches!(cb_match, CbMatch::Multi(_)), + genome_u: n_loci == 1, + genome_m: n_loci > 1, + feature_u, + feature_m: counted_m, + // The region columns split by strand: STAR reports an antisense + // read under `exonicAS`/`intronicAS`, not under `exonic`/`intronic`. + exonic: !class.antisense && class.region == Some(Region::Exonic), + intronic: !class.antisense && class.region == Some(Region::Intronic), + exonic_as: class.antisense && class.region == Some(Region::Exonic), + intronic_as: class.antisense && class.region == Some(Region::Intronic), + mito: !self.mito_chr.is_empty() + && transcripts + .iter() + .any(|t| self.mito_chr.contains(&t.chr_idx)), + counted_u, + counted_m, + }; + let mut stats = stats.lock().unwrap(); + match cb_resolved { + Some(cb) => stats.add_cell(cb, &flag), + None => stats.add_no_cb(&flag), + } + } + /// Process one 5' paired-end solo read (`--soloBarcodeMate 1`): the barcode is /// from mate 1, and both mates align as a pair. Genes are assigned from the /// union of both mates evaluated against the pair's (mate 1's) transcription