Skip to content

--alignEndsType on pairs: extend a mate's two ends while stitching (100M, AS:182 as STAR) #368

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
BenjaminDEMAILLE wants to merge 2 commits into
base: perf/hotpaths
Choose a base branch
from
Open

--alignEndsType on pairs: extend a mate's two ends while stitching (100M, AS:182 as STAR) #368

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
051bc58
--alignEndsType on pairs: extend a mate's two ends while stitching
BenjaminDEMAILLE Nov 2, 2026
f7da8ee
End-to-end extension: check the chromosome edge before the mate spacer
BenjaminDEMAILLE Nov 3, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
4 changes: 0 additions & 4 deletions src/align/read_align.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
};
Expand All @@ -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;
};
Expand Down Expand Up @@ -981,7 +979,6 @@ pub fn align_paired_read(
orig_rev,
false,
false,
0, // mate1
) {
t.is_reverse = stitch_is_reverse;
any_transcript = true;
Expand Down Expand Up @@ -1021,7 +1018,6 @@ pub fn align_paired_read(
orig_rev,
false,
false,
1, // mate2
) {
t.is_reverse = !stitch_is_reverse;
any_transcript = true;
Expand Down
91 changes: 52 additions & 39 deletions src/align/stitch.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand Down Expand Up @@ -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 {
Expand Down Expand Up @@ -1271,6 +1273,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.
Expand Down Expand Up @@ -1350,7 +1358,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(
Expand All @@ -1366,8 +1375,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;
Expand Down Expand Up @@ -1427,8 +1439,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;
Expand Down Expand Up @@ -2047,24 +2062,12 @@ pub(crate) fn finalize_transcript(
original_is_reverse: bool,
no_left_ext: bool,
no_right_ext: bool,
imate: u8,
) -> Option<Transcript> {
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.
Expand Down Expand Up @@ -2102,7 +2105,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()
Expand All @@ -2121,7 +2124,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()
Expand All @@ -2142,7 +2145,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()
Expand All @@ -2161,20 +2164,14 @@ 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
};
(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
Expand Down Expand Up @@ -2509,7 +2506,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,
}
Expand Down Expand Up @@ -2764,6 +2763,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 {
Expand All @@ -2780,7 +2788,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()
Expand All @@ -2801,7 +2809,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()
Expand All @@ -2828,7 +2836,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
Expand All @@ -2850,13 +2858,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());
Expand Down Expand Up @@ -3369,7 +3385,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 {
Expand Down Expand Up @@ -4534,10 +4549,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);
}
Expand Down
Loading
Loading