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 5f9994bc..13dc6466 100644 --- a/src/mapper.rs +++ b/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,6 +1379,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 +1387,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 +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,