diff --git a/CHANGELOG.md b/CHANGELOG.md index 200160f..d15db6e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -213,6 +213,86 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Bug fixes +- **STARsolo's stats files as STAR writes them.** `Solo.out/Barcodes.stats`, + `/Features.stats` and `/UMIperCellSorted.txt` were not + written, SJ and Velocyto had no `Summary.csv`, and the Gene/GeneFull + `Summary.csv` lacked the Q30 and "Unique+Multiple" lines, printed fixed + decimals, and computed several values its own way. They now come from + STAR's counters (`SoloReadBarcodeStats`, `SoloReadFeatureStats`) and its + `statsOutput`: + - each read counts once in `Barcodes.stats`, under its final `cbMatch` + code (a UMI rejection replaces the barcode's); a missing adapter or UMI + is counted as STAR counts it; + - the per-feature counts come from STAR's three stages: per read + (`record`), at counting (oneExact and 1MM_multi rejections) and at UMI + collapse; Velocyto keeps only the UMIs and cells, as `countVelocyto` + never sums the rest; + - "Reads With Valid Barcodes" subtracts both barcode and feature + rejections; "Reads Mapped to : Unique" is the counted reads; + means are integer, medians the element at `nCells/2`, the UMI median + from every barcode's UMIs; + - doubles, here and in `UniqueAndMult-*.mtx` (which printed `0.50000` + for STAR's `0.5`), are written as C++'s `ostream` does; + - `CB_samTagOut` writes `Barcodes.stats`, and so does a run with + `--soloCBwhitelist None` (all zeros, as STAR's: it bins reads only + against a whitelist). That mode still writes no matrices. + + Also: a barcode's exact-match count (used by oneExact and the 1MM_multi + posterior) now counts only reads whose UMI passed too, as STAR's + `addStats` does; it counted exact barcodes with a rejected UMI. Across 30 + STARsolo configurations of the 10x v2 mouse chr19 fixture, 453 of 457 + `Solo.out` files are byte-identical to STARsolo's; the other four are the + `CB_UMI_Complex` homopolymer UMIs of DIVERGENCE.md §1.6. + + STAR's per-line SJ counts are ported too: a junction that is not a matrix + row counts nowhere, a read with no junction is `noNoFeature` (a + 1MM_multi barcode with junctions was too), and 1MM_multi barcodes on SJ + records, which were dropped, are resolved by the posterior as for Gene. `CB_UMI_Complex` reads rejected + twice can land in a different `Barcodes.stats` bin from STAR's + (DIVERGENCE.md §3.4). + + `CellRanger.summary.csv`, a rustar-aligner extension written whenever + Gene and GeneFull both ran, is now opt-in: `--soloCellRangerSummary Yes` + (DIVERGENCE.md §1.8). By default `Solo.out` holds only STARsolo's files. + +- **`--soloCellFilter TopCells N` calls cells as STARsolo.** It took the N + barcodes with the most UMIs. STAR keeps every barcode with at least the + UMIs of the barcode at index N of the descending list, the (N+1)th, so + ties are kept (`SoloFeature_cellFiltering.cpp:46-59`). With `TopCells 100` + on the 10x v2 mouse chr19 fixture STAR calls 106 cells and rustar-aligner + called 100; the filtered matrices are now identical. + +- **`--soloFeatures Velocyto` counts as STARsolo.** It used its own rule (a + junction means spliced, an intronic block unspliced, anything else + ambiguous), read by read. It now ports STAR's: + - each uniquely mapped read is compared with every transcript that + contains it on the `--soloStrand` strand (`alignToTranscriptMinOverlap`, + velocyto's 5-base flank, already ported for bulk use in + `quant/splice_status.rs`); + - a molecule (cell, UMI, taken from the Gene feature's readInfo) keeps the + transcripts all its reads share; it is spliced, unspliced or ambiguous by + STAR's model flags, and dropped if its transcripts span two genes + (`SoloFeature_countVelocyto.cpp`); + - all three matrices hold the same (cell, gene) rows, zeros included; + - `Velocyto/filtered/` is written, for the Gene feature's called cells that + Velocyto detected (`SoloFeature_cellFiltering.cpp:16-30`), unless there + are none (`:106-108`). + + On the 10x v2 mouse chr19 fixture the raw and filtered matrices are now + byte-identical to STARsolo's (spliced/unspliced/ambiguous UMIs 369/57/32, + against 52/100/362 before), in 10 configurations: Unstranded, Reverse, + `1MM_multi`, GeneFull listed first, with `CB`/`UB` tags, `1MM_Directional` + and `Exact` dedup, EmptyDrops_CR, and paired-end cDNA. STAR keeps a readInfo + per feature (`CB`/`UB` read the first feature's, Velocyto the Gene + feature's), and so does rustar-aligner now. `--soloFeatures Velocyto` + without `Gene`, on which STAR crashes, is refused; with `--soloType + SmartSeq` it is refused with STAR's message. + + `--soloVelocytoAmbiguous` is removed. It was a rustar-aligner extension + that folded ambiguous molecules into spliced, on the premise that exon-only + reads land in ambiguous; under STAR's model they are spliced, and ambiguous + means exonic in one isoform and intronic in another. + - **Paired-end reads with mates of 50 bases or fewer align as STAR.** STAR caps the spacing of its seed search starts by the whole read, both mates and the spacer (`ReadAlign_mapOneRead.cpp:41-42`); we capped each mate by diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 69bd33e..2b9ce43 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -130,7 +130,7 @@ So when `umiL` is still 0 at construction, the loop body never runs and all four Poly-A being special in RNA does not apply here: the UMI comes from the barcode read, not the cDNA. In oligo-dT chemistries the UMI sits next to the dT, so a UMI that runs into it reads as poly-T, the one homopolymer STAR's complex path lets through. Matching the accident would count poly-T UMIs as molecules. -**Impact.** Under `--soloType CB_UMI_Complex`, reads whose UMI is a non-A homopolymer are dropped here and counted by STAR. On real data these are rare and are sequencing artefacts. +**Impact.** Under `--soloType CB_UMI_Complex`, reads whose UMI is a non-A homopolymer are dropped here and counted by STAR. On real data these are rare and are sequencing artefacts. The stats files show it too: such reads count as `noUMIhomopolymer` in `Barcodes.stats` here, and under their barcode's match (or `noNoWLmatch`) in STAR's, which also moves `Features.stats` and the "Reads With Valid Barcodes" line of `Summary.csv`. On the 10x v2 mouse chr19 fixture run as `CB_UMI_Complex`, 43 reads; the count matrices are unchanged. **Source.** `src/solo/whitelist.rs` (`check_umi`, `umi_valid_rejects_every_homopolymer`). STAR: `SoloReadBarcode.cpp:16-21`, `SoloReadBarcode_getCBandUMI.cpp:139` and `:353-354`. @@ -159,6 +159,18 @@ Every STAR record matches the alignment STAR computed only once SEQ is moved by **Source.** `src/clip/mod.rs` (`clip_mate_cr4`). STAR: `ClipMate_clip.cpp:14-24, 50-57, 64-70`, `readLoad.cpp:60-61`, `ReadAlign_alignBAM.cpp:224-230`, `ReadAlign_calcCIGAR.cpp:16-22`. +### 1.8 `--soloCellRangerSummary`: a CellRanger-style mapping summary + +A rustar-aligner extension with no STAR equivalent, off by default (`No`). With the default no output changes. + +**What STAR does.** STARsolo writes one `Summary.csv` per feature (`SoloFeature_statsOutput.cpp`), with reads mapped to the genome and to the feature, but no split of the reads by region. + +**What rustar-aligner does.** With `--soloCellRangerSummary Yes` and both `Gene` and `GeneFull` in `--soloFeatures`, it also writes `Solo.out//CellRanger.summary.csv`: of all input reads, the fractions that map uniquely and fall in exonic, intronic or intergenic regions, or antisense to a gene, under the names CellRanger's `metrics_summary.csv` uses ("Reads Mapped Confidently to Exonic Regions" and so on). The definitions are rustar-aligner's own (exonic and intronic come from the two features' overlaps) and are not claimed to equal CellRanger's. `Summary.csv` itself is never changed. + +**Why.** CellRanger users look at these fractions, and the two overlaps needed are computed anyway when both features run. + +**Source.** `src/solo/count.rs` (`write_cellranger_summary`), `src/solo/gene.rs` (`Region`). + --- ## 2. Cases where rustar-aligner outperforms STAR @@ -248,6 +260,22 @@ On the 10k yeast PE benchmark, 4 reads differ in alignment score (AS) because ST --- +### 3.4 `CB_UMI_Complex`: the `Barcodes.stats` bin of two kinds of rejected read + +**What STAR does.** Under `CB_UMI_Complex`, STAR checks the UMI first and then matches each cell-barcode segment to its own whitelist, in order (`SoloReadBarcode_getCBandUMI.cpp:356-420`). The first failure sets `cbMatch`, which picks the read's `Barcodes.stats` bin. Two consequences: +- a read with one mismatch in each of two segments is `cbMatch` -12, counted as `noTooManyMM`; +- a read with a rejected UMI keeps the UMI's code (`noNinUMI` or `noUMIhomopolymer`) even when a segment cannot be extracted. + +**What rustar-aligner does.** It matches the concatenated segments against the product of the segment whitelists. A read with mismatches in two segments is two mismatches from every entry, so it is counted as `noNoWLmatch`. A barcode read too short for a segment is counted as `noNoCB` before its UMI is checked. + +**Why.** Not chosen: per-segment matching is not ported yet. The product whitelist gives the same match as STAR for a read with at most one mismatched segment, which is every read STAR counts. + +**Impact.** Only the bins of reads that both tools reject: the count matrices, `Features.stats` and `Summary.csv` are unchanged. On the 10x v2 mouse chr19 fixture run as `CB_UMI_Complex`, `Barcodes.stats` differs from STAR's only by the reads of §1.6. + +**Source.** `src/solo/mod.rs` (`SoloBarcodeLayout::extract`), `src/solo/whitelist.rs` (`CbWhitelist::load_complex`). STAR: `SoloReadBarcode_getCBandUMI.cpp:356-420`. + +--- + ## 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/ROADMAP.md b/ROADMAP.md index 672153a..690ee01 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -312,7 +312,7 @@ See [docs-old/phase17_features.md](docs-old/phase17_features.md) for sub-phase t **Prerequisite met**: position agreement >99% (SE 99.815% tie-adj, PE 99.883%). Phase unblocked 2026-06-10. -**Status**: A STARsolo-faithful single-cell pipeline, Gene count matrix **byte-identical to STARsolo's**, with GeneFull/SJ/Velocyto features, CB_UMI_Simple/Complex/SmartSeq chemistries, EmptyDrops_CR + CellRanger2.2 cell calling, multi-mapper resolution, `Summary.csv`, and the per-read STARsolo SAM tags. Validated in a native three-way benchmark against STARsolo and CellRanger (see below). Remaining: the Velocyto feature's counts differ from STARsolo's, and `Barcodes.stats`/`Features.stats` are not written. +**Status**: A STARsolo-faithful single-cell pipeline, Gene count matrix **byte-identical to STARsolo's**, with GeneFull/SJ/Velocyto features, CB_UMI_Simple/Complex/SmartSeq chemistries, EmptyDrops_CR + CellRanger2.2 cell calling, multi-mapper resolution, `Summary.csv`, and the per-read STARsolo SAM tags. Validated in a native three-way benchmark against STARsolo and CellRanger (see below). Single-cell quantification layered around the existing aligner: the cDNA read aligns through the normal SE path; a paired **barcode read** (R1 = cell barcode + UMI) is parsed, corrected against a whitelist, assigned to a gene, UMI-deduplicated, and emitted as a sparse per-cell count matrix. Target: faithful port of STARsolo (all features). See [docs-old/phase14_starsolo.md](docs-old/phase14_starsolo.md) for the full design and sub-phase tracking. diff --git a/docs/src/content/docs/reference/cli-parameters.md b/docs/src/content/docs/reference/cli-parameters.md index 0ea6005..ed87d8c 100644 --- a/docs/src/content/docs/reference/cli-parameters.md +++ b/docs/src/content/docs/reference/cli-parameters.md @@ -162,6 +162,7 @@ Run `rustar-aligner --help` for the full machine-generated listing. | `--quantGeneCoverage` | `No` | rustar-aligner extension. `Yes` writes per-gene exonic coverage to `GeneCoverage.out.bin` for DegNorm. See [DegNorm](/rustar-aligner/guides/degnorm/). | | `--degNormSampleId` | basename of `--outFileNamePrefix` | rustar-aligner extension. Sample name stored in `GeneCoverage.out.bin`. | | `--degNormCoverageFiles` | — | rustar-aligner extension. Two or more `GeneCoverage.out.bin` files: runs the DegNorm merge (no reads, no `--runMode`) and writes `DegNorm.out/`. Tuning: `--degNormIter`, `--degNormNmfIter`, `--degNormDownsampleRate`, `--degNormMinimaxCoverage`, `--degNormSkipBaselineSelection`, `--degNormBins`, `--degNormMinHighCoverage` (see the guide). | +| `--soloCellRangerSummary` | `No` | rustar-aligner extension. `Yes` writes `Solo.out//CellRanger.summary.csv` for `Gene` and `GeneFull` (both needed): the fractions of all reads that map uniquely in exonic, intronic or intergenic regions, or antisense to a gene, under CellRanger's metric names. STARsolo's `Summary.csv` is unchanged. | ## Two-pass mode diff --git a/src/lib.rs b/src/lib.rs index ecd87e8..1d0c993 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -616,18 +616,15 @@ fn write_solo_output( sj_stats: &std::sync::Arc, index: &std::sync::Arc, ) -> anyhow::Result<()> { - use std::sync::atomic::Ordering; - let s = &sctx.stats; + let barcode_counts = sctx.stats.values(); info!( - "STARsolo barcode stats: exact={} 1MM={} multiMM={} noMatch={} N-in-CB={} multReject={} N-in-UMI={} UMIhomopolymer={}", - s.yes_exact.load(Ordering::Relaxed), - s.yes_one_mm.load(Ordering::Relaxed), - s.yes_mult_mm.load(Ordering::Relaxed), - s.no_match.load(Ordering::Relaxed), - s.n_in_cb.load(Ordering::Relaxed), - s.mult_rejected.load(Ordering::Relaxed), - s.n_in_umi.load(Ordering::Relaxed), - s.umi_homopolymer.load(Ordering::Relaxed), + "STARsolo barcode stats: {}", + crate::solo::BarcodeStats::NAMES + .iter() + .zip(barcode_counts) + .map(|(n, v)| format!("{n}={v}")) + .collect::>() + .join(" ") ); for (feature, recorder) in sctx.features.iter().zip(&sctx.recorders) { info!( @@ -1048,14 +1045,12 @@ fn run_single_pass( // `CB_samTagOut` produces no Solo.out matrices at all (`Solo.cpp:13`). if params.solo_type != params::SoloType::CbSamTagOut { write_solo_output(sctx, params, &stats, &sj_stats, index)?; + } else { + // STAR still writes Barcodes.stats (`Solo::processAndOutput`). + crate::solo::write_barcode_stats(sctx, params)?; } - if sctx.read_info_enabled() { - let info = sctx - .read_info - .as_ref() - .expect("readInfo enabled") - .lock() - .unwrap(); + if let Some(info) = sctx.tag_read_info() { + let info = info.lock().unwrap(); let umi_len = params.solo_umi_len as usize; writer.finish_with(&mut |record| { crate::io::sam::apply_solo_read_info(record, &info, &sctx.whitelist, umi_len); @@ -2512,7 +2507,10 @@ fn align_reads_solo( } else { crate::params::SamAttributes::empty() }; - let track_read_info = emit_sam && solo_ctx.read_info_enabled(); + // Read indices for STAR's readInfo (the CB/UB tags, Velocyto), and the + // private tag that carries them to the sorted BAM's CB/UB pass. + let track_read_info = solo_ctx.read_info_enabled(); + let stamp_read_index = emit_sam && solo_ctx.tag_read_info().is_some(); // `--soloType CB_samTagOut`: barcodes become SAM tags, nothing is counted. let tag_out_only = params.solo_type == params::SoloType::CbSamTagOut; // Shared, 'static parameters for the per-batch aligner tasks spawned below. @@ -2522,7 +2520,7 @@ fn align_reads_solo( struct SoloReadProduct { sam_records: BufferedSamRecords, per_feature: Vec, - sj: Vec, + sj: crate::solo::SjRecords, velocyto: Option, } @@ -2570,10 +2568,13 @@ fn align_reads_solo( (0..n_feat).map(|_| Vec::new()).collect(); let mut feat_multi_gene: Vec> = (0..n_feat).map(|_| Vec::new()).collect(); - let mut sj_batch: Vec = Vec::new(); + let mut sj_batch = crate::solo::SjRecords::default(); let mut velo_batch: Vec = Vec::new(); - let mut no_feature: Vec = Vec::new(); - let mut no_feature_multi: Vec = Vec::new(); + // STAR's no-feature records, per feature (readInfo is per feature). + let mut no_feature: Vec> = + (0..n_feat).map(|_| Vec::new()).collect(); + let mut no_feature_multi: Vec> = + (0..n_feat).map(|_| Vec::new()).collect(); for result in products { let product = result?; writer.write_batch(&product.sam_records.records)?; @@ -2588,13 +2589,13 @@ fn align_reads_solo( feat_multi_gene[fi].push(mg); } if let Some(r) = fo.no_feature { - no_feature.push(r); + no_feature[fi].push(r); } if let Some(m) = fo.no_feature_multi { - no_feature_multi.push(m); + no_feature_multi[fi].push(m); } } - sj_batch.extend(product.sj); + sj_batch.append(product.sj); velo_batch.extend(product.velocyto); } for (fi, recorder) in solo.recorders.iter().enumerate() { @@ -2606,24 +2607,23 @@ fn align_reads_solo( if !mg.is_empty() { recorder.multi_gene.lock().unwrap().extend(mg); } - } - if let Some(recorder) = solo.recorders.first() { - if !no_feature.is_empty() { - recorder.no_feature.lock().unwrap().extend(no_feature); + let nf = std::mem::take(&mut no_feature[fi]); + if !nf.is_empty() { + recorder.no_feature.lock().unwrap().extend(nf); } - if !no_feature_multi.is_empty() { - recorder - .no_feature_multi - .lock() - .unwrap() - .extend(no_feature_multi); + let nfm = std::mem::take(&mut no_feature_multi[fi]); + if !nfm.is_empty() { + recorder.no_feature_multi.lock().unwrap().extend(nfm); } } if !sj_batch.is_empty() { - solo.sj_records.lock().unwrap().extend(sj_batch); + solo.sj_records.lock().unwrap().append(sj_batch); } if !velo_batch.is_empty() { - solo.velocyto_records.lock().unwrap().extend(velo_batch); + let mut store = solo.velocyto_records.lock().unwrap(); + for r in &velo_batch { + store.push(r); + } } Ok(true) }; @@ -2693,13 +2693,18 @@ fn align_reads_solo( ); let mut buffer = BufferedSamRecords::new(params.out_sam_attributes); stats.record_read_bases(clipped_seq.len() as u64); + solo.record_qualities(sread.barcode_q30, &[&clipped_qual]); if clipped_seq.is_empty() { stats.record_alignment(0, max_multimaps); stats.record_unmapped_reason(crate::stats::UnmappedReason::Other); // No alignment → barcode still counts toward stats (unmapped → no gene). - let outcome = - solo.process_read(&[], 0, sread.barcode.as_ref(), &[]); + let outcome = solo.process_read( + &[], + 0, + sread.barcode.as_ref().map_err(|m| *m), + &[], + ); return Ok(SoloReadProduct { sam_records: buffer, per_feature: outcome.per_feature, @@ -2746,18 +2751,23 @@ fn align_reads_solo( // counted. let (mut outcome, cb_corrected) = if tag_out_only { let mut outcome = crate::solo::SoloReadOutcome::default(); - let corrected = sread.barcode.as_ref().map(|bc| { - let (tags, corrected) = solo.tag_barcode(bc); - outcome.barcode = Some(tags); - corrected - }); + let corrected = sread + .barcode + .as_ref() + .map_err(|m| solo.record_barcode_miss(*m)) + .ok() + .map(|bc| { + let (tags, corrected) = solo.tag_barcode(bc); + outcome.barcode = Some(tags); + corrected + }); (outcome, corrected) } else { ( solo.process_read( &transcripts, transcripts.len(), - sread.barcode.as_ref(), + sread.barcode.as_ref().map_err(|m| *m), &junctions, ), None, @@ -2832,7 +2842,7 @@ fn align_reads_solo( // the read, mapped or unmapped, as in STAR. if !solo_tags.is_empty() { let values = crate::solo::SoloTagStrings::build( - sread.barcode.as_ref(), + sread.barcode.as_ref().ok(), sread.barcode_read.as_ref(), outcome.barcode, cb_corrected, @@ -2851,7 +2861,7 @@ fn align_reads_solo( ); } } - if track_read_info { + if stamp_read_index { crate::io::sam::add_solo_read_index( &mut buffer.records, read_index, @@ -2948,7 +2958,10 @@ fn align_reads_solo_pe( } else { crate::params::SamAttributes::empty() }; - let track_read_info = emit_sam && solo_ctx.read_info_enabled(); + // Read indices for STAR's readInfo (the CB/UB tags, Velocyto), and the + // private tag that carries them to the sorted BAM's CB/UB pass. + let track_read_info = solo_ctx.read_info_enabled(); + let stamp_read_index = emit_sam && solo_ctx.tag_read_info().is_some(); let tag_out_only = params.solo_type == params::SoloType::CbSamTagOut; let keep_barcode_read = params.solo_keeps_barcode_read(); let params_arc = Arc::new(params.clone()); @@ -2956,7 +2969,7 @@ fn align_reads_solo_pe( struct SoloReadProduct { sam_records: BufferedSamRecords, per_feature: Vec, - sj: Vec, + sj: crate::solo::SjRecords, velocyto: Option, } @@ -2992,10 +3005,13 @@ fn align_reads_solo_pe( (0..n_feat).map(|_| Vec::new()).collect(); let mut feat_multi_gene: Vec> = (0..n_feat).map(|_| Vec::new()).collect(); - let mut sj_batch: Vec = Vec::new(); + let mut sj_batch = crate::solo::SjRecords::default(); let mut velo_batch: Vec = Vec::new(); - let mut no_feature: Vec = Vec::new(); - let mut no_feature_multi: Vec = Vec::new(); + // STAR's no-feature records, per feature (readInfo is per feature). + let mut no_feature: Vec> = + (0..n_feat).map(|_| Vec::new()).collect(); + let mut no_feature_multi: Vec> = + (0..n_feat).map(|_| Vec::new()).collect(); for result in products { let product = result?; writer.write_batch(&product.sam_records.records)?; @@ -3010,13 +3026,13 @@ fn align_reads_solo_pe( feat_multi_gene[fi].push(mg); } if let Some(r) = fo.no_feature { - no_feature.push(r); + no_feature[fi].push(r); } if let Some(m) = fo.no_feature_multi { - no_feature_multi.push(m); + no_feature_multi[fi].push(m); } } - sj_batch.extend(product.sj); + sj_batch.append(product.sj); velo_batch.extend(product.velocyto); } for (fi, recorder) in solo.recorders.iter().enumerate() { @@ -3028,24 +3044,23 @@ fn align_reads_solo_pe( if !mg.is_empty() { recorder.multi_gene.lock().unwrap().extend(mg); } - } - if let Some(recorder) = solo.recorders.first() { - if !no_feature.is_empty() { - recorder.no_feature.lock().unwrap().extend(no_feature); + let nf = std::mem::take(&mut no_feature[fi]); + if !nf.is_empty() { + recorder.no_feature.lock().unwrap().extend(nf); } - if !no_feature_multi.is_empty() { - recorder - .no_feature_multi - .lock() - .unwrap() - .extend(no_feature_multi); + let nfm = std::mem::take(&mut no_feature_multi[fi]); + if !nfm.is_empty() { + recorder.no_feature_multi.lock().unwrap().extend(nfm); } } if !sj_batch.is_empty() { - solo.sj_records.lock().unwrap().extend(sj_batch); + solo.sj_records.lock().unwrap().append(sj_batch); } if !velo_batch.is_empty() { - solo.velocyto_records.lock().unwrap().extend(velo_batch); + let mut store = solo.velocyto_records.lock().unwrap(); + for r in &velo_batch { + store.push(r); + } } Ok(true) }; @@ -3095,6 +3110,7 @@ fn align_reads_solo_pe( ); let mut buffer = BufferedSamRecords::new(params.out_sam_attributes); stats.record_read_bases((m1_seq.len() + m2_seq.len()) as u64); + solo.record_qualities(pread.barcode_q30, &[&m1_qual, &m2_qual]); let (results, _pe_chimeric, n_for_mapq, unmapped_reason, best_tr) = align_paired_read( @@ -3187,11 +3203,16 @@ fn align_reads_solo_pe( // counts nothing, exactly as in the single-end loop. let (mut outcome, cb_corrected) = if tag_out_only { let mut outcome = crate::solo::SoloReadOutcome::default(); - let corrected = pread.barcode.as_ref().map(|bc| { - let (tags, corrected) = solo.tag_barcode(bc); - outcome.barcode = Some(tags); - corrected - }); + let corrected = pread + .barcode + .as_ref() + .map_err(|m| solo.record_barcode_miss(*m)) + .ok() + .map(|bc| { + let (tags, corrected) = solo.tag_barcode(bc); + outcome.barcode = Some(tags); + corrected + }); (outcome, corrected) } else if !both_mapped.is_empty() { let pairs: Vec<_> = both_mapped @@ -3201,7 +3222,7 @@ fn align_reads_solo_pe( ( solo.process_read_pe( &pairs, - pread.barcode.as_ref(), + pread.barcode.as_ref().map_err(|m| *m), &junctions, ), None, @@ -3215,13 +3236,21 @@ fn align_reads_solo_pe( solo.process_read( std::slice::from_ref(mapped_transcript), 1, - pread.barcode.as_ref(), + pread.barcode.as_ref().map_err(|m| *m), &junctions, ), None, ) } else { - (solo.process_read(&[], 0, pread.barcode.as_ref(), &[]), None) + ( + solo.process_read( + &[], + 0, + pread.barcode.as_ref().map_err(|m| *m), + &[], + ), + None, + ) }; let read_index = (base + pair_idx as u64) as u32; if track_read_info { @@ -3338,7 +3367,7 @@ fn align_reads_solo_pe( .as_ref() .or_else(|| keep_barcode_read.then_some(&pread.mate1)); let values = crate::solo::SoloTagStrings::build( - pread.barcode.as_ref(), + pread.barcode.as_ref().ok(), barcode_read, outcome.barcode, cb_corrected, @@ -3357,7 +3386,7 @@ fn align_reads_solo_pe( ); } } - if track_read_info { + if stamp_read_index { crate::io::sam::add_solo_read_index( &mut buffer.records, read_index, diff --git a/src/params/mod.rs b/src/params/mod.rs index aac6bdb..47e4ecf 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1188,6 +1188,14 @@ pub struct Parameters { #[arg(long = "outSAMsplicingStatus", default_value = "No")] pub out_sam_splicing_status: String, + /// rustar-aligner extension (not in STAR): `Yes` writes + /// `/CellRanger.summary.csv` for the `Gene` and `GeneFull` + /// features: the fractions of all reads that map uniquely in exonic, + /// intronic or intergenic regions, or antisense to a gene, under + /// CellRanger's metric names. Needs both features. + #[arg(long = "soloCellRangerSummary", default_value = "No")] + pub solo_cellranger_summary: String, + // ── DegNorm (rustar-aligner extension, not STAR) ──────────────────── /// `Yes` writes per-gene, per-exonic-base coverage to /// `GeneCoverage.out.bin` for the DegNorm merge (`--degNormCoverageFiles`). @@ -1449,15 +1457,6 @@ pub struct Parameters { #[arg(long = "soloOutGzip", default_value = "no")] pub solo_out_gzip: String, - /// Velocyto ambiguous-molecule handling (rustar extension beyond STARsolo). - /// `yes` (default) writes the three `spliced`/`unspliced`/`ambiguous` matrices - /// like STARsolo — exon-only molecules with no junction/intron evidence stay in - /// `ambiguous`. `no` resolves those molecules to `spliced` (an exon-only read is - /// most likely mature mRNA; cf. He, Soneson & Patro 2023) and writes only - /// `spliced`/`unspliced`, with no `ambiguous.mtx`. - #[arg(long = "soloVelocytoAmbiguous", default_value = "yes")] - pub solo_velocyto_ambiguous: String, - /// Strand of the read relative to the gene for counting: Forward, Reverse, Unstranded. #[arg(long = "soloStrand", default_value = "Forward")] pub solo_strand: String, @@ -2188,6 +2187,7 @@ impl Parameters { ("--quantGeneSplicing", ¶ms.quant_gene_splicing), ("--outSAMsplicingStatus", ¶ms.out_sam_splicing_status), ("--quantGeneCoverage", ¶ms.quant_gene_coverage), + ("--soloCellRangerSummary", ¶ms.solo_cellranger_summary), ] { if !matches!(value.as_str(), "Yes" | "No") { return Err(command.error( @@ -2418,6 +2418,26 @@ impl Parameters { )); } } + let velocyto = params.solo_features.iter().any(|f| f == "Velocyto"); + if velocyto && params.solo_type == SoloType::SmartSeq { + // `ParametersSolo.cpp:238-241`. + return Err(command.error( + ErrorKind::InvalidValue, + "EXITING because of fatal PARAMETERS error: --soloFeatures Velocyto is \ + presently not compatible with --soloType SmartSeq .\n\ + SOLUTION: re-run without --soloFeatures Velocyto .", + )); + } + // Velocyto takes its cells and UMIs from the Gene feature's readInfo + // (`SoloFeature_countVelocyto.cpp:33-34`). STAR does not check that + // Gene is listed and crashes without it; refused here instead. + if velocyto && !params.solo_features.iter().any(|f| f == "Gene") { + return Err(command.error( + ErrorKind::InvalidValue, + "--soloFeatures Velocyto requires Gene in --soloFeatures: Velocyto takes \ + its cells and UMIs from the Gene feature (STAR crashes without it)", + )); + } // In STAR's order: the whitelist match type (`init_CBmatchWL`, // `ParametersSolo.cpp:672-683`), then CB/UB (`:403-435`), then the // multimappers (`:478-484`). STAR's default match type, 1MM_multi, @@ -2775,6 +2795,11 @@ impl Parameters { self.quant_mode.iter().any(|m| m == "TranscriptomeSAM") } + /// `--soloCellRangerSummary Yes` (rustar-aligner extension). + pub fn solo_cellranger_summary(&self) -> bool { + self.solo_cellranger_summary == "Yes" + } + /// Returns true if `--quantGeneSplicing Yes` (rustar-aligner extension: /// bulk spliced / unspliced / ambiguous gene counts) was requested. pub fn quant_gene_splicing(&self) -> bool { @@ -4217,6 +4242,46 @@ mod tests { assert!(try_parse(&with_ub).is_err()); } + /// Velocyto needs the Gene feature's readInfo, and STAR refuses it with + /// SmartSeq (`ParametersSolo.cpp:238-241`). + #[test] + fn velocyto_refusals() { + let base = |features: &[&'static str], solo_type: &'static str| { + let mut v = vec![ + "--readFilesIn", + "cdna.fq", + "bc.fq", + "--soloType", + solo_type, + "--soloCBwhitelist", + "wl.txt", + "--soloFeatures", + ]; + v.extend_from_slice(features); + v + }; + assert!(try_parse(&base(&["Gene", "Velocyto"], "CB_UMI_Simple")).is_ok()); + let err = try_parse(&base(&["GeneFull", "Velocyto"], "CB_UMI_Simple")).unwrap_err(); + assert!(err.to_string().contains("requires Gene"), "{err}"); + let mut smart_seq = base(&["Gene", "Velocyto"], "SmartSeq"); + smart_seq.extend_from_slice(&["--readFilesManifest", "manifest.tsv"]); + let err = try_parse(&smart_seq).unwrap_err(); + assert!( + err.to_string() + .contains("presently not compatible with --soloType SmartSeq"), + "{err}" + ); + } + + #[test] + fn solo_cellranger_summary_is_yes_or_no() { + let p = try_parse(&["--readFilesIn", "r.fq"]).unwrap(); + assert!(!p.solo_cellranger_summary()); + let p = try_parse(&["--readFilesIn", "r.fq", "--soloCellRangerSummary", "Yes"]).unwrap(); + assert!(p.solo_cellranger_summary()); + assert!(try_parse(&["--readFilesIn", "r.fq", "--soloCellRangerSummary", "yes"]).is_err()); + } + #[test] fn out_sam_order_accepts_star_values_rejects_others() { assert!(try_parse(&["--readFilesIn", "r.fq", "--outSAMorder", "Paired"]).is_ok()); diff --git a/src/quant/splice_status.rs b/src/quant/splice_status.rs index e86680a..4c3afe1 100644 --- a/src/quant/splice_status.rs +++ b/src/quant/splice_status.rs @@ -379,6 +379,17 @@ pub fn collapse_gene_category( } } +/// STAR's `countVelocyto` call for one molecule from the type bits of its +/// transcripts (all of one gene; the caller checks that): +/// `SoloFeature_countVelocyto.cpp:102-133`. +pub fn velocyto_status(types: impl IntoIterator) -> SpliceStatus { + let mut models = Models::default(); + for ty in types { + models.add(ty); + } + models.status() +} + /// The model flags of STAR's `countVelocyto` collapse. #[allow(clippy::struct_excessive_bools)] // mirrors STAR's four flags struct Models { diff --git a/src/solo/count.rs b/src/solo/count.rs index 4e7ba43..2acfb8a 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -483,8 +483,13 @@ pub struct MatrixStats { pub nnz: usize, /// One entry per barcode that received ≥1 UMI (the raw, unfiltered set). pub cells: Vec, - /// Distinct genes with a nonzero count anywhere in the raw matrix. - pub genes_detected: u32, + /// Counted unique-gene reads, all barcodes (STAR's `nReadPerCBunique` + /// summed: `yessubWLmatch_UniqueFeature`). + pub unique_reads: u64, + /// Reads dropped for a 1MM barcode never seen exactly (`noMMtoWLwithoutExact`). + pub dropped_one_exact: u64, + /// 1MM_multi reads whose posterior picked no barcode (`noTooManyWLmatches`). + pub dropped_posterior: u64, } /// Stream the per-cell deduplicated counts into a plain temporary MatrixMarket @@ -500,16 +505,17 @@ fn build_matrix_body( umi_len: usize, pseudocount: f64, dir: &Path, - n_features: usize, - record_read_info: bool, + read_info_slot: Option<&std::sync::Mutex>>, ) -> Result<(tempfile::NamedTempFile, MatrixStats), Error> { + let record_read_info = read_info_slot.is_some(); let mut body_tmp = tempfile::Builder::new() .prefix(".matrix_body") .tempfile_in(dir) .map_err(|e| Error::io(e, dir))?; let mut nnz = 0usize; let mut cell_stats: Vec = Vec::new(); - let mut gene_seen = vec![false; n_features]; + let (unique_reads, dropped_one_exact); + let mut dropped_posterior = 0u64; { let mut body = std::io::BufWriter::new(body_tmp.as_file_mut()); @@ -519,7 +525,9 @@ fn build_matrix_body( let exact_counts = ctx.whitelist.exact_count_snapshot(); // `oneExact`: drop single-1MM corrections to a barcode never seen exactly. let one_exact = ctx.match_type.one_exact(); + let before = records.len(); records.retain(|r| cb_passes_one_exact(one_exact, &exact_counts, r.cb)); + dropped_one_exact = (before - records.len()) as u64; let multi = std::mem::take(&mut *recorder.multi_records.lock().unwrap()); for m in &multi { if let Some(cb) = resolve_multi_cb(&m.candidates, &exact_counts, pseudocount) { @@ -529,9 +537,12 @@ fn build_matrix_body( gene: m.gene, read_index: m.read_index, }); + } else { + dropped_posterior += 1; } } drop(multi); + unique_reads = records.len() as u64; // Group each cell's reads together (parallel sort — the record vec is large). use rayon::prelude::*; @@ -658,9 +669,7 @@ fn build_matrix_body( .collect(); // Sequential merge: byte order preserved (CB-ascending, gene-ascending). - let mut info_guard = record_read_info - .then(|| ctx.read_info.as_ref().map(|m| m.lock().unwrap())) - .flatten(); + let mut info_guard = read_info_slot.map(|m| m.lock().unwrap()); // No-feature reads keep the cell their barcode resolves to and their // raw UMI (`SoloReadFeature_inputRecords.cpp:35-65, 68-103`): the same // oneExact and 1MM_multi resolution as counted reads. @@ -684,9 +693,6 @@ fn build_matrix_body( for co in cell_outs { body.write_all(&co.body).map_err(|e| Error::io(e, dir))?; nnz += co.genes.len(); - for g in co.genes { - gene_seen[g as usize] = true; - } if let Some(s) = co.stat { cell_stats.push(s); } @@ -702,13 +708,14 @@ fn build_matrix_body( body.flush().map_err(|e| Error::io(e, dir))?; } - let genes_detected = gene_seen.iter().filter(|&&s| s).count() as u32; Ok(( body_tmp, MatrixStats { nnz, cells: cell_stats, - genes_detected, + unique_reads, + dropped_one_exact, + dropped_posterior, }, )) } @@ -906,13 +913,10 @@ fn distribute_multi( out } -/// Format a real matrix value compactly (integers without a decimal point). +/// Format a real matrix value as STAR's `ostream << double` does +/// (`SoloFeature_outputResults.cpp:181-223`): [`fmt_double`]. fn fmt_real(v: f64) -> String { - if v.fract().abs() < 1e-9 { - format!("{}", v.round() as i64) - } else { - format!("{v:.5}") - } + fmt_double(v) } /// Write the `UniqueAndMult-.mtx` matrices (real-valued) for the @@ -1212,11 +1216,23 @@ fn called_cells(cells: &[CellStat], filter: &[String]) -> Option> { let arg = |i: usize, d: f64| filter.get(i).and_then(|s| s.parse().ok()).unwrap_or(d); let mut cbs: Vec = match method { "None" => return None, + // STAR's threshold, not a cut at N (`SoloFeature_cellFiltering.cpp:46-54`): + // the UMIs of the barcode at index N of the descending list (the + // N+1th), at least 1, and every barcode with at least that many, ties + // included. "TopCells" => { let n = arg(1, 0.0) as usize; - let mut idx: Vec<&CellStat> = cells.iter().collect(); - idx.sort_by(|a, b| b.n_umis.cmp(&a.n_umis).then(a.cb.cmp(&b.cb))); - idx.into_iter().take(n).map(|c| c.cb).collect() + let mut umis: Vec = cells.iter().map(|c| c.n_umis).collect(); + umis.sort_unstable_by(|a, b| b.cmp(a)); + let Some(&at_n) = umis.get(n.min(umis.len().saturating_sub(1))) else { + return Some(Vec::new()); + }; + let min = at_n.max(1); + cells + .iter() + .filter(|c| c.n_umis >= min) + .map(|c| c.cb) + .collect() } // EmptyDrops_CR is handled by `emptydrops_called`; the knee here is the // fallback / guaranteed-cell base. @@ -1482,20 +1498,8 @@ fn emptydrops_called( Ok(called) } -/// Median of an ascending-sorted slice (0 if empty). -fn median_sorted(sorted: &[u64]) -> u64 { - let n = sorted.len(); - if n == 0 { - 0 - } else if n % 2 == 1 { - sorted[n / 2] - } else { - u64::midpoint(sorted[n / 2 - 1], sorted[n / 2]) - } -} - /// Write the raw gene-count matrix + `Summary.csv` for a finished solo run. -/// No-op (with a warning) when there is no explicit whitelist. +/// Only `Barcodes.stats` (with a warning) when there is no explicit whitelist. pub fn write_gene_matrix( ctx: &SoloContext, params: &crate::params::Parameters, @@ -1503,6 +1507,8 @@ pub fn write_gene_matrix( sj_stats: Option<&crate::junction::SpliceJunctionStats>, genome: &crate::genome::Genome, ) -> Result<(), Error> { + // STAR writes it before any feature (`Solo.cpp:53-64`). + write_barcode_stats(ctx, params)?; let CbWhitelist::List { sorted, .. } = &ctx.whitelist else { log::warn!( "STARsolo: --soloCBwhitelist None matrix output is not yet supported (Phase 14.4); skipping matrix" @@ -1558,23 +1564,31 @@ pub fn write_gene_matrix( let total_reads = align_stats.total_reads.load(Ordering::Relaxed); let mapped_unique = align_stats.uniquely_mapped.load(Ordering::Relaxed); let mapped_multi = align_stats.multi_mapped.load(Ordering::Relaxed); - let valid_barcodes = ctx.stats.yes_exact.load(Ordering::Relaxed) - + ctx.stats.yes_one_mm.load(Ordering::Relaxed) - + ctx.stats.yes_mult_mm.load(Ordering::Relaxed); - let reads_of = |f: crate::solo::SoloFeature| -> u64 { - ctx.features - .iter() - .position(|&x| x == f) - .map_or(0, |i| ctx.feature_reads[i].load(Ordering::Relaxed)) + let summary_run = SummaryContext { + n_reads: total_reads, + mapped_unique, + mapped_multi, + barcodes: &ctx.stats, + rna_q30: ( + ctx.rna_q30[0].load(Ordering::Relaxed), + ctx.rna_q30[1].load(Ordering::Relaxed), + ), + multimappers: !MultiMethod::parse_list(¶ms.solo_multi_mappers).is_empty(), }; + // `--soloCellRangerSummary Yes` only (a rustar-aligner extension): the + // funnel needs both the Gene and GeneFull overlaps. let have_funnel = ctx.features.contains(&crate::solo::SoloFeature::Gene) && ctx.features.contains(&crate::solo::SoloFeature::GeneFull); - let region = have_funnel.then(|| RegionFunnel { - exonic: ctx.region_stats.exonic.load(Ordering::Relaxed), - intronic: ctx.region_stats.intronic.load(Ordering::Relaxed), - intergenic: ctx.region_stats.intergenic.load(Ordering::Relaxed), - antisense: ctx.region_stats.antisense.load(Ordering::Relaxed), - }); + let region = ctx + .region_stats + .as_ref() + .filter(|_| have_funnel) + .map(|r| RegionFunnel { + exonic: r.exonic.load(Ordering::Relaxed), + intronic: r.intronic.load(Ordering::Relaxed), + intergenic: r.intergenic.load(Ordering::Relaxed), + antisense: r.antisense.load(Ordering::Relaxed), + }); let gzip = matches!(params.solo_out_gzip.as_str(), "yes" | "Yes" | "true"); let n_genes = ctx.gene_ann.gene_ids.len(); @@ -1584,6 +1598,9 @@ pub fn write_gene_matrix( let exact_counts = ctx.whitelist.exact_count_snapshot(); let multi_methods = MultiMethod::parse_list(¶ms.solo_multi_mappers); + // The Gene feature's called cells, which Velocyto's filtered matrices use. + let mut gene_called: Option> = None; + // One {prefix}{soloOutFileNames[0]}/{raw,filtered}/ per feature. for (fi, (feature, recorder)) in ctx.features.iter().zip(&ctx.recorders).enumerate() { let feature_dir = params.output_path(&format!("{solo_dir}{}/", feature.dir_name())); @@ -1600,10 +1617,9 @@ pub fn write_gene_matrix( umi_len, pseudocount, &raw_dir, - n_genes, - // STAR fills readInfo from one feature only: the first on the - // --soloFeatures list (`ParametersSolo.cpp:423-434`). - fi == 0 && ctx.read_info_enabled(), + // The feature's readInfo, if it keeps one: the first feature for + // the CB/UB tags, Gene for Velocyto (`ParametersSolo.cpp:418-436`). + ctx.read_info[fi].as_ref(), )?; write_features( &raw_dir.join(&features_name), @@ -1673,6 +1689,10 @@ pub fn write_gene_matrix( } else { called_cells(&mstats.cells, ¶ms.solo_cell_filter) }; + if *feature == crate::solo::SoloFeature::Gene { + gene_called.clone_from(&called); + } + let called_for_summary = called.clone(); if let Some(cbs) = called && !cbs.is_empty() { @@ -1709,9 +1729,13 @@ pub fn write_gene_matrix( } // --soloMultiMappers: UniqueAndMult-.mtx alongside raw. + let (mut multi_counted, mut multi_dropped) = (0u64, 0u64); if !multi_methods.is_empty() { let mut mg = recorder.multi_gene.lock().unwrap(); + let before = mg.len(); mg.retain(|r| cb_passes_one_exact(one_exact, &exact_counts, r.cb)); + multi_counted = mg.len() as u64; + multi_dropped = (before - mg.len()) as u64; build_multi_matrices( &body, &mg, @@ -1730,16 +1754,46 @@ pub fn write_gene_matrix( ); } + // Features.stats: the counting-stage figures (`inputRecords`, + // `collapseUMIall`) join the per-read ones from `process_read`. + { + use crate::solo::FeatureStats as F; + let st = &ctx.feature_stats[fi]; + st.add( + F::NO_MM_TO_WL_WITHOUT_EXACT, + mstats.dropped_one_exact + multi_dropped, + ); + st.add(F::NO_TOO_MANY_WL_MATCHES, mstats.dropped_posterior); + st.add(F::YES_WL_MATCH, mstats.unique_reads + multi_counted); + st.add(F::YES_SUB_WL_MATCH_UNIQUE_FEATURE, mstats.unique_reads); + st.add(F::YES_UMIS, mstats.cells.iter().map(|c| c.n_umis).sum()); + st.add( + F::YES_CELL_BARCODES, + mstats.cells.iter().filter(|c| c.n_genes > 0).count() as u64, + ); + let values: Vec = (0..F::NAMES.len()).map(|i| st.get(i)).collect(); + write_stats_file(&feature_dir.join("Features.stats"), &F::NAMES, &values)?; + } + // Summary.csv, with the called-cell block and UMIperCellSorted.txt + // unless `--soloCellFilter None`. + let cell_summary = match &called_for_summary { + Some(cbs) => Some(CellSummary::new( + &mstats.cells, + cbs, + genes_in_cells(&body, cbs)?, + )), + None => None, + }; write_summary( &feature_dir.join("Summary.csv"), feature.dir_name(), - &mstats, - total_reads, - valid_barcodes, - mapped_unique, - mapped_multi, - reads_of(*feature), + &summary_run, + &ctx.feature_stats[fi], + cell_summary.as_ref(), )?; + if called_for_summary.is_some() { + write_umi_per_cell_sorted(&feature_dir.join("UMIperCellSorted.txt"), &mstats.cells)?; + } log::info!("STARsolo: wrote {}/Summary.csv", feature.dir_name()); // CellRanger-style mapping funnel goes in a SEPARATE additional file so the // faithful Summary.csv is never altered (PR #90 review: keep this release a @@ -1822,9 +1876,34 @@ pub fn write_gene_matrix( gzip, )?; let umi_len = params.solo_umi_len as usize; - let mut sj_records = ctx.sj_records.lock().unwrap(); - sj_records.retain(|r| cb_passes_one_exact(one_exact, &exact_counts, r.cb)); - let nnz = build_sj_matrix( + // STAR's per-line filters (`SoloReadFeature_inputRecords.cpp:23-104`), + // counted only for a junction with a row: STAR skips the others first. + let has_row = |start: u64, end: u64| row.contains_key(&(start, end)); + let crate::solo::SjRecords { + unique: mut sj_records, + multi, + } = std::mem::take(&mut *ctx.sj_records.lock().unwrap()); + let mut dropped_one_exact = 0u64; + sj_records.retain(|r| { + let keep = cb_passes_one_exact(one_exact, &exact_counts, r.cb); + dropped_one_exact += u64::from(!keep && has_row(r.intron_start, r.intron_end)); + keep + }); + let mut dropped_posterior = 0u64; + for m in multi { + if let Some(cb) = resolve_multi_cb(&m.candidates, &exact_counts, pseudocount) { + sj_records.push(crate::solo::SjCountRecord { + cb, + umi: m.umi, + intron_start: m.intron_start, + intron_end: m.intron_end, + exact: false, + }); + } else if has_row(m.intron_start, m.intron_end) { + dropped_posterior += 1; + } + } + let sj_counts = build_sj_matrix( &sj_records, &row, method, @@ -1838,8 +1917,21 @@ pub fn write_gene_matrix( "STARsolo: wrote SJ/raw matrix ({} junctions × {} barcodes, {} entries)", order.len(), sorted.len(), - nnz, + sj_counts.nnz, ); + { + use crate::solo::FeatureStats as F; + let st = &ctx.sj_stats; + st.add(F::NO_MM_TO_WL_WITHOUT_EXACT, dropped_one_exact); + st.add(F::NO_TOO_MANY_WL_MATCHES, dropped_posterior); + st.add(F::YES_WL_MATCH, sj_counts.reads); + st.add(F::YES_SUB_WL_MATCH_EXACT, sj_counts.exact); + st.add(F::YES_SUB_WL_MATCH_UNIQUE_FEATURE, sj_counts.reads); + st.add(F::YES_UMIS, sj_counts.umis); + st.add(F::YES_CELL_BARCODES, sj_counts.cells); + let sj_root = params.output_path(&format!("{solo_dir}SJ/")); + write_feature_stats_and_summary(&sj_root, "SJ", &summary_run, st)?; + } } // Velocyto feature: spliced / unspliced / ambiguous gene×cell matrices. @@ -1859,37 +1951,91 @@ pub fn write_gene_matrix( sorted.len(), gzip, )?; - let umi_len = params.solo_umi_len as usize; - // `--soloVelocytoAmbiguous no` folds exon-only molecules into spliced and - // omits ambiguous.mtx (rustar extension); default `yes` = STARsolo 3-matrix. - let keep_ambiguous = !matches!( - params.solo_velocyto_ambiguous.as_str(), - "no" | "No" | "false" - ); - let mut velocyto_records = ctx.velocyto_records.lock().unwrap(); - velocyto_records.retain(|r| cb_passes_one_exact(one_exact, &exact_counts, r.cb)); - let nnz = build_velocyto_matrices( - &velocyto_records, - method, - umi_len, + // Cells and UMIs come from the Gene feature's readInfo, which UMI + // collapsing of that feature filled above. + let records = std::mem::take(&mut *ctx.velocyto_records.lock().unwrap()); + let counts = { + let gene_info = ctx + .gene_read_info() + .expect("Velocyto keeps the Gene feature's readInfo") + .lock() + .unwrap(); + count_velocyto(&records, &gene_info, &ctx.gene_model) + }; + let nnz = write_velocyto_matrices( + &counts, + |cb| Some(cb + 1), &velo_dir, n_genes, sorted.len(), gzip, - keep_ambiguous, )?; - if keep_ambiguous { - log::info!( - "STARsolo: wrote Velocyto/raw matrices (spliced={} unspliced={} ambiguous={} entries)", - nnz[0], - nnz[1], - nnz[2], + log::info!( + "STARsolo: wrote Velocyto/raw matrices ({} entries each)", + nnz[0] + ); + { + // `countVelocyto` adds only the UMIs and cells (`:144-145`). + use crate::solo::FeatureStats as F; + let st = &ctx.velocyto_stats; + st.add( + F::YES_UMIS, + counts + .iter() + .map(|&(_, _, c)| u64::from(c.iter().sum::())) + .sum(), ); - } else { + let mut cells: Vec = counts.iter().map(|&(cb, _, _)| cb).collect(); + cells.dedup(); + st.add(F::YES_CELL_BARCODES, cells.len() as u64); + let velo_root = params.output_path(&format!("{solo_dir}Velocyto/")); + write_feature_stats_and_summary(&velo_root, "Velocyto", &summary_run, st)?; + } + + // Filtered: the Gene feature's called cells that Velocyto detected, + // i.e. that a read with a Velocyto record matched + // (`SoloFeature_cellFiltering.cpp:16-30`). + if let Some(gene_cells) = &gene_called + && !gene_cells.is_empty() + { + let mut cells: Vec = gene_cells + .iter() + .copied() + .filter(|&cb| records.detected(cb)) + .collect(); + cells.sort_unstable(); + // No detected cell: STAR writes no filtered matrices + // (`SoloFeature_cellFiltering.cpp:106-108`). + if cells.is_empty() { + return Ok(()); + } + let column: HashMap = cells + .iter() + .enumerate() + .map(|(i, &cb)| (cb, i as u32 + 1)) + .collect(); + let filt_dir = params.output_path(&format!("{solo_dir}Velocyto/filtered/")); + std::fs::create_dir_all(&filt_dir).map_err(|e| Error::io(e, &filt_dir))?; + write_features( + &filt_dir.join(&features_name), + &ctx.gene_ann.gene_ids, + &ctx.gene_ann.gene_names, + ¶ms.solo_out_format_features_gene_field3, + gzip, + )?; + write_barcodes_subset(&filt_dir.join(&barcodes_name), &ctx.whitelist, &cells, gzip)?; + let fnnz = write_velocyto_matrices( + &counts, + |cb| column.get(&cb).copied(), + &filt_dir, + n_genes, + cells.len(), + gzip, + )?; log::info!( - "STARsolo: wrote Velocyto/raw matrices, ambiguous folded into spliced (spliced={} unspliced={} entries)", - nnz[0], - nnz[1], + "STARsolo: wrote Velocyto/filtered matrices ({} cells, {} entries each)", + cells.len(), + fnnz[0], ); } } @@ -1910,7 +2056,7 @@ fn build_sj_matrix( n_junctions: usize, n_barcodes: usize, gzip: bool, -) -> Result { +) -> Result { // Group by cell barcode (ascending column order). use rayon::prelude::*; let mut recs: Vec<&crate::solo::SjCountRecord> = records.iter().collect(); @@ -1922,6 +2068,7 @@ fn build_sj_matrix( .tempfile_in(dir) .map_err(|e| Error::io(e, dir))?; let mut nnz = 0usize; + let mut counts = SjCounts::default(); { let mut body = std::io::BufWriter::new(body_tmp.as_file_mut()); let mut i = 0; @@ -1933,6 +2080,8 @@ fn build_sj_matrix( let r = recs[i]; if let Some(&rw) = row.get(&(r.intron_start, r.intron_end)) { *sj_umis.entry(rw).or_default().entry(r.umi).or_insert(0) += 1; + counts.reads += 1; + counts.exact += u64::from(r.exact); } i += 1; } @@ -1942,7 +2091,11 @@ fn build_sj_matrix( .filter(|&(_, c)| c > 0) .collect(); entries.sort_unstable_by_key(|&(rw, _)| rw); + if !entries.is_empty() { + counts.cells += 1; + } for (rw, c) in entries { + counts.umis += c; writeln!(body, "{} {} {}", rw + 1, cb + 1, c).map_err(|e| Error::io(e, dir))?; nnz += 1; } @@ -1960,152 +2113,148 @@ fn build_sj_matrix( std::io::copy(&mut r, w).map_err(|e| Error::io(e, matrix_path))?; Ok(()) })?; - Ok(nnz) + counts.nnz = nnz; + Ok(counts) } -/// Build the `Velocyto` matrices from (cell, UMI, gene, category) records. Per -/// (cell, gene) each UMI is resolved to one category (priority unspliced over -/// spliced over ambiguous — any intron evidence makes the molecule nascent), then -/// UMI-deduplicated per category. Genes are rows, cells columns — same layout as -/// the Gene matrix, written as files scVelo/dynamo ingest directly. -/// -/// With `keep_ambiguous` (default, STARsolo-faithful) three matrices are written: -/// `spliced`/`unspliced`/`ambiguous`. With `keep_ambiguous = false` the exon-only -/// `ambiguous` molecules are folded into `spliced` (an exon-only read is most -/// likely mature mRNA; cf. He, Soneson & Patro 2023) and only `spliced`/`unspliced` -/// are written — no `ambiguous.mtx`. The returned `[usize; 3]` always reports -/// `[spliced, unspliced, ambiguous]` nnz (ambiguous is 0 when folded). -#[allow(clippy::too_many_arguments)] -fn build_velocyto_matrices( - records: &[crate::solo::VelocytoRecord], - method: UmiDedup, - umi_len: usize, - dir: &Path, - n_genes: usize, - n_barcodes: usize, - gzip: bool, - keep_ambiguous: bool, -) -> Result<[usize; 3], Error> { - use crate::solo::VelocytoCategory; - // Category → matrix index (file order) and resolution priority. - let cat_idx = |c: VelocytoCategory| match c { - VelocytoCategory::Spliced => 0usize, - VelocytoCategory::Unspliced => 1, - VelocytoCategory::Ambiguous => 2, - }; - let priority = |c: VelocytoCategory| match c { - VelocytoCategory::Unspliced => 2u8, - VelocytoCategory::Spliced => 1, - VelocytoCategory::Ambiguous => 0, - }; - let names = ["spliced.mtx", "unspliced.mtx", "ambiguous.mtx"]; +/// What the SJ matrix counted: its entries, the records whose junction is a +/// row (STAR's counted SJ "reads", one per junction), their UMIs, and the +/// cells with any. +#[derive(Debug, Default)] +struct SjCounts { + nnz: usize, + reads: u64, + /// Of `reads`, those with an exactly matched barcode. + exact: u64, + umis: u64, + cells: u64, +} - let mut recs: Vec<&crate::solo::VelocytoRecord> = records.iter().collect(); - use rayon::prelude::*; - recs.par_sort_unstable_by_key(|r| r.cb); +/// `Features.stats` and `Summary.csv` of a feature without cell calling (SJ, +/// Velocyto), with its counting-stage figures added. +fn write_feature_stats_and_summary( + dir: &Path, + name: &str, + run: &SummaryContext<'_>, + stats: &crate::solo::FeatureStats, +) -> Result<(), Error> { + use crate::solo::FeatureStats as F; + let values: Vec = (0..F::NAMES.len()).map(|i| stats.get(i)).collect(); + write_stats_file(&dir.join("Features.stats"), &F::NAMES, &values)?; + write_summary(&dir.join("Summary.csv"), name, run, stats, None) +} - // One contiguous [start, end) slice per CB. - let mut bounds: Vec<(usize, usize)> = Vec::new(); - let mut i = 0; - while i < recs.len() { - let cb = recs[i].cb; - let mut j = i + 1; - while j < recs.len() && recs[j].cb == cb { - j += 1; +/// STAR's `countVelocyto` (`SoloFeature_countVelocyto.cpp`): the counts of +/// each (cell, gene), `[spliced, unspliced, ambiguous]`, sorted by cell then +/// gene. +/// +/// - Each record's cell and UMI are the Gene feature's readInfo (`:33-38`); +/// a read with either undefined is skipped. +/// - A molecule (cell, UMI) keeps the transcripts all its reads are +/// compatible with, their type bits OR-ed (`:43-77`); once that set is +/// empty, the molecule's later reads are skipped. +/// - A molecule whose transcripts span more than one gene is dropped; the +/// others are spliced, unspliced or ambiguous by STAR's model flags +/// (`:95-136`, [`velocyto_status`]). +fn count_velocyto( + records: &crate::solo::VelocytoRecords, + read_info: &[crate::solo::ReadInfo], + model: &crate::solo::GeneModel, +) -> Vec<(u32, u32, [u32; 3])> { + use crate::quant::splice_status::velocyto_status; + + // cell -> UMI -> the molecule's transcripts so far. + let mut molecules: HashMap>> = HashMap::default(); + for (read_index, types) in records.iter() { + let Some(info) = read_info.get(read_index as usize) else { + continue; + }; + if info.cb == u32::MAX || info.umi == u64::MAX { + continue; } - bounds.push((i, j)); - i = j; - } - - // Per-cell dedup is independent across cells → run in parallel, each cell - // producing the three matrices' lines (gene-ascending). Merge sequentially - // in CB order so the three .mtx files stay byte-identical to the serial path. - type VeloCellOut = ([Vec; 3], [usize; 3]); - let cell_outs: Vec = bounds - .par_iter() - .map(|&(lo, hi)| { - let cb = recs[lo].cb; - // gene → umi → (resolved category, read count) - let mut gene_umi: HashMap> = - HashMap::default(); - for &r in &recs[lo..hi] { - let e = gene_umi - .entry(r.gene) - .or_default() - .entry(r.umi) - .or_insert((r.category, 0)); - e.1 += 1; - if priority(r.category) > priority(e.0) { - e.0 = r.category; - } + let umis = molecules.entry(info.cb).or_default(); + match umis.get_mut(&info.umi) { + None => { + umis.insert(info.umi, types.collect()); } - // Per gene, dedup UMIs within each resolved category, emit entries. - let mut genes: Vec<&u32> = gene_umi.keys().collect(); - genes.sort_unstable(); - let mut bufs: [Vec; 3] = [Vec::new(), Vec::new(), Vec::new()]; - let mut cnt = [0usize; 3]; - for &g in &genes { - let umis = &gene_umi[g]; - let mut by_cat: [HashMap; 3] = - [HashMap::default(), HashMap::default(), HashMap::default()]; - for (&umi, &(cat, rc)) in umis { - by_cat[cat_idx(cat)].insert(umi, rc); - } - // Fold ambiguous (exon-only) molecules into spliced. A UMI resolves - // to exactly one category, so spliced/ambiguous keys are disjoint. - if !keep_ambiguous { - let amb = std::mem::take(&mut by_cat[2]); - by_cat[0].extend(amb); - } - for (k, buf) in bufs.iter_mut().enumerate() { - let c = dedup_count(&by_cat[k], method, umi_len); - if c > 0 { - let _ = writeln!(buf, "{} {} {}", g + 1, cb + 1, c); - cnt[k] += 1; + Some(kept) if kept.is_empty() => {} + Some(kept) => { + // Both lists are sorted by transcript (`:66-76`). + let mut types = types.peekable(); + let mut both = Vec::with_capacity(kept.len()); + for &(tr, ty) in kept.iter() { + while types.next_if(|&(new, _)| tr > new).is_some() {} + let Some(&(new, new_ty)) = types.peek() else { + break; + }; + if tr == new { + both.push((tr, ty | new_ty)); } } + *kept = both; } - (bufs, cnt) - }) - .collect(); - - let mut bodies: Vec = Vec::new(); - for _ in 0..3 { - bodies.push( - tempfile::Builder::new() - .prefix(".velo_body") - .tempfile_in(dir) - .map_err(|e| Error::io(e, dir))?, - ); + } } - let mut nnz = [0usize; 3]; - { - let mut writers: Vec> = bodies - .iter_mut() - .map(|t| std::io::BufWriter::new(t.as_file_mut())) - .collect(); - for (bufs, cnt) in &cell_outs { - for (k, w) in writers.iter_mut().enumerate() { - w.write_all(&bufs[k]).map_err(|e| Error::io(e, dir))?; - nnz[k] += cnt[k]; + + let mut out = Vec::new(); + let mut cells: Vec = molecules.keys().copied().collect(); + cells.sort_unstable(); + for cb in cells { + let mut gene_counts: std::collections::BTreeMap = + std::collections::BTreeMap::new(); + for trs in molecules[&cb].values() { + let Some(&(tr0, _)) = trs.first() else { + continue; // no transcript in the intersection + }; + let gene = model.tr_gene(tr0 as usize); + if gene == u32::MAX + || trs + .iter() + .any(|&(tr, _)| model.tr_gene(tr as usize) != gene) + { + continue; // multi-gene } + let status = velocyto_status(trs.iter().map(|&(_, ty)| ty)); + gene_counts.entry(gene).or_insert([0; 3])[status as usize] += 1; } - for w in &mut writers { - w.flush().map_err(|e| Error::io(e, dir))?; - } + out.extend(gene_counts.into_iter().map(|(gene, c)| (cb, gene, c))); } + out +} - // Three files by default; only spliced/unspliced when ambiguous is folded. - let n_out = if keep_ambiguous { 3 } else { 2 }; - for (k, body) in bodies.iter().take(n_out).enumerate() { +/// Write the three Velocyto matrices for the cells `column` numbers (1-based +/// column, `None` to leave a cell out), as STAR's `outputResults` does: each +/// (cell, gene) row appears once in every matrix, zeros included, since the +/// three share one `countCellGeneUMI` row per gene. Returns the entries per +/// matrix. +fn write_velocyto_matrices( + counts: &[(u32, u32, [u32; 3])], + column: impl Fn(u32) -> Option, + dir: &Path, + n_genes: usize, + n_columns: usize, + gzip: bool, +) -> Result<[usize; 3], Error> { + let names = ["spliced.mtx", "unspliced.mtx", "ambiguous.mtx"]; + let mut bodies: [Vec; 3] = [Vec::new(), Vec::new(), Vec::new()]; + let mut nnz = [0usize; 3]; + for &(cb, gene, c) in counts { + let Some(col) = column(cb) else { + continue; + }; + for k in 0..3 { + let _ = writeln!(bodies[k], "{} {} {}", gene + 1, col, c[k]); + nnz[k] += 1; + } + } + for (k, body) in bodies.iter().enumerate() { let path = dir.join(names[k]); write_file(&path, gzip, |w| { writeln!(w, "%%MatrixMarket matrix coordinate integer general") .map_err(|e| Error::io(e, &path))?; writeln!(w, "%").map_err(|e| Error::io(e, &path))?; - writeln!(w, "{n_genes} {n_barcodes} {}", nnz[k]).map_err(|e| Error::io(e, &path))?; - let mut r = std::fs::File::open(body.path()).map_err(|e| Error::io(e, body.path()))?; - std::io::copy(&mut r, w).map_err(|e| Error::io(e, &path))?; + writeln!(w, "{n_genes} {n_columns} {}", nnz[k]).map_err(|e| Error::io(e, &path))?; + w.write_all(body).map_err(|e| Error::io(e, &path))?; Ok(()) })?; } @@ -2121,130 +2270,263 @@ struct RegionFunnel { antisense: u64, } -/// Write the STARsolo-faithful `Summary.csv` for one feature: the sequencing / -/// genome-mapping rows plus per-cell UMI/gene statistics over the CR2.2-knee-called -/// cells. The CellRanger-style exonic/intronic/intergenic/antisense funnel is a -/// rustar extension kept in a SEPARATE file (see `write_cellranger_summary`) so -/// this file is never altered relative to STARsolo's own Summary.csv. -#[allow(clippy::too_many_arguments)] -fn write_summary( - path: &Path, - feature_name: &str, - mstats: &MatrixStats, - total_reads: u64, - valid_barcodes: u64, - mapped_unique: u64, - mapped_multi: u64, - feature_mapped: u64, -) -> Result<(), Error> { - let frac = |num: u64| -> f64 { - if total_reads == 0 { - 0.0 +/// A double as C++ `ostream << double` prints it with the default precision: +/// six significant digits, `%g` style, trailing zeros dropped (`0.0357143`, +/// `1`, `5e-05`), and `-inf` / `-nan` for STAR's divisions by zero. +pub(crate) fn fmt_double(v: f64) -> String { + if v.is_nan() { + return "-nan".to_string(); + } + if v.is_infinite() { + return if v > 0.0 { "inf" } else { "-inf" }.to_string(); + } + if v == 0.0 { + return "0".to_string(); + } + let trim = |t: String| { + if t.contains('.') { + t.trim_end_matches('0').trim_end_matches('.').to_string() } else { - num as f64 / total_reads as f64 + t } }; - - // Cell calling: CR2.2 knee on per-barcode UMI totals. - let mut umis_desc: Vec = mstats.cells.iter().map(|c| c.n_umis).collect(); - umis_desc.sort_unstable_by(|a, b| b.cmp(a)); - let thr = knee_cr22(&umis_desc, 3000, 0.99, 10.0); - let cells: Vec<&CellStat> = mstats.cells.iter().filter(|c| c.n_umis >= thr).collect(); - let n_cells = cells.len(); - - // Totals across all barcodes (for sequencing saturation + fraction-in-cells). - let total_reads_counted: u64 = mstats.cells.iter().map(|c| c.n_reads).sum(); - let total_umis_all: u64 = mstats.cells.iter().map(|c| c.n_umis).sum(); - let saturation = if total_reads_counted > 0 { - 1.0 - total_umis_all as f64 / total_reads_counted as f64 + // The exponent after rounding to six significant digits decides the form. + let sci = format!("{v:.5e}"); + let (mantissa, exp) = sci.split_once('e').expect("{:e} has an exponent"); + let exp: i32 = exp.parse().expect("integer exponent"); + if !(-4..6).contains(&exp) { + let sign = if exp < 0 { '-' } else { '+' }; + format!("{}e{sign}{:02}", trim(mantissa.to_string()), exp.abs()) } else { - 0.0 - }; + trim(format!("{v:.*}", (5 - exp) as usize)) + } +} - // Per-cell aggregates over called cells. - let reads_in_cells: u64 = cells.iter().map(|c| c.n_reads).sum(); - let umis_in_cells: u64 = cells.iter().map(|c| c.n_umis).sum(); - let mut reads_sorted: Vec = cells.iter().map(|c| c.n_reads).collect(); - let mut umis_sorted: Vec = cells.iter().map(|c| c.n_umis).collect(); - let mut genes_sorted: Vec = cells.iter().map(|c| c.n_genes as u64).collect(); - reads_sorted.sort_unstable(); - umis_sorted.sort_unstable(); - genes_sorted.sort_unstable(); - let mean = |sum: u64| -> u64 { - if n_cells == 0 { - 0 - } else { - sum / n_cells as u64 +/// Write a STAR stats file: `Barcodes.stats` or `/Features.stats` +/// (`statsOut`, one `setw(50)` name and `setw(15)` value per line). +fn write_stats_file(path: &Path, names: &[&str], values: &[u64]) -> Result<(), Error> { + let text: String = names + .iter() + .zip(values) + .map(|(n, &v)| crate::solo::whitelist::stats_line(n, v)) + .collect(); + std::fs::write(path, text).map_err(|e| Error::io(e, path)) +} + +/// `Solo.out/Barcodes.stats` (`Solo::processAndOutput`). +pub fn write_barcode_stats( + ctx: &SoloContext, + params: &crate::params::Parameters, +) -> Result<(), Error> { + let solo_dir = params + .solo_out_file_names + .first() + .cloned() + .unwrap_or_else(|| "Solo.out/".to_string()); + let dir = params.output_path(&solo_dir); + std::fs::create_dir_all(&dir).map_err(|e| Error::io(e, &dir))?; + write_stats_file( + &dir.join("Barcodes.stats"), + &crate::solo::BarcodeStats::NAMES, + &ctx.stats.values(), + ) +} + +/// The called-cell figures of one gene feature (`SoloFeature_cellFiltering.cpp:75-125`). +struct CellSummary { + n_cells: u64, + reads_in_cells: u64, + mean_reads: u64, + median_reads: u64, + umis_in_cells: u64, + mean_umis: u64, + median_umis: u64, + mean_genes: u64, + median_genes: u64, + genes_detected: u64, +} + +impl CellSummary { + /// STAR's arithmetic: integer means, and the element at `nCells/2` for the + /// medians, of the called cells' reads and genes sorted ascending, and of + /// every barcode's UMIs sorted descending for the UMI median. + fn new(cells: &[CellStat], called: &[u32], genes_detected: u64) -> Self { + let called: std::collections::HashSet = called.iter().copied().collect(); + let inc: Vec<&CellStat> = cells.iter().filter(|c| called.contains(&c.cb)).collect(); + let n = inc.len() as u64; + if n == 0 { + return Self { + n_cells: 0, + reads_in_cells: 0, + mean_reads: 0, + median_reads: 0, + umis_in_cells: 0, + mean_umis: 0, + median_umis: 0, + mean_genes: 0, + median_genes: 0, + genes_detected: 0, + }; } - }; + let mut reads: Vec = inc.iter().map(|c| c.n_reads).collect(); + let mut genes: Vec = inc.iter().map(|c| u64::from(c.n_genes)).collect(); + reads.sort_unstable(); + genes.sort_unstable(); + let mut umis_desc: Vec = cells.iter().map(|c| c.n_umis).collect(); + umis_desc.sort_unstable_by(|a, b| b.cmp(a)); + let reads_in_cells: u64 = reads.iter().sum(); + let umis_in_cells: u64 = inc.iter().map(|c| c.n_umis).sum(); + let mid = (n / 2) as usize; + Self { + n_cells: n, + reads_in_cells, + mean_reads: reads_in_cells / n, + median_reads: reads[mid], + umis_in_cells, + mean_umis: umis_in_cells / n, + median_umis: umis_desc.get(mid).copied().unwrap_or(0), + mean_genes: genes.iter().sum::() / n, + median_genes: genes[mid], + genes_detected, + } + } +} + +/// Run-wide inputs of every feature's `Summary.csv`. +struct SummaryContext<'a> { + n_reads: u64, + mapped_unique: u64, + mapped_multi: u64, + barcodes: &'a crate::solo::BarcodeStats, + rna_q30: (u64, u64), + multimappers: bool, +} +/// `/Summary.csv`, STAR's `SoloFeature::statsOutput` +/// (`SoloFeature_statsOutput.cpp:6-85`). +fn write_summary( + path: &Path, + feature: &str, + run: &SummaryContext<'_>, + stats: &crate::solo::FeatureStats, + cells: Option<&CellSummary>, +) -> Result<(), Error> { + use crate::solo::FeatureStats as F; use std::fmt::Write as _; + let n = run.n_reads as f64; + let ratio = |num: u64, den: u64| num as f64 / den as f64; let mut out = String::new(); let mut row = |k: &str, v: String| { let _ = writeln!(out, "{k},{v}"); }; - row("Number of Reads", total_reads.to_string()); + row("Number of Reads", run.n_reads.to_string()); row( "Reads With Valid Barcodes", - format!("{:.6}", frac(valid_barcodes)), - ); - row("Sequencing Saturation", format!("{saturation:.6}")); - row( - "Reads Mapped to Genome: Unique+Multiple", - format!("{:.6}", frac(mapped_unique + mapped_multi)), + fmt_double(1.0 - (run.barcodes.invalid() + stats.invalid()) as f64 / n), ); row( - "Reads Mapped to Genome: Unique", - format!("{:.6}", frac(mapped_unique)), - ); - row( - &format!("Reads Mapped to {feature_name}: Unique {feature_name}"), - format!("{:.6}", frac(feature_mapped)), - ); - row("Estimated Number of Cells", n_cells.to_string()); - row( - &format!("Unique Reads in Cells Mapped to {feature_name}"), - reads_in_cells.to_string(), - ); - row( - "Fraction of Unique Reads in Cells", - format!( - "{:.6}", - if total_reads_counted > 0 { - reads_in_cells as f64 / total_reads_counted as f64 - } else { - 0.0 - } + "Sequencing Saturation", + fmt_double( + 1.0 - ratio( + stats.get(F::YES_UMIS), + stats.get(F::YES_SUB_WL_MATCH_UNIQUE_FEATURE), + ), ), ); - row("Mean Reads per Cell", mean(reads_in_cells).to_string()); + let (q30, total) = run.barcodes.q30(); + row("Q30 Bases in CB+UMI", fmt_double(ratio(q30, total))); row( - "Median Reads per Cell", - median_sorted(&reads_sorted).to_string(), + "Q30 Bases in RNA read", + fmt_double(ratio(run.rna_q30.0, run.rna_q30.1)), ); - row("UMIs in Cells", umis_in_cells.to_string()); - row("Mean UMI per Cell", mean(umis_in_cells).to_string()); row( - "Median UMI per Cell", - median_sorted(&umis_sorted).to_string(), + "Reads Mapped to Genome: Unique+Multiple", + fmt_double((run.mapped_unique + run.mapped_multi) as f64 / n), ); row( - &format!("Mean {feature_name} per Cell"), - mean(genes_sorted.iter().sum()).to_string(), + "Reads Mapped to Genome: Unique", + fmt_double(run.mapped_unique as f64 / n), ); row( - &format!("Median {feature_name} per Cell"), - median_sorted(&genes_sorted).to_string(), + &format!("Reads Mapped to {feature}: Unique+Multiple {feature}"), + if run.multimappers { + fmt_double(stats.get(F::YES_WL_MATCH) as f64 / n) + } else { + "NoMulti".to_string() + }, ); row( - &format!("Total {feature_name} Detected"), - mstats.genes_detected.to_string(), + &format!("Reads Mapped to {feature}: Unique {feature}"), + fmt_double(stats.get(F::YES_SUB_WL_MATCH_UNIQUE_FEATURE) as f64 / n), ); - + if let Some(c) = cells { + row("Estimated Number of Cells", c.n_cells.to_string()); + row( + &format!("Unique Reads in Cells Mapped to {feature}"), + c.reads_in_cells.to_string(), + ); + row( + "Fraction of Unique Reads in Cells", + fmt_double(ratio( + c.reads_in_cells, + stats.get(F::YES_SUB_WL_MATCH_UNIQUE_FEATURE), + )), + ); + row("Mean Reads per Cell", c.mean_reads.to_string()); + row("Median Reads per Cell", c.median_reads.to_string()); + row("UMIs in Cells", c.umis_in_cells.to_string()); + row("Mean UMI per Cell", c.mean_umis.to_string()); + row("Median UMI per Cell", c.median_umis.to_string()); + row( + &format!("Mean {feature} per Cell"), + c.mean_genes.to_string(), + ); + row( + &format!("Median {feature} per Cell"), + c.median_genes.to_string(), + ); + row( + &format!("Total {feature} Detected"), + c.genes_detected.to_string(), + ); + } std::fs::write(path, out).map_err(|e| Error::io(e, path))?; Ok(()) } +/// `/UMIperCellSorted.txt`: every barcode's UMIs, descending, up to +/// the first zero (`SoloFeature_statsOutput.cpp:74-82`). +fn write_umi_per_cell_sorted(path: &Path, cells: &[CellStat]) -> Result<(), Error> { + use std::fmt::Write as _; + let mut umis: Vec = cells.iter().map(|c| c.n_umis).filter(|&u| u > 0).collect(); + umis.sort_unstable_by(|a, b| b.cmp(a)); + let mut out = String::new(); + for u in umis { + let _ = writeln!(out, "{u}"); + } + std::fs::write(path, out).map_err(|e| Error::io(e, path)) +} + +/// Genes with a count in any of the `called` cells, read off the streamed +/// matrix body (`gene+1 cb+1 count` lines). +fn genes_in_cells(body: &tempfile::NamedTempFile, called: &[u32]) -> Result { + use std::io::BufRead; + let called: std::collections::HashSet = called.iter().map(|&c| c + 1).collect(); + let f = std::fs::File::open(body.path()).map_err(|e| Error::io(e, body.path()))?; + let mut genes = std::collections::HashSet::new(); + for line in std::io::BufReader::new(f).lines() { + let line = line.map_err(|e| Error::io(e, body.path()))?; + let mut it = line.split_ascii_whitespace(); + if let (Some(g), Some(cb), Some(c)) = (it.next(), it.next(), it.next()) + && c != "0" + && cb.parse::().is_ok_and(|cb| called.contains(&cb)) + { + genes.insert(g.to_string()); + } + } + Ok(genes.len() as u64) +} + /// Write the CellRanger-style positional mapping funnel (exonic / intronic / /// intergenic + antisense, as fractions of total reads) to a SEPARATE additional /// file, so the faithful STARsolo `Summary.csv` is never altered by this rustar @@ -2568,17 +2850,131 @@ fn open_maybe_gz(path: &Path) -> Result, Error> { #[cfg(test)] mod tests { use super::*; - use crate::io::fastq::encode_base; - use crate::solo::whitelist::pack_barcode; + /// `fmt_double` against what a C++ program printed for the same values + /// with the default `ostream` (g++, libstdc++). + #[test] + fn fmt_double_matches_cpp_ostream() { + let cases: [(f64, &str); 16] = [ + (1.0 - 432.0 / 0.0, "-inf"), + (f64::NAN, "-nan"), + (0.05773, "0.05773"), + (0.035_714_285_7, "0.0357143"), + (1.0, "1"), + (0.0, "0"), + (0.00448, "0.00448"), + (448.0 / 448.0, "1"), + (0.000_012_345, "1.2345e-05"), + (1_234_567.0, "1.23457e+06"), + (123_456.0, "123456"), + (0.1591, "0.1591"), + (0.962_140_18, "0.96214"), + (5.0e-5, "5e-05"), + (0.923_089_123, "0.923089"), + (0.5, "0.5"), + ]; + for (v, want) in cases { + assert_eq!(fmt_double(v), want, "{v}"); + } + } + + /// STAR's `Barcodes.stats` bins: one per final `cbMatch` code, codes it + /// does not list counted as multiple matches (`addStats`' `default`). #[test] - fn median_sorted_odd_even_empty() { - assert_eq!(median_sorted(&[]), 0); - assert_eq!(median_sorted(&[5]), 5); - assert_eq!(median_sorted(&[1, 2, 3]), 2); - assert_eq!(median_sorted(&[10, 20, 30, 40]), 25); // midpoint(20,30) + fn barcode_stats_follow_stars_switch() { + let st = crate::solo::BarcodeStats::new(); + for code in [0, 0, 1, 3, -1, -2, -3, -11, -12, -23, -24, -99] { + st.record(code); + } + st.no_adapter(); + st.no_umi(); + assert_eq!(st.values(), [1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 2]); + // numInvalidBarcodes: the first nine. + assert_eq!(st.invalid(), 9); } + /// STAR's called-cell figures: integer means; medians are the element at + /// `nCells/2` (no averaging), and the UMI median indexes every barcode's + /// UMIs sorted descending. + #[test] + fn cell_summary_uses_stars_arithmetic() { + let cell = |cb, n_reads, n_umis, n_genes| CellStat { + cb, + n_reads, + n_umis, + n_genes, + }; + let cells = [ + cell(0, 10, 9, 3), + cell(1, 4, 4, 2), + cell(2, 7, 5, 1), + cell(3, 2, 20, 1), + ]; + // Called: 0, 1, 2 (cell 3 has the most UMIs but is not called). + let c = CellSummary::new(&cells, &[0, 1, 2], 4); + assert_eq!((c.n_cells, c.reads_in_cells, c.mean_reads), (3, 21, 7)); + assert_eq!(c.median_reads, 7); // [4, 7, 10][1] + assert_eq!((c.umis_in_cells, c.mean_umis), (18, 6)); + assert_eq!(c.median_umis, 9); // all barcodes descending [20, 9, 5, 4][1] + assert_eq!((c.mean_genes, c.median_genes, c.genes_detected), (2, 2, 4)); + } + + /// STAR's `countVelocyto` collapse: a molecule keeps the transcripts all + /// its reads share, their types OR-ed; an emptied molecule stays empty; a + /// multi-gene molecule and a read with no readInfo are dropped. + #[test] + fn count_velocyto_collapses_molecules_as_star() { + use crate::quant::splice_status::AlignVsTranscript as A; + use crate::solo::{GeneModel, ReadInfo, VelocytoRecord, VelocytoRecords}; + let exon = A::Concordant.type_bits(); + let intron = A::Intron.type_bits(); + // Transcripts 0 and 1 belong to gene 0, transcript 2 to gene 1. + let model = GeneModel::with_tr_genes(vec![0, 0, 1]); + let rec = |read_index: u32, types: &[(u32, u8)]| VelocytoRecord { + read_index, + types: types.to_vec(), + cbs: if read_index == 7 { + vec![7, 70] + } else { + vec![5] + }, + }; + let info = |cb: u32, umi: u64| ReadInfo { cb, umi }; + let read_info = vec![ + info(5, 1), // 0, 1: molecule A, reads intersect on transcript 0 + info(5, 1), + info(5, 2), // 2, 3: molecule B, no shared transcript + info(5, 2), + info(5, 2), // 4: B again, ignored once B is empty + info(5, 3), // 5: molecule C, two genes + ReadInfo::default(), // 6: no cell or UMI + info(7, 1), // 7: another cell, exonic in both isoforms + ]; + let mut records = VelocytoRecords::default(); + for r in [ + rec(0, &[(0, exon), (1, intron)]), + rec(1, &[(0, intron)]), + rec(2, &[(0, exon)]), + rec(3, &[(1, exon)]), + rec(4, &[(0, exon)]), + rec(5, &[(0, exon), (2, exon)]), + rec(6, &[(0, exon)]), + rec(7, &[(0, exon), (1, exon)]), + ] { + records.push(&r); + } + let counts = count_velocyto(&records, &read_info, &model); + // A: transcript 0, exon | intron, with no span → mixed → unspliced. + // B dropped (empty), C dropped (two genes), read 6 skipped. + assert_eq!(counts, vec![(5, 0, [0, 1, 0]), (7, 0, [1, 0, 0])]); + // Detected cells: every barcode a record matched, a 1MM_multi + // read's candidates included (70 is in no readInfo). + let detected: Vec = (0..128).filter(|&cb| records.detected(cb)).collect(); + assert_eq!(detected, vec![5, 7, 70]); + } + use crate::io::fastq::encode_base; + use crate::solo::whitelist::pack_barcode; + #[test] fn distribute_multi_methods() { // Unique counts: gene 0 has 4, gene 1 has none. One ambiguous molecule @@ -2618,10 +3014,23 @@ mod tests { let cells = vec![mk(5, 1000), mk(2, 900), mk(8, 50), mk(1, 40)]; let s = |v: &[&str]| v.iter().map(ToString::to_string).collect::>(); - // TopCells 2: the two highest-UMI cells (cb 5, 2), returned cb-ascending. + // TopCells 2: every cell with at least the third barcode's 50 UMIs + // (cb 5, 2, 8), returned cb-ascending. assert_eq!( called_cells(&cells, &s(&["TopCells", "2"])).unwrap(), - vec![2, 5] + vec![2, 5, 8] + ); + // TopCells is STAR's threshold, not a cut: the UMIs of the N+1th + // barcode, ties included (`SoloFeature_cellFiltering.cpp:46-54`). + let tied = vec![mk(5, 1000), mk(2, 900), mk(8, 900), mk(1, 40)]; + assert_eq!( + called_cells(&tied, &s(&["TopCells", "1"])).unwrap(), + vec![2, 5, 8] + ); + // N past the end: the last barcode's UMIs, i.e. every cell. + assert_eq!( + called_cells(&cells, &s(&["TopCells", "10"])).unwrap(), + vec![1, 2, 5, 8] ); // None: no filtered output. assert!(called_cells(&cells, &s(&["None"])).is_none()); diff --git a/src/solo/gene.rs b/src/solo/gene.rs index ff05ce5..5b702bb 100644 --- a/src/solo/gene.rs +++ b/src/solo/gene.rs @@ -107,6 +107,21 @@ pub struct GeneModel { } impl GeneModel { + /// Solo gene index of transcript `tr`; `u32::MAX` if its gene is not listed. + pub fn tr_gene(&self, tr: usize) -> u32 { + self.tr_gene[tr] + } + + /// A model that knows only each transcript's gene, for tests of the + /// per-molecule steps. + #[cfg(test)] + pub(crate) fn with_tr_genes(tr_gene: Vec) -> Self { + Self { + transcripts: TranscriptomeIndex::default(), + tr_gene, + } + } + pub fn new(transcripts: TranscriptomeIndex, gene_ann: &GeneAnnotation) -> Self { let pos: std::collections::HashMap<&str, u32> = gene_ann .gene_ids @@ -329,41 +344,91 @@ pub fn gene_full_sets( sets } -/// RNA-velocity read category (Sullivan et al. 2025 mature/nascent/ambiguous, -/// reported as scVelo's spliced/unspliced/ambiguous). -#[derive(Debug, Clone, Copy, PartialEq, Eq)] -pub enum VelocytoCategory { - /// Spans an exon–exon junction → processed (mature) mRNA. - Spliced, - /// No junction, but a block extends into an intron → nascent mRNA. - Unspliced, - /// No junction, all blocks wholly within exons → origin indistinguishable. - Ambiguous, -} - -/// Classify a uniquely-mapped read (assigned to gene `g` by gene-body overlap) -/// into its velocity category from the alignment: a splice in the CIGAR means -/// the read is mature; otherwise an aligned block that leaves the exons (into an -/// intron) means nascent; a wholly-exonic block is ambiguous. -pub fn velocyto_category( - transcripts: &[Transcript], - gene_ann: &GeneAnnotation, - g: u32, -) -> VelocytoCategory { - if transcripts.iter().any(|t| t.n_junction > 0) { - return VelocytoCategory::Spliced; - } - let g = g as usize; - let all_exonic = transcripts.iter().all(|t| { - t.exons - .iter() - .all(|e| gene_ann.block_is_exonic(g, e.genome_start, e.genome_end)) - }); - if all_exonic { - VelocytoCategory::Ambiguous - } else { - VelocytoCategory::Unspliced +/// The Velocyto annotation of one uniquely mapped alignment (STAR's +/// `trVelocytoType`, `Transcriptome_classifyAlign.cpp:227-258`; strand +/// filter `:207-209`): for every transcript that contains it on a strand +/// `--soloStrand` keeps, the transcript (an index into `model.transcripts`) +/// and STAR's `reAnn1` type bits from `alignToTranscriptMinOverlap` with +/// velocyto's 5-base flank. +/// A spliced alignment leaves out a transcript whose intron one of its blocks +/// runs into (its junctions are not compared with the transcript's). Sorted +/// by transcript, as `SoloReadFeature_record.cpp:164-165` sorts it. +/// +/// The alignment's blocks are its exons merged across indels, mate by mate, +/// in transcript order (STAR expands a block over `canonSJ` -1/-2, not the +/// mate gap -3); `sjYes` is whether any mate crosses a junction. +pub fn velocyto_transcript_types( + aln: &Transcript, + model: &GeneModel, + strand: SoloStrand, +) -> Vec<(u32, u8)> { + use crate::quant::splice_status::{MIN_OVERLAP_MINUS_ONE, align_to_transcript_min_overlap}; + let idx = &model.transcripts; + let mut out = Vec::new(); + let (Some(first), Some(last)) = (aln.exons.first(), aln.exons.last()) else { + return out; + }; + if idx.tr_starts_sorted.is_empty() { + return out; + } + // Blocks: exons merged across indels within a mate. + let gaps = block_gaps(aln); + let mut blocks: Vec<(u64, u64)> = Vec::with_capacity(aln.exons.len()); + let mut has_junction = false; + for (i, e) in aln.exons.iter().enumerate() { + let joins = i > 0 + && e.i_frag == aln.exons[i - 1].i_frag + && gaps.get(i - 1) != Some(&BlockGap::Junction); + if i > 0 && e.i_frag == aln.exons[i - 1].i_frag && !joins { + has_junction = true; + } + match blocks.last_mut() { + Some(b) if joins => b.1 = e.genome_end - 1, + _ => blocks.push((e.genome_start, e.genome_end - 1)), + } + } + // STAR's transcript walk (`binarySearch1a`, then back while `trEmax` covers + // the alignment end), with the `--soloStrand` filter applied before the + // Velocyto call, as in `classifyAlign`. + let pp = idx + .tr_starts_sorted + .partition_point(|&s| s <= first.genome_start); + if pp == 0 { + return out; + } + let a_g_end = last.genome_end - 1; + let read_str = u32::from(aln.is_reverse); + let mut k = pp - 1; + loop { + let tr = idx.tr_order[k]; + let rel_str = if idx.tr_strand[tr] == 1 { + read_str + } else { + 1 - read_str + }; + let strand_ok = match strand { + SoloStrand::Unstranded => true, + SoloStrand::Forward => rel_str == 0, + SoloStrand::Reverse => rel_str == 1, + }; + if a_g_end < idx.tr_end[tr] + && strand_ok + && let Some(status) = align_to_transcript_min_overlap( + &blocks, + has_junction, + &idx.tr_exons[tr], + MIN_OVERLAP_MINUS_ONE, + ) + { + out.push((tr as u32, status.type_bits())); + } + if !(idx.tr_end_max_sorted[k] > a_g_end && k > 0) { + break; + } + k -= 1; } + out.sort_unstable_by_key(|&(tr, _)| tr); + out } /// CellRanger-style positional region of a uniquely-mapped read (independent of @@ -878,6 +943,94 @@ mod tests { assert_eq!(sets.per_align, vec![Vec::::new()]); } + /// Velocyto's per-read transcript types: only transcripts on the strand + /// `--soloStrand` keeps, each with STAR's type bits. + #[test] + fn velocyto_types_follow_the_strand_and_the_exons() { + use crate::quant::splice_status::AlignVsTranscript as A; + let ann = annotation(); + let model = model_of(&annotation_exons(), &ann); + // G1 (+) is a single exon [100, 200): a read inside it is concordant. + let sense = read_at(120, 180, false); + let types = velocyto_transcript_types(&sense, &model, SoloStrand::Forward); + assert_eq!(types.len(), 1); + assert_eq!(types[0].1, A::Concordant.type_bits()); + // The same read on the other strand: kept only when unstranded. + let anti = read_at(120, 180, true); + assert_eq!( + velocyto_transcript_types(&anti, &model, SoloStrand::Forward), + [] + ); + assert_eq!( + velocyto_transcript_types(&anti, &model, SoloStrand::Unstranded).len(), + 1 + ); + // Outside every transcript: nothing. + let far = read_at(600, 650, false); + assert_eq!( + velocyto_transcript_types(&far, &model, SoloStrand::Unstranded), + [] + ); + } + + /// A spliced alignment leaves out a transcript only if one of its blocks + /// runs into that transcript's intron: STAR does not ask its junctions to + /// be the transcript's (`Transcriptome_classifyAlign.cpp:103-158`). A mate + /// pair's gap is not a junction. + #[test] + fn velocyto_types_drop_transcripts_whose_intron_a_spliced_block_enters() { + use crate::quant::splice_status::AlignVsTranscript as A; + // G3 (+): T1 spliced, exons [500,600) and [800,900); T2 retains the + // intron, [500,900). + let g = genome(); + let mut retained = gtf_exon(501, 900, '+', "G3"); + retained + .attributes + .insert("transcript_id".to_string(), "G3_t2".to_string()); + let exons = vec![ + gtf_exon(501, 600, '+', "G3"), + gtf_exon(801, 900, '+', "G3"), + retained, + ]; + let ann = GeneAnnotation::from_gtf_exons(&exons, &g); + let model = model_of(&exons, &ann); + let bits = |t: &Transcript| { + let mut b: Vec = velocyto_transcript_types(t, &model, SoloStrand::Forward) + .into_iter() + .map(|(_, ty)| ty) + .collect(); + b.sort_unstable(); + b + }; + let exon = A::Concordant.type_bits(); + let intron = A::Intron.type_bits(); + // T1's junction, and one neither transcript has: both blocks are + // exonic in both transcripts, so both are kept. + assert_eq!(bits(&spliced_read((550, 600), (800, 850))), [exon, exon]); + assert_eq!(bits(&spliced_read((550, 580), (820, 850))), [exon, exon]); + // A spliced read whose first block runs into T1's intron: T1 is out. + assert_eq!(bits(&spliced_read((550, 620), (820, 850))), [exon]); + // Unspliced, the same block spans T1's exon and intron. + let span = A::ExonIntronSpan.type_bits(); + let mut both = vec![exon, span]; + both.sort_unstable(); + assert_eq!(bits(&read_at(550, 620, false)), both); + // In the intron: intronic in T1, exonic in T2. + let mut both = vec![exon, intron]; + both.sort_unstable(); + assert_eq!(bits(&read_at(650, 700, false)), both); + // A pair with a mate in each exon: the mates' gap is no junction, so + // T2 keeps it, and both transcripts hold it. + use noodles::sam::alignment::record::cigar::{Op, op::Kind}; + let mate = |start: u64, end: u64, is_reverse: bool| { + let mut t = read_at(start, end, is_reverse); + t.cigar = vec![Op::new(Kind::Match, (end - start) as usize)]; + t + }; + let pair = crate::solo::pair_alignment(&mate(550, 590, false), &mate(810, 850, true)); + assert_eq!(bits(&pair), [exon, exon]); + } + /// `GeneFull` gene sets: per alignment, the sense genes whose body it /// overlaps. #[test] diff --git a/src/solo/mod.rs b/src/solo/mod.rs index 66f7d0a..b91d1fb 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -18,13 +18,13 @@ pub mod smartseq; pub mod transcript3p; pub mod whitelist; -pub use count::{UmiDedup, UmiFiltering, write_gene_matrix}; +pub use count::{UmiDedup, UmiFiltering, write_barcode_stats, write_gene_matrix}; pub use gene::{ - GeneAssignment, GeneModel, GeneSets, Region, SoloFeature, SoloStrand, VelocytoCategory, - assign_gene_se, classify_align, classify_read, gene_full_sets, velocyto_category, + GeneAssignment, GeneModel, GeneSets, ReadClass, Region, SoloFeature, SoloStrand, + assign_gene_se, classify_align, classify_read, gene_full_sets, velocyto_transcript_types, }; pub use whitelist::{ - CbCandidate, CbMatch, CbMatchStats, CbMatchType, CbWhitelist, UmiCheck, check_umi, pack_barcode, + BarcodeStats, CbCandidate, CbMatch, CbMatchType, CbWhitelist, UmiCheck, check_umi, pack_barcode, }; use crate::align::transcript::Transcript; @@ -269,25 +269,28 @@ impl SoloBarcodeLayout { } /// Extract the CB (concatenating segments for `Complex`) and UMI from one - /// barcode read. `None` if the read is shorter than [`Self::min_read_len`]. - pub fn extract(&self, barcode_read: &EncodedRead) -> Option { + /// barcode read, or why it cannot be (STAR's `getCBandUMI`: no adapter, + /// then no UMI, then no CB). + pub fn extract(&self, barcode_read: &EncodedRead) -> Result { let seq = &barcode_read.sequence; let qual = &barcode_read.quality; - if seq.len() < self.min_read_len() { - return None; - } match self { Self::Simple { cb_start, cb_len, umi_start, umi_len, - } => Some(CellBarcode { - cb_seq: seq[*cb_start..cb_start + cb_len].to_vec(), - cb_qual: slice_or_empty(qual, *cb_start, *cb_len), - umi_seq: seq[*umi_start..umi_start + umi_len].to_vec(), - umi_qual: slice_or_empty(qual, *umi_start, *umi_len), - }), + } => { + if seq.len() < self.min_read_len() { + return Err(BarcodeMiss::NoCb); + } + Ok(CellBarcode { + cb_seq: seq[*cb_start..cb_start + cb_len].to_vec(), + cb_qual: slice_or_empty(qual, *cb_start, *cb_len), + umi_seq: seq[*umi_start..umi_start + umi_len].to_vec(), + umi_qual: slice_or_empty(qual, *umi_start, *umi_len), + }) + } Self::Anchored { cb_segments, umi, @@ -299,7 +302,7 @@ impl SoloBarcodeLayout { // whatever the anchors; then the UMI, then each CB segment. let (dist, adapter_start) = local_align_hamming_dist(seq, adapter); if dist > *max_mismatches { - return None; + return Err(BarcodeMiss::NoAdapter); } let span = |(a, b): &(Anchored, Anchored)| { barcode_span( @@ -308,15 +311,15 @@ impl SoloBarcodeLayout { seq.len(), ) }; - let (us, ul) = span(umi)?; + let (us, ul) = span(umi).ok_or(BarcodeMiss::NoUmi)?; let mut cb_seq = Vec::new(); let mut cb_qual = Vec::new(); for segment in cb_segments { - let (s0, l) = span(segment)?; + let (s0, l) = span(segment).ok_or(BarcodeMiss::NoCb)?; cb_seq.extend_from_slice(&seq[s0..s0 + l]); cb_qual.extend_from_slice(&slice_or_empty(qual, s0, l)); } - Some(CellBarcode { + Ok(CellBarcode { cb_seq, cb_qual, umi_seq: seq[us..us + ul].to_vec(), @@ -324,13 +327,19 @@ impl SoloBarcodeLayout { }) } Self::Complex { cb_segments, umi } => { + if seq.len() < umi.0 + umi.1 { + return Err(BarcodeMiss::NoUmi); + } + if seq.len() < self.min_read_len() { + return Err(BarcodeMiss::NoCb); + } let mut cb_seq = Vec::new(); let mut cb_qual = Vec::new(); for &(s, l) in cb_segments { cb_seq.extend_from_slice(&seq[s..s + l]); cb_qual.extend_from_slice(&slice_or_empty(qual, s, l)); } - Some(CellBarcode { + Ok(CellBarcode { cb_seq, cb_qual, umi_seq: seq[umi.0..umi.0 + umi.1].to_vec(), @@ -341,6 +350,49 @@ impl SoloBarcodeLayout { } } +/// Why no barcode could be read, as STAR counts it in `Barcodes.stats`. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum BarcodeMiss { + /// No adapter within `--soloAdapterMismatchesNmax` (`noNoAdapter`). + NoAdapter, + /// The UMI cannot be extracted (`noNoUMI`). + NoUmi, + /// A cell barcode cannot be extracted (`noNoCB`, `cbMatch` -11). + NoCb, +} + +/// The `(Q30 bases, bases)` of a barcode read that STAR adds to its CB+UMI +/// quality histogram (`SoloReadBarcode_getCBandUMI.cpp:243-266`): the CB and +/// UMI qualities for `CB_UMI_Simple`, else the whole barcode read (its first +/// `--soloBarcodeReadLength` bases when that is set), extracted or not. +pub fn barcode_q30( + layout: &SoloBarcodeLayout, + read: &EncodedRead, + bc: Result<&CellBarcode, BarcodeMiss>, + read_len: Option, +) -> (u64, u64) { + let q30 = |q: &[u8]| { + ( + q.iter().filter(|&&b| b >= 33 + 30).count() as u64, + q.len() as u64, + ) + }; + match (layout, bc) { + (SoloBarcodeLayout::Simple { .. }, Ok(bc)) => { + let (a, b) = q30(&bc.cb_qual); + let (c, d) = q30(&bc.umi_qual); + (a + c, b + d) + } + (SoloBarcodeLayout::Simple { .. }, Err(_)) => (0, 0), + _ => { + let n = read_len + .unwrap_or(read.quality.len()) + .min(read.quality.len()); + q30(&read.quality[..n]) + } + } +} + fn slice_or_empty(data: &[u8], start: usize, len: usize) -> Vec { if start + len <= data.len() { data[start..start + len].to_vec() @@ -534,8 +586,10 @@ pub struct SoloReadReader { /// One cDNA read paired with its (optional) extracted barcode. pub struct SoloRead { pub cdna: EncodedRead, - /// `None` when the barcode read was too short to extract CB+UMI. - pub barcode: Option, + /// The extracted barcode, or why there is none. + pub barcode: Result, + /// The barcode read's `(Q30 bases, bases)` for `Summary.csv` ([`barcode_q30`]). + pub barcode_q30: (u64, u64), /// The whole barcode read, kept only for the `sS`/`sQ` SAM tags. pub barcode_read: Option, } @@ -568,10 +622,17 @@ impl SoloReadReader { (Some(cdna), Some(mut bc)) => { self.barcode_len.apply(&mut bc)?; let barcode = self.layout.extract(&bc); + let barcode_q30 = barcode_q30( + &self.layout, + &bc, + barcode.as_ref().map_err(|m| *m), + self.barcode_len.expected, + ); let barcode_read = self.keep_barcode_read.then_some(bc); Ok(Some(SoloRead { cdna, barcode, + barcode_q30, barcode_read, })) } @@ -636,8 +697,10 @@ pub fn open_reader(params: &Parameters) -> Result { pub struct SoloPairedRead { pub mate1: EncodedRead, pub mate2: EncodedRead, - /// `None` when the barcode read (or mate 1) was too short to extract CB+UMI. - pub barcode: Option, + /// The extracted barcode, or why there is none. + pub barcode: Result, + /// The barcode read's `(Q30 bases, bases)` for `Summary.csv` ([`barcode_q30`]). + pub barcode_q30: (u64, u64), /// The whole barcode read, kept only for the `sS`/`sQ` SAM tags. `None` for /// a `--soloBarcodeMate 1` run, where mate 1 itself is the barcode read. pub barcode_read: Option, @@ -682,26 +745,39 @@ impl SoloPairedReader { pub fn next_read(&mut self) -> Result, Error> { match (self.mate1.next_encoded()?, self.mate2.next_encoded()?) { (Some(mate1), Some(mate2)) => { - let (barcode, barcode_read) = match &mut self.barcode { - // Separate barcode read: must stay in lockstep with the mates. - Some(reader) => { - let mut bc = reader.next_encoded()?.ok_or_else(|| { - Error::from(std::io::Error::new( - std::io::ErrorKind::UnexpectedEof, - "solo: barcode read file has fewer reads than the cDNA mate files", - )) - })?; - self.barcode_len.apply(&mut bc)?; - let extracted = self.layout.extract(&bc); - (extracted, self.keep_barcode_read.then_some(bc)) - } - // `--soloBarcodeMate 1`: the barcode is a prefix of mate 1. - None => (self.layout.extract(&mate1), None), + // A separate barcode read must stay in lockstep with the mates; + // with `--soloBarcodeMate 1` the barcode is a prefix of mate 1. + let (barcode, barcode_q30, barcode_read) = if let Some(reader) = &mut self.barcode { + let mut bc = reader.next_encoded()?.ok_or_else(|| { + Error::from(std::io::Error::new( + std::io::ErrorKind::UnexpectedEof, + "solo: barcode read file has fewer reads than the cDNA mate files", + )) + })?; + self.barcode_len.apply(&mut bc)?; + let extracted = self.layout.extract(&bc); + let q = barcode_q30( + &self.layout, + &bc, + extracted.as_ref().map_err(|m| *m), + self.barcode_len.expected, + ); + (extracted, q, self.keep_barcode_read.then_some(bc)) + } else { + let extracted = self.layout.extract(&mate1); + let q = barcode_q30( + &self.layout, + &mate1, + extracted.as_ref().map_err(|m| *m), + None, + ); + (extracted, q, None) }; Ok(Some(SoloPairedRead { mate1, mate2, barcode, + barcode_q30, barcode_read, })) } @@ -990,13 +1066,45 @@ pub struct SoloCountRecord { pub read_index: u32, } -/// `read_index` of a record from a run that does not need STAR's readInfo. +/// `read_index` of a record from a run that does not need STAR's readInfo +/// (before [`SoloReadOutcome::set_read_index`] stamps it). pub const NO_READ_INDEX: u32 = u32::MAX; /// `gene` of a no-feature record ([`FeatureOutcome::no_feature`]), STAR's /// feature `-1`. pub const NO_FEATURE: u32 = u32::MAX; +/// The SJ feature's records: one per junction a uniquely mapped read crosses, +/// as STAR writes one line per junction (`outputReadCB`, +/// `SoloReadFeature_record.cpp:253-262`). +#[derive(Debug, Default)] +pub struct SjRecords { + pub unique: Vec, + /// Records of 1MM_multi barcodes, resolved by the posterior when the matrix + /// is built (`SoloReadFeature_inputRecords.cpp:66-104`). + pub multi: Vec, +} + +impl SjRecords { + pub fn is_empty(&self) -> bool { + self.unique.is_empty() && self.multi.is_empty() + } + + pub fn append(&mut self, other: SjRecords) { + self.unique.extend(other.unique); + self.multi.extend(other.multi); + } +} + +/// An SJ record whose barcode matched several whitelist entries by 1MM. +#[derive(Debug, Clone, PartialEq, Eq)] +pub struct SjMultiRecord { + pub candidates: Vec, + pub umi: u64, + pub intron_start: u64, + pub intron_end: u64, +} + /// One (cell, UMI, splice-junction) observation for the `SJ` feature. The /// junction is identified by its absolute intron coordinates; it is mapped to a /// matrix row (the `SJ.out.tab` order) at output time. @@ -1006,16 +1114,78 @@ pub struct SjCountRecord { pub umi: u64, pub intron_start: u64, pub intron_end: u64, + /// The read's barcode matched the whitelist exactly (`cbMatch` 0), for + /// `yessubWLmatchExact`. + pub exact: bool, } -/// One (cell, UMI, gene) observation for the `Velocyto` feature, tagged with the -/// read's spliced/unspliced/ambiguous category. -#[derive(Debug, Clone, Copy, PartialEq, Eq)] +/// One read's record for the `Velocyto` feature (STAR writes `iRead nTr tr +/// type ...`, `SoloReadFeature_record.cpp:160-171`): the read, and the +/// transcripts its unique alignment is compatible with, each with STAR's type +/// bits ([`velocyto_transcript_types`]). Its cell and UMI come later, from the +/// Gene feature's readInfo. +#[derive(Debug, Clone, PartialEq, Eq)] pub struct VelocytoRecord { - pub cb: u32, - pub umi: u64, - pub gene: u32, - pub category: VelocytoCategory, + pub read_index: u32, + pub types: Vec<(u32, u8)>, + /// The whitelist entries the read's barcode matched (several for an + /// unresolved 1MM_multi barcode): STAR counts each as a cell this feature + /// detected (`cbReadCount`, `SoloReadFeature_record.cpp:187-189`), which + /// decides the cells of the filtered matrices. + pub cbs: Vec, +} + +/// The Velocyto feature's records for the whole run, kept compact: STAR +/// writes them to a temporary file, and a run can have hundreds of millions. +/// Each read is its index and transcript count, and each transcript its index +/// and type bits packed in one `u32` (`tr << 4 | bits`; STAR's transcript +/// index is a `uint32`, and the type bits are the four of +/// [`AlignVsTranscript`](crate::quant::splice_status::AlignVsTranscript)). +/// The cells the records' barcodes matched are a bitset over the barcode +/// index. +#[derive(Debug, Default)] +pub struct VelocytoRecords { + reads: Vec<(u32, u32)>, + types: Vec, + cells: Vec, +} + +impl VelocytoRecords { + pub fn push(&mut self, r: &VelocytoRecord) { + self.reads.push((r.read_index, r.types.len() as u32)); + self.types.extend(r.types.iter().map(|&(tr, ty)| { + debug_assert!(tr < 1 << 28 && ty < 1 << 4); + tr << 4 | u32::from(ty) + })); + for &cb in &r.cbs { + let word = cb as usize / 64; + if word >= self.cells.len() { + self.cells.resize(word + 1, 0); + } + self.cells[word] |= 1 << (cb % 64); + } + } + + /// Each read's index and its `(transcript, type bits)`, sorted by + /// transcript. + pub fn iter(&self) -> impl Iterator + '_)> { + let mut start = 0; + self.reads.iter().map(move |&(read_index, n)| { + let packed = &self.types[start..start + n as usize]; + start += n as usize; + ( + read_index, + packed.iter().map(|&p| (p >> 4, (p & 0xf) as u8)), + ) + }) + } + + /// Whether a record's barcode matched cell `cb`. + pub fn detected(&self, cb: u32) -> bool { + self.cells + .get(cb as usize / 64) + .is_some_and(|w| w & (1 << (cb % 64)) != 0) + } } /// A read whose cell barcode matched multiple whitelist entries by 1MM @@ -1089,27 +1259,32 @@ pub struct SoloContext { pub gene_ann: GeneAnnotation, /// Transcript models for the `Gene` feature and GX/GN (`classifyAlign`). pub gene_model: GeneModel, - pub stats: CbMatchStats, + /// `Barcodes.stats` and the CB+UMI Q30 tally. + pub stats: BarcodeStats, + /// `Features.stats` per quantified feature, parallel to `features`. + pub feature_stats: Vec, + /// `Features.stats` of the SJ and Velocyto features. + pub sj_stats: FeatureStats, + pub velocyto_stats: FeatureStats, + /// `(Q30 bases, bases)` of the cDNA reads' clipped sequence, for + /// "Q30 Bases in RNA read" (`ReadAlign_oneRead.cpp:27-33`). + pub rna_q30: [AtomicU64; 2], /// Quantified features (`Gene`, `GeneFull`, …), each with its own recorder /// and `Solo.out//raw/` output. Parallel to `recorders`. pub features: Vec, pub recorders: Vec, - /// Reads uniquely assigned to a gene per feature (parallel to `features`), - /// among valid-barcode reads — the STARsolo "Reads Mapped to : - /// Unique" metric. - pub feature_reads: Vec, /// CellRanger-style positional mapping funnel over uniquely-mapped reads - /// (independent of barcode), populated only when both `Gene` and `GeneFull` - /// features run. - pub region_stats: RegionStats, + /// (independent of barcode), kept only for `--soloCellRangerSummary Yes` + /// and written only when both `Gene` and `GeneFull` features run. + pub region_stats: Option, /// `--soloFeatures SJ`: collect per-cell splice-junction counts. pub sj_enabled: bool, /// (cell, UMI, junction) observations for the SJ feature. - pub sj_records: Mutex>, + pub sj_records: Mutex, /// `--soloFeatures Velocyto`: collect spliced/unspliced/ambiguous counts. pub velocyto_enabled: bool, - /// (cell, UMI, gene, category) observations for the Velocyto feature. - pub velocyto_records: Mutex>, + /// The Velocyto feature's per-read records. + pub velocyto_records: Mutex, /// `--soloMultiMappers` includes a non-`Unique` method → capture gene- /// ambiguous reads for distribution into `UniqueAndMult-*.mtx`. pub want_multi: bool, @@ -1129,7 +1304,14 @@ pub struct SoloContext { /// UMI each read ended up counted under, filled by UMI collapsing and read /// back when the sorted BAM is written to fill `CB`/`UB`. `None` unless one /// of those tags was requested. - pub read_info: Option>>, + /// + /// One slot per feature, as STAR keeps a readInfo per feature: the first + /// feature's for `CB`/`UB` (`samAttrFeature`), the Gene feature's for + /// Velocyto (`ParametersSolo.cpp:420-422`). + pub read_info: Vec>>>, + /// The slot of `read_info` the `CB`/`UB` SAM tags read (the first + /// feature's), when those tags were requested. + pub tag_read_info_slot: Option, } /// One `readInfo` entry: what a read was counted as, once collapsing has run. @@ -1150,6 +1332,54 @@ impl Default for ReadInfo { } } +/// STAR's `SoloReadFeatureStats` (`SoloReadFeatureStats.h`) for one feature, +/// written to `Solo.out//Features.stats` and behind `Summary.csv`. +#[derive(Debug, Default)] +pub struct FeatureStats { + counts: [AtomicU64; FeatureStats::NAMES.len()], +} + +impl FeatureStats { + /// STAR's names, in its order. + pub const NAMES: [&'static str; 11] = [ + "noUnmapped", + "noNoFeature", + "MultiFeature", + "subMultiFeatureMultiGenomic", + "noTooManyWLmatches", + "noMMtoWLwithoutExact", + "yesWLmatch", + "yessubWLmatchExact", + "yessubWLmatch_UniqueFeature", + "yesCellBarcodes", + "yesUMIs", + ]; + pub const NO_UNMAPPED: usize = 0; + pub const NO_NO_FEATURE: usize = 1; + pub const MULTI_FEATURE: usize = 2; + pub const SUB_MULTI_FEATURE_MULTI_GENOMIC: usize = 3; + pub const NO_TOO_MANY_WL_MATCHES: usize = 4; + pub const NO_MM_TO_WL_WITHOUT_EXACT: usize = 5; + pub const YES_WL_MATCH: usize = 6; + pub const YES_SUB_WL_MATCH_EXACT: usize = 7; + pub const YES_SUB_WL_MATCH_UNIQUE_FEATURE: usize = 8; + pub const YES_CELL_BARCODES: usize = 9; + pub const YES_UMIS: usize = 10; + + pub fn add(&self, i: usize, n: u64) { + self.counts[i].fetch_add(n, Ordering::Relaxed); + } + + pub fn get(&self, i: usize) -> u64 { + self.counts[i].load(Ordering::Relaxed) + } + + /// STAR's `numInvalidBarcodes`: reads dropped for their barcode at counting. + pub fn invalid(&self) -> u64 { + self.get(Self::NO_TOO_MANY_WL_MATCHES) + self.get(Self::NO_MM_TO_WL_WITHOUT_EXACT) + } +} + /// Per-region read tallies for the `Summary.csv` mapping funnel (uniquely-mapped /// reads, mirroring CellRanger's "confidently mapped to ... regions"). #[derive(Default)] @@ -1166,8 +1396,9 @@ pub struct RegionStats { pub struct SoloReadOutcome { pub per_feature: Vec, /// SJ-feature records for this read (one per crossed junction); empty unless - /// `--soloFeatures SJ` and the read is uniquely mapped with a resolved CB. - pub sj: Vec, + /// `--soloFeatures SJ` and the read is uniquely mapped with a whitelisted + /// barcode and a valid UMI. + pub sj: SjRecords, /// Velocyto record for this read (resolved CB, gene-assigned), if enabled. pub velocyto: Option, /// Barcode facts for the per-read SAM tags. `None` when the barcode read was @@ -1177,7 +1408,8 @@ pub struct SoloReadOutcome { impl SoloReadOutcome { /// Stamp the input read index onto every count record this read produced, so - /// UMI collapsing can fill STAR's readInfo (the `CB`/`UB` SAM tags). + /// UMI collapsing can fill STAR's readInfo (for the `CB`/`UB` SAM tags and + /// Velocyto). pub fn set_read_index(&mut self, read_index: u32) { for fo in &mut self.per_feature { if let Some(r) = &mut fo.record { @@ -1193,6 +1425,9 @@ impl SoloReadOutcome { m.read_index = read_index; } } + if let Some(v) = &mut self.velocyto { + v.read_index = read_index; + } } } @@ -1278,8 +1513,9 @@ pub struct FeatureOutcome { pub multi_gene: Option, /// STAR's no-feature record (`SoloReadFeature_record.cpp:180-182`): a read /// with a valid barcode that this feature does not count to a single gene, - /// kept only to fill readInfo (the `CB`/`UB` tags). `gene` is [`NO_FEATURE`]; a 1MM_multi - /// barcode goes in `no_feature_multi` instead. + /// kept only to fill readInfo (for the `CB`/`UB` tags or Velocyto). `gene` + /// is [`NO_FEATURE`]; a 1MM_multi barcode goes in `no_feature_multi` + /// instead. pub no_feature: Option, pub no_feature_multi: Option, } @@ -1369,12 +1605,22 @@ impl SoloContext { } }; let recorders = features.iter().map(|_| SoloRecorder::new()).collect(); - 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"); + // STAR's readInfo, per feature; sized once the read count is known + // (`reserve_read_info`). + let read_info = features + .iter() + .enumerate() + .map(|(fi, &f)| { + let tags = fi == 0 && params.solo_read_info_needed(); + let velocyto = velocyto_enabled && f == SoloFeature::Gene; + (tags || velocyto).then(|| Mutex::new(Vec::new())) + }) + .collect(); Ok(Self { layout: SoloBarcodeLayout::from_params(params), whitelist, @@ -1382,15 +1628,18 @@ impl SoloContext { strand, gene_ann, gene_model, - stats: CbMatchStats::new(), + stats: BarcodeStats::new(), + feature_stats: features.iter().map(|_| FeatureStats::default()).collect(), + sj_stats: FeatureStats::default(), + velocyto_stats: FeatureStats::default(), + rna_q30: [AtomicU64::new(0), AtomicU64::new(0)], features, recorders, - feature_reads, - region_stats: RegionStats::default(), + region_stats: params.solo_cellranger_summary().then(RegionStats::default), sj_enabled, - sj_records: Mutex::new(Vec::new()), + sj_records: Mutex::new(SjRecords::default()), velocyto_enabled, - velocyto_records: Mutex::new(Vec::new()), + velocyto_records: Mutex::new(VelocytoRecords::default()), want_multi, cell_read_stats: (params.solo_cell_read_stats == "CB") .then(|| Mutex::new(crate::solo::cell_reads::CellReadStats::new())), @@ -1420,16 +1669,28 @@ impl SoloContext { .transpose()?, transcript3p: transcript3p .then(|| Mutex::new(crate::solo::transcript3p::Transcript3pAcc::new())), - // Sized once the read count is known (`reserve_read_info`). - read_info: params - .solo_read_info_needed() - .then(|| Mutex::new(Vec::new())), + read_info, + tag_read_info_slot: params.solo_read_info_needed().then_some(0), }) } - /// Whether this run has to track STAR's readInfo (the `CB`/`UB` SAM tags). + /// Whether this run has to track STAR's readInfo (for the `CB`/`UB` SAM + /// tags or Velocyto). pub fn read_info_enabled(&self) -> bool { - self.read_info.is_some() + self.read_info.iter().any(Option::is_some) + } + + /// The readInfo the `CB`/`UB` tags are filled from (the first feature's), + /// when those tags were requested. + pub fn tag_read_info(&self) -> Option<&Mutex>> { + self.read_info.get(self.tag_read_info_slot?)?.as_ref() + } + + /// The Gene feature's readInfo, which Velocyto takes its cells and UMIs + /// from. + pub fn gene_read_info(&self) -> Option<&Mutex>> { + let fi = self.features.iter().position(|&f| f == SoloFeature::Gene)?; + self.read_info[fi].as_ref() } /// `--soloType CB_samTagOut`: match the barcode to the whitelist and stop @@ -1441,7 +1702,7 @@ impl SoloContext { let cb_match = self .whitelist .match_cb(&bc.cb_seq, &bc.cb_qual, self.match_type); - self.stats.record_cb(&cb_match); + self.record_barcode_code(cb_match.star_code(), &cb_match); let corrected = match cb_match.resolved_index() { Some(idx) => self .whitelist @@ -1462,7 +1723,7 @@ impl SoloContext { /// Size the readInfo array for `n_reads` input reads, once the alignment /// pass has counted them. No-op unless `CB`/`UB` were requested. pub fn reserve_read_info(&self, n_reads: usize) { - if let Some(info) = &self.read_info { + for info in self.read_info.iter().flatten() { let mut info = info.lock().unwrap(); info.clear(); info.resize(n_reads, ReadInfo::default()); @@ -1531,6 +1792,113 @@ impl SoloContext { .collect() } + /// Count a read with no barcode in `Barcodes.stats`: a missing adapter or + /// UMI directly (`SoloReadBarcode_getCBandUMI.cpp:340-350`), a missing CB + /// as `cbMatch` -11. + pub fn record_barcode_miss(&self, miss: BarcodeMiss) { + match miss { + BarcodeMiss::NoAdapter => self.stats.no_adapter(), + BarcodeMiss::NoUmi => self.stats.no_umi(), + BarcodeMiss::NoCb => self.add_stats(-11), + } + } + + /// STAR's `addStats`: a read's bin by its final `cbMatch`, none without a + /// whitelist (`SoloReadBarcode_getCBandUMI.cpp:95-96`). + fn add_stats(&self, code: i32) { + if !matches!(self.whitelist, CbWhitelist::NoWhitelist { .. }) { + self.stats.record(code); + } + } + + /// Add one read's quality tallies: its barcode read's (`barcode_q30`) and + /// its cDNA mates' clipped qualities (STAR's `qualHist`, bases at or above + /// `readQualityScoreBase + 30`). + pub fn record_qualities(&self, barcode_q30: (u64, u64), mates: &[&[u8]]) { + self.stats.add_q30(barcode_q30.0, barcode_q30.1); + for q in mates { + let q30 = q.iter().filter(|&&b| b >= 33 + 30).count() as u64; + self.rna_q30[0].fetch_add(q30, Ordering::Relaxed); + self.rna_q30[1].fetch_add(q.len() as u64, Ordering::Relaxed); + } + } + + /// Count a read by its final `cbMatch` code (`addStats`), and for an exact + /// match with a valid UMI (code 0), its barcode's exact count. + fn record_barcode_code(&self, code: i32, cb_match: &CbMatch) { + self.add_stats(code); + if code == 0 + && let CbMatch::Exact(idx) = cb_match + { + self.whitelist.bump_exact(*idx); + } + } + + /// STAR's per-read feature counts for a read with a valid barcode and UMI + /// (`SoloReadFeature_record.cpp:75-176`), plus `yessubWLmatchExact` for + /// its counted records (`SoloReadFeature_inputRecords.cpp:111-113`; an + /// exact barcode always survives the later filters). + fn record_feature_stats( + &self, + n_loci: usize, + class: &ReadClass, + out: &SoloReadOutcome, + exact: bool, + junctions: &[(u64, u64)], + ) { + use FeatureStats as F; + for (fi, (&feature, st)) in self.features.iter().zip(&self.feature_stats).enumerate() { + if n_loci == 0 { + st.add(F::NO_UNMAPPED, 1); + continue; + } + let assignment = match feature { + SoloFeature::Gene => class.gene, + SoloFeature::GeneFull => class.gene_full, + }; + match assignment { + GeneAssignment::Gene(_) => { + if exact { + st.add(F::YES_SUB_WL_MATCH_EXACT, 1); + } + } + GeneAssignment::Ambiguous => { + st.add(F::MULTI_FEATURE, 1); + if n_loci > 1 { + st.add(F::SUB_MULTI_FEATURE_MULTI_GENOMIC, 1); + } + let counted = out + .per_feature + .get(fi) + .is_some_and(|fo| fo.multi_gene.is_some()); + if exact && counted { + st.add(F::YES_SUB_WL_MATCH_EXACT, 1); + } + } + GeneAssignment::NoFeature | GeneAssignment::Unmapped => { + st.add(F::NO_NO_FEATURE, 1); + } + } + } + if self.sj_enabled { + let st = &self.sj_stats; + if n_loci == 0 { + st.add(F::NO_UNMAPPED, 1); + } else if n_loci > 1 { + st.add(F::SUB_MULTI_FEATURE_MULTI_GENOMIC, 1); + st.add(F::MULTI_FEATURE, 1); + } else if junctions.is_empty() { + st.add(F::NO_NO_FEATURE, 1); + } + // The rest are counted per junction record that has a matrix row, + // when the matrix is built (`SoloReadFeature_inputRecords.cpp:23-27`). + } + // Velocyto: STAR's `record` counts its reads too, but `countVelocyto` + // never sums the per-thread stats into the feature's (only + // `countCBgeneUMI` calls `addStats`, `SoloFeature_countCBgeneUMI.cpp:47`), + // so Velocyto's Features.stats carries only the UMIs and cells. + } + /// Process one solo read: match the cell barcode, validate the UMI, assign /// a gene, and (on success) produce a count record. Stats are recorded /// here; the returned records are appended to the recorder by the caller. @@ -1538,7 +1906,7 @@ impl SoloContext { &self, cdna_transcripts: &[Transcript], n_loci: usize, - barcode: Option<&CellBarcode>, + barcode: Result<&CellBarcode, BarcodeMiss>, junctions: &[(u64, u64)], ) -> SoloReadOutcome { let mut out = SoloReadOutcome::default(); @@ -1547,8 +1915,7 @@ impl SoloContext { // per-feature gene assignment and the CellRanger-style mapping funnel, so // this is no more work than the old per-feature `assign_gene_se` calls. let want_exon = self.features.contains(&SoloFeature::Gene); - // Velocyto assigns its gene by gene-body overlap, so it needs `want_body`. - let want_body = self.features.contains(&SoloFeature::GeneFull) || self.velocyto_enabled; + let want_body = self.features.contains(&SoloFeature::GeneFull); let class = classify_read( cdna_transcripts, &self.gene_ann, @@ -1563,43 +1930,50 @@ impl SoloContext { // "confidently mapped" = MAPQ 255 ≈ a single alignment), independent of // barcode validity. `n_loci` is the number of genomic loci the read (or // pair) maps to — for SE this equals `cdna_transcripts.len()`. - if n_loci == 1 { + if let Some(funnel) = &self.region_stats + && n_loci == 1 + { match class.region { Some(Region::Exonic) => { - self.region_stats.exonic.fetch_add(1, Ordering::Relaxed); + funnel.exonic.fetch_add(1, Ordering::Relaxed); } Some(Region::Intronic) => { - self.region_stats.intronic.fetch_add(1, Ordering::Relaxed); + funnel.intronic.fetch_add(1, Ordering::Relaxed); } Some(Region::Intergenic) => { - self.region_stats.intergenic.fetch_add(1, Ordering::Relaxed); + funnel.intergenic.fetch_add(1, Ordering::Relaxed); } None => {} } if class.antisense { - self.region_stats.antisense.fetch_add(1, Ordering::Relaxed); + funnel.antisense.fetch_add(1, Ordering::Relaxed); } } - // No barcode read (too short) → nothing to count (region already tallied). - let Some(bc) = barcode else { - return out; + // No barcode → nothing to count. + let bc = match barcode { + Ok(bc) => bc, + Err(miss) => { + self.record_barcode_miss(miss); + return out; + } }; // Cell-barcode match. let cb_match = self .whitelist .match_cb(&bc.cb_seq, &bc.cb_qual, self.match_type); - self.stats.record_cb(&cb_match); // Per-read SAM-tag facts, recorded before any early return so that a // rejected barcode still tags its alignments. STAR reports a UMI // rejection in place of the CB code (`getCBandUMI.cpp:304`). let umi_status = check_umi(&bc.umi_seq); + let code = umi_status + .star_code() + .unwrap_or_else(|| cb_match.star_code()); + self.record_barcode_code(code, &cb_match); out.barcode = Some(SoloBarcodeTags { - cb_match: umi_status - .star_code() - .unwrap_or_else(|| cb_match.star_code()), + cb_match: code, cb_index: cb_match.resolved_index(), }); @@ -1610,46 +1984,63 @@ impl SoloContext { }; // UMI validity. - let umi = match umi_status { - UmiCheck::Ok(packed) => { - self.stats.record_umi(&UmiCheck::Ok(packed)); - packed - } - rejected => { - self.stats.record_umi(&rejected); - return out; - } + let UmiCheck::Ok(umi) = umi_status else { + return out; }; - // SJ feature: record (cell, UMI, junction) for each crossed junction. - // Only for resolved CBs (1MM_multi deferral is not applied to SJ). - if self.sj_enabled - && !junctions.is_empty() - && let Some(cb) = cb_resolved - { - out.sj = junctions - .iter() - .map(|&(intron_start, intron_end)| SjCountRecord { - cb, - umi, - intron_start, - intron_end, - }) - .collect(); + // SJ feature: record (cell, UMI, junction) for each crossed junction; + // a 1MM_multi barcode's cell is resolved when the matrix is built. + if self.sj_enabled && !junctions.is_empty() { + match (cb_resolved, &cb_match) { + (Some(cb), _) => { + out.sj.unique = junctions + .iter() + .map(|&(intron_start, intron_end)| SjCountRecord { + cb, + umi, + intron_start, + intron_end, + exact: code == 0, + }) + .collect(); + } + (None, CbMatch::Multi(cands)) => { + out.sj.multi = junctions + .iter() + .map(|&(intron_start, intron_end)| SjMultiRecord { + candidates: cands.clone(), + umi, + intron_start, + intron_end, + }) + .collect(); + } + (None, _) => {} + } } - // Velocyto feature: gene from gene-body overlap, then classify the read - // spliced/unspliced/ambiguous. Resolved CB only. + // Velocyto: the transcripts a uniquely mapped read is compatible with + // (`classifyAlign` computes them only for `nAlignG==1`), recorded for + // any read whose barcode and UMI pass (`record` returns before + // `cbMatch<0`), counted or not; its cell and UMI come from the Gene + // feature's readInfo when the molecules are collapsed. if self.velocyto_enabled - && let Some(cb) = cb_resolved - && let GeneAssignment::Gene(gene) = class.gene_full + && n_loci == 1 + && let [aln] = cdna_transcripts { - out.velocyto = Some(VelocytoRecord { - cb, - umi, - gene, - category: velocyto_category(cdna_transcripts, &self.gene_ann, gene), - }); + let types = velocyto_transcript_types(aln, &self.gene_model, self.strand); + if !types.is_empty() { + let cbs = match &cb_match { + CbMatch::Exact(i) | CbMatch::Corrected(i) => vec![*i], + CbMatch::Multi(cands) => cands.iter().map(|c| c.wl_index).collect(), + CbMatch::NoMatch | CbMatch::NinCb | CbMatch::MultMatchRejected => Vec::new(), + }; + out.velocyto = Some(VelocytoRecord { + read_index: NO_READ_INDEX, + types, + cbs, + }); + } } // The CB match + UMI are shared across features; reuse the cached @@ -1664,10 +2055,10 @@ impl SoloContext { SoloFeature::Gene => class.gene, SoloFeature::GeneFull => class.gene_full, }; - // STAR's no-feature record, feature 0 only: the readInfo of - // the `CB`/`UB` tags comes from the tag feature, the first one. + // STAR's no-feature record, for the features whose readInfo + // is kept (the first for `CB`/`UB`, Gene for Velocyto). let no_feature = |fo: &mut FeatureOutcome| { - if fi != 0 || self.read_info.is_none() { + if self.read_info.get(fi).is_none_or(Option::is_none) { return; } match (cb_resolved, &cb_match) { @@ -1719,9 +2110,6 @@ impl SoloContext { return fo; } }; - // Reads uniquely mapped to a gene under this feature, among - // valid-barcode reads (STARsolo "Reads Mapped to "). - self.feature_reads[fi].fetch_add(1, Ordering::Relaxed); match (cb_resolved, &cb_match) { (Some(cb), _) => { fo.record = Some(SoloCountRecord { @@ -1745,6 +2133,9 @@ impl SoloContext { }) .collect(); + // The read's `cbMatch` is >= 0 here: STAR's `SoloReadFeature::record` + // counts it for every feature. + self.record_feature_stats(n_loci, &class, &out, code == 0, junctions); self.record_cell_read( cb_resolved, &cb_match, @@ -1834,7 +2225,7 @@ impl SoloContext { pub fn process_read_pe( &self, pairs: &[(&Transcript, &Transcript)], - barcode: Option<&CellBarcode>, + barcode: Result<&CellBarcode, BarcodeMiss>, junctions: &[(u64, u64)], ) -> SoloReadOutcome { // A mate pair is one alignment in STAR (one `Transcript`, both mates' @@ -2308,10 +2699,10 @@ mod tests { } #[test] - fn extract_too_short_returns_none() { + fn extract_too_short_is_no_cb() { let layout = v2_layout(); let read = encoded_read("short", "AAAAAAAACCCC", "IIIIIIIIIIII"); - assert!(layout.extract(&read).is_none()); + assert_eq!(layout.extract(&read), Err(BarcodeMiss::NoCb)); } #[test] @@ -2522,7 +2913,7 @@ mod tests { quality: vec![b'I'; 12], file_index: 0, }; - assert!(layout.extract(&read).is_some()); + assert!(layout.extract(&read).is_ok()); // No adapter at all: no barcode. A barcode read at an unknown offset // is not salvageable, so guessing would be worse than declining. @@ -2532,7 +2923,7 @@ mod tests { quality: vec![b'I'; 12], file_index: 0, }; - assert!(layout.extract(&read).is_none()); + assert_eq!(layout.extract(&read), Err(BarcodeMiss::NoAdapter)); } /// The discriminating case: an early offset that fits the budget against a @@ -2609,7 +3000,7 @@ mod tests { quality: vec![b'I'; 12], file_index: 0, }; - assert!(layout.extract(&read).is_none()); + assert_eq!(layout.extract(&read), Err(BarcodeMiss::NoAdapter)); } #[test] @@ -2639,6 +3030,6 @@ mod tests { quality: vec![b'I'; 7], file_index: 0, }; - assert!(layout.extract(&read).is_none()); + assert_eq!(layout.extract(&read), Err(BarcodeMiss::NoCb)); } } diff --git a/src/solo/whitelist.rs b/src/solo/whitelist.rs index 40a6fa1..a5287ed 100644 --- a/src/solo/whitelist.rs +++ b/src/solo/whitelist.rs @@ -466,7 +466,14 @@ impl CbWhitelist { } /// Increment the exact-match count for sorted whitelist index `idx`. - fn bump_exact(&self, idx: u32) { + /// + /// STAR counts a read here only when its final `cbMatch` is 0, i.e. the + /// barcode matched exactly *and* the UMI passed + /// (`SoloReadBarcode::addStats`, `SoloReadBarcode_getCBandUMI.cpp:93-100`; + /// only a `MATCH_CellRanger` build also counts reads whose UMI fails), so + /// the caller bumps it once the read's code is known. These counts feed + /// `oneExact` and the 1MM_multi posterior. + pub fn bump_exact(&self, idx: u32) { if let Self::List { exact_counts, .. } = self { exact_counts[idx as usize].fetch_add(1, Ordering::Relaxed); } @@ -502,7 +509,6 @@ impl CbWhitelist { PackResult::ManyN => CbMatch::NinCb, PackResult::NoN(packed) => { if let Some(idx) = self.search(packed) { - self.bump_exact(idx); return CbMatch::Exact(idx); } if !match_type.mm1 { @@ -590,49 +596,100 @@ fn open_maybe_gzip(path: &Path) -> Result, Error> { // Stats (STAR cbMatch categories) // --------------------------------------------------------------------------- -/// Per-run barcode-matching statistics, mirroring STAR's `SoloReadBarcodeStats`. +/// STAR's `SoloReadBarcodeStats` (`SoloReadBarcodeStats.h`), written to +/// `Solo.out/Barcodes.stats`, plus the CB+UMI base-quality tally behind +/// `Summary.csv`'s "Q30 Bases in CB+UMI". #[derive(Debug, Default)] -pub struct CbMatchStats { - pub yes_exact: AtomicU64, - pub yes_one_mm: AtomicU64, - pub yes_mult_mm: AtomicU64, - pub no_match: AtomicU64, - pub n_in_cb: AtomicU64, - pub mult_rejected: AtomicU64, - pub n_in_umi: AtomicU64, - pub umi_homopolymer: AtomicU64, +pub struct BarcodeStats { + counts: [AtomicU64; BarcodeStats::NAMES.len()], + q30_bases: AtomicU64, + bases: AtomicU64, } -impl CbMatchStats { +impl BarcodeStats { + /// STAR's names, in its order. + pub const NAMES: [&'static str; 12] = [ + "noNoAdapter", + "noNoUMI", + "noNoCB", + "noNinCB", + "noNinUMI", + "noUMIhomopolymer", + "noNoWLmatch", + "noTooManyMM", + "noTooManyWLmatches", + "yesWLmatchExact", + "yesOneWLmatchWithMM", + "yesMultWLmatchWithMM", + ]; + const NO_NO_ADAPTER: usize = 0; + const NO_NO_UMI: usize = 1; + pub fn new() -> Self { Self::default() } - /// Record one CB match outcome. - pub fn record_cb(&self, m: &CbMatch) { - let c = match m { - CbMatch::Exact(_) => &self.yes_exact, - CbMatch::Corrected(_) => &self.yes_one_mm, - CbMatch::Multi(_) => &self.yes_mult_mm, - CbMatch::NoMatch => &self.no_match, - CbMatch::NinCb => &self.n_in_cb, - CbMatch::MultMatchRejected => &self.mult_rejected, + /// Count one read by its final `cbMatch` code, STAR's + /// `SoloReadBarcode::addStats` switch (`SoloReadBarcode_getCBandUMI.cpp:93-130`): + /// a code it does not list falls to `default`, which counts as a + /// multiple-match read. + pub fn record(&self, code: i32) { + let i = match code { + 0 => 9, + 1 => 10, + -1 => 6, + -2 => 3, + -3 => 8, + -11 => 2, + -12 => 7, + -23 => 4, + -24 => 5, + _ => 11, }; - c.fetch_add(1, Ordering::Relaxed); + self.counts[i].fetch_add(1, Ordering::Relaxed); } - /// Record one UMI check outcome (only the rejection cases are counted). - pub fn record_umi(&self, u: &UmiCheck) { - match u { - UmiCheck::NinUmi => { - self.n_in_umi.fetch_add(1, Ordering::Relaxed); - } - UmiCheck::Homopolymer => { - self.umi_homopolymer.fetch_add(1, Ordering::Relaxed); - } - UmiCheck::Ok(_) => {} - } + /// `CB_UMI_Complex` with no adapter within the mismatch budget: STAR counts + /// it directly, outside `addStats` (`:340-343`). + pub fn no_adapter(&self) { + self.counts[Self::NO_NO_ADAPTER].fetch_add(1, Ordering::Relaxed); + } + + /// `CB_UMI_Complex` with no extractable UMI (`:347-350`). + pub fn no_umi(&self) { + self.counts[Self::NO_NO_UMI].fetch_add(1, Ordering::Relaxed); + } + + /// Add one read's CB+UMI `(Q30 bases, bases)` (`qualHist`, counted + /// against `readQualityScoreBase + 30`). + pub fn add_q30(&self, q30: u64, bases: u64) { + self.q30_bases.fetch_add(q30, Ordering::Relaxed); + self.bases.fetch_add(bases, Ordering::Relaxed); } + + /// The counts, in [`Self::NAMES`] order. + pub fn values(&self) -> [u64; 12] { + std::array::from_fn(|i| self.counts[i].load(Ordering::Relaxed)) + } + + /// STAR's `numInvalidBarcodes`: the first nine counts. + pub fn invalid(&self) -> u64 { + self.values()[..9].iter().sum() + } + + /// `(Q30 bases, all bases)` of the CB+UMI qualities. + pub fn q30(&self) -> (u64, u64) { + ( + self.q30_bases.load(Ordering::Relaxed), + self.bases.load(Ordering::Relaxed), + ) + } +} + +/// One line of a STAR stats file: the name right-aligned in 50 columns, the +/// value in 15 (`statsOut`: `setw(50) << name << setw(15) << value`). +pub fn stats_line(name: &str, value: u64) -> String { + format!("{name:>50}{value:>15}\n") } #[cfg(test)] @@ -682,8 +739,13 @@ mod tests { CbMatch::Exact(idx) => assert_eq!(wl.barcode_string(idx).unwrap(), "ACGT"), other => panic!("expected exact, got {other:?}"), } - let counts = wl.exact_count_snapshot(); - assert_eq!(counts.iter().sum::(), 1); + // Matching does not count: STAR counts an exact read only once its + // UMI has passed too, so the caller bumps it then. + assert_eq!(wl.exact_count_snapshot().iter().sum::(), 0); + if let CbMatch::Exact(idx) = m { + wl.bump_exact(idx); + } + assert_eq!(wl.exact_count_snapshot().iter().sum::(), 1); } #[test] diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index 8294975..7e71c04 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -1216,7 +1216,7 @@ fn test_starsolo_gene_matrix() { "expected 8 resolved solo records in log, stderr was:\n{stderr}" ); assert!( - stderr.contains("exact=8"), + stderr.contains("yesWLmatchExact=8"), "expected 8 exact CB matches in log, stderr was:\n{stderr}" ); @@ -1275,7 +1275,8 @@ fn test_starsolo_gene_matrix() { // --------------------------------------------------------------------------- // Test 9a' — Summary.csv stays STARsolo-faithful; the CellRanger mapping funnel // (exonic/intronic/intergenic/antisense) is split out into a separate -// CellRanger.summary.csv (PR #90 review: keep the faithful Summary.csv unaltered). +// CellRanger.summary.csv, written only with `--soloCellRangerSummary Yes` (a +// rustar-aligner extension; the faithful Summary.csv is never altered). // Needs Gene + GeneFull — the exonic/intronic split requires both queries. // --------------------------------------------------------------------------- #[test] @@ -1304,42 +1305,56 @@ fn test_starsolo_summary_split() { fs::write(&wl_path, format!("{cb}\nCCCCGGGGTTTTAAAA\n")).unwrap(); } - let output_dir = tmpdir.path().join("out_split"); - fs::create_dir_all(&output_dir).unwrap(); - let prefix = format!("{}/", output_dir.display()); - cargo_bin_cmd!("rustar-aligner") - .args([ - "--runMode", - "alignReads", - "--genomeDir", - genome_dir.to_str().unwrap(), - "--readFilesIn", - cdna_path.to_str().unwrap(), - bc_path.to_str().unwrap(), - "--soloType", - "CB_UMI_Simple", - "--soloCBwhitelist", - wl_path.to_str().unwrap(), - "--soloCBstart", - "1", - "--soloCBlen", - "16", - "--soloUMIstart", - "17", - "--soloUMIlen", - "12", - "--soloFeatures", - "Gene", - "GeneFull", - "--soloStrand", - "Forward", - "--sjdbGTFfile", - gtf.to_str().unwrap(), - "--outFileNamePrefix", - &prefix, - ]) - .assert() - .success(); + // Run once with the defaults and once with the extension on. + let run = |name: &str, extra: &[&str]| { + let output_dir = tmpdir.path().join(name); + fs::create_dir_all(&output_dir).unwrap(); + let prefix = format!("{}/", output_dir.display()); + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + cdna_path.to_str().unwrap(), + bc_path.to_str().unwrap(), + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + wl_path.to_str().unwrap(), + "--soloCBstart", + "1", + "--soloCBlen", + "16", + "--soloUMIstart", + "17", + "--soloUMIlen", + "12", + "--soloFeatures", + "Gene", + "GeneFull", + "--soloStrand", + "Forward", + "--sjdbGTFfile", + gtf.to_str().unwrap(), + "--outFileNamePrefix", + &prefix, + ]) + .args(extra) + .assert() + .success(); + output_dir + }; + // A rustar-aligner extension: nothing beyond STARsolo's files by default. + let default_dir = run("out_split_default", &[]); + assert!( + !default_dir + .join("Solo.out/Gene/CellRanger.summary.csv") + .exists(), + "CellRanger.summary.csv must be opt-in" + ); + let output_dir = run("out_split", &["--soloCellRangerSummary", "Yes"]); let gene = output_dir.join("Solo.out").join("Gene"); let summary = fs::read_to_string(gene.join("Summary.csv")).unwrap(); @@ -1447,6 +1462,135 @@ fn test_starsolo_sj_feature() { assert_eq!(matrix.lines().last().unwrap(), "1 1 1"); } +// --------------------------------------------------------------------------- +// Test 9b2 — STARsolo stats files: Barcodes.stats by final cbMatch code, and +// SJ Features.stats with a 1MM_multi read resolved by the posterior. +// +// The posterior needs the competing barcode's exact count to stay 0: a read +// that matches it exactly but has an N in its UMI is counted as noNinUMI and +// does not count as an exact match (`addStats`). If it did, the multi read +// would fail cbMinP (5/6 < 0.975) and be dropped. +// --------------------------------------------------------------------------- +#[test] +fn test_starsolo_stats_files() { + let tmpdir = TempDir::new().unwrap(); + let genome = build_genome(); + let fasta = write_fasta(&tmpdir, &genome); + let gtf = write_gtf(&tmpdir); + let genome_dir = tmpdir.path().join("genome"); + build_index(&fasta, &genome_dir, "7", Some(>f)); + + let cdna_path = tmpdir.path().join("cdna.fq"); + let barcode_path = tmpdir.path().join("barcode.fq"); + let wl_path = tmpdir.path().join("whitelist.txt"); + let cb = "AAAACCCCGGGGTTTT"; + let cb2 = "AAAACCCCGGGGTTTA"; + // Spliced across the intron [10050,10250); unspliced inside exon 1. + let mut spliced = genome[10025..10050].to_vec(); + spliced.extend_from_slice(&genome[10250..10275]); + let exonic = genome[10000..10050].to_vec(); + let reads: Vec<(&[u8], &str, &str)> = vec![ + (&spliced, cb, "ACGTACGTAC"), // 4 exact reads, one molecule + (&spliced, cb, "ACGTACGTAC"), + (&spliced, cb, "ACGTACGTAC"), + (&spliced, cb, "ACGTACGTAC"), + (&spliced, "AAAACCCCGGGGTTTC", "TGCATGCATG"), // 1MM to cb and cb2 + (&spliced, cb2, "ACGTNCGTAC"), // exact, N in UMI + (&spliced, "ACACACACACACACAC", "GATCGATCGA"), // no whitelist match + (&exonic, cb, "CATGCATGCA"), // exact, no junction + ]; + { + let mut cf = fs::File::create(&cdna_path).unwrap(); + let mut bf = fs::File::create(&barcode_path).unwrap(); + for (i, (seq, bc, umi)) in reads.iter().enumerate() { + writeln!(cf, "@r{i}").unwrap(); + cf.write_all(seq).unwrap(); + writeln!(cf, "\n+\n{}", "I".repeat(seq.len())).unwrap(); + writeln!(bf, "@r{i}\n{bc}{umi}\n+\n{}", "I".repeat(26)).unwrap(); + } + fs::write(&wl_path, format!("{cb}\n{cb2}\nCCCCGGGGTTTTAAAA\n")).unwrap(); + } + + let output_dir = tmpdir.path().join("out_stats"); + fs::create_dir_all(&output_dir).unwrap(); + let prefix = format!("{}/", output_dir.display()); + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + cdna_path.to_str().unwrap(), + barcode_path.to_str().unwrap(), + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + wl_path.to_str().unwrap(), + "--soloFeatures", + "Gene", + "SJ", + "--soloStrand", + "Forward", + "--sjdbGTFfile", + gtf.to_str().unwrap(), + "--outFileNamePrefix", + &prefix, + ]) + .assert() + .success(); + + let solo = output_dir.join("Solo.out"); + let stats = |path: PathBuf| -> Vec<(String, u64)> { + fs::read_to_string(&path) + .unwrap() + .lines() + .map(|l| { + let mut f = l.split_whitespace(); + let name = f.next().unwrap().to_string(); + (name, f.next().unwrap().parse().unwrap()) + }) + .collect() + }; + let nonzero = |v: Vec<(String, u64)>| -> Vec<(String, u64)> { + v.into_iter().filter(|(_, n)| *n > 0).collect() + }; + let named = |v: &[(&str, u64)]| -> Vec<(String, u64)> { + v.iter().map(|&(n, c)| (n.to_string(), c)).collect() + }; + + let barcodes = stats(solo.join("Barcodes.stats")); + assert_eq!(barcodes.len(), 12, "{barcodes:?}"); + assert_eq!( + nonzero(barcodes), + named(&[ + ("noNinUMI", 1), + ("noNoWLmatch", 1), + ("yesWLmatchExact", 5), + ("yesMultWLmatchWithMM", 1), + ]) + ); + + // SJ: the five spliced reads with a whitelisted barcode, the multi one + // resolved to `cb`; the exonic read has no junction. + let sj = stats(solo.join("SJ").join("Features.stats")); + assert_eq!(sj.len(), 11, "{sj:?}"); + assert_eq!( + nonzero(sj), + named(&[ + ("noNoFeature", 1), + ("yesWLmatch", 5), + ("yessubWLmatchExact", 4), + ("yessubWLmatch_UniqueFeature", 5), + ("yesCellBarcodes", 1), + ("yesUMIs", 2), + ]) + ); + let sj_matrix = fs::read_to_string(solo.join("SJ").join("raw").join("matrix.mtx")).unwrap(); + // Whitelist order: cb2 (…TTTA) is column 1, cb column 2. + assert_eq!(sj_matrix.lines().last().unwrap(), "1 2 2", "{sj_matrix}"); +} + // --------------------------------------------------------------------------- // Test 9c — STARsolo --soloMultiMappers (gene-ambiguous distribution) // @@ -1546,8 +1690,9 @@ fn test_starsolo_multimappers() { let um = fs::read_to_string(raw.join("UniqueAndMult-Uniform.mtx")).unwrap(); assert!(um.contains("coordinate real general"), "um header:\n{um}"); let rows: Vec<&str> = um.lines().filter(|l| !l.starts_with('%')).skip(1).collect(); - assert!(rows.contains(&"1 1 0.50000"), "expected G1 0.5, got:\n{um}"); - assert!(rows.contains(&"2 1 0.50000"), "expected G3 0.5, got:\n{um}"); + // Written as STAR's `ostream << double` does. + assert!(rows.contains(&"1 1 0.5"), "expected G1 0.5, got:\n{um}"); + assert!(rows.contains(&"2 1 0.5"), "expected G3 0.5, got:\n{um}"); assert!(rows.contains(&"3 1 1"), "expected G2 1, got:\n{um}"); } @@ -1700,9 +1845,9 @@ fn test_starsolo_smartseq_paired() { // --------------------------------------------------------------------------- // Test 9f — STARsolo Velocyto (spliced / unspliced / ambiguous) // -// Three reads on gene G1, one per category: a junction-spanning read (spliced), -// a purely intronic read (unspliced), and a wholly-exonic read with no junction -// (ambiguous, per Sullivan 2025). Distinct UMIs → one molecule in each matrix. +// Three reads on gene G1 (one transcript), each its own molecule: across the +// junction and inside exon 1 are both spliced, inside the intron unspliced, as +// STAR's model calls them (`SoloFeature_countVelocyto.cpp`). // --------------------------------------------------------------------------- #[test] fn test_starsolo_velocyto() { @@ -1720,10 +1865,13 @@ fn test_starsolo_velocyto() { // category → cDNA read + a distinct (non-homopolymer) 12 bp UMI. let mut spliced = genome[10025..10050].to_vec(); // Exon1 end ... spliced.extend_from_slice(&genome[10250..10275]); // ... + Exon2 start → junction + // STAR's model (`SoloFeature_countVelocyto.cpp`) against T1 alone: a read + // inside exons is spliced whether or not it crosses the junction, and a + // read inside the intron is unspliced. let reads: [(Vec, &str); 3] = [ - (spliced, "ACGTACGTACGT"), // spliced + (spliced, "ACGTACGTACGT"), // across the junction → spliced (genome[10100..10150].to_vec(), "TGCATGCATGCA"), // intronic → unspliced - (genome[10000..10050].to_vec(), "GATCGATCGATC"), // exonic, no junction → ambiguous + (genome[10000..10050].to_vec(), "GATCGATCGATC"), // inside exon 1 → spliced ]; { let mut cf = fs::File::create(&cdna_path).unwrap(); @@ -1737,11 +1885,20 @@ fn test_starsolo_velocyto() { fs::write(&wl_path, format!("{cb}\nCCCCGGGGTTTTAAAA\n")).unwrap(); } - let output_dir = tmpdir.path().join("out_velo"); - fs::create_dir_all(&output_dir).unwrap(); - let prefix = format!("{}/", output_dir.display()); - cargo_bin_cmd!("rustar-aligner") - .args([ + // Gene alone, then GeneFull first with CB/UB: the tags then read + // GeneFull's readInfo and Velocyto the Gene feature's, as in STAR. + for (tag, features, tags) in [ + ("gene", &["Gene", "Velocyto"][..], false), + ( + "genefull_first", + &["GeneFull", "Gene", "Velocyto"][..], + true, + ), + ] { + let output_dir = tmpdir.path().join(format!("out_velo_{tag}")); + fs::create_dir_all(&output_dir).unwrap(); + let prefix = format!("{}/", output_dir.display()); + let mut args = vec![ "--runMode", "alignReads", "--genomeDir", @@ -1761,38 +1918,66 @@ fn test_starsolo_velocyto() { "17", "--soloUMIlen", "12", - "--soloFeatures", - "Velocyto", "--soloStrand", "Forward", "--sjdbGTFfile", gtf.to_str().unwrap(), "--outFileNamePrefix", &prefix, - ]) - .assert() - .success(); + "--soloFeatures", + ]; + args.extend_from_slice(features); + if tags { + args.extend_from_slice(&[ + "--outSAMtype", + "BAM", + "SortedByCoordinate", + "--outSAMattributes", + "NH", + "CB", + "UB", + ]); + } + cargo_bin_cmd!("rustar-aligner") + .args(&args) + .assert() + .success(); - let raw = output_dir.join("Solo.out").join("Velocyto").join("raw"); - // Each category matrix holds exactly its one molecule for G1 (row 1, col 1). - for name in ["spliced", "unspliced", "ambiguous"] { - let m = fs::read_to_string(raw.join(format!("{name}.mtx"))).unwrap(); - assert_eq!( - m.lines().last().unwrap(), - "1 1 1", - "{name}.mtx should have G1=1:\n{m}" - ); + // One (G1, cell) row in every matrix, zeros included, as STAR writes + // them; the cell is called, so the filtered matrices hold it too. + let velo = output_dir.join("Solo.out").join("Velocyto"); + for (dir, n_cells) in [("raw", 2), ("filtered", 1)] { + for (name, want) in [ + ("spliced", "1 1 2"), + ("unspliced", "1 1 1"), + ("ambiguous", "1 1 0"), + ] { + let m = fs::read_to_string(velo.join(dir).join(format!("{name}.mtx"))).unwrap(); + assert_eq!( + m.lines().nth(2).unwrap(), + format!("1 {n_cells} 1"), + "{tag} {dir}/{name}.mtx header:\n{m}" + ); + assert_eq!( + m.lines().last().unwrap(), + want, + "{tag} {dir}/{name}.mtx:\n{m}" + ); + } + } + let filtered_barcodes = + fs::read_to_string(velo.join("filtered").join("barcodes.tsv")).unwrap(); + assert_eq!(filtered_barcodes, format!("{cb}\n"), "{tag}"); } } // --------------------------------------------------------------------------- -// Test 9f — Velocyto fold mode (`--soloVelocytoAmbiguous no`): exon-only -// molecules are folded into spliced and no `ambiguous.mtx` is written -// (cf. He, Soneson & Patro 2023). Same fixture as Test 9 but the exon-only -// molecule (distinct UMI) lands in spliced, giving G1 = 2 there. +// Test 9f2 — Velocyto ambiguous molecules: with a retained-intron isoform, a +// read in the intron is exonic in one transcript and intronic in the other, +// which STAR calls ambiguous (`SoloFeature_countVelocyto.cpp:125-133`). // --------------------------------------------------------------------------- #[test] -fn test_starsolo_velocyto_fold_ambiguous() { +fn test_starsolo_velocyto_retained_intron_is_ambiguous() { let tmpdir = TempDir::new().unwrap(); let genome = build_genome(); let fasta = write_fasta(&tmpdir, &genome); @@ -1806,10 +1991,12 @@ fn test_starsolo_velocyto_fold_ambiguous() { let cb = "AAAACCCCGGGGTTTT"; let mut spliced = genome[10025..10050].to_vec(); spliced.extend_from_slice(&genome[10250..10275]); + // With T2 retaining the intron, the intronic read is exonic in T2 and + // intronic in T1: ambiguous. The other two are spliced in both. let reads: [(Vec, &str); 3] = [ (spliced, "ACGTACGTACGT"), // spliced - (genome[10100..10150].to_vec(), "TGCATGCATGCA"), // intronic → unspliced - (genome[10000..10050].to_vec(), "GATCGATCGATC"), // exonic, no junction → ambiguous + (genome[10100..10150].to_vec(), "TGCATGCATGCA"), // T1 intron, T2 exon → ambiguous + (genome[10000..10050].to_vec(), "GATCGATCGATC"), // spliced ]; { let mut cf = fs::File::create(&cdna_path).unwrap(); @@ -1823,7 +2010,16 @@ fn test_starsolo_velocyto_fold_ambiguous() { fs::write(&wl_path, format!("{cb}\nCCCCGGGGTTTTAAAA\n")).unwrap(); } - let output_dir = tmpdir.path().join("out_velo_fold"); + // T2: the same gene with the intron retained. + { + let mut f = fs::OpenOptions::new().append(true).open(>f).unwrap(); + writeln!( + f, + "chr1\ttest\texon\t10001\t10300\t.\t+\t.\tgene_id \"G1\"; transcript_id \"T2\";" + ) + .unwrap(); + } + let output_dir = tmpdir.path().join("out_velo_retained"); fs::create_dir_all(&output_dir).unwrap(); let prefix = format!("{}/", output_dir.display()); cargo_bin_cmd!("rustar-aligner") @@ -1848,9 +2044,8 @@ fn test_starsolo_velocyto_fold_ambiguous() { "--soloUMIlen", "12", "--soloFeatures", + "Gene", "Velocyto", - "--soloVelocytoAmbiguous", - "no", "--soloStrand", "Forward", "--sjdbGTFfile", @@ -1862,22 +2057,14 @@ fn test_starsolo_velocyto_fold_ambiguous() { .success(); let raw = output_dir.join("Solo.out").join("Velocyto").join("raw"); - assert!( - !raw.join("ambiguous.mtx").exists(), - "ambiguous.mtx must not be written in fold mode" - ); - let unspliced = fs::read_to_string(raw.join("unspliced.mtx")).unwrap(); - assert_eq!( - unspliced.lines().last().unwrap(), - "1 1 1", - "unspliced unchanged by folding:\n{unspliced}" - ); - let spliced_m = fs::read_to_string(raw.join("spliced.mtx")).unwrap(); - assert_eq!( - spliced_m.lines().last().unwrap(), - "1 1 2", - "spliced should gain the folded exon-only molecule:\n{spliced_m}" - ); + for (name, want) in [ + ("spliced", "1 1 2"), + ("unspliced", "1 1 0"), + ("ambiguous", "1 1 1"), + ] { + let m = fs::read_to_string(raw.join(format!("{name}.mtx"))).unwrap(); + assert_eq!(m.lines().last().unwrap(), want, "{name}.mtx:\n{m}"); + } } // --------------------------------------------------------------------------- @@ -2571,9 +2758,12 @@ fn test_run_mode_solo_cell_filtering_calls_cells_from_a_raw_matrix() { "soloCellFiltering", raw.to_str().unwrap(), out.to_str().unwrap(), + // STAR's TopCells N keeps every barcode with at least the UMIs of + // the barcode at index N (`SoloFeature_cellFiltering.cpp:46-54`): + // 4 sets the threshold at the fifth deep barcode. "--soloCellFilter", "TopCells", - "5", + "4", ]) .assert() .success(); @@ -3361,9 +3551,11 @@ fn test_solo_cb_sam_tag_out() { assert_eq!(by_cr["TTTTGGGGCCCCAAAA"]["sM"], "-1"); // No UB (there is no UMI collapsing at all) and no count matrices. assert!(tags.iter().all(|t| !t.contains_key("UB"))); - assert!( - !output_dir.join("Solo.out").exists(), - "CB_samTagOut must not write Solo.out" + // No count matrices: STAR writes only Barcodes.stats (`Solo::processAndOutput`). + assert_eq!( + solo_out_files(&output_dir), + ["Barcodes.stats"], + "CB_samTagOut writes only Barcodes.stats" ); } @@ -3457,5 +3649,27 @@ fn test_solo_paired_cdna_with_separate_barcode_read() { // sS is the whole barcode read, not either cDNA mate. assert_eq!(t.get("sS").map(|s| s.len()), Some(26)); } - assert!(!output_dir.join("Solo.out").exists()); + assert_eq!(solo_out_files(&output_dir), ["Barcodes.stats"]); +} + +/// The files under `Solo.out/`, relative and sorted. +fn solo_out_files(output_dir: &Path) -> Vec { + let root = output_dir.join("Solo.out"); + let mut out = Vec::new(); + let mut stack = vec![root.clone()]; + while let Some(dir) = stack.pop() { + let Ok(entries) = fs::read_dir(&dir) else { + continue; + }; + for e in entries.flatten() { + let p = e.path(); + if p.is_dir() { + stack.push(p); + } else { + out.push(p.strip_prefix(&root).unwrap().display().to_string()); + } + } + } + out.sort(); + out }