From 3016b851c63b2c797480fa9cf5575fa4fcb8ca63 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 29 Sep 2026 19:27:14 +0200 Subject: [PATCH] fix(align): pair-subset dedup only compares mates on one chromosome and strand STAR runs blocksOverlap inside a single window, so both transcripts always share chromosome and strand. Require the same in dedup_pair_subsets so a diagonal coincidence across strands or chromosomes can never count as a subset. Adds tests for the strand guard and blocksOverlap diagonal logic. Refs #31 Co-Authored-By: Claude Opus 5.5 --- src/align/read_align.rs | 30 +++++++++++++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/src/align/read_align.rs b/src/align/read_align.rs index 3811b88..8439481 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -144,8 +144,13 @@ fn dedup_pair_subsets(pairs: &mut Vec) { .map(|e| (e.read_end - e.read_start) as u32) .sum() }; + // STAR compares transcripts within one window, so both share chromosome + // and strand; require the same here so a diagonal coincidence across + // strands or chromosomes can never count as overlap. let subset_of = |a: &Transcript, b: &Transcript| -> bool { - mapped(a).saturating_sub(blocks_overlap_transcripts(a, b)) == 0 + a.chr_idx == b.chr_idx + && a.is_reverse == b.is_reverse + && mapped(a).saturating_sub(blocks_overlap_transcripts(a, b)) == 0 }; let mut kept: Vec = Vec::with_capacity(pairs.len()); @@ -1673,6 +1678,29 @@ mod tests { assert_eq!(pairs.len(), 2); } + /// A mate on the opposite strand is never a subset, even when its blocks + /// sit on the same read-to-genome diagonal. + #[test] + fn a_lower_scoring_pair_on_the_other_strand_is_kept() { + let full = pair_for_dedup(0, 1_000, 0, 100, 5_000, 174); + let mut clipped = pair_for_dedup(0, 1_009, 9, 100, 5_000, 173); + clipped.mate1_transcript.is_reverse = true; + let mut pairs = vec![full, clipped]; + dedup_pair_subsets(&mut pairs); + assert_eq!(pairs.len(), 2); + } + + /// STAR's `blocksOverlap`: only bases on the same diagonal count. + #[test] + fn blocks_overlap_counts_only_the_shared_diagonal() { + let a = pair_for_dedup(0, 1_000, 0, 100, 5_000, 0).mate1_transcript; + let b = pair_for_dedup(0, 1_009, 9, 100, 5_000, 0).mate1_transcript; + let c = pair_for_dedup(0, 1_010, 9, 100, 5_000, 0).mate1_transcript; + assert_eq!(blocks_overlap_transcripts(&a, &b), 91); + assert_eq!(blocks_overlap_transcripts(&b, &a), 91); + assert_eq!(blocks_overlap_transcripts(&a, &c), 0); + } + fn pair_for_dedup( chr_idx: usize, m1_genome_start: u64,