Skip to content

alignMatesGapMax limits the gap between mates, not the fragment span - #373

Open
BenjaminDEMAILLE wants to merge 1 commit into
perf/hotpathsfrom
fix/mates-gap-max
Open

BenjaminDEMAILLE wants to merge 1 commit into
perf/hotpathsfrom
fix/mates-gap-max

Conversation

@BenjaminDEMAILLE

Copy link
Copy Markdown
Contributor

Part of #354 (--alignMatesGapMax limits the gap between mates; --outMultimapperOrder Random checked).

Base: perf/hotpaths. Merges cleanly with #367 (fix/half-mapped-multimap) and #368 (fix/end-to-end-pairs), in that order.

Problem

try_pair_transcripts rejected a pair when the fragment span (right.genome_end - left.genome_start) was above --alignMatesGapMax. STAR limits the gap between the mates. So under --alignMatesGapMax 300, a pair with 100-base mates and a 149-base gap (span 349) was reported unmapped "too short". STAR maps it:

STAR    r0  163  chr1  29360  255  100M  =  29609   349  AS:i:196 NH:i:1
STAR    r0   83  chr1  29609  255  100M  =  29360  -349  AS:i:196 NH:i:1
rustar  r0  141  *  0  0  *  *  0  0  AS:i:196 NH:i:0      (unmapped)
rustar  r0   77  *  0  0  *  *  0  0  AS:i:196 NH:i:0

check_proper_pair also used the span against --alignMatesGapMax for the 0x2 flag.

What STAR does (2.7.11b source)

  • The only per-pair use of alignMatesGapMax is in stitchAlignToTranscript.cpp:354: stitching mate B onto mate A fails when gBstart > end of mate A's last exon + alignMatesGapMax. That is the gap, not the span. It is skipped when the value is 0.
  • Its other uses are window sizing (Genome_genomeLoad.cpp:382-405: winBinNbits, winFlankNbins) and the old chimeric detection (ReadAlign_chimericDetectionOld.cpp:299). rustar already ports both (Parameters::redefine_window_params, chimeric/detect.rs); I checked them and changed nothing.
  • The 0x2 flag (ReadAlign_outputTranscriptSAM.cpp:79-83) is geometric: it fails only for mates that protrude past each other. With the default --alignEndsProtrude 0 (the only value rustar takes) stitching already rejects those pairs, so every pair is proper. alignMatesGapMax plays no part.

Change

  • try_pair_transcripts (src/align/read_align.rs): reject when right.genome_start > left.genome_end + alignMatesGapMax, as STAR's stitcher does. The stitcher's own check (stitch.rs, forward-gap check) was already right.
  • check_proper_pair returns true. I kept the function and its call sites so the branch merges cleanly with Half-mapped pairs go through the same selection and filters as pairs #367, which rewrites those call sites.

--outMultimapperOrder Random: only DIVERGENCE.md §1.1, nothing changed

STAR's multMapSelect (ReadAlign_multMapSelect.cpp:58-87) moves the best-scoring alignments to the top, shuffles them and the rest separately with the per-thread mt19937, and marks trMult[0] primary. rustar keeps the same candidates (same score range), puts the best-scoring ones first, shuffles them with a per-read seed, and marks the first one primary. The remaining difference is the random stream itself (§1.1), plus rustar leaves the non-best alignments sorted rather than shuffled. Both are order only.

Evidence, --outMultimapperOrder Random --runThreadN 1, yeast 10k, this branch:

SE PE
reads identical record for record 8,565 8,115
reads that differ only in order / 0x100 / HI 362 276
reads whose set of alignments differs 0 0
reads where the primary is a best-AS alignment (STAR / rustar) 8,927 / 8,927 8,391 / 8,391
reads with one best-AS alignment: same primary 8,299 / 8,299 7,886 / 7,886

On the 4,000-pair simulation below with #367 and #368 merged, 551 reads differ under Random. Set the unmapped-mate records aside and 405 differ, all only in order and primary, with the same set of alignments; the primary is a best-AS alignment in all 4,000 reads in both tools. The other differences are in the unmapped-mate record: with #367, a half-mapped multimapper writes the unmapped record for the mate opposite the primary, so a different random primary also changes that record. That is still the primary choice.

Measured against STAR 2.7.11b

Each tool with its own index, --runThreadN 1. "Reads differ" counts reads whose records (fields 2-9 and tags) are not all identical.

Simulation: 4,000 pairs on a random 300 kb chr1 with two 2 kb repeats, 100-base mates, fragments 80-900, 1% mismatches. Options --alignMatesGapMax 300 --outFilterScoreMinOverLread 0.3 --outFilterMatchNminOverLread 0.3 --outSAMunmapped Within.

Branch Simulation Yeast 10k PE, --alignMatesGapMax 300
perf/hotpaths 1,689 1,571
this PR 700 0
perf/hotpaths + #367 + #368 989 1,571
this PR + #367 + #368 0 0

The 700 left on this branch alone are half-mapped multimappers, which #367 fixes.

CONTRIBUTING.md check, yeast 10k, defaults: SE 9,724 records and PE 18,600 records identical to STAR, SJ.out.tab identical (test/compare_sam.py: 8,927/8,927 and 8,391/8,391).

Tests

test_mates_gap_max_limits_gap_not_span (replaces test_check_proper_pair_too_far): mates at 1000-1100 and 1300-1400 (gap 200, span 400) pair under a limit of 300 and are proper; under 199 they do not pair. It fails with the old span check. cargo test --release, cargo clippy --all-targets -- -D warnings and cargo fmt --check are clean on macOS.

Not checked

  • parity_suite.sh (chr21, fusion, STARsolo): the default alignMatesGapMax is 0 and both changed checks are skipped at 0.
  • --outMultimapperOrder Random with --outSAMmultNmax or --outSAMprimaryFlag AllBestScore.
  • The reviewer's own 4,000-pair set: I used my own simulation with the same options.

🤖 Generated with Claude Code

STAR rejects a pair only when mate B starts more than alignMatesGapMax
after mate A's end (stitchAlignToTranscript.cpp:354). rustar compared the
whole fragment span, so pairs with a small gap but a span above the limit
were reported unmapped. The proper-pair flag no longer depends on
alignMatesGapMax either: STAR's 0x2 test is geometric.

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