From 80775fbfdde39a190667d42969b4d23c51c0800f Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 3 Nov 2026 10:59:02 +0100 Subject: [PATCH] peOverlap EndToEnd: keep the merged read's chromosome-edge transcript With --peOverlapNbasesMin and --alignEndsType EndToEnd, a pair whose overlapping mates hang over a chromosome edge stayed unmapped; STAR maps it soft-clipped at the edge (1S99M / 100M). STAR does not drop a transcript whose end-to-end extension hits the chromosome edge: it adds -999999999 and nMMmax+1 to it and goes on (stitchWindowAligns.cpp:44-55), then floors the score at 0 when it applies the genomic-length penalty (:224), so the window keeps it, unextended. For an ordinary read it is never reported. For the merged mates, once any window scores above 0, peOverlapSEtoPE splits every window transcript back into mates and rescores it from its exons (ReadAlign_peOverlapMergeMap.cpp:276-292), and the kept edge transcript wins. The merged read now keeps such transcripts at score 0, and all of its transcripts are split and rescored, not only those in the score range. A read whose transcripts all score 0 has no window, as in STAR (ReadAlign_stitchPieces.cpp:344). Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 6 +++ src/align/read_align.rs | 16 ++++-- src/align/score.rs | 21 ++++++++ src/align/stitch.rs | 38 +++++++++++---- tests/end_to_end_pairs.rs | 100 ++++++++++++++++++++++++++++++++++++-- 5 files changed, 163 insertions(+), 18 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9dddbb4..a2a1a52 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -52,6 +52,12 @@ Sections commonly used: Features, Bug fixes, Other changes. soft-clipped (`96M4S`, `AS:190` where STAR has `100M`, `AS:182`), and a mate that cannot reach its end could be clipped into a pair. A single-end read that runs past a chromosome end now counts as unmapped "other", as in STAR. +- `--peOverlapNbasesMin` with `--alignEndsType EndToEnd`: a pair whose + overlapping mates hang over a chromosome edge maps soft-clipped at the edge + (`1S99M` / `100M`), as in STAR, instead of staying unmapped. The merged + read's transcript that hits the edge is kept at score 0, unextended, and + every transcript of the merged read is split back into mates and rescored, + not only the best ones. ### Other changes diff --git a/src/align/read_align.rs b/src/align/read_align.rs index ddaf83d..acda14c 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -409,7 +409,9 @@ fn align_read_inner( // does NOT filter clusters by seed coverage. Removed to match STAR. // Step 3: Stitch seeds within each cluster - let scorer = AlignmentScorer::from_params(params); + let mut scorer = AlignmentScorer::from_params(params); + // Only the merged mates run without the mapped filter. + scorer.keep_end_to_end_kills = !mapped_filter; let mut transcripts = Vec::new(); // STAR's `trAll`: every window's transcripts, before any cross-window // dedup or filtering, which is what its chimeric detection reads. @@ -515,6 +517,13 @@ fn align_read_inner( transcripts.retain(|_| keep.next().unwrap_or(true)); } + // No transcript above score 0 leaves STAR's `trBest` at `trInit`: no window + // (`ReadAlign_stitchPieces.cpp:344`). The merged mates are then not split + // (`ReadAlign_peOverlapMergeMap.cpp:32`). + if transcripts.iter().all(|t| t.score <= 0) { + transcripts.clear(); + } + // STAR's `trBest` (ReadAlign_stitchPieces.cpp), over the transcripts in window // order: greater score, or equal score and shorter genomic length. Its score // and mismatches go on the unmapped records. @@ -534,8 +543,9 @@ fn align_read_inner( }); // STAR's multMapSelect (ReadAlign_multMapSelect.cpp): score range, then the - // output order and the primary. - if !transcripts.is_empty() { + // output order and the primary. The merged mates skip the range: STAR splits + // every window transcript back into mates (`ReadAlign_peOverlapMergeMap.cpp:276`). + if mapped_filter && !transcripts.is_empty() { let max_score = transcripts.iter().map(|t| t.score).max().unwrap_or(0); let score_threshold = max_score - params.out_filter_multimap_score_range; transcripts.retain(|t| t.score >= score_threshold); diff --git a/src/align/score.rs b/src/align/score.rs index 36d3606..e03d9e4 100644 --- a/src/align/score.rs +++ b/src/align/score.rs @@ -59,6 +59,13 @@ pub struct AlignmentScorer { /// `--chimSegmentMin > 0`: STAR then records every finalized transcript in /// its window, skipping the score gate (`stitchWindowAligns.cpp:247`). pub record_all_window_transcripts: bool, + /// Keep a transcript whose end-to-end extension hit a chromosome edge, + /// unextended and at window score 0, as STAR does. Only the merged read of + /// `--peOverlapNbasesMin` sets it: STAR rescores every window transcript of + /// that read once split into mates (`ReadAlign_peOverlapMergeMap.cpp:276-292`), + /// where a kept one can win. shortcut: other reads drop it, as STAR never + /// reports a score-0 transcript; upgrade if chimeric detection must see it. + pub keep_end_to_end_kills: bool, } /// The three intron checks STAR applies when a transcript is finalized @@ -151,6 +158,7 @@ impl AlignmentScorer { intron_filter: IntronFilter::default(), multimap_score_range: 1, record_all_window_transcripts: false, + keep_end_to_end_kills: false, } } @@ -198,6 +206,7 @@ impl AlignmentScorer { }, multimap_score_range: params.out_filter_multimap_score_range, record_all_window_transcripts: params.chim_segment_min > 0, + keep_end_to_end_kills: false, } } @@ -1257,6 +1266,7 @@ mod tests { intron_filter: IntronFilter::default(), multimap_score_range: 1, record_all_window_transcripts: false, + keep_end_to_end_kills: false, }; // Intron from position 2, length 12 (spans positions 2-13 inclusive) @@ -1305,6 +1315,7 @@ mod tests { intron_filter: IntronFilter::default(), multimap_score_range: 1, record_all_window_transcripts: false, + keep_end_to_end_kills: false, }; let motif = scorer.detect_splice_motif(2, 12, &genome); @@ -1352,6 +1363,7 @@ mod tests { intron_filter: IntronFilter::default(), multimap_score_range: 1, record_all_window_transcripts: false, + keep_end_to_end_kills: false, }; let motif = scorer.detect_splice_motif(2, 12, &genome); @@ -1397,6 +1409,7 @@ mod tests { intron_filter: IntronFilter::default(), multimap_score_range: 1, record_all_window_transcripts: false, + keep_end_to_end_kills: false, }; let motif = scorer.detect_splice_motif(2, 12, &genome); @@ -1435,6 +1448,7 @@ mod tests { intron_filter: IntronFilter::default(), multimap_score_range: 1, record_all_window_transcripts: false, + keep_end_to_end_kills: false, }; let (score, gap_type) = scorer.score_gap(0, 5, 0, &genome); @@ -1471,6 +1485,7 @@ mod tests { intron_filter: IntronFilter::default(), multimap_score_range: 1, record_all_window_transcripts: false, + keep_end_to_end_kills: false, }; // Small gap (< align_intron_min) is deletion @@ -1515,6 +1530,7 @@ mod tests { intron_filter: IntronFilter::default(), multimap_score_range: 1, record_all_window_transcripts: false, + keep_end_to_end_kills: false, }; // Gap starting at position 2 (GT), length 26 (>= 21) is splice junction @@ -1557,6 +1573,7 @@ mod tests { intron_filter: IntronFilter::default(), multimap_score_range: 1, record_all_window_transcripts: false, + keep_end_to_end_kills: false, }; let annotated_score = scorer.score_annotated_junction(0, true); @@ -1600,6 +1617,7 @@ mod tests { intron_filter: IntronFilter::default(), multimap_score_range: 1, record_all_window_transcripts: false, + keep_end_to_end_kills: false, }; // CT-AC motif: (1,3,0,1) — reverse complement of GT-AG @@ -1711,6 +1729,7 @@ mod tests { intron_filter: IntronFilter::default(), multimap_score_range: 1, record_all_window_transcripts: false, + keep_end_to_end_kills: false, }; // Gap of exactly 589824 starting at position 100 should be splice junction @@ -1794,6 +1813,7 @@ mod tests { intron_filter: IntronFilter::default(), multimap_score_range: 1, record_all_window_transcripts: false, + keep_end_to_end_kills: false, }; // Gap of 1001 (> 1000 max) should be deletion, not splice junction @@ -1849,6 +1869,7 @@ mod tests { intron_filter: IntronFilter::default(), multimap_score_range: 1, record_all_window_transcripts: false, + keep_end_to_end_kills: false, } } diff --git a/src/align/stitch.rs b/src/align/stitch.rs index b06ad0b..a84d9f9 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -10,9 +10,9 @@ use noodles::sam::alignment::record::cigar; pub(crate) const PE_SPACER_BASE: u8 = 11; /// Sentinel score returned by end-to-end extension (`--alignEndsType`) when the -/// forced extension runs into a chromosome boundary. STAR sets a large negative -/// score here; we use it to reject the transcript in `finalize_transcript`. -const EXTEND_TO_END_KILL: i32 = i32::MIN / 2; +/// forced extension runs into a chromosome boundary: STAR's own value +/// (`extendAlign.cpp:27`), which it adds to the transcript's score. +const EXTEND_TO_END_KILL: i32 = -999_999_999; /// Count mismatches in an alignment by comparing read sequence to genome sequence. /// @@ -2506,9 +2506,15 @@ fn wt_final_score(wt: &WorkingTranscript, scorer: &AlignmentScorer) -> i32 { match (wt.exons.first(), wt.exons.last()) { (Some(first), Some(last)) => { let span = last.genome_end - first.genome_start; - // Not floored: an end-to-end extension can take a transcript below zero, - // and then the window's score gate must reject it. - wt.score + scorer.genomic_length_penalty(span) + // STAR floors at 0 when it applies the penalty (`stitchWindowAligns.cpp:224`), + // which keeps a killed transcript at score 0. Only the merged mates keep + // kills (see `keep_end_to_end_kills`); elsewhere the gate must reject them. + let score = wt.score + scorer.genomic_length_penalty(span); + if scorer.keep_end_to_end_kills && scorer.score_genomic_length_log2_scale != 0.0 { + score.max(0) + } else { + score + } } _ => wt.score, } @@ -2867,16 +2873,25 @@ fn stitch_recurse( // An end-to-end extension that ran into a chromosome boundary kills the // transcript: it cannot reach the read's end, and it may not be clipped. - if first_ext.max_score == EXTEND_TO_END_KILL - || second_ext.max_score == EXTEND_TO_END_KILL - { + let killed = [&first_ext, &second_ext].map(|e| e.max_score == EXTEND_TO_END_KILL); + if killed.contains(&true) && !scorer.keep_end_to_end_kills { break 'leaf; } // Apply the extensions to wt - if first_ext.extend_len > 0 || second_ext.extend_len > 0 { + if first_ext.extend_len > 0 || second_ext.extend_len > 0 || killed.contains(&true) { mark = Some(wt.mark_ends()); } + // STAR adds a killed extension and goes on (`stitchWindowAligns.cpp:44-55`): + // no bases, the kill score and nMMmax+1 mismatches. + for (e, _) in [&first_ext, &second_ext] + .iter() + .zip(killed) + .filter(|(_, k)| *k) + { + wt.score += e.max_score; + wt.n_mismatch += e.n_mismatch; + } if do_left_first { if first_ext.extend_len > 0 { let first = wt.exons.first_mut().unwrap(); @@ -4067,6 +4082,7 @@ mod tests { intron_filter: crate::align::score::IntronFilter::default(), multimap_score_range: 1, record_all_window_transcripts: false, + keep_end_to_end_kills: false, }; let cluster = SeedCluster { alignments: Vec::new(), @@ -4360,6 +4376,7 @@ mod tests { intron_filter: crate::align::score::IntronFilter::default(), multimap_score_range: 1, record_all_window_transcripts: false, + keep_end_to_end_kills: false, }; // Left overhang (prev.length) = 3, below min of 5 @@ -4407,6 +4424,7 @@ mod tests { intron_filter: crate::align::score::IntronFilter::default(), multimap_score_range: 1, record_all_window_transcripts: false, + keep_end_to_end_kills: false, }; // Both overhangs >= 5 diff --git a/tests/end_to_end_pairs.rs b/tests/end_to_end_pairs.rs index 47f5728..c0509ab 100644 --- a/tests/end_to_end_pairs.rs +++ b/tests/end_to_end_pairs.rs @@ -57,13 +57,15 @@ type Record = (u32, u32, String, i32, u32); /// The alignment records of a run, and its `Log.final.out`. fn run(reads: &[&[u8]], options: &[&str]) -> (Vec, String) { + let fasta = format!(">c1\n{}\n", String::from_utf8_lossy(&genome())); + run_on(&fasta, reads, options) +} + +/// [`run`] on the genome in `fasta_text`. +fn run_on(fasta_text: &str, reads: &[&[u8]], options: &[&str]) -> (Vec, String) { let dir = TempDir::new().unwrap(); let fasta = dir.path().join("g.fa"); - fs::write( - &fasta, - format!(">c1\n{}\n", String::from_utf8_lossy(&genome())), - ) - .unwrap(); + fs::write(&fasta, fasta_text).unwrap(); let index = dir.path().join("i"); fs::create_dir_all(&index).unwrap(); cargo_bin_cmd!("rustar-aligner") @@ -233,3 +235,91 @@ fn single_end_read_past_the_chromosome_start_has_no_alignment() { assert_eq!(count("Number of reads unmapped: other"), "1"); assert_eq!(count("Number of reads unmapped: too short"), "0"); } + +/// `--peOverlapNbasesMin` maps overlapping mates as one merged read. When that read +/// hangs over a chromosome edge, STAR keeps its end-to-end transcript at score 0 +/// instead of dropping it (`stitchWindowAligns.cpp:44-55, 224`), and once another +/// window scores above 0 it splits every window transcript back into mates and +/// rescores them (`ReadAlign_peOverlapMergeMap.cpp:276-292`). The kept one then wins, +/// soft-clipped at the edge. Here the other window is a diverged copy on c2 of the +/// merged read's edge end. +#[test] +fn overlapping_mates_over_a_chromosome_edge_are_clipped_like_star() { + let g = genome(); + let diverged = |s: &[u8]| -> Vec { + let positions: Vec = (5..s.len()).step_by(10).collect(); + mutated(s, &positions) + }; + let options = [ + "--alignEndsType", + "EndToEnd", + "--peOverlapNbasesMin", + "10", + "--outFilterScoreMinOverLread", + "0.3", + "--outFilterMatchNminOverLread", + "0.3", + ]; + for over in [&b"T"[..], b"GA"] { + let k = over.len() as u32; + let decoy_start = [over, &g[..40], &diverged(&g[40..150])].concat(); + let decoy_end = [ + &diverged(&g[g.len() - 150..g.len() - 40])[..], + &g[g.len() - 40..], + over, + ] + .concat(); + let c2 = [ + &g[5000..6000], + &decoy_start, + &g[6000..7000], + &decoy_end, + &g[7000..8000], + ] + .concat(); + let fasta = format!( + ">c1\n{}\n>c2\n{}\n", + String::from_utf8_lossy(&g), + String::from_utf8_lossy(&c2) + ); + for fragment in [130u32, 169] { + let f = fragment as usize; + // Over the start: mate 2 forward from the overhang, mate 1 reverse. + let frag = [over, &g[..f - over.len()]].concat(); + let (records, _) = run_on( + &fasta, + &[&reverse_complement(&frag[f - 100..]), &frag[..100]], + &options, + ); + assert_eq!( + records, + [ + record(163, 1, &format!("{k}S{}M", 100 - k), 198 - k as i32, 0), + record(83, fragment - 99 - k, "100M", 198 - k as i32, 0), + ], + "start, overhang {k}, fragment {fragment}" + ); + // Over the end: mate 1 forward, mate 2 reverse into the overhang. + let frag = [&g[g.len() - f + over.len()..], over].concat(); + let (records, _) = run_on( + &fasta, + &[&frag[..100], &reverse_complement(&frag[f - 100..])], + &options, + ); + assert_eq!( + records, + [ + record(99, 20_001 - fragment + k, "100M", 198 - k as i32, 0), + record( + 147, + 19_901 + k, + &format!("{}M{k}S", 100 - k), + 198 - k as i32, + 0 + ), + ], + "end, overhang {k}, fragment {fragment}" + ); + } + } +}