Fixes bug in assigning MAPQ 60 to reads that are not proper pairs.

CHANGES.md

@@ -19,6 +19,8 @@
   select it. This is equivalent to `-k 16 -s 12 -l 2 -u 2 -m 100`.
   We found these settings to increase accuracy on error-prone reads at the cost
   of runtime.
+* #596: Fix not properly paired reads getting MAPQ 60 even if there were
+  alternative equally good mapping locations.
 
 ## v0.17.0 (2025-12-18)
 

src/mapper.rs

@@ -717,6 +717,32 @@ pub fn align_paired_end_read(
                 alignment_pairs.swap(0, random_index);
             }
 
+            // Compute per-read MAPQ.
+            let mapqs = {
+                // When a single non-proper pair survives (the a_indv_max fallback used
+                // when the two reads map to different chromosomes), the joint
+                // alignment_quality would return 60 even for reads with
+                // multiple equally-good alternative chains. Use chain-based quality
+                // instead, exactly as the Proper fast path already does.
+                if alignment_pairs.len() == 1
+                    && is_proper_pair_opt(
+                        alignment_pairs[0].alignment1.as_ref(),
+                        alignment_pairs[0].alignment2.as_ref(),
+                        insert_size_distribution.mu,
+                        insert_size_distribution.sigma,
+                    ) == PairStatus::NotProper
+                {
+                    [
+                        mapping_quality(&nams_pair[0]),
+                        mapping_quality(&nams_pair[1]),
+                    ]
+                } else {
+                    let joint = alignment_quality(&alignment_pairs, |ap| ap.score as f32);
+
+                    [joint, joint]
+                }
+            };
+
             let secondary_dropoff = 2 * aligner.scores.mismatch + aligner.scores.gap_open;
             sam_records.extend(aligned_pairs_to_sam(
                 &alignment_pairs,
@@ -729,6 +755,7 @@ pub fn align_paired_end_read(
                 insert_size_distribution.mu,
                 insert_size_distribution.sigma,
                 &details,
+                mapqs,
             ));
         }
     }
@@ -1352,14 +1379,14 @@ fn aligned_pairs_to_sam(
     mu: f32,
     sigma: f32,
     details: &[Details; 2],
+    mapqs: [u8; 2],
 ) -> Vec<SamRecord> {
     let mut records = vec![];
     if high_scores.is_empty() {
         records.extend(sam_output.make_unmapped_pair([record1, record2], details));
         return records;
     }
 
-    let mapq = alignment_quality(high_scores, |ap| ap.score as f32);
     let best_aln_pair = &high_scores[0];
 
     if max_secondary == 0 {
@@ -1371,7 +1398,7 @@ fn aligned_pairs_to_sam(
             [alignment1.as_ref(), alignment2.as_ref()],
             references,
             [record1, record2],
-            [mapq, mapq],
+            mapqs,
             details,
             AlignmentStatus::Primary,
             pair_status,
@@ -1386,16 +1413,16 @@ fn aligned_pairs_to_sam(
             if s_max - s_score < secondary_dropoff {
                 let pair_status =
                     is_proper_pair_opt(alignment1.as_ref(), alignment2.as_ref(), mu, sigma);
-                let mapq = if alignment_status == AlignmentStatus::Primary {
-                    mapq
+                let effective_mapqs = if alignment_status == AlignmentStatus::Primary {
+                    mapqs
                 } else {
-                    0
+                    [0, 0]
                 };
                 records.extend(sam_output.make_paired_records(
                     [alignment1.as_ref(), alignment2.as_ref()],
                     references,
                     [record1, record2],
-                    [mapq, mapq],
+                    effective_mapqs,
                     details,
                     alignment_status,
                     pair_status,