From 5b6783ae0152b3ce8885b2d1a43d438a77cd7aff Mon Sep 17 00:00:00 2001 From: Kristoffer Sahlin Date: Tue, 9 Jun 2026 22:45:34 +0200 Subject: [PATCH 1/2] Fix not-proper pairs getting MAPQ 60 Previously the individual best non-proper pair would get a mapq of 60, regardless if it had multiple best scoring chains or not. This commit fixes the bug. If the best pair is individual: check for each read if second best chain has the same score, if so, assign mapq 0, otherwise compute mapq from our formula. If there is no second chain the read gets 60. --- src/mapper.rs | 42 ++++++++++++++++++++++++++++++++++++------ 1 file changed, 36 insertions(+), 6 deletions(-) diff --git a/src/mapper.rs b/src/mapper.rs index 5f9994bc..d34dcd1b 100644 --- a/src/mapper.rs +++ b/src/mapper.rs @@ -717,6 +717,35 @@ pub fn align_paired_end_read( alignment_pairs.swap(0, random_index); } + // Compute per-read MAPQ. + // 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(len=1) would return 60 even for reads with hundreds + // of equally-good alternative NAMs. Use NAM-based per-read quality + // instead, exactly as the Proper fast path already does. + let mapqs = { + let joint = alignment_quality(&alignment_pairs, |ap| ap.score as f32); + if alignment_pairs.len() == 1 { + let p = &alignment_pairs[0]; + let proper = is_proper_pair_opt( + p.alignment1.as_ref(), + p.alignment2.as_ref(), + insert_size_distribution.mu, + insert_size_distribution.sigma, + ) == PairStatus::Proper; + if !proper { + [ + mapping_quality(&nams_pair[0]), + mapping_quality(&nams_pair[1]), + ] + } else { + [joint, joint] + } + } else { + [joint, joint] + } + }; + let secondary_dropoff = 2 * aligner.scores.mismatch + aligner.scores.gap_open; sam_records.extend(aligned_pairs_to_sam( &alignment_pairs, @@ -729,6 +758,7 @@ pub fn align_paired_end_read( insert_size_distribution.mu, insert_size_distribution.sigma, &details, + mapqs, )); } } @@ -1352,6 +1382,7 @@ fn aligned_pairs_to_sam( mu: f32, sigma: f32, details: &[Details; 2], + mapqs: [u8; 2], ) -> Vec { let mut records = vec![]; if high_scores.is_empty() { @@ -1359,7 +1390,6 @@ fn aligned_pairs_to_sam( 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 +1401,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 +1416,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, From e9cb27809ab3f9098d77152c3b7386f7c02a2267 Mon Sep 17 00:00:00 2001 From: Marcel Martin Date: Fri, 12 Jun 2026 11:23:49 +0200 Subject: [PATCH 2/2] Simplify slightly; changelog --- CHANGES.md | 2 ++ src/mapper.rs | 37 +++++++++++++++++-------------------- 2 files changed, 19 insertions(+), 20 deletions(-) diff --git a/CHANGES.md b/CHANGES.md index 1c270d7a..5f108a1c 100644 --- a/CHANGES.md +++ b/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) diff --git a/src/mapper.rs b/src/mapper.rs index d34dcd1b..13dc6466 100644 --- a/src/mapper.rs +++ b/src/mapper.rs @@ -718,30 +718,27 @@ pub fn align_paired_end_read( } // Compute per-read MAPQ. - // 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(len=1) would return 60 even for reads with hundreds - // of equally-good alternative NAMs. Use NAM-based per-read quality - // instead, exactly as the Proper fast path already does. let mapqs = { - let joint = alignment_quality(&alignment_pairs, |ap| ap.score as f32); - if alignment_pairs.len() == 1 { - let p = &alignment_pairs[0]; - let proper = is_proper_pair_opt( - p.alignment1.as_ref(), - p.alignment2.as_ref(), + // 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::Proper; - if !proper { - [ - mapping_quality(&nams_pair[0]), - mapping_quality(&nams_pair[1]), - ] - } else { - [joint, joint] - } + ) == 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] } };