From 5be003f48773afb698fd4c196b22aad89a658e1c Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 27 Aug 2026 21:49:33 +0200 Subject: [PATCH] fix(align): a soft-clipped variant of a pair is not a second locus #31 reported NH up to 20 where STAR caps at 14 on the same data. The extra loci are not different places in the genome: they are the same alignments with the leading exon soft-clipped off, scoring one point lower, kept as separate multimapper entries. STAR removes them in stitchWindowAligns.cpp, where a transcript whose blocks are a subset of another's and which scores lower is dropped. The equivalent check here is gated on the two transcripts having the same number of exons, and these variants have one fewer, so they slip past the very comparison meant to catch them. The rule is applied to whole pairs, after the score-range filter: a pair whose mate1 and mate2 are both covered by a better-scoring pair is dropped. Equal scores keep both, matching STAR's strict comparison. Measured, nf-core/rnaseq test data (50 000 pairs, --outFilterMultimapNmax 20): deepest NH STAR 14 before 20 after 14 uniquely mapped STAR 41691 before 41665 after 41684 multi-mapped STAR 934 before 961 after 942 and on the project's yeast tier (10 000 pairs of ERR12389696), no regression: uniquely mapped STAR 7861 before 7860 after 7860 same chr/pos/CIGAR before 98.4825% after 98.4825% same NH before 99.9644% after 99.9763% Three earlier attempts at this regressed the yeast tier and are not in this diff: applying the rule inside the stitcher, and applying it to single transcripts either before or after the quality filters. All three cost ~39 uniquely mapped reads, because a subset alignment is sometimes the only one that survives the score and overhang gates, and because align_read is reused for the paired mate-overlap merge with those gates disabled. Closes #31. --- CHANGELOG.md | 9 +++ src/align/read_align.rs | 162 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 171 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 14b4346..23d012b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -42,6 +42,15 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Bug fixes +- Paired-end multimapper depth no longer counts soft-clipped variants of an + alignment as separate loci. STAR drops an alignment whose blocks are a subset + of another's and which scores lower (`stitchWindowAligns.cpp`); the same rule + now applies to whole pairs, after the score-range filter. On the + nf-core/rnaseq test data the deepest multimapper falls from NH=20 to NH=14, + which is STAR's own maximum, and the uniquely-mapped count moves from 41 665 + to 41 684 against STAR's 41 691. No change on the yeast tier: 7860 uniquely + mapped either way, 98.48% of mates at the same position, with NH agreement + rising from 99.964% to 99.976%. Closes #31. - Read names are cut at `--readNameSeparator` (default `/`), as STAR does. A read named `foo/1` was previously emitted as `foo/1` where STAR emits `foo`. diff --git a/src/align/read_align.rs b/src/align/read_align.rs index ec66a16..3811b88 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -132,6 +132,75 @@ pub enum PairedAlignmentResult { }, } +/// Drop a pair whose two mates are both covered by another pair scoring at +/// least as well, using STAR's mapped-length overlap test per mate. +fn dedup_pair_subsets(pairs: &mut Vec) { + if pairs.len() < 2 { + return; + } + let mapped = |t: &Transcript| -> u32 { + t.exons + .iter() + .map(|e| (e.read_end - e.read_start) as u32) + .sum() + }; + let subset_of = |a: &Transcript, b: &Transcript| -> bool { + mapped(a).saturating_sub(blocks_overlap_transcripts(a, b)) == 0 + }; + + let mut kept: Vec = Vec::with_capacity(pairs.len()); + for p in pairs.drain(..) { + let dominated = kept.iter().any(|k| { + k.combined_wt_score > p.combined_wt_score + && subset_of(&p.mate1_transcript, &k.mate1_transcript) + && subset_of(&p.mate2_transcript, &k.mate2_transcript) + }); + if dominated { + continue; + } + kept.retain(|k| { + !(p.combined_wt_score > k.combined_wt_score + && subset_of(&k.mate1_transcript, &p.mate1_transcript) + && subset_of(&k.mate2_transcript, &p.mate2_transcript)) + }); + kept.push(p); + } + *pairs = kept; +} + +/// Bases covered by both transcripts, counting only blocks that sit on the +/// same read-to-genome diagonal. Port of STAR's `blocksOverlap`, including its +/// rule of advancing both cursors when two blocks end at the same read +/// position. +fn blocks_overlap_transcripts(t1: &Transcript, t2: &Transcript) -> u32 { + let (mut i1, mut i2, mut overlap) = (0usize, 0usize, 0u32); + while i1 < t1.exons.len() && i2 < t2.exons.len() { + let a = &t1.exons[i1]; + let b = &t2.exons[i2]; + let (rs1, re1) = (a.read_start, a.read_end); + let (rs2, re2) = (b.read_start, b.read_end); + + if rs1 >= re2 { + i2 += 1; + } else if rs2 >= re1 { + i1 += 1; + } else { + let diag1 = a.genome_start as i64 - rs1 as i64; + let diag2 = b.genome_start as i64 - rs2 as i64; + if diag1 == diag2 { + overlap += (re1.min(re2) - rs1.max(rs2)) as u32; + } + if re1 >= re2 { + i2 += 1; + } + if re2 >= re1 { + i1 += 1; + } + } + } + overlap +} + /// Align a read to the genome. /// /// # Algorithm @@ -1142,6 +1211,13 @@ pub fn align_paired_read( .unwrap_or(0); let score_threshold = best_score - params.out_filter_multimap_score_range; joint_pairs.retain(|pa| pa.combined_wt_score >= score_threshold); + + // STAR drops an alignment whose blocks are a subset of another's and + // which scores lower (`stitchWindowAligns.cpp`). Applied here to whole + // pairs: the extra loci #31 reports are pairs whose mate1 is a + // soft-clipped variant of another pair's mate1, with the same mate2, + // scoring one point less. + dedup_pair_subsets(&mut joint_pairs); } // Deterministic primary tie-break (combined score, then a fixed positional @@ -1550,6 +1626,92 @@ pub(crate) fn pe_junctions_consistent(left: &Transcript, right: &Transcript) -> #[cfg(test)] mod tests { use super::*; + + /// The case #31 reports: a pair whose mate1 is a soft-clipped variant of + /// another pair's mate1, with the same mate2 and a lower score, is not a + /// second locus. + #[test] + fn a_pair_that_is_a_subset_of_a_better_one_is_dropped() { + let full = pair_for_dedup(0, 1_000, 0, 100, 5_000, 174); + let clipped = pair_for_dedup(0, 1_009, 9, 100, 5_000, 173); + let mut pairs = vec![full, clipped]; + dedup_pair_subsets(&mut pairs); + assert_eq!(pairs.len(), 1); + assert_eq!(pairs[0].combined_wt_score, 174); + } + + /// Order must not decide the outcome: the subset arriving first is still + /// the one removed. + #[test] + fn the_subset_is_dropped_whichever_order_it_arrives_in() { + let full = pair_for_dedup(0, 1_000, 0, 100, 5_000, 174); + let clipped = pair_for_dedup(0, 1_009, 9, 100, 5_000, 173); + let mut pairs = vec![clipped, full]; + dedup_pair_subsets(&mut pairs); + assert_eq!(pairs.len(), 1); + assert_eq!(pairs[0].combined_wt_score, 174); + } + + /// Two genuinely different loci are both kept, even at different scores: + /// neither mate1 covers the other. + #[test] + fn distinct_loci_are_both_kept() { + let a = pair_for_dedup(0, 1_000, 0, 100, 5_000, 174); + let b = pair_for_dedup(0, 40_000, 0, 100, 44_000, 173); + let mut pairs = vec![a, b]; + dedup_pair_subsets(&mut pairs); + assert_eq!(pairs.len(), 2); + } + + /// Equal scores keep both, as STAR does: its test is a strict `<`. + #[test] + fn an_equal_scoring_subset_is_kept() { + let full = pair_for_dedup(0, 1_000, 0, 100, 5_000, 174); + let clipped = pair_for_dedup(0, 1_009, 9, 100, 5_000, 174); + let mut pairs = vec![full, clipped]; + dedup_pair_subsets(&mut pairs); + assert_eq!(pairs.len(), 2); + } + + fn pair_for_dedup( + chr_idx: usize, + m1_genome_start: u64, + m1_read_start: usize, + m1_read_end: usize, + m2_genome_start: u64, + score: i32, + ) -> PairedAlignment { + let mate = |gs: u64, rs: usize, re: usize| Transcript { + chr_idx, + genome_start: gs, + genome_end: gs + (re - rs) as u64, + is_reverse: false, + exons: vec![crate::align::transcript::Exon { + genome_start: gs, + genome_end: gs + (re - rs) as u64, + read_start: rs, + read_end: re, + i_frag: 0, + }], + cigar: Vec::new(), + score, + n_mismatch: 0, + n_gap: 0, + n_junction: 0, + junction_motifs: Vec::new(), + junction_annotated: Vec::new(), + }; + PairedAlignment { + mate1_transcript: mate(m1_genome_start, m1_read_start, m1_read_end), + mate2_transcript: mate(m2_genome_start, 0, 100), + mate1_region: (m1_read_start, m1_read_end), + mate2_region: (0, 100), + is_proper_pair: true, + insert_size: 0, + combined_wt_score: score, + combined_n_match: 200, + } + } use crate::genome::Genome; use crate::index::packed_array::PackedArray; use crate::index::sa_index::SaIndex;