Skip to content

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

Open
BenjaminDEMAILLE wants to merge 1 commit into
fix/end-to-end-pairsfrom
fix/pe-overlap-chr-edge
Open

BenjaminDEMAILLE wants to merge 1 commit into
fix/end-to-end-pairsfrom
fix/pe-overlap-chr-edge

Conversation

@BenjaminDEMAILLE

Copy link
Copy Markdown
Contributor

Part of #354 (--peOverlapNbasesMin with EndToEnd: mates overlapping at a chromosome edge).

Stacked on #368. The base is fix/end-to-end-pairs; only the last commit is this PR. It merges cleanly with #367 (fix/half-mapped-multimap).

Problem

With --alignEndsType EndToEnd --peOverlapNbasesMin 10, a pair whose mates overlap and hang 1 to 4 bases over a chromosome edge stayed unmapped. STAR maps it, with the overhang soft-clipped on the edge mate:

STAR     q1870  163  chr1  1   1S99M  AS:195 nM:1      rustar  q1870  141  *  (unmapped)
STAR     q1870   83  chr1  69  100M   AS:195 nM:1      rustar  q1870   77  *  (unmapped)
STAR     q3260  163  chr1  6   8S92M  AS:188 nM:1      rustar  q3260  141  *  (unmapped)

A clip at the edge looks odd under EndToEnd. It comes from how STAR handles the merged read.

What STAR does (2.7.11b source)

  • extendAlign.cpp:25-32: an end-to-end extension that reaches a chromosome edge returns true with extendL=0, maxScore=-999999999, nMM=nMMmax+1.
  • stitchWindowAligns.cpp:44-55, 65-77: the caller adds that to the transcript and goes on. Nothing drops it.
  • stitchWindowAligns.cpp:221-225: when scoreGenomicLengthLog2scale != 0 (default -0.25), Score = max(0, Score). The transcript is kept in the window at score 0, unextended.
  • ReadAlign_stitchPieces.cpp:328, 344: a score-0 transcript never becomes trBest. If no window scores above 0, nW=0 and the read is unmapped. So for an ordinary read the kept transcript is never seen.
  • ReadAlign_peOverlapMergeMap.cpp:31-39, 266-306: for the merged mates, once nW>0, peOverlapSEtoPE converts every transcript of every window back into a pair and rescores it from its exons with alignScore (Transcript_alignScore.cpp). The -999999999 and the extra mismatches are gone. The kept edge transcript is now the best one, and its unextended end is the soft clip.
  • The window that scored above 0 is a by-chance alignment of the merged read (in the fuzz case a 57M26I86M, AS:13). Mapping the merged read single-end with permissive filters shows both: the 1S168M at AS:0 nM:171 and the AS:13 one.
  • The clip length is wherever the seed started (8S above for a 3-base overhang), since no extension is applied.

Change

  • src/align/stitch.rs: EXTEND_TO_END_KILL is STAR's -999999999. When the scorer's new keep_end_to_end_kills is set, a killed transcript is not dropped: its score and mismatches get the kill values, and wt_final_score floors at 0 when the length penalty is on, as STAR does.
  • src/align/read_align.rs: align_read_inner sets keep_end_to_end_kills only for the merged read (mapped_filter == false), and skips the score range there, so every transcript is split and rescored. A read whose transcripts all score 0 now has no window ("other", trBest = trInit), as in STAR; for the merged read this means the pair is not replaced.
  • Other reads still drop the killed transcript. STAR keeps it there too but never reports it (it can only matter with a best score of 1, or for chimeric detection). Noted as a shortcut: on the field.

Measured against STAR 2.7.11b

Fuzz: 4000 pairs per seed, 2 random 20 kb chromosomes, fragments 130-175 bp (mates 100 bp, overlap 25-70), 70% at a chromosome start or end with a 0-4 base random overhang, 0-2 mismatches, 30% interior. --alignEndsType EndToEnd --peOverlapNbasesMin 10 --outFilterScoreMinOverLread 0.3 --outFilterMatchNminOverLread 0.3 --outSAMunmapped Within. Reads whose SAM fields 1-9 and tags differ from STAR:

seed #368 (base) this PR this PR + #367
1 2 2 0
2 2 1 0
3 6 4 0
4 5 4 0
5 1 1 0
6 3 2 0
7 3 2 0
8 3 2 0
total (32000 pairs) 25 18 0

The 7 fixed here are all clipped edge pairs like the ones above. The 18 left are half-mapped pairs (multimapper NH, the unmapped mate's record), which #367 fixes; merged together, all 32000 match.

Same fuzz, other modes (seeds 1-2): --alignEndsType Local 0 differences before and after. Extend5pOfReads12: 815 → 221 and 791 → 190; see Not checked.

Tests

tests/end_to_end_pairs.rs: overlapping_mates_over_a_chromosome_edge_are_clipped_like_star. Overhangs of 1 and 2 bases, fragments 130 and 169, over the start and the end of c1. A diverged copy of the merged read's edge end on c2 gives the window above 0, so the case does not depend on chance. Expected records are STAR 2.7.11b's output on the same genome and reads (1S99M/2S98M at 1, 99M1S/98M2S at 19902/19903, AS 197/196, nM 0). It fails on #368 (no records) and passes here. run gains run_on to take a FASTA. cargo test --release, cargo clippy --all-targets -- -D warnings and cargo fmt --check are clean on macOS.

Not checked

  • Extend5pOfReads12 with --peOverlapNbasesMin still differs on ~200 of 4000 fuzz pairs. Cause seen: STAR's peOverlapSEtoPE also converts a merged transcript that lies in one mate only (AS:7, then the pair fails mappedFilter), while convert_merged_transcript_to_pe drops it, so rustar keeps the separate-mate result. Separate issue, not fixed here.
  • Chimeric detection does not see the kept score-0 transcripts (STAR's trAll has them).
  • The parity suite (yeast, chr21) was not run: the change is limited to the merged read of --peOverlapNbasesMin, plus the all-score-0 rule.

🤖 Generated with Claude Code

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 <noreply@anthropic.com>

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant