peOverlap EndToEnd: keep the merged read's chromosome-edge transcript

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
 

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);

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,
         }
     }
 

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

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<Record>, 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<Record>, 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<u8> {
+        let positions: Vec<usize> = (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}"
+            );
+        }
+    }
+}