From 051bc58ec74dfc4ba98eefe473538224693cb073 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 3 Nov 2026 00:58:53 +0100 Subject: [PATCH 1/2] --alignEndsType on pairs: extend a mate's two ends while stitching Both ends of a mate are now extended in the stitcher, with ext[iFrag][iEnd] applied at the mate boundary and at the transcript's outer ends. The extension's mismatches count in the pair's score and limits. finalize_transcript worked on one mate's slice after the pair was already chosen, so inner ends stayed clipped and a mate that cannot reach its end could be clipped into a pair. Transcript scores are not floored at 0 before the window's score gate, so a negative end-to-end transcript loses to an existing one. Part of #354 (finding 3). Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 7 ++ src/align/read_align.rs | 4 - src/align/stitch.rs | 81 ++++++++------ tests/end_to_end_pairs.rs | 216 ++++++++++++++++++++++++++++++++++++++ 4 files changed, 269 insertions(+), 39 deletions(-) create mode 100644 tests/end_to_end_pairs.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index b18a76e..9dddbb4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -45,6 +45,13 @@ Sections commonly used: Features, Bug fixes, Other changes. - `--soloUMIfiltering` applies to the STARsolo `SJ` feature as it does to genes: a read that crosses two junctions is one UMI on two features, a tie `MultiGeneUMI_CR` drops. +- `--alignEndsType EndToEnd`, `Extend5pOfRead1`, `Extend5pOfReads12` and + `Extend3pOfRead1` on paired-end reads: both ends of a mate are extended while + the pair's transcript is stitched, so the extension's mismatches count in the + pair's score and limits, as in STAR. The inner ends of the mates were left + 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. ### Other changes diff --git a/src/align/read_align.rs b/src/align/read_align.rs index 35e450a..ddaf83d 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -857,7 +857,6 @@ pub fn align_paired_read( m1_orig_rev, m1_orig_rev, // no_left_ext = inner for reverse (orig_is_rev=true) !m1_orig_rev, // no_right_ext = inner for forward (orig_is_rev=false) - 0, // mate1 ) else { continue; }; @@ -870,7 +869,6 @@ pub fn align_paired_read( m2_orig_rev, m2_orig_rev, // no_left_ext = inner for reverse (orig_is_rev=true) !m2_orig_rev, // no_right_ext = inner for forward (orig_is_rev=false) - 1, // mate2 ) else { continue; }; @@ -981,7 +979,6 @@ pub fn align_paired_read( orig_rev, false, false, - 0, // mate1 ) { t.is_reverse = stitch_is_reverse; any_transcript = true; @@ -1021,7 +1018,6 @@ pub fn align_paired_read( orig_rev, false, false, - 1, // mate2 ) { t.is_reverse = !stitch_is_reverse; any_transcript = true; diff --git a/src/align/stitch.rs b/src/align/stitch.rs index 4fe809f..3cfec7f 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -1271,6 +1271,12 @@ impl WorkingTranscript { } } +/// A seed's fragment in `--alignEndsType`'s table. Single-end seeds carry mate_id 2; +/// STAR's iFrag is 0 there. +fn ext_frag(mate_id: u8) -> usize { + usize::from(mate_id == 1) +} + /// Fill the gap between the last exon in a WorkingTranscript and the next WA entry. /// Returns None if stitching fails (mismatch limit, overhang too short, etc.). /// Matches STAR's stitchAlignToTranscript.cpp logic. @@ -1350,7 +1356,8 @@ fn stitch_align_in_place( // STAR stitchAlignToTranscript.cpp:374-381: right-extend mate A to fragment boundary. // extendAlign(R, G, rAend+1, gAend+1, 1, 1, DEF_readSeqLengthMax, nMatch, nMM, ...) - // extendToEnd=false (local mode); spacer at position 150 naturally stops the extension. + // The spacer at the mate's end stops the extension, also an end-to-end one. These + // are the mates' inner ends, their 3' ends. // nMatch approximated as score + n_mismatch (exact for seeds-only paths, approx otherwise). let n_match_m1 = (new_wt.score + new_wt.n_mismatch as i32).max(0) as usize; let right_ext = extend_alignment_cached( @@ -1366,8 +1373,11 @@ fn stitch_align_in_place( scorer.p_mm_max, index, cluster.is_reverse, - false, // internal stitch extension: alignEndsType applied at finalize + scorer.align_ends_type.ext[ext_frag(last_mate)][1], ); + if right_ext.max_score == EXTEND_TO_END_KILL { + return false; + } if right_ext.extend_len > 0 { let last = new_wt.exons.last_mut().unwrap(); last.read_end += right_ext.extend_len; @@ -1427,8 +1437,11 @@ fn stitch_align_in_place( scorer.p_mm_max, index, cluster.is_reverse, - false, // internal stitch extension: alignEndsType applied at finalize + scorer.align_ends_type.ext[ext_frag(wa.mate_id)][1], ); + if left_ext.max_score == EXTEND_TO_END_KILL { + return false; + } if left_ext.extend_len > 0 { let last = new_wt.exons.last_mut().unwrap(); last.read_start -= left_ext.extend_len; @@ -2047,24 +2060,12 @@ pub(crate) fn finalize_transcript( original_is_reverse: bool, no_left_ext: bool, no_right_ext: bool, - imate: u8, ) -> Option { use crate::align::transcript::Exon; let alignment_start = wt.read_start; let alignment_end = wt.read_end; - // --alignEndsType: decide whether each end is force-extended (no soft-clip). - // `read_seq` is a single mate's slice, so `original_is_reverse` here equals - // STAR's `(Str != imate)`: when the slice is reverse-complemented, the read's - // 5' end sits on the right. STAR indexes ext[imate][iEnd] with iEnd 0 = 5', - // 1 = 3', so the left/start end maps to `original_is_reverse` and the - // right/end to its complement. - let ext = scorer.align_ends_type.ext; - let mi = (imate as usize).min(1); - let extend_left_to_end = ext[mi][original_is_reverse as usize]; - let extend_right_to_end = ext[mi][(!original_is_reverse) as usize]; - // Guard: exon positions must be within read bounds. Out-of-bounds positions can // arise when a combined PE read's WorkingTranscript has junction shifts that extend // exon read_end past the individual mate boundary. Filter rather than underflow. @@ -2102,7 +2103,7 @@ pub(crate) fn finalize_transcript( scorer.p_mm_max, index, cluster.is_reverse, - extend_left_to_end, + false, // --alignEndsType was applied in the stitcher ) } else { zero_extend.clone() @@ -2121,7 +2122,7 @@ pub(crate) fn finalize_transcript( scorer.p_mm_max, index, cluster.is_reverse, - extend_right_to_end, + false, // --alignEndsType was applied in the stitcher ) } else { zero_extend.clone() @@ -2142,7 +2143,7 @@ pub(crate) fn finalize_transcript( scorer.p_mm_max, index, cluster.is_reverse, - extend_right_to_end, + false, // --alignEndsType was applied in the stitcher ) } else { zero_extend.clone() @@ -2161,7 +2162,7 @@ pub(crate) fn finalize_transcript( scorer.p_mm_max, index, cluster.is_reverse, - extend_left_to_end, + false, // --alignEndsType was applied in the stitcher ) } else { zero_extend @@ -2169,12 +2170,6 @@ pub(crate) fn finalize_transcript( (left, right) }; - // End-to-end extension that ran into a chromosome boundary kills the - // transcript (STAR: no soft-clipping past the read end in EndToEnd mode). - if left_extend.max_score == EXTEND_TO_END_KILL || right_extend.max_score == EXTEND_TO_END_KILL { - return None; - } - // STAR finalization check (`stitchWindowAligns.cpp:96-108`): exon lengths // including repeat lengths (shiftSJ). An unannotated junction needs // `alignSJoverhangMin + shiftSJ` on both sides. An annotated one needs @@ -2509,7 +2504,9 @@ 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; - (wt.score + scorer.genomic_length_penalty(span)).max(0) + // 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) } _ => wt.score, } @@ -2764,6 +2761,15 @@ fn stitch_recurse( }; let do_left_first = !original_is_reverse; + // --alignEndsType, as STAR's `ext[iFrag][iEnd]` with iEnd 0 = 5': the + // left end is a 5' end where the mate is on its own strand. + let ends = scorer.align_ends_type.ext; + let frag_first = ext_frag(wt.exons[0].mate_id); + let frag_last = ext_frag(wt.exons[wt.exons.len() - 1].mate_id); + let extend_left_to_end = + ends[frag_first][usize::from(usize::from(original_is_reverse) != frag_first)]; + let extend_right_to_end = + ends[frag_last][usize::from(usize::from(original_is_reverse) == frag_last)]; let first_ext = if do_left_first { // Left extension (5' of forward read) if wt.read_start > 0 { @@ -2780,7 +2786,7 @@ fn stitch_recurse( scorer.p_mm_max, index, cluster.is_reverse, - false, // internal stitch extension: alignEndsType applied at finalize + extend_left_to_end, ) } else { zero_ext.clone() @@ -2801,7 +2807,7 @@ fn stitch_recurse( scorer.p_mm_max, index, cluster.is_reverse, - false, // internal stitch extension: alignEndsType applied at finalize + extend_right_to_end, ) } else { zero_ext.clone() @@ -2828,7 +2834,7 @@ fn stitch_recurse( scorer.p_mm_max, index, cluster.is_reverse, - false, // internal stitch extension: alignEndsType applied at finalize + extend_right_to_end, ) } else { zero_ext @@ -2850,13 +2856,21 @@ fn stitch_recurse( scorer.p_mm_max, index, cluster.is_reverse, - false, // internal stitch extension: alignEndsType applied at finalize + extend_left_to_end, ) } else { zero_ext } }; + // 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 + { + break 'leaf; + } + // Apply the extensions to wt if first_ext.extend_len > 0 || second_ext.extend_len > 0 { mark = Some(wt.mark_ends()); @@ -3369,7 +3383,6 @@ pub(crate) fn stitch_seeds_with_jdb_debug( stitch_is_reverse, false, false, - 0, // single-end read is mate 0 ) { // Restore original reverse-strand flag and read sequence for SAM output. if stitch_is_reverse { @@ -4534,10 +4547,8 @@ mod tests { wt.genome_end = 8; wt.score = 6; wt.n_mismatch = 1; - let t = finalize_transcript( - &wt, &read, &index, &scorer, &cluster, false, false, false, 0, - ) - .expect("finalizes"); + let t = finalize_transcript(&wt, &read, &index, &scorer, &cluster, false, false, false) + .expect("finalizes"); assert_eq!(t.cigar_string(), "8M"); assert_eq!(t.n_mismatch, 1); } diff --git a/tests/end_to_end_pairs.rs b/tests/end_to_end_pairs.rs new file mode 100644 index 0000000..bb4ceb5 --- /dev/null +++ b/tests/end_to_end_pairs.rs @@ -0,0 +1,216 @@ +//! `--alignEndsType` on the two ends of a mate that meet in a pair. Expected values +//! are STAR 2.7.11b's observed output on this fixture (same genome and reads, +//! `--genomeSAindexNbases 7`). + +use assert_cmd::cargo::cargo_bin_cmd; +use std::fs; +use std::path::Path; +use tempfile::TempDir; + +fn genome() -> Vec { + let mut state = 11u32; + (0..20_000) + .map(|_| { + state = state.wrapping_mul(1_103_515_245).wrapping_add(12345); + b"ACGT"[((state >> 16) & 3) as usize] + }) + .collect() +} + +fn reverse_complement(seq: &[u8]) -> Vec { + seq.iter().rev().map(|&b| complement(b)).collect() +} + +fn complement(base: u8) -> u8 { + match base { + b'A' => b'T', + b'C' => b'G', + b'G' => b'C', + _ => b'A', + } +} + +/// `seq` with a mismatch at each of `positions`. +fn mutated(seq: &[u8], positions: &[usize]) -> Vec { + let mut out = seq.to_vec(); + for &p in positions { + out[p] = complement(out[p]); + } + out +} + +const LAST_FOUR: [usize; 4] = [96, 97, 98, 99]; +const FIRST_FOUR: [usize; 4] = [0, 1, 2, 3]; + +/// A pair 400 bases long at 2001, with 4 mismatches at each of the mates' `first` +/// and `second` positions on the forward strand. +fn pair(first: &[usize], second: &[usize]) -> (Vec, Vec) { + let g = genome(); + ( + mutated(&g[2000..2100], first), + reverse_complement(&mutated(&g[2300..2400], second)), + ) +} + +/// FLAG, POS, CIGAR, AS, nM. +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 dir = TempDir::new().unwrap(); + let fasta = dir.path().join("g.fa"); + fs::write( + &fasta, + format!(">c1\n{}\n", String::from_utf8_lossy(&genome())), + ) + .unwrap(); + let index = dir.path().join("i"); + fs::create_dir_all(&index).unwrap(); + cargo_bin_cmd!("rustar-aligner") + .args(["--runMode", "genomeGenerate", "--genomeSAindexNbases", "7"]) + .args(["--genomeDir", index.to_str().unwrap()]) + .args(["--genomeFastaFiles", fasta.to_str().unwrap()]) + .assert() + .success(); + let files: Vec = reads + .iter() + .enumerate() + .map(|(i, seq)| { + let path = dir.path().join(format!("r{i}.fq")); + let qual = "I".repeat(seq.len()); + fs::write( + &path, + format!("@r\n{}\n+\n{qual}\n", String::from_utf8_lossy(seq)), + ) + .unwrap(); + path.to_str().unwrap().to_string() + }) + .collect(); + let prefix = format!("{}/", dir.path().display()); + cargo_bin_cmd!("rustar-aligner") + .args(["--genomeDir", index.to_str().unwrap()]) + .arg("--readFilesIn") + .args(&files) + .args(["--outFileNamePrefix", &prefix, "--outSAMtype", "SAM"]) + .args(options) + .assert() + .success(); + let sam = fs::read_to_string(Path::new(&prefix).join("Aligned.out.sam")).unwrap(); + let tag = |fields: &[&str], tag: &str| -> i64 { + let found = fields.iter().find(|f| f.starts_with(tag)).unwrap(); + found[tag.len() + 3..].parse().unwrap() + }; + let records = sam + .lines() + .filter(|l| !l.starts_with('@')) + .map(|l| l.split('\t').collect::>()) + .filter(|f| f[1].parse::().unwrap() & 4 == 0) + .map(|f| { + ( + f[1].parse().unwrap(), + f[3].parse().unwrap(), + f[5].to_string(), + tag(&f, "AS") as i32, + tag(&f, "nM") as u32, + ) + }) + .collect(); + let log = fs::read_to_string(Path::new(&prefix).join("Log.final.out")).unwrap(); + (records, log) +} + +fn record(flag: u32, pos: u32, cigar: &str, score: i32, nm: u32) -> Record { + (flag, pos, cigar.to_string(), score, nm) +} + +fn end_to_end_pair(mate1: &[u8], mate2: &[u8]) -> Vec { + run(&[mate1, mate2], &["--alignEndsType", "EndToEnd"]).0 +} + +#[test] +fn inner_ends_are_extended_and_charged() { + let (m1, m2) = pair(&LAST_FOUR, &FIRST_FOUR); + assert_eq!( + end_to_end_pair(&m1, &m2), + [ + record(99, 2001, "100M", 182, 8), + record(147, 2301, "100M", 182, 8), + ] + ); + // The same pair with the mates swapped, on the other strand. + assert_eq!( + end_to_end_pair(&m2, &m1), + [ + record(163, 2001, "100M", 182, 8), + record(83, 2301, "100M", 182, 8), + ] + ); +} + +#[test] +fn outer_ends_are_charged_to_the_pair() { + let (m1, m2) = pair(&FIRST_FOUR, &LAST_FOUR); + assert_eq!( + end_to_end_pair(&m1, &m2), + [ + record(99, 2001, "100M", 182, 8), + record(147, 2301, "100M", 182, 8), + ] + ); +} + +#[test] +fn a_mate_past_the_chromosome_start_cannot_be_clipped_into_a_pair() { + let g = genome(); + let mut over_the_edge = b"ACGTT".to_vec(); + over_the_edge.extend_from_slice(&g[0..95]); + let mate2 = reverse_complement(&g[300..400]); + assert!(end_to_end_pair(&over_the_edge, &mate2).is_empty()); + + // With filters low enough for one mate, only the mate that reaches its ends maps. + let (records, _) = run( + &[&over_the_edge, &mate2], + &[ + "--alignEndsType", + "EndToEnd", + "--outFilterScoreMinOverLread", + "0.3", + "--outFilterMatchNminOverLread", + "0.3", + ], + ); + assert_eq!(records, [record(153, 301, "100M", 98, 0)]); +} + +#[test] +fn one_sided_modes_extend_only_their_inner_end() { + let (m1, m2) = pair(&LAST_FOUR, &FIRST_FOUR); + // Read 1's 3' end is forced, mate 2's 3' end stays clipped. + let (records, _) = run(&[&m1, &m2], &["--alignEndsType", "Extend3pOfRead1"]); + assert_eq!( + records, + [ + record(99, 2001, "100M", 186, 4), + record(147, 2305, "4S96M", 186, 4), + ] + ); +} + +#[test] +fn single_end_read_past_the_chromosome_start_has_no_alignment() { + let g = genome(); + let mut over_the_edge = b"ACGTT".to_vec(); + over_the_edge.extend_from_slice(&g[0..95]); + let (records, log) = run(&[&over_the_edge], &["--alignEndsType", "EndToEnd"]); + assert!(records.is_empty()); + let count = |label: &str| { + log.lines() + .find(|l| l.contains(label)) + .and_then(|l| l.split('|').nth(1)) + .unwrap() + .trim() + .to_string() + }; + assert_eq!(count("Number of reads unmapped: other"), "1"); + assert_eq!(count("Number of reads unmapped: too short"), "0"); +} From f7da8eeee1f2f0447ece0a1d8b894bdc9232cb37 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 3 Nov 2026 09:22:18 +0100 Subject: [PATCH 2/2] End-to-end extension: check the chromosome edge before the mate spacer STAR's extendAlign checks the chromosome boundary first and the fragment spacer second, so a mate whose 3' end sits on a chromosome edge is killed under EndToEnd / Extend3pOfRead1 instead of stopping at the spacer. With mate 1 the last 100 bases of the chromosome and mate 2 its reverse complement, STAR reports mate 2 alone (153 100M AS:98); this branch reported a pair with AS:198. It now matches STAR, and with an unrelated mate 2 both leave the read unmapped. Co-Authored-By: Claude Opus 5.5 --- src/align/stitch.rs | 10 ++++++---- tests/end_to_end_pairs.rs | 19 +++++++++++++++++++ 2 files changed, 25 insertions(+), 4 deletions(-) diff --git a/src/align/stitch.rs b/src/align/stitch.rs index 3cfec7f..b06ad0b 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -306,10 +306,6 @@ fn extend_alignment( break; } let read_base = read_seq[read_pos]; - // No extension through the inter-mate spacer. - if read_base == PE_SPACER_BASE { - break; - } let genome_pos = if direction > 0 { genome_start + i as u64 } else { @@ -339,6 +335,12 @@ fn extend_alignment( n_mismatch: n_mm_max + 1, }; } + // No extension through the inter-mate spacer. Checked after the + // boundary, as STAR does (`extendAlign.cpp:24-31`): a mate ending on a + // chromosome edge is killed, not stopped at the spacer. + if read_base == PE_SPACER_BASE { + break; + } ext_len = i + 1; // Ns (read or genome) carry no score. if read_base == 4 || genome_base == 4 { diff --git a/tests/end_to_end_pairs.rs b/tests/end_to_end_pairs.rs index bb4ceb5..47f5728 100644 --- a/tests/end_to_end_pairs.rs +++ b/tests/end_to_end_pairs.rs @@ -182,6 +182,25 @@ fn a_mate_past_the_chromosome_start_cannot_be_clipped_into_a_pair() { assert_eq!(records, [record(153, 301, "100M", 98, 0)]); } +/// STAR checks the chromosome edge before the mate spacer (`extendAlign.cpp`), +/// so a mate whose 3' end sits on the edge cannot be extended into a pair. +#[test] +fn a_mate_ending_on_the_chromosome_end_cannot_be_extended_into_a_pair() { + let g = genome(); + let at_end = &g[g.len() - 100..]; + let low = [ + "--outFilterScoreMinOverLread", + "0.3", + "--outFilterMatchNminOverLread", + "0.3", + ]; + for mode in ["EndToEnd", "Extend3pOfRead1"] { + let options = [&["--alignEndsType", mode][..], &low[..]].concat(); + let (records, _) = run(&[at_end, &reverse_complement(at_end)], &options); + assert_eq!(records, [record(153, 19_901, "100M", 98, 0)], "{mode}"); + } +} + #[test] fn one_sided_modes_extend_only_their_inner_end() { let (m1, m2) = pair(&LAST_FOUR, &FIRST_FOUR);