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}" + ); + } + } +}