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;