diff --git a/CHANGELOG.md b/CHANGELOG.md index b6e4893..9a63f97 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -97,8 +97,31 @@ Sections commonly used: Features, Bug fixes, Other changes. glibc's malloc and per-thread heaps that return whole segments to the OS when abandoned, so allocator cache size stays bounded. +- **`--soloCellReadStats CB`** writes `Solo.out//CellReads.stats`: + one row per cell barcode with fifteen counters describing what + happened to its reads — barcode match quality, unique or multi + genomic mapping, feature assignment, exonic/intronic and their + antisense counterparts, mitochondrial, and whether the read reached + the matrix — plus the per-cell UMI and gene totals. Reads whose + barcode never resolved are summed into a `CBnotInPasslist` row rather + than dropped, so the columns account for the whole input. + `--genomeChrSetMitochondrial` names the chromosomes behind the `mito` + column. +- **`--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 +- `--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 @@ -131,3 +154,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. + diff --git a/DIVERGENCE.md b/DIVERGENCE.md index bd957ed..f1d9f98 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -65,6 +65,20 @@ On the 10k yeast PE benchmark, 4 reads differ in alignment score (AS) because ST --- +### 3.2 `CellReads.stats` row order + +**What STAR does.** `--soloCellReadStats CB` emits its rows by iterating a libc++ `std::unordered_map`, so the order is a hash-table walk rather than a sort. At the map sizes this produces, libc++ chains new entries at the head of their bucket and walks buckets in order, which comes out as the reverse of each barcode's first appearance in read order. + +**What rustar-aligner does.** Emits that same reverse-first-appearance order, including across threads: the per-read accumulator merges in read order, so a threaded run writes the same file as a serial one. + +**Why.** Reproducing the order where it is reproducible costs nothing and keeps a byte-comparison against STAR usable on the sizes where it can work at all. + +**Impact.** Past libc++'s load factor the map rehashes, and the order then depends on the bucket count, which depends on how many distinct barcodes were seen; beyond that size the order diverges. The **values never do** — only which line they appear on. Reading the file by barcode rather than by position is unaffected either way. + +**Source.** `src/solo/cell_reads.rs`, locked by `rows_are_emitted_in_reverse_first_appearance_order` and `merging_partials_preserves_order_and_sums`. STAR: `SoloFeature_statsOutput.cpp`. + +--- + ## 4. Implementation divergences (no intended output difference) These differ in *how* a result is produced, not *what* is produced. They are documented so a reviewer chasing a discrepancy knows the mechanism differs by design. diff --git a/src/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/lib.rs b/src/lib.rs index 6fce173..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** @@ -68,17 +68,20 @@ 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(); - } - - match params.run_mode { + // + // 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), 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 0536a85..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())] @@ -1113,6 +1119,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, @@ -1162,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)) @@ -1353,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", @@ -1371,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 { @@ -1408,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( @@ -1431,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, @@ -1541,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() { @@ -1552,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() { @@ -1710,6 +1758,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(), @@ -1884,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")); @@ -1978,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, @@ -2018,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/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/cell_reads.rs b/src/solo/cell_reads.rs new file mode 100644 index 0000000..5d140a5 --- /dev/null +++ b/src/solo/cell_reads.rs @@ -0,0 +1,288 @@ +//! `--soloCellReadStats CB`: the per-cell-barcode read summary STARsolo writes +//! as `Solo.out//CellReads.stats`. +//! +//! One row per cell barcode, fifteen counters describing what happened to the +//! reads carrying it — how the barcode matched, whether the read mapped +//! uniquely, whether it landed on a feature, where in the gene, and whether it +//! reached the matrix — plus the per-cell UMI and gene totals. The reads whose +//! barcode never resolved are not dropped; they are summed into a single +//! `CBnotInPasslist` row, so the columns account for every read rather than +//! only the ones that succeeded. +//! +//! # D24: row order +//! +//! STAR iterates a libc++ `std::unordered_map` to emit these rows, so the order +//! is a hash-table walk, not a sort. For the small maps this produces, libc++ +//! chains new entries at the head of their bucket and walks buckets in order, +//! which comes out as the reverse of first appearance in read order. That is +//! what this reproduces. +//! +//! It is not reproducible in general: past the load factor libc++ rehashes, and +//! the order after a rehash depends on the bucket count, which depends on how +//! many distinct barcodes were seen. At that size the order diverges. The +//! **values never do** — only which line they appear on. A consumer that reads +//! this file by barcode rather than by position is unaffected either way, and +//! sorting by barcode is the only stable thing to do with it. +//! +//! Reads are folded in under a mutex held by `SoloContext`, so the order is the +//! order reads were processed in regardless of thread count. + +use std::collections::{BTreeMap, HashSet}; +use std::fmt::Write as _; + +/// What happened to one read, as the fourteen optional flags STAR tracks. +/// +/// `cbMatch` is not here: every read that reaches the accumulator matched a +/// barcode well enough to be attributed somewhere, so it is always set. +/// +/// Fourteen bools rather than a bitfield because they are written once per read +/// and read once per fold, and the names are the column names of the file. +#[allow(clippy::struct_excessive_bools)] +#[derive(Debug, Clone, Copy, Default)] +pub struct CellReadFlag { + /// The barcode was an exact whitelist hit. + pub cb_perfect: bool, + /// Corrected via a single one-mismatch neighbour. + pub cb_mm_unique: bool, + /// Corrected via several one-mismatch neighbours, resolved by the posterior. + pub cb_mm_multiple: bool, + /// The read mapped to exactly one genomic locus. + pub genome_u: bool, + /// The read mapped to several genomic loci. + pub genome_m: bool, + /// It fell on exactly one feature (gene). + pub feature_u: bool, + /// It fell on several features. + pub feature_m: bool, + /// Exonic, on the annotated strand. + pub exonic: bool, + /// Intronic, on the annotated strand. + pub intronic: bool, + /// Exonic, antisense to the annotation. + pub exonic_as: bool, + /// Intronic, antisense to the annotation. + pub intronic_as: bool, + /// On a chromosome named by `--genomeChrSetMitochondrial`. + pub mito: bool, + /// Counted into the unique-gene matrix. + pub counted_u: bool, + /// Counted through the multi-gene distribution. + pub counted_m: bool, +} + +/// The accumulator behind `CellReads.stats`. +#[derive(Debug, Default, Clone)] +pub struct CellReadStats { + /// Whitelist index to its fifteen counters. + cells: BTreeMap, + /// The single bucket for reads whose barcode did not resolve. + no_cb: [u64; 15], + /// Whitelist indices in first-appearance order; emitted reversed (D24). + order: Vec, + seen: HashSet, +} + +impl CellReadStats { + pub fn new() -> Self { + Self::default() + } + + fn fold(v: &mut [u64; 15], f: &CellReadFlag) { + v[0] += 1; // cbMatch: every read that got this far + for (i, set) in [ + f.cb_perfect, + f.cb_mm_unique, + f.cb_mm_multiple, + f.genome_u, + f.genome_m, + f.feature_u, + f.feature_m, + f.exonic, + f.intronic, + f.exonic_as, + f.intronic_as, + f.mito, + f.counted_u, + f.counted_m, + ] + .into_iter() + .enumerate() + { + if set { + v[i + 1] += 1; + } + } + } + + /// Record a read that resolved to whitelist cell `cb`. + pub fn add_cell(&mut self, cb: u32, flag: &CellReadFlag) { + if self.seen.insert(cb) { + self.order.push(cb); + } + Self::fold(self.cells.entry(cb).or_insert([0; 15]), flag); + } + + /// Record a read whose barcode did not resolve, or whose UMI was rejected. + pub fn add_no_cb(&mut self, flag: &CellReadFlag) { + Self::fold(&mut self.no_cb, flag); + } + + /// Render the file. `umi_gene` gives each cell its final + /// `(nUMIunique, nGenesUnique)`; a cell absent from it prints zeros. + /// `barcode_of` renders a whitelist index as its barcode string. + pub fn render( + &self, + barcode_of: impl Fn(u32) -> String, + umi_gene: &BTreeMap, + ) -> String { + let mut s = String::from( + "CB\tcbMatch\tcbPerfect\tcbMMunique\tcbMMmultiple\tgenomeU\tgenomeM\tfeatureU\t\ + featureM\texonic\tintronic\texonicAS\tintronicAS\tmito\tcountedU\tcountedM\t\ + nUMIunique\tnGenesUnique\tnUMImulti\tnGenesMulti\n", + ); + s.push_str("CBnotInPasslist"); + for v in &self.no_cb { + s.push('\t'); + s.push_str(&v.to_string()); + } + s.push_str("\t0\t0\t0\t0\n"); + // Reverse first-appearance order — see D24 in the module docs. + for &cb in self.order.iter().rev() { + s.push_str(&barcode_of(cb)); + for v in &self.cells[&cb] { + s.push('\t'); + s.push_str(&v.to_string()); + } + let (n_umi, n_gene) = umi_gene.get(&cb).copied().unwrap_or((0, 0)); + // nUMImulti and nGenesMulti stay zero: multi-gene UMIs are not + // collapsed into per-cell totals here. + let _ = writeln!(s, "\t{n_umi}\t{n_gene}\t0\t0"); + } + s + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn flag_perfect_counted() -> CellReadFlag { + CellReadFlag { + cb_perfect: true, + genome_u: true, + feature_u: true, + exonic: true, + counted_u: true, + ..Default::default() + } + } + + fn header_and_rows(text: &str) -> (Vec<&str>, Vec<&str>) { + let mut lines = text.lines(); + let header: Vec<&str> = lines.next().unwrap().split('\t').collect(); + (header, lines.collect()) + } + + /// Every column has a name, and every row has a value under each of them. + #[test] + fn each_row_fills_the_header() { + let mut st = CellReadStats::new(); + st.add_cell(0, &flag_perfect_counted()); + st.add_no_cb(&CellReadFlag::default()); + let out = st.render(|i| format!("CB{i}"), &BTreeMap::new()); + let (header, rows) = header_and_rows(&out); + assert_eq!(header.len(), 20); + assert_eq!(rows.len(), 2, "the passlist-miss row plus one cell"); + for row in rows { + assert_eq!(row.split('\t').count(), header.len()); + } + } + + /// A read is counted once under `cbMatch` and once under each flag it sets, + /// so the flag columns can exceed neither `cbMatch` nor each other's logic. + #[test] + fn flags_accumulate_per_read() { + let mut st = CellReadStats::new(); + for _ in 0..3 { + st.add_cell(7, &flag_perfect_counted()); + } + st.add_cell( + 7, + &CellReadFlag { + cb_mm_unique: true, + genome_m: true, + ..Default::default() + }, + ); + let out = st.render(|i| format!("CB{i}"), &BTreeMap::new()); + let row: Vec<&str> = out.lines().nth(2).unwrap().split('\t').collect(); + assert_eq!(row[0], "CB7"); + assert_eq!(row[1], "4", "cbMatch counts every read"); + assert_eq!(row[2], "3", "cbPerfect"); + assert_eq!(row[3], "1", "cbMMunique"); + assert_eq!(row[5], "3", "genomeU"); + assert_eq!(row[6], "1", "genomeM"); + assert_eq!(row[14], "3", "countedU"); + } + + /// Reads whose barcode never resolved are summed rather than dropped, so + /// the file accounts for the whole input. + #[test] + fn unresolved_reads_land_in_the_passlist_miss_row() { + let mut st = CellReadStats::new(); + st.add_cell(0, &flag_perfect_counted()); + for _ in 0..5 { + st.add_no_cb(&CellReadFlag { + genome_u: true, + ..Default::default() + }); + } + let out = st.render(|i| format!("CB{i}"), &BTreeMap::new()); + let row: Vec<&str> = out.lines().nth(1).unwrap().split('\t').collect(); + assert_eq!(row[0], "CBnotInPasslist"); + assert_eq!(row[1], "5"); + assert_eq!(row[5], "5", "genomeU"); + assert_eq!(&row[16..], ["0", "0", "0", "0"], "no UMI columns for it"); + } + + /// D24: rows come out in reverse first-appearance order, which is what + /// STAR's libc++ hash-map walk produces at these sizes. + #[test] + fn rows_are_emitted_in_reverse_first_appearance_order() { + let mut st = CellReadStats::new(); + for cb in [4u32, 1, 9] { + st.add_cell(cb, &flag_perfect_counted()); + } + st.add_cell(1, &flag_perfect_counted()); // seen again: order unchanged + let out = st.render(|i| format!("CB{i}"), &BTreeMap::new()); + let cbs: Vec<&str> = out + .lines() + .skip(2) + .map(|l| l.split('\t').next().unwrap()) + .collect(); + assert_eq!(cbs, ["CB9", "CB1", "CB4"]); + } + + /// The UMI and gene totals come from the final matrix, not from the read + /// counters, so a cell missing from that map prints zeros rather than + /// inheriting a read count. + #[test] + fn umi_and_gene_totals_come_from_the_matrix() { + let mut st = CellReadStats::new(); + st.add_cell(2, &flag_perfect_counted()); + st.add_cell(3, &flag_perfect_counted()); + let mut umi_gene = BTreeMap::new(); + umi_gene.insert(2u32, (17u32, 5u32)); + let out = st.render(|i| format!("CB{i}"), &umi_gene); + let rows: Vec> = out + .lines() + .skip(2) + .map(|l| l.split('\t').collect()) + .collect(); + // Reverse order: CB3 first, then CB2. + assert_eq!(rows[0][0], "CB3"); + assert_eq!(&rows[0][16..18], ["0", "0"]); + assert_eq!(rows[1][0], "CB2"); + assert_eq!(&rows[1][16..18], ["17", "5"]); + } +} diff --git a/src/solo/count.rs b/src/solo/count.rs index 9af3edf..58f8524 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 @@ -1881,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/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/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 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![], } } 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 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(); +}