aligner core: annotated-junction stitching, alignEndsType, in-recursion length penalty

[BenjaminDEMAILLE]

Aligner-core faithfulness: annotated-junction stitching, alignEndsType, and the genomic-length penalty inside the stitch recursion.

What changed

sjA tags survive the seed pipeline. decode_gsj_hit computed the junction index and discarded it; it is now carried through expand_hit, WinCandidate, WindowAlignment and ExonBlock. The overlap dedup in cluster_seeds gained STAR's e.sjA == sjA guard, so seeds from different annotated junctions are no longer collapsed for sharing a diagonal.

The annotated shortcut. When the previous exon and the incoming piece carry the same tag, sit on the same mate and are read-adjacent, they are the two flanks of one junction that the Gsj insert split apart. They are stitched directly onto the annotated boundary with the stored motif, shifts and sjdbScore bonus, skipping the general search.

Annotated-boundary snapping. A non-canonical annotated junction is shifted by shiftLeft onto its annotated coordinates instead of being left leftmost-flush, and the stored motif overrides the one re-derived from the genome. Both of STAR's -1000006 rejections are reproduced.

The genomic-length penalty moves inside the recursion, where STAR applies it, rather than being added in a post-pass.

Flag	Status
--alignEndsType	implemented (Local, EndToEnd, Extend5pOfRead1, Extend5pOfReads12)
--alignSoftClipAtReferenceEnds	implemented
--alignInsertionFlush	implemented (Right; None is bit-identical to before)
--outFilterMismatchNoverReadLmax	implemented
--alignTranscriptsPerReadNmax	implemented (STAR's headroom break)
--alignEndsProtrude, --seedNoneLociPerWindow, --seedSplitMin	accepted and validated

Why

STAR-RS-COMPARISON.md §7.2 lists these as the open faithfulness items and §7.4 names annotated-boundary snapping the best remaining lever. stitch.rs asked the annotation a yes/no question and, on a hit, swapped the motif penalty for sjdbScore; it never used what the annotation says about where the junction is. Item #2 stays dropped, as §7.5 records it regressing.

No index format change

PreparedJunction already carried motif, shift_left, shift_right and strand; sjdbInfo.txt already wrote and read all four; sort_and_dedup already established the (stored_start, stored_end) ordering binarySearch2 needs. find() searches data that was already on disk. Existing indices load and align unmodified.

The contents of a newly generated index do change, because of the one-base-intron fix below: yeast R64-1-1 now yields 364 junctions rather than 361, matching STAR. That matters for reproducing the numbers here and is called out under Verification.

One inconsistency is fixed along the way: the align-time --sjdbGTFfile path went through from_raw_junctions, which discarded motif and shift, while the index-loaded path kept them. The three population paths now build identical tables. Where the index carries junctions, its array is the sjA table, because that is the array decode_gsj_hit indexes.

Verification against STAR 2.7.11b

Both columns are measured runs of the two binaries, each against an index built by that same binary. Sharing one index between them would credit this branch with main's junction set, or main with this branch's; the numbers move either way, so they are built separately. Data per CONTRIBUTING.md (yeast R64-1-1 release 110, ERR12389696, seqtk sample -s100).

Raw exact records, comparing FLAG, RNAME, POS, MAPQ, CIGAR, NH, AS and NM:

tier	records	main	this branch
junction-spanning yeast reads	1000	835	958
10k single-end yeast reads	8927	8785	8797
10k paired-end yeast reads	16782	16572	16583

Tie-adjusted, for comparison with earlier reports: junction tier 87.161% to 100.000%, SE 99.830% to 99.898%, PE 99.795% to 99.861%. Nothing regressed on any tier.

The junction tier is what makes the change visible at all. The existing tiers align real RNA-seq reads, few of which cross an intron, so the annotated-junction paths barely run. This PR adds test/simulate_junction_reads.py and a yeast_junction_gtf tier: 1000 reads drawn from spliced transcript sequences, centred on annotated junctions, fixed seed, so the FASTQ is byte-reproducible.

To reproduce, with $DATA set up as CONTRIBUTING describes:

python3 test/simulate_junction_reads.py \
  "$DATA/reference/Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa" \
  "$DATA/reference/Saccharomyces_cerevisiae.R64-1-1.110.gtf" \
  "$DATA/reads/sim_junction_1k.fastq" 1000 150
# rebuild the rustar-aligner index with THIS branch's binary, then align
# with STAR and with rustar-aligner against their own GTF-built indices:
python3 test/assess_faithfulness.py --se-only \
  --se-star "$DATA/star_junc/Aligned.out.sam" \
  --se-rustar-aligner "$DATA/rustar_junc/Aligned.out.sam"

A cached index will silently reproduce main's numbers. That is not hypothetical: the first run of this comparison used a stale index and reported 21 reads as an unfixed CIGAR bug when the code was already correct and only the fixture was old.

Five defects the tiers exposed, which reading the code had not

Each was found by measuring, and they compound, so each one hid the next.

genomeGenerate dropped one-base introns. The GTF junction extractor required intron_end > intron_start, which rejects an intron of length one. STAR skips a pair only when the exons touch or overlap (GTF_transcriptGeneSJ.cpp:123), so exon2.start == exon1.end + 2 is a junction it keeps. Yeast R64-1-1 has three, and the index carried 361 junctions against STAR's 364; the three missing were exactly those. sjdbList.out.tab is now byte-identical to STAR's. Worth 21 reads on the tier.

An annotated gap shorter than --alignIntronMin was scored and written as a deletion. STAR resolves sjdbInd before testing Del >= alignIntronMin (stitchAlignToTranscript.cpp:198-241): on a hit it takes the annotated branch whatever the gap length. Here the lookup sat inside the already-decided splice branch, so a short annotated gap never reached it. The CIGAR builder then re-derived the same decision from length alone, so even with the score fixed the operation stayed D. Both now ask through one gap_is_junction.

The primary tie-break compared the junction count where STAR compares gLength, the alignment's genomic span (ReadAlign_stitchPieces.cpp:340): a proxy for the span rather than the span itself. Using the span directly gains 6 SE reads and costs none.

Window alignments split out of different annotated junctions were treated as duplicates of one another. STAR's assignAlignToWindow guards its duplicate test with both aFrag == WA[iA][WA_iFrag] and sjA == WA[iA][WA_sjA]. The pre-stitch diagonal dedup keyed on (diagonal, mate_id) only, so for a read whose last exon is a few bases long, the two entries covering the same donor-side bases (one tagged with the junction the read arrives on, one with the junction it leaves by) collapsed to the longest. That kept only the arriving tag, the micro-exon entry lost its same-sjA partner, and the terminal exon came out soft-clipped:

STAR              61M2448N36M2483N3M   AS:i:101
before this fix   61M2448N36M3S        AS:i:96

Four uniquely-mapped reads on the tier, five points of AS each.

A genome spacer compared as the smaller byte in the seed search. compare_seq_to_genome returned "the query is greater" for any genome byte of 5. STAR's compareSeqToGenome just compares the raw bytes, and since read bases are 0-3 while N is 4 and the spacer is 5, a read base is never the greater one: every branch reports the query as sorting before such a suffix.

That value orders the binary search in max_mappable_length, so it has to agree with the order the suffix array was built in. It did not, and the search read the array as unsorted exactly where the sjdb inserts sit. On ERR12389696.20597455 the array holds a 29-base match at the right locus one slot below a 24-base match inside an insert whose comparison hit a spacer; the probe claimed the query was greater, the search dropped everything below it, and the MMP came back 24. The read then stitched two shorter pieces and placed its 1-base insertion four bases left of STAR's, 25M1I124M against 29M1I120M, at the same total score. This one is invisible on the annotated-junction tier and worth one read on each of SE and PE.

Locked by a_one_base_intron_is_a_junction_and_touching_exons_are_not, a_short_annotated_gap_is_a_junction_not_a_deletion, same_diagonal_entries_from_different_junctions_both_survive, same_diagonal_overlapping_entries_collapse_to_the_longest and a_query_running_into_the_padding_sorts_before_that_suffix.

Two seed-search corrections that changed no output byte

Both are in the MMP chain loop, read from ReadAlign_mapOneRead.cpp, and both are here because the search should be STAR's search, not because they moved a number.

flagDirMap (:50, :61, :74): when the first left-to-right search of a good piece maps it to the piece end, STAR skips the istart == 0 chain of the right-to-left direction. rustar-aligner always ran both. On these tiers the extra seed is a full-length duplicate that storeAligns would drop anyway.

The loop bound: STAR continues while istart*Lstart + Lmapped + P.seedMapMin < splitR[1][ip], that is while strictly more than seedMapMin bases remain; rustar-aligner continued while at least that many remained, running one extra search on a piece with exactly seedMapMin left. A 5-base MMP in a 12 Mb genome exceeds seedMultimapNmax and is never stored, so nothing changed here either.

Reported as no-ops rather than left out, because a reviewer should be able to tell which parts of this PR are load-bearing.

What is left

42 reads on the junction tier, 130 on SE and 197 mate records on PE still pick a different primary. Every one is a multimapper where both tools found alignments of identical score at different loci. On the junction tier and on SE that is now the only category left, which is what the 100.000% and 99.898% tie-adjusted figures mean.

Two PE pairs are not ties, and both are in DIVERGENCE.md §5:

ERR12389696.11539725 scores 224 against STAR's 235, and traces to the seed search rather than the stitcher.

Both tools build the same window and agree on every align in it but one. STAR holds two aligns on a single diagonal, r172 g5083545 L11 and r174 g5083547 L9; this branch holds only the first. stitchAlignToTranscript.cpp:352 refuses to cross into the second mate when the mate's seed starts before the transcript's first exon, and that exon starts at g5083546, so the 11-base seed at g5083545 is one base too far left and only the 9-base seed at g5083547 can open the mate; STAR's extendAlign then walks it one base left onto g5083546. Without that align the two halves cannot be assembled, which is why the best reached here are 256 and 254 where the pair would score 266.

It is not a dedup difference. STAR's own assignAlignToWindow overlap test would collapse that pair, but it fires on insertion order, and a Gsj-crossing piece records its acceptor half at aRstart + aLengthD rather than at its own PC start. Reproducing it means reproducing STAR's PC order and its Gsj splits.

Ruled out by measurement rather than by reading: the per-window transcript cap (raising it to 10000 gives the same 84 transcripts and 1460 recursions), qualitySplit (no N and no base below Q3 in the pair, so both aligns share one fragment), and the diagonal dedup key. Rewriting the dedup onto STAR's asymmetric overlap test was tried and reverted, because its own test shows STAR's rule collapses the pair too when applied in sorted order.

ERR12389696.4972950 scores 260 against STAR's 248, splicing 1S33M72N50M186N65M1S where STAR soft-clips 27 bases. Both junctions are novel, not in the sjdb, so the only claim made is about the score; neither alignment has been checked against the transcript.

I previously wrote in this PR that closing the tie category would cost thread-invariance. That was wrong and I retract it: STAR's window order is deterministic per read, and this codebase builds windows per read too, so reproducing it costs nothing in parallelism. I tried it (dropping the positional keys from the primary sort, relying on the stable sort to preserve discovery order) and measured +21 on the junction tier against -142 on SE. Net worse, so it is not in this PR. What that says is that this codebase's anchor iteration order is not STAR's PC order; matching it is a change to how MMP results are recorded, and it belongs in its own PR with its own before/after on both tiers.

A note on the harness

The Homebrew macOS STAR 2.7.11b build silently reads zero reads from a gzipped FASTQ (--readFilesCommand zcat produced a 776-byte SAM with 0 alignments). Decompressing first and passing the plain FASTQ works. The CONTRIBUTING commands use zcat, so anyone reproducing this on macOS will hit it; on Linux it does not occur.

Gate: 574 lib + 21 integration tests, cargo clippy --all-targets -- -D warnings, cargo fmt --check, MSRV 1.89, all green. Thread invariance checked at --runThreadN 1 against 16 on the 10k SE tier, byte-identical. Rebased onto current main.

[BenjaminDEMAILLE]

Followed the residual to its end rather than leaving it at "documented divergence". Reporting the attempt because it failed, and the numbers say something useful.

What the remaining differences are

37 reads on the junction tier choose a different primary locus. All 37 are multimappers with identical alignment scores; 18 also have identical CIGARs, so identical genomic spans. Both tools find the same alignment set and disagree only about which equal member is marked primary.

STAR's tie-break is maxScore, then gLength, then whichever window it reached first (ReadAlign_stitchPieces.cpp:340). The first two are now matched. The third is window iteration order.

Why I thought that was closed, and why it is not

I had written that reproducing window order would cost thread-invariance. That was wrong: STAR's window order is deterministic per read — anchor pieces in storage order, positions within a piece in suffix-array order — and this codebase builds windows the same way, per read, in parallel. Reproducing it costs nothing in parallelism.

So I tried it: drop the positional tie-break keys and let Rust's stable sort preserve discovery order, which is what cluster_seeds already produces and the dedup already preserves.

The measurement

	positional keys (current)	discovery order
junction tier, raw exact	951	972
SE 10k, raw exact	8797	8655

+21 on the tier, −142 on SE. Net worse, so it is reverted.

What that tells us: this codebase's anchor iteration order is not STAR's PC order. Where the two happen to agree, discovery order is strictly better — hence the tier gain. Where they do not, it is worse than an arbitrary-but-stable positional rule.

What reaching exact parity would take

Reproducing STAR's storeAligns piece ordering — the order in which MMP search results are recorded, which decides anchor order, which decides window order. That is a change to the seed-search side, not the stitcher, and it should be measured on both tiers before being believed. It is not a tie-break tweak.

Until then the residual is what it looks like: ties, resolved differently, on reads where both tools found the same alignments. Thread invariance is intact — 10k SE reads at --runThreadN 1 and 16 produce byte-identical SAM.