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

CHANGELOG.md

@@ -13,6 +13,65 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 ### Features
 
+- **Annotated-junction stitching, `--alignEndsType`, and the in-recursion
+  genomic-length penalty** — closes items #4b, #7, #8 and part of #9 of
+  `STAR-RS-COMPARISON.md` §7.2.
+
+  - `SpliceJunctionDb` gains STAR's `sjdb*` table and a `binarySearch2`-style
+    `find()`, so an annotated junction can be addressed by index (`sjA`) and
+    not merely tested for existence. No index format change: `sjdbInfo.txt`
+    already carried motif, shiftLeft, shiftRight and strand.
+  - Seeds derived from the sjdb genome insert carry their `sjA` tag through the
+    window into the stitcher, which takes STAR's annotated-junction shortcut
+    when two pieces turn out to be the two flanks of one junction.
+  - Non-canonical annotated junctions are snapped onto their annotated
+    coordinates instead of being left wherever the leftmost-flush scan put
+    them, and the stored motif overrides the one re-derived from the genome.
+  - `--alignSoftClipAtReferenceEnds`, `--alignInsertionFlush Right` and
+    `--outFilterMismatchNoverReadLmax` are implemented; `--alignEndsProtrude`,
+    `--alignTranscriptsPerReadNmax`, `--seedNoneLociPerWindow` and
+    `--seedSplitMin` are accepted. (`--alignEndsType` itself came from #145.)
+  - The genomic-length penalty is applied inside the stitch recursion, where it
+    can affect which transcripts survive, rather than only at finalisation.
+  - A one-base intron is a junction. `GTF_transcriptGeneSJ.cpp:123` skips an
+    exon pair only when `exon2.start <= exon1.end + 1`; rustar-aligner also
+    skipped `exon2.start == exon1.end + 2`, so three yeast junctions were
+    missing from the index and reads crossing them were written with `1D`
+    rather than `1N`.
+  - An annotated gap is never a deletion. `--alignIntronMin` decides between
+    `N` and `D` only for unannotated gaps, so the sjdb is consulted before the
+    length test at both CIGAR decision sites.
+  - The primary transcript is chosen by score, then by genomic span, which is
+    the comparison `ReadAlign_stitchPieces.cpp:340` makes
+    (`maxScore`, then `gLength`).
+  - Window alignments split out of different annotated junctions are no longer
+    treated as duplicates of one another. STAR's `assignAlignToWindow` compares
+    `sjA` as well as fragment and diagonal; keying the pre-stitch dedup on the
+    diagonal alone unpaired the two halves of a terminal micro-exon, which then
+    came out soft-clipped instead of spliced.
+  - A genome spacer compares as the larger byte in the seed search, as it does
+    in STAR's `compareSeqToGenome`. It was hardcoded to the opposite, which
+    made the binary search in `max_mappable_length` read the suffix array as
+    unsorted wherever the sjdb inserts sit and drop the half holding the real
+    maximum, returning a short MMP.
+
+  Measured against native STAR 2.7.11b, raw exact records (FLAG, RNAME, POS,
+  MAPQ, CIGAR, NH, AS, NM), with the index rebuilt from this branch:
+
+  | 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 |
+
+  No read regressed. Every remaining difference on all three tiers is a read
+  where both tools found alignments of identical score at different loci.
+  Output against an unannotated index is unchanged.
+
+  Reproducing this needs an index built by this branch: the one-base-intron
+  fix changes `genomeGenerate` output (364 yeast junctions, not 361), and a
+  cached index built before it silently reproduces the old numbers.
+
 - **STARsolo single-cell quantification (`--soloType`)** — the 10x
   Chromium / plate-based count-matrix pipeline, ported from STAR and
   verified against real STARsolo (#90).

DIVERGENCE.md

@@ -24,9 +24,11 @@ Divergences are grouped by kind:
 
 **Why.** Determinism and thread-count invariance: the same read produces the same primary regardless of `--runThreadN`. STAR's exact mt19937 stream cannot be reproduced under per-read parallelism, and matching it would forfeit reproducibility.
 
-**Impact.** With the default `--outMultimapperOrder Old_2.4`, **no RNG is consulted at all** — the primary is the deterministic best alignment (max score → smaller genomic length → earliest discovered), which is STAR-faithful. The divergence is observable only under `--outMultimapperOrder Random`, and only in *which* equal-scoring locus is marked primary — never in the set of reported alignments.
+**Impact.** With the default `--outMultimapperOrder Old_2.4`, **no RNG is consulted at all**: the primary is the deterministic best alignment. The first two keys are STAR's own (`ReadAlign_stitchPieces.cpp:340` compares `maxScore`, then `gLength`). Where those two tie, STAR takes the earliest window in its iteration order and rustar-aligner takes the smallest genomic position; that last key is a divergence, and it is the one described below. The RNG divergence proper is observable only under `--outMultimapperOrder Random`, and in both cases only in *which* equal-scoring locus is marked primary, never in the set of reported alignments.
 
-This is the reason faithfulness is reported **tie-adjusted**. On the 10k yeast benchmark, 299 SE and 475 PE primary-selection differences are all genuine ties: both tools find the identical alignment set, and differ only in which equal-scoring member is primary (from SA-iteration order or the RNG-seed difference above). Excluding those ties, SE is 99.815% and PE 99.883% exact.
+**On the residual ties.** STAR's window order is deterministic per read (anchor pieces in `PC` order, positions in suffix-array order) and rustar-aligner also builds windows per read, so reproducing it would not cost thread-invariance. It has been measured, not assumed: substituting seed discovery order for the positional key gained 21 reads on the annotated-junction tier and lost 142 on the 10k single-end tier, because this codebase's anchor iteration order is not STAR's `PC` order. Closing the gap means matching how MMP results are recorded, which has not been done. Until it is, the positional key stays, because it is total, cheap and thread-invariant.
+
+This is the reason faithfulness is reported **tie-adjusted**. On the 10k yeast benchmark, all 130 differing single-end records and 197 of the 201 differing paired-end mate records are genuine ties: both tools find the identical alignment set and differ only in which equal-scoring member is primary. On the annotated-junction tier all 42 remaining differences are of this kind, which is what its 100.000% tie-adjusted figure means. The handful that are not ties are in [§5](#5-known-residual-single-read-differences). Raw counts are reported alongside, never the tie-adjusted figure alone.
 
 **Source.** `src/rng.rs`, `src/align/read_align.rs` (`per_read_seed`, `shuffle_tied_prefix`), `src/params/mod.rs` (`MultimapperOrder`). STAR: `ReadAlign_multMapSelect.cpp`, `ReadAlignChunk` RNG seeding.
 
@@ -85,9 +87,14 @@ rustar-aligner uses an in-tree splitmix64 (`src/rng.rs`) rather than the `rand`
 
 These are **not** deliberate divergences — they are tracked residual diffs on the 10k yeast benchmark, kept here for completeness. Each is a single read; none is a systematic behaviour difference.
 
-- **1 SE CIGAR-only diff** — `ERR12389696.13573895`: both tools align to XV:218357, MAPQ 255, identical score (AS=133), but place a 1-base insertion differently (`100M1I45M4S` vs STAR's `108M1I37M4S`). The 71-base seed is found at a different position within a long homopolymer run (a seed-level tie); resolving it requires reproducing STAR's exact Lmapped chain path.
-- **1 STAR-only PE mate** — `ERR12389696.18919121`: an SA-level difference.
-- **1 rustar-aligner-only PE mate** — `ERR12389696.6302610`: a pre-existing false positive.
+- **1 PE pair where rustar-aligner scores lower**: `ERR12389696.11539725`, AS 224 against STAR's 235, emitting `31S95M459N24M` for mate 1 where STAR emits `22S10M468N94M459N24M`.
+
+  Traced end to end. Both tools build the same window and agree on every seed in it but one. STAR's window holds two aligns on a single diagonal, `r172 g5083545 L11` and `r174 g5083547 L9`; rustar-aligner 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 the first exon here starts at `g5083546`: the 11-base seed at `g5083545` is one base too far left, so only the 9-base seed at `g5083547` can open the mate, after which STAR's `extendAlign` walks it one base left to `g5083546`. Without that align rustar-aligner cannot assemble the two halves, and the best it reaches are 256 and 254 where the pair would score 266.
+
+  The missing align is not a dedup difference. STAR's `assignAlignToWindow` overlap test would itself collapse the 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 that means reproducing STAR's `PC` order and its Gsj splits, which is the seed-search side, not the stitcher. Ruled out along the way: the per-window transcript cap (raising it changes nothing), `qualitySplit` (the pair has no `N` and no base below Q3, so both aligns share one fragment), and the diagonal dedup key.
+- **1 PE pair where rustar-aligner scores higher**: `ERR12389696.4972950`, AS 260 against STAR's 248; mate 2 is spliced `1S33M72N50M186N65M1S` where STAR soft-clips 27 bases to `27S122M1S`. Both junctions are novel, not in the sjdb, so "higher-scoring" is all that is claimed here; neither alignment has been verified against the transcript. Recorded here rather than in [§2](#2-cases-where-rustar-aligner-outperforms-star) for that reason.
+
+Counts are from the 10k yeast SE and PE tiers with an index built by the same binary. The annotated-junction tier has no non-tie residual, and neither does the single-end tier.
 
 See `CLAUDE.md` ("Known Issues" / "PE Status") for the current status of these.
 

docs-old/dev/divergences.md

@@ -0,0 +1,58 @@
+# Deliberate divergences from STAR 2.7.11b
+
+rustar-aligner aims to be a faithful port: given the same inputs and flags, it should
+produce the same bytes as STAR 2.7.11b. This file records the places where it
+deliberately does not, and why.
+
+Every entry here is a case where STAR's behaviour is undefined, plainly wrong, or
+impossible to reproduce, and where reproducing it would mean shipping a known
+defect. Each is covered by a test that asserts the **correct** result, never
+STAR's.
+
+Divergences that arise from bugs on our side are not listed here. They are bugs,
+and they get fixed.
+
+---
+
+## Format
+
+Each entry gives: what STAR does, what rustar-aligner does instead, why, and the
+test that locks the behaviour in.
+
+---
+
+## D-01 · Non-canonical annotated junctions carry no strand
+
+**STAR.** Tracks a per-junction strand (`sjStr`) alongside the motif, and for an
+annotated junction takes that strand from `sjdbStrand`. SJ.out.tab column 4 then
+reports the annotated strand even when the motif is non-canonical.
+
+**rustar-aligner.** Has no per-junction strand on the working transcript, so
+column 4 reports `0` (undefined) for a non-canonical annotated junction where
+STAR would report `1` or `2`. Canonical junctions are unaffected: their strand
+is implied by the motif, which both derive identically.
+
+**Why.** Not a deliberate improvement, just not yet ported. Recorded here so the
+difference is not mistaken for a bug in the annotated-junction path, which is
+otherwise faithful. Adding `junction_strands` to `WorkingTranscript` and
+`Transcript` would close it.
+
+**Test.** None yet; this entry is the marker.
+
+---
+
+## Divergences considered and rejected
+
+Cases where STAR-rs, the other Rust port, made a different choice that
+rustar-aligner deliberately does not adopt.
+
+### Transcriptome-BAM primary flag
+
+STAR picks the primary alignment for the transcriptome BAM with a per-thread
+RNG, so the output depends on `--runThreadN`. STAR-rs resolves this by always
+taking the first alignment (`j == 0`).
+
+rustar-aligner already picks by a per-read seed, which is thread-count invariant
+and therefore fixes the same defect. The two ports resolve it differently and
+both diverge from STAR; there is no fidelity argument for switching to `j == 0`,
+so rustar-aligner keeps its own rule.

src/align/read_align.rs

@@ -299,6 +299,16 @@ pub fn align_read(
     };
 
     for (ci, cluster) in clusters.iter().enumerate() {
+        // STAR's `alignTranscriptsPerReadNmax` headroom break
+        // (`ReadAlign_stitchPieces.cpp`): stop collecting windows once one more
+        // window's worth of transcripts could overrun the per-read cap. The
+        // test is on the headroom, not on the running total, so the cap is
+        // never exceeded rather than merely noticed afterwards.
+        if transcripts.len() + params.align_transcripts_per_window_nmax
+            >= params.align_transcripts_per_read_nmax
+        {
+            break;
+        }
         let debug_name = if debug_read { read_name } else { "" };
         let cluster_transcripts = stitch_seeds_with_jdb_debug(
             cluster,
@@ -364,11 +374,16 @@ pub fn align_read(
         });
     }
 
-    // Deterministic primary tie-break (score, then a fixed positional order).
+    // Deterministic primary tie-break. STAR compares `maxScore`, then
+    // `gLength` — the alignment's genomic span — and leaves anything still
+    // tied to whichever window it reached first
+    // (`ReadAlign_stitchPieces.cpp:340`). The span is reproducible here; the
+    // window order is not, since windows are built per read in parallel, so
+    // the remaining keys are positional and fixed (see DIVERGENCE.md §1.1).
     transcripts.sort_by(|a, b| {
         b.score
             .cmp(&a.score)
-            .then_with(|| a.n_junction.cmp(&b.n_junction))
+            .then_with(|| (a.genome_end - a.genome_start).cmp(&(b.genome_end - b.genome_start)))
             .then_with(|| a.chr_idx.cmp(&b.chr_idx))
             .then_with(|| a.genome_start.cmp(&b.genome_start))
             .then_with(|| a.is_reverse.cmp(&b.is_reverse))
@@ -731,15 +746,20 @@ pub fn align_paired_read(
     // → Nstart=7, starts={0,43,...,129,...}). Using combined length creates a spurious
     // start at position 129 (between mates) that can produce anchors widening windows
     // beyond STAR's range, causing window overflow and eviction of valid 7M exon seeds.
-    let mut combined_seeds = Seed::find_seeds(
+    // The piece offsets are STAR's `splitR[0][ip]`: mate1 starts the concatenated
+    // read, mate2 starts one base past the spacer. Only the `flagDirMap`
+    // shortcut consults them.
+    let mut combined_seeds = Seed::find_seeds_at(
         &combined_read[..len1],
+        0,
         index,
         params.seed_map_min,
         params,
         debug_name,
     )?;
-    let mut m2_seeds = Seed::find_seeds(
+    let mut m2_seeds = Seed::find_seeds_at(
         &combined_read[len1 + 1..],
+        len1 + 1,
         index,
         params.seed_map_min,
         params,