Batch PRs #166, #170, #171, #157, #159

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/<feature>/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 <raw dir> <output prefix>`** 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.
+

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.

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();

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

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);
         }

src/align/transcript.rs

@@ -29,8 +29,6 @@ pub struct Transcript {
     pub junction_motifs: Vec<crate::align::score::SpliceMotif>,
     /// Whether each junction is annotated in the GTF (for jM +20 offset)
     pub junction_annotated: Vec<bool>,
-    /// Original read sequence
-    pub read_seq: Vec<u8>,
 }
 
 /// 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();

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],
         }
     }
 

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";