more MAPQ fixes

Author: Psy-Fer

CLAUDE.md

@@ -32,7 +32,7 @@ Always run `cargo clippy`, `cargo fmt --check`, and `cargo test` before consider
 
 ## Current Status
 
-**268 tests passing, 0 clippy warnings.** SE: 99.7% position agreement (adjusted), 99.9% CIGAR (pos-agreeing), 2.2% splice rate (STAR: 2.2%), 66 shared junctions, **100.0% MAPQ agreement, MAPQ inflation: 0, deflation: 0**. 127 position disagreements (ALL verified as genuine ties). 1 CIGAR-only disagree (ERR12389696.13573895, insertion placement, seed-level tie). **1 truly actionable SE issue remains** (CIGAR insertion placement only). PE: **8390/8390 both-mapped (0 gap, exact STAR match)**, 0 half-mapped, **6 MAPQ inflations** (was 24), **2 MAPQ deflations** (new, spurious subset WT), **98.879% PE faithfulness** (Phase 16.46: per-position PE dedup removal). See [ROADMAP.md](ROADMAP.md) for detailed phase tracking and [docs/](docs/) for per-phase notes.
+**268 tests passing, 0 clippy warnings.** SE: 8796/8926 compare_sam.py (98.5%), 2.2% splice rate (STAR: 2.2%), 66 shared junctions, **100.0% MAPQ agreement, MAPQ inflation: 0, deflation: 0**. 127 position disagreements (ALL verified as genuine ties). 1 CIGAR-only disagree (ERR12389696.13573895, insertion placement, seed-level tie). **0 STAR-only / 0 ruSTAR-only SE reads**. PE: **8390/8390 both-mapped (0 gap, exact STAR match)**, 0 half-mapped, **4 MAPQ inflations** (rDNA/repeat secondary loci), **98.903% PE faithfulness** (Phase 16.50). See [ROADMAP.md](ROADMAP.md) for detailed phase tracking and [docs/](docs/) for per-phase notes.
 
 ## Source Layout
 

ROADMAP.md

@@ -51,7 +51,7 @@ Paired-end (Phase 8) builds on threaded infrastructure. GTF/junctions (Phase 7)
 | 12 | Chimeric Detection | ✅ | 170 | SE chimeric, Chimeric.out.junction |
 | [13](docs/phase13_accuracy.md) | Performance + Accuracy | ✅ | 205 | 94.5% pos, 97.8% CIGAR, 2.1% splice |
 | [15](docs/phase15_sam_tags.md) | SAM Tags + PE Fix | ✅ | 235 | NH/HI/AS/NM/nM/XS/jM/jI/MD, PE fix |
-| [16](docs/phase16_algorithm.md) | Algorithm Parity | ✅* | 268 | SE: 99.7% pos, 2.2% splice, **MAPQ 100%**; PE: **8390/8390 (0 gap)**, 99.0% per-mate pos, 98.9% CIGAR, **6 MAPQ inflations**, 0 deflations; faithfulness: SE 98.57%, PE 98.89%, SJ 96.97% (Phase 16.48) |
+| [16](docs/phase16_algorithm.md) | Algorithm Parity | ✅* | 268 | SE: **8796/8926 (0 STAR-only)**, 2.2% splice, **MAPQ 100%**; PE: **8390/8390 (0 gap)**, 99.0% per-mate pos, 98.9% CIGAR, **4 MAPQ inflations**, 0 deflations; faithfulness: SE 98.5%+, PE 98.903%, SJ 96.97% (Phase 16.50) |
 | [17](docs/phase17_features.md) | Features + Polish | ✅* | 268 | Log.final.out, clippy cleanup, sorted BAM planned |
 | 14 | STARsolo | DEFERRED | — | Waiting for accuracy parity |
 

src/align/read_align.rs

@@ -564,7 +564,6 @@ fn build_combined_read(mate1: &[u8], mate2: &[u8]) -> Vec<u8> {
     v
 }
 
-
 /// Align paired-end reads using STAR's joint combined-read approach (Phase 16.11).
 ///
 /// # Algorithm (STAR-faithful)

src/align/seed.rs

@@ -390,8 +390,7 @@ fn find_seed_at_position(
     // matched_level directly as the MMP. This causes chains to advance by shorter
     // amounts, ensuring intermediate positions (missed if advancing by the true MMP)
     // are not skipped.
-    let (match_length, narrowed_start, narrowed_end) = if bounds_tight
-        && matched_level < sa_nbases
+    let (match_length, narrowed_start, narrowed_end) = if bounds_tight && matched_level < sa_nbases
     {
         // STAR short-circuit (maxMappableLength2strands.cpp):
         // "if (Lind < gSAindexNbases && iSA1noN && iSA2good) { maxL=Lind; }"

src/align/stitch.rs

@@ -81,7 +81,9 @@ fn score_region(
             break;
         }
         let read_base = read_seq[read_pos];
-        let Some(genome_base) = index.genome.get_base(genome_start + i as u64 + genome_offset)
+        let Some(genome_base) = index
+            .genome
+            .get_base(genome_start + i as u64 + genome_offset)
         else {
             break;
         };
@@ -580,8 +582,32 @@ pub fn cluster_seeds(
     // STAR resets nWA=0 after window creation, then re-assigns ALL seeds (including
     // anchors) through assignAlignToWindow. We match this by not pre-loading anchors
     // in Phase 2 and processing all seeds here through the same capacity logic.
-    let mut too_many_anchors = false; // STAR: MARKER_TOO_MANY_ANCHORS_PER_WINDOW
-    'outer: for (seed_idx, seed) in seeds.iter().enumerate() {
+    //
+    // STAR processes all pieces in a sorted pass (rStart asc, length desc for ties).
+    // This sorting ensures that when two seeds overlap on the same diagonal, the longer
+    // one enters the window first and blocks the shorter one via overlap detection —
+    // preventing window capacity overflow from many short copies on the same diagonal.
+    //
+    // We replicate this by: (1) collecting all candidates per window, (2) pre-deduping
+    // overlapping entries per diagonal (longest wins, shorter blocked before processing),
+    // then (3) processing survivors in original discovery order. This achieves the same
+    // overlap suppression without globally reordering seeds across windows.
+
+    // Per-window candidate entries collected before sorting.
+    struct WinCandidate {
+        seed_idx: usize,
+        sa_pos: u64,
+        forward_pos: u64,
+        length: usize,
+        n_loci: usize,
+        is_anchor: bool,
+        ps_rstart: usize, // positive-strand read start (sort key)
+    }
+
+    let mut win_candidates: Vec<Vec<WinCandidate>> =
+        (0..windows.len()).map(|_| Vec::new()).collect();
+
+    for (seed_idx, seed) in seeds.iter().enumerate() {
         let n_loci = seed.sa_end - seed.sa_start;
         if n_loci == 0 {
             continue;
@@ -608,6 +634,97 @@ pub fn cluster_seeds(
                 _ => continue,
             };
 
+            let ps_rstart = if windows[win_idx].is_reverse {
+                read_len - (length + seed.read_pos)
+            } else {
+                seed.read_pos
+            };
+
+            win_candidates[win_idx].push(WinCandidate {
+                seed_idx,
+                sa_pos,
+                forward_pos,
+                length,
+                n_loci,
+                is_anchor: is_anchor_seed,
+                ps_rstart,
+            });
+        }
+    }
+
+    // Pre-dedup overlapping diagonals (order-independent, length wins):
+    // Two candidate entries on the same diagonal with overlapping read ranges represent
+    // the same alignment position — keep only the longest. This is what STAR achieves
+    // via sorted-order overlap detection (longer seed enters first, shorter blocked).
+    // Pre-computing this dedup ensures correct results regardless of discovery order,
+    // fixing window capacity overflow without changing capacity-eviction behavior
+    // for other reads (seeds on unique diagonals are unaffected).
+    //
+    // After pre-dedup, Phase 4 processes survivors in original discovery order.
+    // The overlap detection in the main Phase 4 loop below is still present but will
+    // never find an overlap (since pre-dedup already resolved all of them). It serves
+    // as a safety net for any edge cases not covered by pre-dedup.
+    let win_n = windows.len();
+    let mut win_blocked: Vec<Vec<bool>> = (0..win_n)
+        .map(|i| vec![false; win_candidates[i].len()])
+        .collect();
+
+    for win_idx in 0..win_n {
+        let candidates = &win_candidates[win_idx];
+        if candidates.is_empty() {
+            continue;
+        }
+        // Sort indices by length descending so that the longest entry per diagonal
+        // is processed first (and blocks shorter overlapping entries on the same diagonal).
+        let mut by_len: Vec<usize> = (0..candidates.len()).collect();
+        by_len.sort_by(|&a, &b| candidates[b].length.cmp(&candidates[a].length));
+
+        // For each diagonal, track accepted [ps_rstart, ps_rend) ranges.
+        let mut diag_ranges: HashMap<i64, Vec<(usize, usize)>> = HashMap::new();
+        for &ci in &by_len {
+            let cand = &candidates[ci];
+            let diag = cand.forward_pos as i64 - cand.ps_rstart as i64;
+            let ps_rend = cand.ps_rstart + cand.length;
+
+            let blocked = diag_ranges.get(&diag).is_some_and(|ranges| {
+                ranges.iter().any(|&(rs, re)| {
+                    (cand.ps_rstart >= rs && cand.ps_rstart < re) || (ps_rend >= rs && ps_rend < re)
+                })
+            });
+
+            if blocked {
+                win_blocked[win_idx][ci] = true;
+            } else {
+                diag_ranges
+                    .entry(diag)
+                    .or_default()
+                    .push((cand.ps_rstart, ps_rend));
+            }
+        }
+    }
+
+    // Process each window's candidates in original discovery order,
+    // skipping pre-dedup-blocked entries.
+    let mut too_many_anchors = false; // STAR: MARKER_TOO_MANY_ANCHORS_PER_WINDOW
+    'outer: for win_idx in 0..win_n {
+        if !windows[win_idx].alive {
+            continue;
+        }
+        for (ci, cand) in win_candidates[win_idx].iter().enumerate() {
+            if win_blocked[win_idx][ci] {
+                continue; // blocked by a longer overlapping entry on same diagonal
+            }
+
+            let seed_idx = cand.seed_idx;
+            let seed = &seeds[seed_idx];
+            let length = cand.length;
+            let forward_pos = cand.forward_pos;
+            let sa_pos = cand.sa_pos;
+            let n_loci = cand.n_loci;
+            let is_anchor_seed = cand.is_anchor;
+            let new_ps_rstart = cand.ps_rstart;
+            let new_ps_rend = new_ps_rstart + length;
+
             let window = &mut windows[win_idx];
 
             // STAR's WALrec early rejection (assignAlignToWindow line 20):
@@ -617,51 +734,29 @@ pub fn cluster_seeds(
                 continue;
             }
 
-            // STAR's overlap detection (assignAlignToWindow lines 28-84):
-            // Check if new entry overlaps an existing entry on the same diagonal.
-            // Same diagonal means (genome_pos - read_pos) is equal. If overlap
-            // found, keep the longer entry. This deduplicates shorter seeds
-            // subsumed by longer seeds at the same alignment position.
-            //
-            // CRITICAL: STAR uses positive-strand read coordinates for overlap
-            // detection (stitchPieces.cpp line 151: aRstart = Lread - (aLength+aRstart)
-            // for reverse-strand). Without this, reverse-strand seeds covering the
-            // same genomic exon have different diagonals → no overlap → anchor
-            // proliferation → window capacity overflow → splice-target seed eviction.
-            let (new_ps_rstart, new_ps_rend) = if window.is_reverse {
-                let ps = read_len - (length + seed.read_pos);
-                (ps, ps + length)
-            } else {
-                (seed.read_pos, seed.read_pos + length)
-            };
+            // Safety-net overlap detection: after pre-dedup, no overlapping entries
+            // should remain. This check handles any edge cases and matches STAR's
+            // assignAlignToWindow overlap logic for correctness.
             let new_diag = forward_pos as i64 - new_ps_rstart as i64;
             let mut overlap_idx = None;
             for (i, wa) in window.alignments.iter().enumerate() {
-                let (wa_ps_rstart, wa_ps_rend) = if window.is_reverse {
-                    let ps = read_len - (wa.length + wa.read_pos);
-                    (ps, ps + wa.length)
+                let wa_ps_rstart = if window.is_reverse {
+                    read_len - (wa.length + wa.read_pos)
                 } else {
-                    (wa.read_pos, wa.read_pos + wa.length)
+                    wa.read_pos
                 };
+                let wa_ps_rend = wa_ps_rstart + wa.length;
                 let wa_diag = wa.genome_pos as i64 - wa_ps_rstart as i64;
-                if new_diag == wa_diag {
-                    // Same diagonal — check for overlapping read ranges
-                    // in positive-strand coordinates (matches STAR's exact
-                    // overlap condition from assignAlignToWindow)
-                    if (new_ps_rstart >= wa_ps_rstart && new_ps_rstart < wa_ps_rend)
-                        || (new_ps_rend >= wa_ps_rstart && new_ps_rend < wa_ps_rend)
-                    {
-                        overlap_idx = Some(i);
-                        break;
-                    }
+                if new_diag == wa_diag
+                    && ((new_ps_rstart >= wa_ps_rstart && new_ps_rstart < wa_ps_rend)
+                        || (new_ps_rend >= wa_ps_rstart && new_ps_rend < wa_ps_rend))
+                {
+                    overlap_idx = Some(i);
+                    break;
                 }
             }
             if let Some(oi) = overlap_idx {
                 if length > window.alignments[oi].length {
-                    // New entry is longer — remove old and re-insert at correct
-                    // sorted position (matches STAR assignAlignToWindow lines 26-62).
-                    // The new entry may have a different length → different ps_rstart
-                    // → different sort position.
                     window.alignments.remove(oi);
                     let insert_pos = window.alignments.partition_point(|wa| {
                         let wa_ps = if window.is_reverse {
@@ -685,7 +780,6 @@ pub fn cluster_seeds(
                         },
                     );
                 }
-                // Either way (replaced or not), don't add as new entry
                 continue;
             }