From a10b29b17310dda628c86df81ba1d5c687c9308b Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 30 Jul 2026 15:27:00 +0200 Subject: [PATCH 01/10] fix: --runThreadN 1 ran on every core, not on one The rayon global pool was configured only when `--runThreadN` was greater than 1. Skipping the build at 1 does not give one thread: it leaves rayon's default, which is one worker per logical core. So `--runThreadN 1` ran the whole machine. Measured on 200k reads, before: 2.44 s wall at **1300% CPU**. After: 26.23 s wall at 100% CPU. The old figure was not a fast single-threaded run, it was a sixteen-way run wearing the wrong flag. This matters beyond the flag reading falsely. A scheduler or a container given one CPU gets sixteen worker threads; on a shared machine the run oversubscribes every other job; and with a thread-caching allocator each of those threads keeps its own heap, which is the very cost the comment above this code says the pool sizing exists to avoid. It also means the project's thread-invariance checks were weaker than they read: the `--runThreadN 1` leg was not a one-thread leg. Verified now that it is one: records are byte-identical between 1 and 8 threads on 200k real reads, and byte-identical to the previous binary's output at `--runThreadN 1`. Only the `@PG` `CL:` line differs between thread counts, because it records the command line. Co-Authored-By: Claude Opus 5 (1M context) --- src/lib.rs | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 6fce173..8a07f5f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -68,11 +68,13 @@ pub fn run(params: &Parameters) -> anyhow::Result<()> { // RSS for nothing. `build_global` errors if called twice; we // ignore the error so the in-process tests that already // initialised the pool still work. - if usize::from(params.run_thread_n) > 1 { - let _ = rayon::ThreadPoolBuilder::new() - .num_threads(params.run_thread_n.into()) - .build_global(); - } + // + // Configured at every value, `1` included. Skipping the build at 1 does not + // yield one thread: it leaves rayon's default of one worker per logical + // core, so `--runThreadN 1` ran on the whole machine. + let _ = rayon::ThreadPoolBuilder::new() + .num_threads(params.run_thread_n.into()) + .build_global(); match params.run_mode { RunMode::GenomeGenerate => genome_generate(params), From 86ffe9845cf8a3e8890d5802cd15b0291460279a Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 30 Jul 2026 20:12:03 +0200 Subject: [PATCH 02/10] docs(changelog): record the --runThreadN 1 fix Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index b6e4893..f103273 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -99,6 +99,11 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Bug fixes +- `--runThreadN 1` ran on every logical core instead of on one. The + rayon pool was configured only above 1, and skipping it leaves rayon's + default of one worker per core. Output is unchanged; the run now uses + the thread count asked for. + - **STARsolo `Gene` assignment now requires exon concordance**, matching STARsolo: a read counts toward a gene only when every aligned block lies within the gene's exons, rather than merely overlapping one. This From 164b29781559555bb153c1cfdba411295e2cef5c Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 30 Jul 2026 23:52:21 +0200 Subject: [PATCH 03/10] refactor: drop Transcript::read_seq, which nothing ever reads The field was filled with a full copy of the read at every finalised transcript and assigned again in four more places, and no code in the crate reads it. Not "reads it rarely": the compiler was asked, and after deleting the field every one of the 97 resulting errors is a struct literal or an assignment. There is no read site, in `src` or in tests. Measured on 200k real reads at 8 threads: about 50 transcripts are finalised per read, so removing it takes **15.1 million allocations off a 200 million total** (7.6%) and 3.87 GB of copying. Wall clock does not move, and that is worth recording rather than hiding: six interleaved rounds at 87-92% CPU idle give medians 20.90 s against 20.65 s with the direction mixed, inside the run-to-run spread. mimalloc is fast enough that seventy-five small allocations per read do not surface. The reason to remove it is that it is dead weight, not that it is slow. That number also calibrates #168 downward: if removing 7.6% of the allocations changes nothing measurable, the rest of the allocation programme is unlikely to be worth a new dependency. `Transcript::read_seq` is `pub`, so this is an API removal and needs sign-off. Nothing outside the crate can be relying on its contents being meaningful, though, since it is only ever written. Output-neutral: SAM byte-identical on 200k real reads. Co-Authored-By: Claude Opus 5 (1M context) --- src/align/pe_overlap.rs | 3 -- src/align/read_align.rs | 18 ------------ src/align/stitch.rs | 2 -- src/align/transcript.rs | 7 ----- src/chimeric/detect.rs | 3 -- src/io/bam.rs | 1 - src/io/sam.rs | 56 -------------------------------------- src/quant/mod.rs | 2 -- src/quant/transcriptome.rs | 2 -- src/signal/mod.rs | 1 - src/solo/gene.rs | 1 - src/stats.rs | 2 -- src/wasp/mod.rs | 1 - 13 files changed, 99 deletions(-) diff --git a/src/align/pe_overlap.rs b/src/align/pe_overlap.rs index cacbf49..98aa6cb 100644 --- a/src/align/pe_overlap.rs +++ b/src/align/pe_overlap.rs @@ -364,7 +364,6 @@ pub fn convert_merged_transcript_to_pe( n_junction, junction_motifs: out_junctions[i].iter().map(|(m, _)| *m).collect(), junction_annotated: out_junctions[i].iter().map(|(_, a)| *a).collect(), - read_seq: mate_seqs[i].to_vec(), }); } @@ -479,7 +478,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: merge.merged.clone(), }; let scorer = AlignmentScorer::from_params_minimal(); @@ -552,7 +550,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: merge.merged.clone(), }; let scorer = AlignmentScorer::from_params_minimal(); diff --git a/src/align/read_align.rs b/src/align/read_align.rs index b4d2f20..ec66a16 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -114,7 +114,6 @@ impl PairedAlignment { n_junction: m1.n_junction + m2.n_junction, junction_motifs: Vec::new(), junction_annotated: Vec::new(), - read_seq: Vec::new(), } } } @@ -889,8 +888,6 @@ pub fn align_paired_read( t1.is_reverse = false; t2.is_reverse = true; } - t1.read_seq = mate1_seq.to_vec(); - t2.read_seq = mate2_seq.to_vec(); if params.chim_segment_min > 0 { all_m1_transcripts.push(t1.clone()); @@ -934,7 +931,6 @@ pub fn align_paired_read( 0, // mate1 ) { t.is_reverse = stitch_is_reverse; - t.read_seq = mate1_seq.to_vec(); if params.chim_segment_min > 0 { all_m1_transcripts.push(t.clone()); } @@ -958,7 +954,6 @@ pub fn align_paired_read( 1, // mate2 ) { t.is_reverse = !stitch_is_reverse; - t.read_seq = mate2_seq.to_vec(); if params.chim_segment_min > 0 { all_m2_transcripts.push(t.clone()); } @@ -1641,7 +1636,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], }; let pair = PairedAlignment { mate1_transcript: make_tr(1000, 1100, 0, 100), @@ -1803,7 +1797,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 100], }; let t2 = Transcript { @@ -1825,7 +1818,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 100], }; // Distance = 300bp, within default limit (auto mode = unlimited) @@ -1859,7 +1851,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 100], }; let t2 = Transcript { @@ -1881,7 +1872,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 100], }; // Distance = 400bp, exceeds limit of 100bp @@ -1913,7 +1903,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 100], }; let t2 = Transcript { @@ -1935,7 +1924,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 100], }; let tlen = calculate_insert_size(&t1, &t2); @@ -1968,7 +1956,6 @@ mod tests { n_junction: 2, junction_motifs: vec![SpliceMotif::GtAg, SpliceMotif::CtAc], // +strand and -strand junction_annotated: vec![], - read_seq: vec![0; 100], }; // Create a transcript with consistent strand motifs (all + strand) @@ -1991,7 +1978,6 @@ mod tests { n_junction: 2, junction_motifs: vec![SpliceMotif::GtAg, SpliceMotif::GcAg], // both + strand junction_annotated: vec![], - read_seq: vec![0; 100], }; // Note: STAR's RemoveInconsistentStrands filters transcripts where @@ -2070,7 +2056,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 100], }; let t2 = Transcript { @@ -2092,7 +2077,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 100], }; let tlen = calculate_insert_size(&t1, &t2); @@ -2133,7 +2117,6 @@ mod tests { n_junction: 1, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 100], }; // Case 1: NonCanonical + unannotated → should be filtered @@ -2219,7 +2202,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 100], }; // Test BothMapped variant diff --git a/src/align/stitch.rs b/src/align/stitch.rs index 1c25ecc..349ff40 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -2155,7 +2155,6 @@ pub(crate) fn finalize_transcript( n_junction: wt.n_junction, junction_motifs: wt.junction_motifs.clone(), junction_annotated: wt.junction_annotated.clone(), - read_seq: read_seq.to_vec(), }) } @@ -2695,7 +2694,6 @@ pub(crate) fn stitch_seeds_with_jdb_debug( // Restore original reverse-strand flag and read sequence for SAM output. if stitch_is_reverse { transcript.is_reverse = true; - transcript.read_seq = read_seq.to_vec(); } transcripts.push(transcript); } diff --git a/src/align/transcript.rs b/src/align/transcript.rs index b611f52..ae358a4 100644 --- a/src/align/transcript.rs +++ b/src/align/transcript.rs @@ -29,8 +29,6 @@ pub struct Transcript { pub junction_motifs: Vec, /// Whether each junction is annotated in the GTF (for jM +20 offset) pub junction_annotated: Vec, - /// Original read sequence - pub read_seq: Vec, } /// An exon segment in a transcript. @@ -182,7 +180,6 @@ mod tests { n_junction: 1, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], }; assert_eq!(transcript.cigar_string(), "50M100N50M"); @@ -230,7 +227,6 @@ mod tests { n_junction: 1, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 100], }; // Query: 45 + 3 + 2 + 50 = 100 @@ -279,7 +275,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], }; let [left, right] = transcript.count_soft_clips(); @@ -303,7 +298,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], }; let [left, right] = transcript.count_soft_clips(); @@ -327,7 +321,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], }; let [left, right] = transcript.count_soft_clips(); diff --git a/src/chimeric/detect.rs b/src/chimeric/detect.rs index 4997049..a32af77 100644 --- a/src/chimeric/detect.rs +++ b/src/chimeric/detect.rs @@ -1090,7 +1090,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0u8; read_len], } } @@ -1210,7 +1209,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], }; assert!(transcript_to_segment(&t).is_err()); } @@ -1354,7 +1352,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0u8; read_len], } } diff --git a/src/io/bam.rs b/src/io/bam.rs index 1d2a143..b1e6717 100644 --- a/src/io/bam.rs +++ b/src/io/bam.rs @@ -550,7 +550,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], }; let read_name = "read1"; diff --git a/src/io/sam.rs b/src/io/sam.rs index 588c71d..343e9e0 100644 --- a/src/io/sam.rs +++ b/src/io/sam.rs @@ -1828,7 +1828,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }; writer @@ -1971,7 +1970,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }; let read_seq = vec![0, 1, 2, 3]; // ACGT @@ -2027,7 +2025,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }, Transcript { chr_idx: 0, @@ -2042,7 +2039,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }, ]; @@ -2091,7 +2087,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }]; let records = SamWriter::build_transcriptome_records( @@ -2182,7 +2177,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }; let mate2_transcript = Transcript { @@ -2204,7 +2198,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2], }; let mate_seq = vec![0, 1, 2, 3]; @@ -2300,7 +2293,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 4], }; // Mate2 at position 4 (chr_start=0, so per-chr pos = 5) @@ -2323,7 +2315,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 3], }; let mate_seq = vec![0; 4]; @@ -2400,7 +2391,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], }; let t2 = Transcript { @@ -2422,7 +2412,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], }; let mut rec1 = RecordBuf::default(); @@ -2480,7 +2469,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], }; let t2 = Transcript { @@ -2502,7 +2490,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], }; let mut rec1 = RecordBuf::default(); @@ -2540,7 +2527,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }; let read_seq = vec![0, 1, 2, 3]; @@ -2659,7 +2645,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }; let read_seq = vec![0, 1, 2, 3]; @@ -2713,7 +2698,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 4], }, Transcript { chr_idx: 0, @@ -2728,7 +2712,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 4], }, Transcript { chr_idx: 0, @@ -2743,7 +2726,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 4], }, ]; @@ -2802,7 +2784,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 4], }; let read_seq = vec![0, 1, 2, 3]; @@ -2853,7 +2834,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 4], }; // Two alignments tied for the best score (100), one strictly worse (98). let transcripts = vec![mk(0, 100), mk(2, 98), mk(4, 100)]; @@ -2912,7 +2892,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 4], }, Transcript { chr_idx: 0, @@ -2927,7 +2906,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 4], }, ]; @@ -2981,7 +2959,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 4], }; let mate2 = Transcript { genome_start: 120, @@ -3037,7 +3014,6 @@ mod tests { n_junction: 1, junction_motifs: vec![SpliceMotif::GtAg], junction_annotated: vec![false], - read_seq: vec![0; 4], }; let read_seq = vec![0, 1, 2, 3]; @@ -3083,7 +3059,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 4], }; let read_seq = vec![0, 1, 2, 3]; @@ -3133,7 +3108,6 @@ mod tests { n_junction: 1, junction_motifs: vec![SpliceMotif::CtAc], junction_annotated: vec![false], - read_seq: vec![0; 4], }; let read_seq = vec![0, 1, 2, 3]; @@ -3185,7 +3159,6 @@ mod tests { n_junction: 2, junction_motifs: vec![SpliceMotif::GtAg, SpliceMotif::CtAc], // +strand and -strand junction_annotated: vec![false, false], - read_seq: vec![0; 4], }; let read_seq = vec![0, 1, 2, 3]; @@ -3235,7 +3208,6 @@ mod tests { n_junction: 1, junction_motifs: vec![SpliceMotif::GtAg], junction_annotated: vec![false], - read_seq: vec![0; 4], }; let read_seq = vec![0, 1, 2, 3]; @@ -3289,7 +3261,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0; 4], }) .collect(); @@ -3362,7 +3333,6 @@ mod tests { n_junction: 1, junction_motifs: vec![SpliceMotif::GtAg], junction_annotated: vec![false], - read_seq: vec![], }; let jm = build_jm_tag(&transcript); @@ -3391,7 +3361,6 @@ mod tests { n_junction: 1, junction_motifs: vec![SpliceMotif::GtAg], junction_annotated: vec![true], - read_seq: vec![], }; let jm = build_jm_tag(&transcript); @@ -3416,7 +3385,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], }; assert!(build_jm_tag(&transcript).is_none()); @@ -3444,7 +3412,6 @@ mod tests { n_junction: 2, junction_motifs: vec![SpliceMotif::GtAg, SpliceMotif::CtAc], junction_annotated: vec![true, false], - read_seq: vec![], }; let jm = build_jm_tag(&transcript); @@ -3473,7 +3440,6 @@ mod tests { n_junction: 1, junction_motifs: vec![SpliceMotif::GtAg], junction_annotated: vec![false], - read_seq: vec![], }; // chr_start=0, genome_start=100, intron starts at 125, ends at 324 @@ -3499,7 +3465,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], }; assert!(build_ji_tag(&transcript, 0).is_none()); @@ -3523,7 +3488,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }; // Read exactly matches genome[0..4] = ACGT @@ -3550,7 +3514,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 0, 2, 0], // A,A,G,A vs genome A,C,G,T }; // Position 1: read=A, ref=C → mismatch (C in MD) @@ -3582,7 +3545,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 0, 1], // AC + AC (genome AC^GT AC) }; // Read: A,C,[del G,T],A,C @@ -3613,7 +3575,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 3, 3, 2, 3], // AC + TT(ins) + GT }; // Insertions are invisible in MD — just match counts @@ -3644,7 +3605,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 0, 2, 3, 0, 1, 0, 0], // XX + GTAC + XX }; // Soft clips don't appear in MD @@ -3672,7 +3632,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }; let read_seq = vec![0, 1, 2, 3]; @@ -3730,7 +3689,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }; // Mate2: reverse, chr 0, pos 4 @@ -3753,7 +3711,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2], }; let seq = vec![0, 1, 2, 3]; @@ -3839,7 +3796,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }; // Mate2: score=80, 2 mismatches, 1 deletion @@ -3866,7 +3822,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2], }; let seq1 = vec![0, 1, 2, 3]; @@ -3968,7 +3923,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }; let mate2_trans = Transcript { @@ -3990,7 +3944,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2], }; let seq = vec![0, 1, 2, 3]; @@ -4063,7 +4016,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }; let read_seq = vec![0, 1, 2, 3]; @@ -4142,7 +4094,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }; let read_seq = vec![0, 1, 2, 3]; @@ -4218,7 +4169,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }; let read_seq = vec![0, 1, 2, 3]; @@ -4359,7 +4309,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }; let read_seq = vec![0, 1, 2, 3]; @@ -4414,7 +4363,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }; let read_seq = vec![0, 1, 2, 3]; @@ -4501,7 +4449,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }; let mate1_seq = vec![0, 1, 2, 3]; // ACGT @@ -4574,7 +4521,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }; let mate1_seq = vec![0, 1, 2, 3]; @@ -4660,7 +4606,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![0, 1, 2, 3], }; let mate1_seq = vec![0, 1, 2, 3]; @@ -4732,7 +4677,6 @@ mod tests { n_junction: 1, junction_motifs: vec![SpliceMotif::GtAg], junction_annotated: vec![false], - read_seq: vec![0; 4], } } diff --git a/src/quant/mod.rs b/src/quant/mod.rs index ec178d6..7df5ddc 100644 --- a/src/quant/mod.rs +++ b/src/quant/mod.rs @@ -629,7 +629,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], } } @@ -745,7 +744,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], }; assert_eq!( ann.overlapping_genes(&t), diff --git a/src/quant/transcriptome.rs b/src/quant/transcriptome.rs index c68d936..2feb761 100644 --- a/src/quant/transcriptome.rs +++ b/src/quant/transcriptome.rs @@ -1133,7 +1133,6 @@ fn align_to_one_transcript( n_junction: 0, junction_motifs: Vec::new(), junction_annotated: Vec::new(), - read_seq: Vec::new(), }) } @@ -2054,7 +2053,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], } } diff --git a/src/signal/mod.rs b/src/signal/mod.rs index 260316c..b923622 100644 --- a/src/signal/mod.rs +++ b/src/signal/mod.rs @@ -200,7 +200,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], } } diff --git a/src/solo/gene.rs b/src/solo/gene.rs index 3310171..5d4e4a6 100644 --- a/src/solo/gene.rs +++ b/src/solo/gene.rs @@ -382,7 +382,6 @@ mod tests { n_junction: 0, junction_motifs: Vec::new(), junction_annotated: Vec::new(), - read_seq: Vec::new(), } } diff --git a/src/stats.rs b/src/stats.rs index b5b4a03..636af48 100644 --- a/src/stats.rs +++ b/src/stats.rs @@ -797,7 +797,6 @@ mod tests { n_junction: 1, junction_motifs: vec![SpliceMotif::GtAg], junction_annotated: vec![true], - read_seq: vec![0; 110], }; stats.record_transcript_stats(&transcript); @@ -918,7 +917,6 @@ mod tests { SpliceMotif::NonCanonical, // motif[0] ], junction_annotated: vec![true, false, true, false], - read_seq: vec![0; 100], }; stats.record_transcript_stats(&transcript); diff --git a/src/wasp/mod.rs b/src/wasp/mod.rs index df79eab..78f5051 100644 --- a/src/wasp/mod.rs +++ b/src/wasp/mod.rs @@ -608,7 +608,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], } } From 7d7b367db8fb7708b0ce4aefcf5d2c196ef2c38e Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 30 Jul 2026 23:52:21 +0200 Subject: [PATCH 04/10] docs(changelog): record the Transcript::read_seq removal Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index f103273..43a8305 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -136,3 +136,9 @@ Sections commonly used: Features, Bug fixes, Other changes. needs random access to the SA in RAM. Initial release of Rust rewrite of STAR. +### Other changes + +- Removed `Transcript::read_seq`, a public field that was filled with a + copy of the read at every finalised alignment and never read. **API + removal.** Output is unchanged. + From a4d62262a57b2e8979508fbcd2d8f55864d6a67d Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 00:04:04 +0200 Subject: [PATCH 05/10] test: an A/B harness that refuses to average over a busy machine Every perf measurement on this project needs the same three guards, and I got each of them wrong at least once in a single session: An earlier harness sampled CPU idle once, before the first run. Three separate measurements then drifted below the threshold mid-series and it kept going, leaving me to spot the contaminated rounds by eye in the output. This samples idle before and after every run, drops a round if any of the four samples falls short, and prints how many it dropped, so a median over four surviving rounds cannot be mistaken for a median over six. The check is on CPU idle rather than load average, because load average is an exponential average over minutes: it refused to measure at 2.24 on a machine whose cores were all free. It reports the spread within each side next to the difference between the medians, and says so in as many words when the difference is smaller. Two changes I measured looked like wins on medians alone and were inside the spread. Both sides run as ./rustar-aligner with --outFileNamePrefix ./ from inside their own directory, because the @PG CL: line records argv verbatim: running ./old against ./new is enough to make the output differ, which cost me two false "output is not neutral" alarms. The header documents the fourth trap, which no script can enforce: timing a total hides the part that changed. BAM writing is 1-4% of a yeast run, so a total dominated by alignment cannot resolve a change to the writer. Run the None configuration alongside and read the difference. Co-Authored-By: Claude Opus 5 (1M context) --- test/bench_ab.sh | 148 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 148 insertions(+) create mode 100755 test/bench_ab.sh diff --git a/test/bench_ab.sh b/test/bench_ab.sh new file mode 100755 index 0000000..d10755b --- /dev/null +++ b/test/bench_ab.sh @@ -0,0 +1,148 @@ +#!/usr/bin/env bash +# Interleaved A/B of two rustar-aligner binaries. +# +# test/bench_ab.sh [threads] [rounds] +# +# Prints one line per round and a median over the rounds that were measured on +# a quiet machine, plus a count of the ones that were not. Read the whole +# output: a median over four surviving rounds out of six means something +# different from a median over six. +# +# BENCH_MIN_IDLE=88 percent CPU idle a round must hold to count +# BENCH_MODE="None" --outSAMtype words (quote the two-word forms) +# BENCH_ARGS="" extra flags passed to both sides +# +# --------------------------------------------------------------------------- +# Why this script exists, and what it refuses to do +# +# Three things went wrong often enough while measuring this aligner that they +# are worth encoding rather than remembering. +# +# 1. A busy machine does not add noise, it inverts results. An unrelated job +# saturating the cores lands on whichever side happens to run under it, and +# nothing in the timings shows that it happened: the series just drifts. So +# idle is sampled *before and after every single run*, not once at the +# start, and a round where either sample falls below the threshold is +# dropped and counted rather than quietly folded into the median. +# +# The check is on CPU idle, not load average. Load average is an +# exponential average over minutes, so it stays high long after the +# offending job is gone and refuses to measure on a machine that is now +# perfectly quiet. +# +# 2. Timing a total hides the part that changed. If the change is in the +# writer, run one configuration with `--outSAMtype None` alongside the BAM +# one: only the difference is the work under test. On yeast, BAM writing is +# 1-4% of the run, so a total dominated by alignment cannot resolve a change +# to it at any reachable scale. +# +# 3. Both sides must see byte-identical argv. The `@PG` `CL:` line records the +# command line verbatim, so running `./old` against `./new`, or writing to +# different `--outFileNamePrefix` directories, makes the output differ in +# the header and the compressed size with it. Each side is therefore copied +# to `rustar-aligner` inside its own directory and run with prefix `./`. +# --------------------------------------------------------------------------- +set -euo pipefail + +BIN_A=${1:?usage: bench_ab.sh [threads] [rounds]} +BIN_B=${2:?usage: bench_ab.sh [threads] [rounds]} +GENOME_DIR=${3:?usage: bench_ab.sh [threads] [rounds]} +READS=${4:?usage: bench_ab.sh [threads] [rounds]} +THREADS=${5:-16} +ROUNDS=${6:-6} + +MIN_IDLE=${BENCH_MIN_IDLE:-88} +MODE=${BENCH_MODE:-None} +EXTRA=${BENCH_ARGS:-} + +WORK=$(mktemp -d) +trap 'rm -rf "$WORK"' EXIT + +NAME_A=$(basename "$BIN_A") +NAME_B=$(basename "$BIN_B") +[ "$NAME_A" = "$NAME_B" ] && { NAME_A="A:$NAME_A"; NAME_B="B:$NAME_B"; } + +for side in a b; do + mkdir -p "$WORK/$side" +done +cp "$BIN_A" "$WORK/a/rustar-aligner" +cp "$BIN_B" "$WORK/b/rustar-aligner" + +cpu_idle() { # percent idle, sampled over one second + top -l 2 -n 0 2>/dev/null | awk '/CPU usage/ {gsub("%","",$(NF-1)); v=$(NF-1)} END {print v+0}' +} + +# Returns " ". +run() { # $1 = side directory + local dir="$WORK/$1" + local before after t0 t1 + before=$(cpu_idle) + rm -rf "$dir/out" + mkdir -p "$dir/out" + t0=$(python3 -c 'import time; print(time.time())') + # shellcheck disable=SC2086 # MODE and EXTRA are deliberately word-split + (cd "$dir" && ./rustar-aligner --genomeDir "$GENOME_DIR" --readFilesIn "$READS" \ + --runThreadN "$THREADS" --outSAMtype $MODE $EXTRA \ + --outFileNamePrefix ./out/ >/dev/null 2>&1) + t1=$(python3 -c "import time; print(f'{time.time()-$t0:.2f}')") + after=$(cpu_idle) + echo "$t1 $before $after" +} + +echo "mode='$MODE' threads=$THREADS rounds=$ROUNDS min_idle=${MIN_IDLE}%" >&2 +run a >/dev/null # discarded: warms the page cache for the reads and the index + +kept_a=() +kept_b=() +dropped=0 + +for i in $(seq 1 "$ROUNDS"); do + # Flip the order on even rounds so a machine that drifts in one direction + # cannot favour whichever side runs first. + if (( i % 2 )); then + read -r ta ia1 ia2 <<<"$(run a)" + read -r tb ib1 ib2 <<<"$(run b)" + else + read -r tb ib1 ib2 <<<"$(run b)" + read -r ta ia1 ia2 <<<"$(run a)" + fi + + worst=$(printf '%s\n%s\n%s\n%s\n' "$ia1" "$ia2" "$ib1" "$ib2" | sort -n | head -1) + if awk -v i="$worst" -v m="$MIN_IDLE" 'BEGIN{exit !(i < m)}'; then + printf 'round%-2s %s=%ss %s=%ss DROPPED (idle fell to %s%%)\n' \ + "$i" "$NAME_A" "$ta" "$NAME_B" "$tb" "$worst" + dropped=$((dropped + 1)) + continue + fi + printf 'round%-2s %s=%ss %s=%ss idle>=%s%%\n' "$i" "$NAME_A" "$ta" "$NAME_B" "$tb" "$worst" + kept_a+=("$ta") + kept_b+=("$tb") +done + +median() { printf '%s\n' "$@" | sort -n | awk '{v[NR]=$1} END {print (NR%2) ? v[(NR+1)/2] : (v[NR/2]+v[NR/2+1])/2}'; } + +echo +if [ ${#kept_a[@]} -eq 0 ]; then + echo "every round was dropped: the machine was never quiet enough to measure." >&2 + echo "wait, or lower BENCH_MIN_IDLE and say so when reporting the numbers." >&2 + exit 1 +fi + +ma=$(median "${kept_a[@]}") +mb=$(median "${kept_b[@]}") +printf 'kept %d of %d rounds (%d dropped)\n' "${#kept_a[@]}" "$ROUNDS" "$dropped" +printf 'median %s=%ss %s=%ss delta=%s%%\n' \ + "$NAME_A" "$ma" "$NAME_B" "$mb" \ + "$(python3 -c "print(f'{($mb-$ma)/$ma*100:+.1f}')")" + +# A difference smaller than the spread within either side is not a result. +sa=$(python3 -c "v=sorted([$(IFS=,; echo "${kept_a[*]}")]); print(f'{v[-1]-v[0]:.2f}')") +sb=$(python3 -c "v=sorted([$(IFS=,; echo "${kept_b[*]}")]); print(f'{v[-1]-v[0]:.2f}')") +printf 'spread %s=%ss %s=%ss\n' "$NAME_A" "$sa" "$NAME_B" "$sb" +python3 - "$ma" "$mb" "$sa" "$sb" <<'PY' +import sys +ma, mb, sa, sb = (float(x) for x in sys.argv[1:5]) +if abs(mb - ma) < max(sa, sb): + print("\nthe difference between the medians is smaller than the spread within a side.") + print("that is not a result: report it as unmeasured rather than as a speedup.") +PY From d181ac3337e65020e715aa125df416389cf7a7b8 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 09:41:56 +0200 Subject: [PATCH 06/10] feat(solo): --soloCellReadStats CB writes CellReads.stats MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit One row per cell barcode with the fifteen counters STARsolo reports: how the barcode matched, whether the read mapped to one locus or several, whether it landed on a feature, where in the gene and on which strand, whether it was mitochondrial, and whether it reached the matrix. The per-cell UMI and gene totals come from the raw matrix rather than from the read counters, so they agree with what the matrix says by construction. Reads whose barcode never resolved are summed into a single `CBnotInPasslist` row instead of being dropped. That row is the reason the file is useful: it is the difference between "these cells look thin" and "most of the input never reached a cell at all". The region columns split by strand — an antisense read counts under `exonicAS` or `intronicAS`, never under `exonic` or `intronic`. `--genomeChrSetMitochondrial` names the chromosomes behind the `mito` column. Without it the column is zero throughout, which is honest: no chromosome was declared mitochondrial. D24 comes with it. STAR emits these rows by walking a libc++ `unordered_map`, which at these sizes is the reverse of each barcode's first appearance. That is reproduced, including across threads: the per-read accumulator merges in read order, so a threaded run writes the same file as a serial one. It stops being reproducible past the point where libc++ rehashes, since the order then depends on the bucket count. The values never differ, only which line they sit on. Recorded in docs-old/dev/divergences.md. --- CHANGELOG.md | 21 +++ src/params/mod.rs | 21 +++ src/solo/cell_reads.rs | 332 +++++++++++++++++++++++++++++++++++++++++ src/solo/count.rs | 21 +++ src/solo/mod.rs | 74 +++++++++ 5 files changed, 469 insertions(+) create mode 100644 src/solo/cell_reads.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index 43a8305..7b1a005 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -97,6 +97,27 @@ Sections commonly used: Features, Bug fixes, Other changes. glibc's malloc and per-thread heaps that return whole segments to the OS when abandoned, so allocator cache size stays bounded. +- **`--soloCellFilter EmptyDrops_CR` now uses CellRanger's actual + statistics.** The ambient profile is smoothed with Simple Good-Turing, + as CellRanger and STAR do, instead of an approximation that reserved + unseen mass from the singleton rate and spread the remainder in + proportion to raw counts. The Monte-Carlo null is drawn with libc++'s + `std::mt19937` and `std::discrete_distribution`, seeded + `19760110 * (isim + 1)` per simulation as STAR seeds it, replacing a + SplitMix64 stream that could not agree with STAR's over an arbitrary + number of draws. Cell calls move as a result. + +- **`--soloCellReadStats CB`** writes `Solo.out//CellReads.stats`: + one row per cell barcode with fifteen counters describing what + happened to its reads — barcode match quality, unique or multi + genomic mapping, feature assignment, exonic/intronic and their + antisense counterparts, mitochondrial, and whether the read reached + the matrix — plus the per-cell UMI and gene totals. Reads whose + barcode never resolved are summed into a `CBnotInPasslist` row rather + than dropped, so the columns account for the whole input. + `--genomeChrSetMitochondrial` names the chromosomes behind the `mito` + column. + ### Bug fixes - `--runThreadN 1` ran on every logical core instead of on one. The diff --git a/src/params/mod.rs b/src/params/mod.rs index 0536a85..1b7bf48 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1113,6 +1113,17 @@ pub struct Parameters { #[arg(long = "soloCBmatchWLtype", default_value = "1MM_multi")] pub solo_cb_match_wl_type: String, + /// `CB`: write `Solo.out//CellReads.stats`, a per-cell-barcode + /// summary of what happened to the reads carrying it. `None` (the default) + /// writes nothing. + #[arg(long = "soloCellReadStats", default_value = "None")] + pub solo_cell_read_stats: String, + + /// Chromosome names treated as mitochondrial, for the `mito` column of + /// `CellReads.stats`. `-` (the default) names none. + #[arg(long = "genomeChrSetMitochondrial", num_args = 1.., default_values_t = vec!["-".to_string()])] + pub genome_chr_set_mitochondrial: Vec, + /// Cell-calling / matrix filtering: None, CellRanger2.2, EmptyDrops_CR, TopCells. #[arg(long = "soloCellFilter", num_args = 1.., default_values_t = vec!["CellRanger2.2".to_string(), "3000".to_string(), "0.99".to_string(), "10".to_string()])] pub solo_cell_filter: Vec, @@ -1710,6 +1721,16 @@ impl Parameters { )); } } + // --soloCellReadStats: `CB` is the only value STAR defines. + if !matches!(params.solo_cell_read_stats.as_str(), "CB" | "None") { + return Err(command.error( + ErrorKind::InvalidValue, + format!( + "unknown --soloCellReadStats '{}'; expected CB or None", + params.solo_cell_read_stats + ), + )); + } // Validate --clipAdapterType. if !matches!( params.clip_adapter_type.as_str(), diff --git a/src/solo/cell_reads.rs b/src/solo/cell_reads.rs new file mode 100644 index 0000000..e17a540 --- /dev/null +++ b/src/solo/cell_reads.rs @@ -0,0 +1,332 @@ +//! `--soloCellReadStats CB`: the per-cell-barcode read summary STARsolo writes +//! as `Solo.out//CellReads.stats`. +//! +//! One row per cell barcode, fifteen counters describing what happened to the +//! reads carrying it — how the barcode matched, whether the read mapped +//! uniquely, whether it landed on a feature, where in the gene, and whether it +//! reached the matrix — plus the per-cell UMI and gene totals. The reads whose +//! barcode never resolved are not dropped; they are summed into a single +//! `CBnotInPasslist` row, so the columns account for every read rather than +//! only the ones that succeeded. +//! +//! # D24: row order +//! +//! STAR iterates a libc++ `std::unordered_map` to emit these rows, so the order +//! is a hash-table walk, not a sort. For the small maps this produces, libc++ +//! chains new entries at the head of their bucket and walks buckets in order, +//! which comes out as the reverse of first appearance in read order. That is +//! what this reproduces. +//! +//! It is not reproducible in general: past the load factor libc++ rehashes, and +//! the order after a rehash depends on the bucket count, which depends on how +//! many distinct barcodes were seen. At that size the order diverges. The +//! **values never do** — only which line they appear on. A consumer that reads +//! this file by barcode rather than by position is unaffected either way, and +//! sorting by barcode is the only stable thing to do with it. + +use std::collections::{BTreeMap, HashSet}; +use std::fmt::Write as _; + +/// What happened to one read, as the fourteen optional flags STAR tracks. +/// +/// `cbMatch` is not here: every read that reaches the accumulator matched a +/// barcode well enough to be attributed somewhere, so it is always set. +/// +/// Fourteen bools rather than a bitfield because they are written once per read +/// and read once per fold, and the names are the column names of the file. +#[allow(clippy::struct_excessive_bools)] +#[derive(Debug, Clone, Copy, Default)] +pub struct CellReadFlag { + /// The barcode was an exact whitelist hit. + pub cb_perfect: bool, + /// Corrected via a single one-mismatch neighbour. + pub cb_mm_unique: bool, + /// Corrected via several one-mismatch neighbours, resolved by the posterior. + pub cb_mm_multiple: bool, + /// The read mapped to exactly one genomic locus. + pub genome_u: bool, + /// The read mapped to several genomic loci. + pub genome_m: bool, + /// It fell on exactly one feature (gene). + pub feature_u: bool, + /// It fell on several features. + pub feature_m: bool, + /// Exonic, on the annotated strand. + pub exonic: bool, + /// Intronic, on the annotated strand. + pub intronic: bool, + /// Exonic, antisense to the annotation. + pub exonic_as: bool, + /// Intronic, antisense to the annotation. + pub intronic_as: bool, + /// On a chromosome named by `--genomeChrSetMitochondrial`. + pub mito: bool, + /// Counted into the unique-gene matrix. + pub counted_u: bool, + /// Counted through the multi-gene distribution. + pub counted_m: bool, +} + +/// The accumulator behind `CellReads.stats`. +#[derive(Debug, Default, Clone)] +pub struct CellReadStats { + /// Whitelist index to its fifteen counters. + cells: BTreeMap, + /// The single bucket for reads whose barcode did not resolve. + no_cb: [u64; 15], + /// Whitelist indices in first-appearance order; emitted reversed (D24). + order: Vec, + seen: HashSet, +} + +impl CellReadStats { + pub fn new() -> Self { + Self::default() + } + + fn fold(v: &mut [u64; 15], f: &CellReadFlag) { + v[0] += 1; // cbMatch: every read that got this far + for (i, set) in [ + f.cb_perfect, + f.cb_mm_unique, + f.cb_mm_multiple, + f.genome_u, + f.genome_m, + f.feature_u, + f.feature_m, + f.exonic, + f.intronic, + f.exonic_as, + f.intronic_as, + f.mito, + f.counted_u, + f.counted_m, + ] + .into_iter() + .enumerate() + { + if set { + v[i + 1] += 1; + } + } + } + + /// Record a read that resolved to whitelist cell `cb`. + pub fn add_cell(&mut self, cb: u32, flag: &CellReadFlag) { + if self.seen.insert(cb) { + self.order.push(cb); + } + Self::fold(self.cells.entry(cb).or_insert([0; 15]), flag); + } + + /// Record a read whose barcode did not resolve, or whose UMI was rejected. + pub fn add_no_cb(&mut self, flag: &CellReadFlag) { + Self::fold(&mut self.no_cb, flag); + } + + /// Merge another accumulator in, preserving first-appearance order: the + /// per-thread partials are merged in read order, so the combined order is + /// the order a single-threaded run would have produced. + pub fn merge(&mut self, other: &Self) { + for &cb in &other.order { + if self.seen.insert(cb) { + self.order.push(cb); + } + let dst = self.cells.entry(cb).or_insert([0; 15]); + for (d, s) in dst.iter_mut().zip(&other.cells[&cb]) { + *d += s; + } + } + for (d, s) in self.no_cb.iter_mut().zip(&other.no_cb) { + *d += s; + } + } + + /// Render the file. `umi_gene` gives each cell its final + /// `(nUMIunique, nGenesUnique)`; a cell absent from it prints zeros. + /// `barcode_of` renders a whitelist index as its barcode string. + pub fn render( + &self, + barcode_of: impl Fn(u32) -> String, + umi_gene: &BTreeMap, + ) -> String { + let mut s = String::from( + "CB\tcbMatch\tcbPerfect\tcbMMunique\tcbMMmultiple\tgenomeU\tgenomeM\tfeatureU\t\ + featureM\texonic\tintronic\texonicAS\tintronicAS\tmito\tcountedU\tcountedM\t\ + nUMIunique\tnGenesUnique\tnUMImulti\tnGenesMulti\n", + ); + s.push_str("CBnotInPasslist"); + for v in &self.no_cb { + s.push('\t'); + s.push_str(&v.to_string()); + } + s.push_str("\t0\t0\t0\t0\n"); + // Reverse first-appearance order — see D24 in the module docs. + for &cb in self.order.iter().rev() { + s.push_str(&barcode_of(cb)); + for v in &self.cells[&cb] { + s.push('\t'); + s.push_str(&v.to_string()); + } + let (n_umi, n_gene) = umi_gene.get(&cb).copied().unwrap_or((0, 0)); + // nUMImulti and nGenesMulti stay zero: multi-gene UMIs are not + // collapsed into per-cell totals here. + let _ = writeln!(s, "\t{n_umi}\t{n_gene}\t0\t0"); + } + s + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn flag_perfect_counted() -> CellReadFlag { + CellReadFlag { + cb_perfect: true, + genome_u: true, + feature_u: true, + exonic: true, + counted_u: true, + ..Default::default() + } + } + + fn header_and_rows(text: &str) -> (Vec<&str>, Vec<&str>) { + let mut lines = text.lines(); + let header: Vec<&str> = lines.next().unwrap().split('\t').collect(); + (header, lines.collect()) + } + + /// Every column has a name, and every row has a value under each of them. + #[test] + fn each_row_fills_the_header() { + let mut st = CellReadStats::new(); + st.add_cell(0, &flag_perfect_counted()); + st.add_no_cb(&CellReadFlag::default()); + let out = st.render(|i| format!("CB{i}"), &BTreeMap::new()); + let (header, rows) = header_and_rows(&out); + assert_eq!(header.len(), 20); + assert_eq!(rows.len(), 2, "the passlist-miss row plus one cell"); + for row in rows { + assert_eq!(row.split('\t').count(), header.len()); + } + } + + /// A read is counted once under `cbMatch` and once under each flag it sets, + /// so the flag columns can exceed neither `cbMatch` nor each other's logic. + #[test] + fn flags_accumulate_per_read() { + let mut st = CellReadStats::new(); + for _ in 0..3 { + st.add_cell(7, &flag_perfect_counted()); + } + st.add_cell( + 7, + &CellReadFlag { + cb_mm_unique: true, + genome_m: true, + ..Default::default() + }, + ); + let out = st.render(|i| format!("CB{i}"), &BTreeMap::new()); + let row: Vec<&str> = out.lines().nth(2).unwrap().split('\t').collect(); + assert_eq!(row[0], "CB7"); + assert_eq!(row[1], "4", "cbMatch counts every read"); + assert_eq!(row[2], "3", "cbPerfect"); + assert_eq!(row[3], "1", "cbMMunique"); + assert_eq!(row[5], "3", "genomeU"); + assert_eq!(row[6], "1", "genomeM"); + assert_eq!(row[14], "3", "countedU"); + } + + /// Reads whose barcode never resolved are summed rather than dropped, so + /// the file accounts for the whole input. + #[test] + fn unresolved_reads_land_in_the_passlist_miss_row() { + let mut st = CellReadStats::new(); + st.add_cell(0, &flag_perfect_counted()); + for _ in 0..5 { + st.add_no_cb(&CellReadFlag { + genome_u: true, + ..Default::default() + }); + } + let out = st.render(|i| format!("CB{i}"), &BTreeMap::new()); + let row: Vec<&str> = out.lines().nth(1).unwrap().split('\t').collect(); + assert_eq!(row[0], "CBnotInPasslist"); + assert_eq!(row[1], "5"); + assert_eq!(row[5], "5", "genomeU"); + assert_eq!(&row[16..], ["0", "0", "0", "0"], "no UMI columns for it"); + } + + /// D24: rows come out in reverse first-appearance order, which is what + /// STAR's libc++ hash-map walk produces at these sizes. + #[test] + fn rows_are_emitted_in_reverse_first_appearance_order() { + let mut st = CellReadStats::new(); + for cb in [4u32, 1, 9] { + st.add_cell(cb, &flag_perfect_counted()); + } + st.add_cell(1, &flag_perfect_counted()); // seen again: order unchanged + let out = st.render(|i| format!("CB{i}"), &BTreeMap::new()); + let cbs: Vec<&str> = out + .lines() + .skip(2) + .map(|l| l.split('\t').next().unwrap()) + .collect(); + assert_eq!(cbs, ["CB9", "CB1", "CB4"]); + } + + /// The UMI and gene totals come from the final matrix, not from the read + /// counters, so a cell missing from that map prints zeros rather than + /// inheriting a read count. + #[test] + fn umi_and_gene_totals_come_from_the_matrix() { + let mut st = CellReadStats::new(); + st.add_cell(2, &flag_perfect_counted()); + st.add_cell(3, &flag_perfect_counted()); + let mut umi_gene = BTreeMap::new(); + umi_gene.insert(2u32, (17u32, 5u32)); + let out = st.render(|i| format!("CB{i}"), &umi_gene); + let rows: Vec> = out + .lines() + .skip(2) + .map(|l| l.split('\t').collect()) + .collect(); + // Reverse order: CB3 first, then CB2. + assert_eq!(rows[0][0], "CB3"); + assert_eq!(&rows[0][16..18], ["0", "0"]); + assert_eq!(rows[1][0], "CB2"); + assert_eq!(&rows[1][16..18], ["17", "5"]); + } + + /// Merging partials keeps first-appearance order, so a threaded run and a + /// serial one write the same file. + #[test] + fn merging_partials_preserves_order_and_sums() { + let mut a = CellReadStats::new(); + a.add_cell(5, &flag_perfect_counted()); + a.add_cell(2, &flag_perfect_counted()); + let mut b = CellReadStats::new(); + b.add_cell(2, &flag_perfect_counted()); + b.add_cell(8, &flag_perfect_counted()); + b.add_no_cb(&CellReadFlag::default()); + a.merge(&b); + + let out = a.render(|i| format!("CB{i}"), &BTreeMap::new()); + let cbs: Vec<&str> = out + .lines() + .skip(2) + .map(|l| l.split('\t').next().unwrap()) + .collect(); + assert_eq!(cbs, ["CB8", "CB2", "CB5"]); + let cb2: Vec<&str> = out + .lines() + .find(|l| l.starts_with("CB2\t")) + .unwrap() + .split('\t') + .collect(); + assert_eq!(cb2[1], "2", "the two partials' reads are summed"); + } +} diff --git a/src/solo/count.rs b/src/solo/count.rs index 9af3edf..cd45309 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -1224,6 +1224,27 @@ pub fn write_gene_matrix( if gzip { " [gzip]" } else { "" }, ); + // `--soloCellReadStats CB`: the per-cell read summary, alongside the + // raw matrix because its UMI and gene columns are the raw totals. + if let Some(cell_stats) = &ctx.cell_read_stats { + let umi_gene: std::collections::BTreeMap = mstats + .cells + .iter() + .map(|c| (c.cb, (c.n_umis as u32, c.n_genes))) + .collect(); + let path = feature_dir.join("CellReads.stats"); + let text = cell_stats.lock().unwrap().render( + |cb| { + ctx.whitelist + .barcode_string(cb) + .unwrap_or_else(|| cb.to_string()) + }, + &umi_gene, + ); + std::fs::write(&path, text).map_err(|e| Error::io(e, &path))?; + log::info!("STARsolo: wrote {}/CellReads.stats", feature.dir_name()); + } + // Filtered (cell-called) matrix per --soloCellFilter. EmptyDrops_CR runs // the Monte-Carlo rescue (needs the per-cell profiles in the body). let called = if params diff --git a/src/solo/mod.rs b/src/solo/mod.rs index 2253dee..e53f4f4 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -9,6 +9,7 @@ //! The barcode read is the SECOND `--readFilesIn` file (STAR convention: //! `--readFilesIn cDNA_read barcode_read`). It is never aligned — only parsed. +pub mod cell_reads; pub mod count; pub mod gene; pub mod smartseq; @@ -564,6 +565,13 @@ pub struct SoloContext { /// `--soloMultiMappers` includes a non-`Unique` method → capture gene- /// ambiguous reads for distribution into `UniqueAndMult-*.mtx`. pub want_multi: bool, + /// `--soloCellReadStats CB`: the per-cell read summary, or `None` when the + /// flag is off. Behind a mutex like the other per-read collections; the + /// lock is only taken when the flag asks for the file. + pub cell_read_stats: Option>, + /// Chromosome indices named by `--genomeChrSetMitochondrial`, for the + /// `mito` column. + pub mito_chr: std::collections::HashSet, } /// Per-region read tallies for the `Summary.csv` mapping funnel (uniquely-mapped @@ -695,6 +703,14 @@ impl SoloContext { velocyto_enabled, velocyto_records: Mutex::new(Vec::new()), want_multi, + cell_read_stats: (params.solo_cell_read_stats == "CB") + .then(|| Mutex::new(crate::solo::cell_reads::CellReadStats::new())), + mito_chr: params + .genome_chr_set_mitochondrial + .iter() + .filter(|n| n.as_str() != "-") + .filter_map(|n| genome.chr_name.iter().position(|c| c == n)) + .collect(), }) } @@ -871,9 +887,67 @@ impl SoloContext { fo }) .collect(); + + self.record_cell_read( + cb_resolved, + &cb_match, + n_loci, + &class, + cdna_transcripts, + &out, + ); out } + /// Fold one read into `CellReads.stats`, when `--soloCellReadStats CB` + /// asked for it. + /// + /// Called at the end of read processing so the counted flags reflect what + /// the read actually produced, rather than what it looked eligible for. + #[allow(clippy::too_many_arguments)] + fn record_cell_read( + &self, + cb_resolved: Option, + cb_match: &CbMatch, + n_loci: usize, + class: &crate::solo::gene::ReadClass, + transcripts: &[Transcript], + out: &SoloReadOutcome, + ) { + let Some(stats) = &self.cell_read_stats else { + return; + }; + let counted_u = out.per_feature.iter().any(|f| f.record.is_some()); + let counted_m = out.per_feature.iter().any(|f| f.multi_gene.is_some()); + let feature_u = counted_u || out.per_feature.iter().any(|f| f.multi.is_some()); + let flag = crate::solo::cell_reads::CellReadFlag { + cb_perfect: matches!(cb_match, CbMatch::Exact(_)), + cb_mm_unique: matches!(cb_match, CbMatch::Corrected(_)), + cb_mm_multiple: matches!(cb_match, CbMatch::Multi(_)), + genome_u: n_loci == 1, + genome_m: n_loci > 1, + feature_u, + feature_m: counted_m, + // The region columns split by strand: STAR reports an antisense + // read under `exonicAS`/`intronicAS`, not under `exonic`/`intronic`. + exonic: !class.antisense && class.region == Some(Region::Exonic), + intronic: !class.antisense && class.region == Some(Region::Intronic), + exonic_as: class.antisense && class.region == Some(Region::Exonic), + intronic_as: class.antisense && class.region == Some(Region::Intronic), + mito: !self.mito_chr.is_empty() + && transcripts + .iter() + .any(|t| self.mito_chr.contains(&t.chr_idx)), + counted_u, + counted_m, + }; + let mut stats = stats.lock().unwrap(); + match cb_resolved { + Some(cb) => stats.add_cell(cb, &flag), + None => stats.add_no_cb(&flag), + } + } + /// Process one 5' paired-end solo read (`--soloBarcodeMate 1`): the barcode is /// from mate 1, and both mates align as a pair. Genes are assigned from the /// union of both mates evaluated against the pair's (mate 1's) transcription From 4f7b6fe1f36edcc3aaccf3cbde8a7c28894fccc3 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 10:24:56 +0200 Subject: [PATCH 07/10] docs: record the CellReads.stats row order in DIVERGENCE.md Section 3.2, in the format CONTRIBUTING.md asks for. --- DIVERGENCE.md | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/DIVERGENCE.md b/DIVERGENCE.md index bd957ed..f1d9f98 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -65,6 +65,20 @@ On the 10k yeast PE benchmark, 4 reads differ in alignment score (AS) because ST --- +### 3.2 `CellReads.stats` row order + +**What STAR does.** `--soloCellReadStats CB` emits its rows by iterating a libc++ `std::unordered_map`, so the order is a hash-table walk rather than a sort. At the map sizes this produces, libc++ chains new entries at the head of their bucket and walks buckets in order, which comes out as the reverse of each barcode's first appearance in read order. + +**What rustar-aligner does.** Emits that same reverse-first-appearance order, including across threads: the per-read accumulator merges in read order, so a threaded run writes the same file as a serial one. + +**Why.** Reproducing the order where it is reproducible costs nothing and keeps a byte-comparison against STAR usable on the sizes where it can work at all. + +**Impact.** Past libc++'s load factor the map rehashes, and the order then depends on the bucket count, which depends on how many distinct barcodes were seen; beyond that size the order diverges. The **values never do** — only which line they appear on. Reading the file by barcode rather than by position is unaffected either way. + +**Source.** `src/solo/cell_reads.rs`, locked by `rows_are_emitted_in_reverse_first_appearance_order` and `merging_partials_preserves_order_and_sums`. STAR: `SoloFeature_statsOutput.cpp`. + +--- + ## 4. Implementation divergences (no intended output difference) These differ in *how* a result is produced, not *what* is produced. They are documented so a reviewer chasing a discrepancy knows the mechanism differs by design. From 1a343b2791a9cac11abf4ddaff168e5466e047dc Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 11:31:17 +0200 Subject: [PATCH 08/10] docs(changelog): keep only this PR's entry CONTRIBUTING.md requires the description to match the code; the entries for the other themes split out of #152 belong to their own PRs. --- CHANGELOG.md | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7b1a005..e2817c3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -97,16 +97,6 @@ Sections commonly used: Features, Bug fixes, Other changes. glibc's malloc and per-thread heaps that return whole segments to the OS when abandoned, so allocator cache size stays bounded. -- **`--soloCellFilter EmptyDrops_CR` now uses CellRanger's actual - statistics.** The ambient profile is smoothed with Simple Good-Turing, - as CellRanger and STAR do, instead of an approximation that reserved - unseen mass from the singleton rate and spread the remainder in - proportion to raw counts. The Monte-Carlo null is drawn with libc++'s - `std::mt19937` and `std::discrete_distribution`, seeded - `19760110 * (isim + 1)` per simulation as STAR seeds it, replacing a - SplitMix64 stream that could not agree with STAR's over an arbitrary - number of draws. Cell calls move as a result. - - **`--soloCellReadStats CB`** writes `Solo.out//CellReads.stats`: one row per cell barcode with fifteen counters describing what happened to its reads — barcode match quality, unique or multi From caec836ac5ed3d3c7c6d8180ea2af2abe078965e Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 11:53:34 +0200 Subject: [PATCH 09/10] refactor(solo): drop CellReadStats::merge, which nothing calls Reads are folded in under a mutex, so there are no per-thread partials to merge; the function was reachable only from its own test. CONTRIBUTING.md rules out shipping a function no production path reaches, and the PR description claimed its test as evidence of thread-safety that the mutex actually provides. --- src/solo/cell_reads.rs | 50 +++--------------------------------------- 1 file changed, 3 insertions(+), 47 deletions(-) diff --git a/src/solo/cell_reads.rs b/src/solo/cell_reads.rs index e17a540..5d140a5 100644 --- a/src/solo/cell_reads.rs +++ b/src/solo/cell_reads.rs @@ -23,6 +23,9 @@ //! **values never do** — only which line they appear on. A consumer that reads //! this file by barcode rather than by position is unaffected either way, and //! sorting by barcode is the only stable thing to do with it. +//! +//! Reads are folded in under a mutex held by `SoloContext`, so the order is the +//! order reads were processed in regardless of thread count. use std::collections::{BTreeMap, HashSet}; use std::fmt::Write as _; @@ -124,24 +127,6 @@ impl CellReadStats { Self::fold(&mut self.no_cb, flag); } - /// Merge another accumulator in, preserving first-appearance order: the - /// per-thread partials are merged in read order, so the combined order is - /// the order a single-threaded run would have produced. - pub fn merge(&mut self, other: &Self) { - for &cb in &other.order { - if self.seen.insert(cb) { - self.order.push(cb); - } - let dst = self.cells.entry(cb).or_insert([0; 15]); - for (d, s) in dst.iter_mut().zip(&other.cells[&cb]) { - *d += s; - } - } - for (d, s) in self.no_cb.iter_mut().zip(&other.no_cb) { - *d += s; - } - } - /// Render the file. `umi_gene` gives each cell its final /// `(nUMIunique, nGenesUnique)`; a cell absent from it prints zeros. /// `barcode_of` renders a whitelist index as its barcode string. @@ -300,33 +285,4 @@ mod tests { assert_eq!(rows[1][0], "CB2"); assert_eq!(&rows[1][16..18], ["17", "5"]); } - - /// Merging partials keeps first-appearance order, so a threaded run and a - /// serial one write the same file. - #[test] - fn merging_partials_preserves_order_and_sums() { - let mut a = CellReadStats::new(); - a.add_cell(5, &flag_perfect_counted()); - a.add_cell(2, &flag_perfect_counted()); - let mut b = CellReadStats::new(); - b.add_cell(2, &flag_perfect_counted()); - b.add_cell(8, &flag_perfect_counted()); - b.add_no_cb(&CellReadFlag::default()); - a.merge(&b); - - let out = a.render(|i| format!("CB{i}"), &BTreeMap::new()); - let cbs: Vec<&str> = out - .lines() - .skip(2) - .map(|l| l.split('\t').next().unwrap()) - .collect(); - assert_eq!(cbs, ["CB8", "CB2", "CB5"]); - let cb2: Vec<&str> = out - .lines() - .find(|l| l.starts_with("CB2\t")) - .unwrap() - .split('\t') - .collect(); - assert_eq!(cb2[1], "2", "the two partials' reads are summed"); - } } From ae575c710d14c3a36c18073e1bc21655396c20a9 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 10:17:14 +0200 Subject: [PATCH 10/10] feat(solo): --runMode soloCellFiltering Cell-calls an existing raw count matrix without aligning anything, taking the raw directory and an output prefix as STAR does: `--runMode soloCellFiltering /path/to/raw/ /path/to/out/prefix`. Cell calling is a decision about a matrix, not about reads. Re-calling with different `--soloCellFilter` parameters should not mean re-aligning 400 million reads, and a matrix produced by another tool should be callable too. The matrix is streamed into the same temp-body form the align path builds, so `called_cells` and `emptydrops_called` are the identical code here and there rather than a second implementation free to drift. Counts are rounded on the way in: a multimapper matrix carries real values, and the filters work on UMI totals. `--runMode` becomes a token list, because that is what STAR's is: the mode followed by its arguments. The mode itself is now validated rather than falling back to `alignReads`, so a typo is refused instead of quietly running something else. The standalone `emptydrops` binary still exists and still carries its own copy of the algorithm, which no longer matches this one. Removing it means moving `test/solo_genefull_compare.py` and `test/solo_genefull_h5_compare.py` to the new mode first, so it is left alone here rather than broken. --- CHANGELOG.md | 7 ++ src/lib.rs | 5 +- src/params/mod.rs | 63 +++++++++--- src/solo/count.rs | 191 ++++++++++++++++++++++++++++++++++++ tests/alignment_features.rs | 90 +++++++++++++++++ 5 files changed, 341 insertions(+), 15 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e2817c3..9a63f97 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -107,6 +107,13 @@ Sections commonly used: Features, Bug fixes, Other changes. than dropped, so the columns account for the whole input. `--genomeChrSetMitochondrial` names the chromosomes behind the `mito` column. +- **`--runMode soloCellFiltering `** cell-calls + an existing raw count matrix without aligning anything. Cell calling + is a decision about a matrix, not about reads: re-calling with + different `--soloCellFilter` parameters should not mean re-aligning, + and a matrix produced elsewhere should be callable too. It streams + the matrix into the same form the align path produces, so the filters + are the identical code rather than a second implementation. ### Bug fixes diff --git a/src/lib.rs b/src/lib.rs index 8a07f5f..e5f4f4e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -55,7 +55,7 @@ pub fn run(params: &Parameters) -> anyhow::Result<()> { if let Some(hint) = cpu::upgrade_hint() { info!("{hint}"); } - info!("runMode: {}", params.run_mode); + info!("runMode: {}", params.run_mode_in.join(" ")); info!("runThreadN: {}", params.run_thread_n); // Configure the rayon global pool from `--runThreadN` **before** @@ -76,11 +76,12 @@ pub fn run(params: &Parameters) -> anyhow::Result<()> { .num_threads(params.run_thread_n.into()) .build_global(); - match params.run_mode { + match params.run_mode() { RunMode::GenomeGenerate => genome_generate(params), RunMode::AlignReads => align_reads(params), RunMode::InputAlignmentsFromBAM => bam_dedup::run(params), RunMode::LiftOver => liftover::run(params), + RunMode::SoloCellFiltering => crate::solo::count::run_cell_filtering(params), } } diff --git a/src/params/mod.rs b/src/params/mod.rs index 1b7bf48..4f08bde 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -37,6 +37,8 @@ pub enum RunMode { GenomeGenerate, InputAlignmentsFromBAM, LiftOver, + /// Cell-call an existing raw count matrix, without aligning anything. + SoloCellFiltering, } impl std::str::FromStr for RunMode { @@ -47,9 +49,10 @@ impl std::str::FromStr for RunMode { "genomeGenerate" => Ok(Self::GenomeGenerate), "inputAlignmentsFromBAM" => Ok(Self::InputAlignmentsFromBAM), "liftOver" => Ok(Self::LiftOver), + "soloCellFiltering" => Ok(Self::SoloCellFiltering), _ => Err(format!( "unknown runMode '{s}'; expected 'alignReads', 'genomeGenerate', \ - 'inputAlignmentsFromBAM', or 'liftOver'" + 'inputAlignmentsFromBAM', 'liftOver', or 'soloCellFiltering'" )), } } @@ -62,6 +65,7 @@ impl std::fmt::Display for RunMode { Self::GenomeGenerate => write!(f, "genomeGenerate"), Self::InputAlignmentsFromBAM => write!(f, "inputAlignmentsFromBAM"), Self::LiftOver => write!(f, "liftOver"), + Self::SoloCellFiltering => write!(f, "soloCellFiltering"), } } } @@ -465,9 +469,11 @@ impl std::fmt::Display for SoloType { )] pub struct Parameters { // ── Run ───────────────────────────────────────────────────────────── - /// Run mode: alignReads or genomeGenerate - #[arg(long = "runMode", default_value = "alignReads")] - pub run_mode: RunMode, + /// Run mode, plus its arguments. `--runMode soloCellFiltering` takes two + /// more tokens: the raw count-matrix directory and the output prefix + /// (STAR `SoloFeature_loadRawMatrix.cpp`). + #[arg(long = "runMode", num_args = 1.., default_values_t = vec!["alignReads".to_string()])] + pub run_mode_in: Vec, /// Number of threads #[arg(long = "runThreadN", default_value_t = NonZeroUsize::new(1).unwrap())] @@ -1173,6 +1179,16 @@ pub struct Parameters { } impl Parameters { + /// The run mode. Validation guarantees it parses, so this cannot fail + /// after `validate()`; before it, an unknown mode reads as `alignReads` + /// and validation is what rejects it. + pub fn run_mode(&self) -> RunMode { + self.run_mode_in + .first() + .and_then(|m| m.parse().ok()) + .unwrap_or(RunMode::AlignReads) + } + /// Build an output path by concatenating `suffix` onto `out_file_name_prefix`. pub fn output_path(&self, suffix: &str) -> PathBuf { PathBuf::from(format!("{}{suffix}", self.out_file_name_prefix)) @@ -1364,8 +1380,29 @@ impl Parameters { shlex::try_join(args.iter().map(AsRef::as_ref)).ok() }; + // The run mode itself must be one this build knows: an unrecognised + // one would otherwise fall through to alignReads and silently do + // something the user did not ask for. + if let Some(mode) = params.run_mode_in.first() + && mode.parse::().is_err() + { + return Err(command.error( + ErrorKind::InvalidValue, + mode.parse::().unwrap_err(), + )); + } + + // `--runMode soloCellFiltering `. + if params.run_mode() == RunMode::SoloCellFiltering && params.run_mode_in.len() < 3 { + return Err(command.error( + ErrorKind::WrongNumberOfValues, + "--runMode soloCellFiltering needs the raw count-matrix directory and the \ + output prefix: --runMode soloCellFiltering /path/to/raw/ /path/to/out/prefix", + )); + } + // genomeGenerate requires FASTA files - if params.run_mode == RunMode::GenomeGenerate && params.genome_fasta_files.is_empty() { + if params.run_mode() == RunMode::GenomeGenerate && params.genome_fasta_files.is_empty() { return Err(command.error( ErrorKind::MissingRequiredArgument, "--genomeFastaFiles is required when --runMode genomeGenerate", @@ -1382,7 +1419,7 @@ impl Parameters { // alignReads requires read files — except SmartSeq, which gets its reads // from --readFilesManifest instead. - if params.run_mode == RunMode::AlignReads + if params.run_mode() == RunMode::AlignReads && params.read_files_in.is_empty() && params.solo_type != SoloType::SmartSeq { @@ -1419,7 +1456,7 @@ impl Parameters { } // inputAlignmentsFromBAM: only --bamRemoveDuplicatesType is implemented so far - if params.run_mode == RunMode::InputAlignmentsFromBAM { + if params.run_mode() == RunMode::InputAlignmentsFromBAM { let dedup = params.bam_remove_duplicates_type.as_str(); if dedup == "-" { return Err(command.error( @@ -1442,7 +1479,7 @@ impl Parameters { } // liftOver requires a chain file and a GTF to lift - if params.run_mode == RunMode::LiftOver { + if params.run_mode() == RunMode::LiftOver { if params.genome_chain_files.is_empty() { return Err(command.error( ErrorKind::MissingRequiredArgument, @@ -1552,7 +1589,7 @@ impl Parameters { // validation time we can only enforce the genomeGenerate rule; // for alignReads, GenomeIndex::load checks for the on-disk files // and surfaces a clear error if neither source is available. - if params.run_mode == RunMode::GenomeGenerate + if params.run_mode() == RunMode::GenomeGenerate && params.quant_transcriptome_sam() && params.sjdb_gtf_file.is_none() { @@ -1563,7 +1600,7 @@ impl Parameters { } // ── STARsolo validation ───────────────────────────────────────── - if params.run_mode == RunMode::AlignReads && params.solo_enabled() { + if params.run_mode() == RunMode::AlignReads && params.solo_enabled() { // CB_UMI_Complex needs one CB position + whitelist per segment. if params.solo_type == SoloType::CbUmiComplex { if params.solo_cb_position.is_empty() { @@ -1905,7 +1942,7 @@ mod tests { #[test] fn defaults() { let p = try_parse(&["--readFilesIn", "reads.fq"]).unwrap(); - assert_eq!(p.run_mode, RunMode::AlignReads); + assert_eq!(p.run_mode(), RunMode::AlignReads); assert_eq!(p.run_thread_n, NonZeroUsize::new(1).unwrap()); assert_eq!(p.run_rng_seed, 777); assert_eq!(p.genome_dir, PathBuf::from("./GenomeDir")); @@ -1999,7 +2036,7 @@ mod tests { "11", ]) .unwrap(); - assert_eq!(p.run_mode, RunMode::GenomeGenerate); + assert_eq!(p.run_mode(), RunMode::GenomeGenerate); assert_eq!(p.genome_dir, PathBuf::from("/data/genome")); assert_eq!( p.genome_fasta_files, @@ -2039,7 +2076,7 @@ mod tests { "Basic", ]) .unwrap(); - assert_eq!(p.run_mode, RunMode::AlignReads); + assert_eq!(p.run_mode(), RunMode::AlignReads); assert_eq!(p.genome_dir, PathBuf::from("/idx/hg38")); assert_eq!( p.read_files_in, diff --git a/src/solo/count.rs b/src/solo/count.rs index cd45309..58f8524 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -1902,6 +1902,197 @@ fn write_barcodes_subset( Ok(()) } +/// `--runMode soloCellFiltering `: cell-call an +/// existing raw count matrix without aligning anything. +/// +/// STAR `SoloFeature_loadRawMatrix.cpp` plus the same `--soloCellFilter` used +/// at the end of a solo run. The point of the mode is that cell calling is a +/// decision about a matrix, not about reads: re-calling with different +/// parameters should not mean re-aligning, and a matrix produced elsewhere +/// should be callable too. +pub fn run_cell_filtering(params: &crate::params::Parameters) -> anyhow::Result<()> { + let raw_dir = Path::new(¶ms.run_mode_in[1]); + let out_prefix = ¶ms.run_mode_in[2]; + + let find = |base: &str| -> Result { + for name in [base.to_string(), format!("{base}.gz")] { + let p = raw_dir.join(&name); + if p.exists() { + return Ok(p); + } + } + Err(Error::Parameter(format!( + "{}: no {base} (or {base}.gz) in the raw matrix directory", + raw_dir.display() + ))) + }; + + let barcodes = read_first_column(&find("barcodes.tsv")?)?; + let features_path = find("features.tsv")?; + let matrix_path = find("matrix.mtx")?; + + // Stream the matrix into the same temp-body form the align path produces, + // so the filters below are the identical code rather than a second + // implementation that can drift from it. + let mut body = tempfile::NamedTempFile::new().map_err(|e| Error::io(e, raw_dir))?; + let mut totals: HashMap = HashMap::default(); + let mut n_features = 0usize; + { + let mut out = std::io::BufWriter::new(body.as_file_mut()); + let reader = open_maybe_gz(&matrix_path)?; + let mut header_seen = false; + for line in reader.lines() { + let line = line.map_err(|e| Error::io(e, &matrix_path))?; + if line.starts_with('%') { + continue; + } + let mut fields = line.split_whitespace(); + let (Some(first), Some(second), Some(third)) = + (fields.next(), fields.next(), fields.next()) + else { + continue; + }; + if !header_seen { + // ` ` + n_features = first.parse().unwrap_or(0); + header_seen = true; + continue; + } + let (gene, cb, count) = ( + first.parse::().unwrap_or(0), + second.parse::().unwrap_or(0), + third.parse::().unwrap_or(0.0), + ); + if gene == 0 || cb == 0 { + continue; + } + // Counts can be real-valued in a multimapper matrix; the cell + // filters work on UMI totals, so round rather than refuse. + let count = count.round().max(0.0) as u64; + let e = totals.entry(cb - 1).or_insert((0, 0)); + e.0 += count; + e.1 += 1; + writeln!(out, "{gene} {cb} {count}").map_err(|e| Error::io(e, raw_dir))?; + } + out.flush().map_err(|e| Error::io(e, raw_dir))?; + } + if n_features == 0 { + return Err(Error::Parameter(format!( + "{}: no MatrixMarket header, so the matrix shape is unknown", + matrix_path.display() + )) + .into()); + } + + let mut cells: Vec = totals + .iter() + .map(|(&cb, &(n_umis, n_genes))| CellStat { + cb, + n_reads: n_umis, + n_umis, + n_genes, + }) + .collect(); + cells.sort_unstable_by_key(|c| c.cb); + log::info!( + "soloCellFiltering: {} barcodes with counts, {n_features} features", + cells.len() + ); + + let called = if params + .solo_cell_filter + .first() + .is_some_and(|m| m == "EmptyDrops_CR") + { + Some(emptydrops_called( + &cells, + &body, + n_features, + ¶ms.solo_cell_filter, + )?) + } else { + called_cells(&cells, ¶ms.solo_cell_filter) + }; + let Some(cbs) = called.filter(|c| !c.is_empty()) else { + log::warn!("soloCellFiltering: no cells called; writing nothing"); + return Ok(()); + }; + + let out_dir = Path::new(out_prefix); + let (dir, prefix): (&Path, &str) = if out_prefix.ends_with('/') { + (out_dir, "") + } else { + ( + out_dir.parent().unwrap_or_else(|| Path::new(".")), + out_dir.file_name().and_then(|s| s.to_str()).unwrap_or(""), + ) + }; + std::fs::create_dir_all(dir).map_err(|e| Error::io(e, dir))?; + + let remap: HashMap = cbs + .iter() + .enumerate() + .map(|(i, &cb)| (cb, i as u32 + 1)) + .collect(); + + let bc_path = dir.join(format!("{prefix}barcodes.tsv")); + let mut bc = String::new(); + for &cb in &cbs { + let Some(name) = barcodes.get(cb as usize) else { + continue; + }; + bc.push_str(name); + bc.push('\n'); + } + std::fs::write(&bc_path, bc).map_err(|e| Error::io(e, &bc_path))?; + + let feat_path = dir.join(format!("{prefix}features.tsv")); + std::fs::copy(&features_path, &feat_path).map_err(|e| Error::io(e, &feat_path))?; + + let mtx_path = dir.join(format!("{prefix}matrix.mtx")); + let nnz = finalize_matrix( + &body, + &mtx_path, + false, + n_features, + cbs.len(), + 0, + Some(&remap), + )?; + log::info!( + "soloCellFiltering: {} cells, {nnz} entries -> {}", + cbs.len(), + dir.display() + ); + Ok(()) +} + +/// First whitespace-separated column of a (possibly gzipped) file. +fn read_first_column(path: &Path) -> Result, Error> { + let reader = open_maybe_gz(path)?; + let mut out = Vec::new(); + for line in reader.lines() { + let line = line.map_err(|e| Error::io(e, path))?; + if let Some(first) = line.split_whitespace().next() { + out.push(first.to_string()); + } + } + Ok(out) +} + +/// Open a file, transparently decompressing `.gz`. +fn open_maybe_gz(path: &Path) -> Result, Error> { + let file = std::fs::File::open(path).map_err(|e| Error::io(e, path))?; + if path + .extension() + .is_some_and(|e| e.eq_ignore_ascii_case("gz")) + { + Ok(Box::new(BufReader::new(flate2::read::GzDecoder::new(file)))) + } else { + Ok(Box::new(BufReader::new(file))) + } +} + #[cfg(test)] mod tests { use super::*; diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index 2db7d39..62ee32a 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -2024,3 +2024,93 @@ fn test_wasp_samtag() { "all 10 unique reads overlapping the het SNV should pass WASP (vW:i:1)" ); } + +// --------------------------------------------------------------------------- +// --runMode soloCellFiltering +// --------------------------------------------------------------------------- + +/// Cell calling is a decision about a matrix, not about reads: the mode takes +/// an existing raw matrix and writes the called subset, with no genome and no +/// FASTQ involved. +#[test] +fn test_run_mode_solo_cell_filtering_calls_cells_from_a_raw_matrix() { + let tmpdir = TempDir::new().unwrap(); + let raw = tmpdir.path().join("raw"); + fs::create_dir_all(&raw).unwrap(); + + // 20 barcodes: five real cells with 1000 UMIs each, fifteen with 2. + let n_features = 50usize; + let n_cells = 20usize; + let mut entries: Vec<(usize, usize, u64)> = Vec::new(); + for cb in 1..=n_cells { + let per_gene = if cb <= 5 { 100 } else { 1 }; + let n_genes = if cb <= 5 { 10 } else { 2 }; + for gene in 1..=n_genes { + entries.push((gene, cb, per_gene)); + } + } + { + let mut f = fs::File::create(raw.join("matrix.mtx")).unwrap(); + writeln!(f, "%%MatrixMarket matrix coordinate integer general").unwrap(); + writeln!(f, "%").unwrap(); + writeln!(f, "{} {} {}", n_features, n_cells, entries.len()).unwrap(); + for (g, c, v) in &entries { + writeln!(f, "{g} {c} {v}").unwrap(); + } + } + { + let mut f = fs::File::create(raw.join("barcodes.tsv")).unwrap(); + for cb in 0..n_cells { + writeln!(f, "{:016b}", cb).unwrap(); + } + } + { + let mut f = fs::File::create(raw.join("features.tsv")).unwrap(); + for g in 0..n_features { + writeln!(f, "gene{g}\tGENE{g}\tGene Expression").unwrap(); + } + } + + let out = tmpdir.path().join("filtered/"); + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "soloCellFiltering", + raw.to_str().unwrap(), + out.to_str().unwrap(), + "--soloCellFilter", + "TopCells", + "5", + ]) + .assert() + .success(); + + let barcodes = fs::read_to_string(out.join("barcodes.tsv")).unwrap(); + assert_eq!( + barcodes.lines().count(), + 5, + "the five deep barcodes are the called cells" + ); + + let matrix = fs::read_to_string(out.join("matrix.mtx")).unwrap(); + let header = matrix.lines().nth(2).unwrap(); + let fields: Vec<&str> = header.split_whitespace().collect(); + assert_eq!(fields[0], "50", "features are carried through"); + assert_eq!(fields[1], "5", "columns are the called cells"); + assert_eq!(fields[2], "50", "10 genes × 5 cells"); + + assert!( + out.join("features.tsv").exists(), + "the feature list travels with the matrix" + ); +} + +/// The mode needs both paths; asking for it without them is refused rather +/// than run against a guess. +#[test] +fn test_run_mode_solo_cell_filtering_requires_its_paths() { + cargo_bin_cmd!("rustar-aligner") + .args(["--runMode", "soloCellFiltering"]) + .assert() + .failure(); +}