Merge origin/main (PRs #5 run-rng-seed, scverse#6 sam-rg-line) into dev

- shuffle_tied_prefix: deterministic per-read RNG for tied primary selection
- --runRNGseed param (default 777, matches STAR)
- --outSAMattrRGline: @rg header + RG:Z tags
- Various code simplifications and formatting

Author: Psy-Fer

Cargo.lock

@@ -49,6 +49,7 @@ flate2 = "1"
 rayon = "1"
 dashmap = "6"
 chrono = "0.4"
+rand = "0.8"
 
 [dev-dependencies]
 tempfile = "3"

Cargo.toml

@@ -10,6 +10,38 @@ use crate::error::Error;
 use crate::index::GenomeIndex;
 use crate::params::{IntronMotifFilter, IntronStrandFilter, Parameters};
 use crate::stats::UnmappedReason;
+use rand::{SeedableRng, rngs::StdRng, seq::SliceRandom};
+use std::hash::{DefaultHasher, Hash, Hasher};
+
+/// Derive a deterministic per-read RNG seed from `run_rng_seed` + the read name.
+///
+/// STAR seeds `std::mt19937` once per chunk/thread (`runRNGseed*(iChunk+1)`),
+/// then advances the state sequentially per read. ruSTAR parallelises per-read
+/// via rayon, so we instead fold the read name into the seed — this keeps tie
+/// breaks reproducible regardless of thread count while still honoring the
+/// user's `--runRNGseed` value.
+fn per_read_seed(run_rng_seed: u64, read_name: &str) -> u64 {
+    let mut hasher = DefaultHasher::new();
+    read_name.hash(&mut hasher);
+    run_rng_seed.wrapping_mul(hasher.finish().wrapping_add(1))
+}
+
+/// Shuffle the prefix of `items` whose `score_fn` equals the first element's score.
+///
+/// Mirrors STAR's `ReadAlign_multMapSelect` / `funPrimaryAlignMark`: best-scoring
+/// alignments are randomized so primary selection (index 0) is not biased by
+/// upstream sort order. Non-tied elements are left alone.
+fn shuffle_tied_prefix<T>(items: &mut [T], score_fn: impl Fn(&T) -> i32, seed: u64) {
+    let Some(first) = items.first() else {
+        return;
+    };
+    let best = score_fn(first);
+    let tied = items.iter().take_while(|t| score_fn(t) == best).count();
+    if tied < 2 {
+        return;
+    }
+    items[..tied].shuffle(&mut StdRng::seed_from_u64(seed));
+}
 
 /// Result of aligning a single read: (transcripts, chimeric_alignments, n_for_mapq, unmapped_reason)
 pub type AlignReadResult = (
@@ -302,6 +334,13 @@ pub fn align_read(
             .then_with(|| a.is_reverse.cmp(&b.is_reverse))
     });
 
+    // Randomize primary among best-scoring ties (ReadAlign_multMapSelect.cpp:71-79).
+    shuffle_tied_prefix(
+        &mut transcripts,
+        |t| t.score,
+        per_read_seed(params.run_rng_seed, read_name),
+    );
+
     // Score-range filter: keep only alignments within outFilterMultimapScoreRange of the best.
     // (STAR's multMapSelect step — must run before quality filters.)
     if !transcripts.is_empty() {
@@ -902,6 +941,13 @@ pub fn align_paired_read(
             })
     });
 
+    // Randomize primary among best-scoring pairs (STAR's funPrimaryAlignMark).
+    shuffle_tied_prefix(
+        &mut joint_pairs,
+        |pa| pa.combined_wt_score,
+        per_read_seed(params.run_rng_seed, read_name),
+    );
+
     // Step 4: quality filter (mappedFilter).
     filter_paired_transcripts(&mut joint_pairs, params);
 
@@ -1122,9 +1168,13 @@ fn calculate_insert_size(mate1_trans: &Transcript, mate2_trans: &Transcript) ->
 /// Filter paired transcripts by quality thresholds.
 /// STAR's mappedFilter applies ALL quality checks to trBest (the highest-scoring transcript).
 fn filter_paired_transcripts(paired_alns: &mut Vec<PairedAlignment>, params: &Parameters) {
-    let best_pa = paired_alns
-        .iter()
-        .max_by_key(|pa| pa.combined_wt_score);
+    // STAR's mappedFilter (ReadAlign_mappedFilter.cpp) applies ALL quality thresholds to trBest
+    // (the highest-scoring transcript), NOT to each individual transcript. If trBest passes,
+    // all transcripts in the score window are included (they affect NH/MAPQ). If trBest fails,
+    // the read is unmapped.
+    //
+    // Step 1: find the best pair and check quality thresholds on it.
+    let best_pa = paired_alns.iter().max_by_key(|pa| pa.combined_wt_score);
     if let Some(best) = best_pa {
         let mate1_len = (best.mate1_region.1 - best.mate1_region.0) as f64;
         let mate2_len = (best.mate2_region.1 - best.mate2_region.0) as f64;
@@ -1133,17 +1183,15 @@ fn filter_paired_transcripts(paired_alns: &mut Vec<PairedAlignment>, params: &Pa
         let combined_score = best.combined_wt_score;
 
         if combined_score < params.out_filter_score_min
-            || combined_score
-                < (params.out_filter_score_min_over_lread * combined_lread_m1) as i32
+            || combined_score < (params.out_filter_score_min_over_lread * combined_lread_m1) as i32
         {
             paired_alns.clear();
             return;
         }
 
         let combined_match = best.combined_n_match;
         if combined_match < params.out_filter_match_nmin
-            || combined_match
-                < (params.out_filter_match_nmin_over_lread * combined_lread_m1) as u32
+            || combined_match < (params.out_filter_match_nmin_over_lread * combined_lread_m1) as u32
         {
             paired_alns.clear();
             return;
@@ -1814,4 +1862,55 @@ mod tests {
             assert!(mate1_is_mapped);
         }
     }
+
+    #[test]
+    fn shuffle_tied_prefix_is_deterministic() {
+        // Same seed + same input → same permutation on reruns.
+        let items: Vec<(i32, u32)> = (0..8).map(|i| (100, i)).collect();
+        let mut a = items.clone();
+        let mut b = items.clone();
+        shuffle_tied_prefix(&mut a, |t| t.0, 12345);
+        shuffle_tied_prefix(&mut b, |t| t.0, 12345);
+        assert_eq!(a, b);
+    }
+
+    #[test]
+    fn shuffle_tied_prefix_respects_ties() {
+        // Only the top-score prefix gets shuffled; lower-scored tail is left alone.
+        let mut items = vec![(100, 0u32), (100, 1), (100, 2), (50, 3), (40, 4)];
+        shuffle_tied_prefix(&mut items, |t| t.0, 777);
+        // Last two elements (non-tied) stay in place.
+        assert_eq!(items[3], (50, 3));
+        assert_eq!(items[4], (40, 4));
+        // Tied prefix contains the original three items in some order.
+        let mut top: Vec<u32> = items[..3].iter().map(|t| t.1).collect();
+        top.sort();
+        assert_eq!(top, vec![0, 1, 2]);
+    }
+
+    #[test]
+    fn shuffle_tied_prefix_different_seeds_can_diverge() {
+        // Probabilistic: for a tied set of 8, at least two seeds should disagree
+        // on the chosen primary. (Exhaustive over a small seed range is fine.)
+        let base: Vec<(i32, u32)> = (0..8).map(|i| (100, i)).collect();
+        let mut firsts = std::collections::HashSet::new();
+        for seed in 0..32u64 {
+            let mut v = base.clone();
+            shuffle_tied_prefix(&mut v, |t| t.0, seed);
+            firsts.insert(v[0].1);
+        }
+        assert!(
+            firsts.len() >= 2,
+            "expected different seeds to pick different primaries, got {:?}",
+            firsts
+        );
+    }
+
+    #[test]
+    fn shuffle_tied_prefix_noop_when_no_ties() {
+        let mut items = vec![(100, 0u32), (90, 1), (80, 2)];
+        let before = items.clone();
+        shuffle_tied_prefix(&mut items, |t| t.0, 42);
+        assert_eq!(items, before);
+    }
 }

src/align/read_align.rs

@@ -2567,7 +2567,6 @@ pub(crate) fn stitch_seeds_with_jdb_debug(
     Ok(transcripts)
 }
 
-
 /// Shared core: preprocessing + recursive stitcher, returns working transcripts + context.
 #[allow(clippy::too_many_arguments)]
 pub(crate) fn stitch_seeds_core(
@@ -2914,8 +2913,6 @@ pub(crate) fn stitch_seeds_core(
     ))
 }
 
-
-
 #[cfg(test)]
 mod tests {
     use super::*;

src/align/stitch.rs

@@ -125,9 +125,10 @@ mod tests {
         let read_qual = vec![30, 30, 30, 30];
 
         // Build record using SAM builder
-        let record =
-            crate::io::sam::SamWriter::build_unmapped_record(read_name, &read_seq, &read_qual)
-                .unwrap();
+        let record = crate::io::sam::SamWriter::build_unmapped_record(
+            read_name, &read_seq, &read_qual, None,
+        )
+        .unwrap();
 
         let result = writer.write_batch(&[record]);
         assert!(result.is_ok(), "Writing unmapped read should succeed");
@@ -204,12 +205,14 @@ mod tests {
                 "read1",
                 &[0, 1, 2, 3],
                 &[30, 30, 30, 30],
+                None,
             )
             .unwrap(),
             crate::io::sam::SamWriter::build_unmapped_record(
                 "read2",
                 &[0, 1, 2, 3],
                 &[30, 30, 30, 30],
+                None,
             )
             .unwrap(),
         ];