diff --git a/CHANGELOG.md b/CHANGELOG.md
index b6e4893..6605b33 100644
--- a/CHANGELOG.md
+++ b/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).
diff --git a/DIVERGENCE.md b/DIVERGENCE.md
index bd957ed..47e3082 100644
--- a/DIVERGENCE.md
+++ b/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.
 
diff --git a/docs-old/dev/divergences.md b/docs-old/dev/divergences.md
new file mode 100644
index 0000000..6127bd6
--- /dev/null
+++ b/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.
diff --git a/src/align/read_align.rs b/src/align/read_align.rs
index b4d2f20..5dd3caf 100644
--- a/src/align/read_align.rs
+++ b/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,
diff --git a/src/align/score.rs b/src/align/score.rs
index 7a34f55..b906b4d 100644
--- a/src/align/score.rs
+++ b/src/align/score.rs
@@ -52,6 +52,9 @@ pub struct AlignmentScorer {
     /// Read-end extension policy (alignEndsType). `ext[iMate][iEnd]==true` forces
     /// full end-to-end extension (no terminal soft-clip) of that mate/end.
     pub align_ends_type: crate::params::AlignEndsType,
+    pub flush_right: bool,
+    pub soft_clip_at_reference_ends: bool,
+    pub p_mm_max_read: f64,
 }
 
 impl AlignmentScorer {
@@ -80,6 +83,9 @@ impl AlignmentScorer {
             align_spliced_mate_map_lmin_over_lmate: 0.66,
             out_filter_score_min_over_lread: 0.66,
             align_ends_type: crate::params::AlignEndsType::default(),
+            flush_right: false,
+            soft_clip_at_reference_ends: true,
+            p_mm_max_read: 1.0,
         }
     }
 
@@ -120,6 +126,9 @@ impl AlignmentScorer {
             out_filter_score_min_over_lread: params.out_filter_score_min_over_lread,
             // Parsed+validated in Parameters::validate; default to Local if unset.
             align_ends_type: params.align_ends_type.parse().unwrap_or_default(),
+            flush_right: params.align_insertion_flush == "Right",
+            soft_clip_at_reference_ends: params.align_soft_clip_at_reference_ends != "No",
+            p_mm_max_read: params.out_filter_mismatch_nover_read_lmax,
         }
     }
 
@@ -132,6 +141,17 @@ impl AlignmentScorer {
         ((genomic_span as f64).log2() * self.score_genomic_length_log2_scale - 0.5).ceil() as i32
     }
 
+    /// STAR's `outFilterMismatchNmaxTotal`: the tightest of the absolute cap,
+    /// a fraction of the mapped length, and a fraction of the read length.
+    ///
+    /// `mapped_len` is the alignment's span in the read; `read_len` is the
+    /// whole read (both mates summed, for a pair).
+    pub fn mismatch_nmax_total(&self, mapped_len: usize, read_len: usize) -> u32 {
+        let by_mapped = (self.p_mm_max * mapped_len as f64) as u32;
+        let by_read = (self.p_mm_max_read * read_len as f64) as u32;
+        self.n_mm_max.min(by_mapped).min(by_read)
+    }
+
     /// Annotated junctions score `sjdb_score`; unannotated junctions score `motif_score`.
     pub fn score_annotated_junction(&self, motif_score: i32, annotated: bool) -> i32 {
         if annotated {
@@ -693,6 +713,9 @@ mod tests {
             align_spliced_mate_map_lmin_over_lmate: 0.66,
             out_filter_score_min_over_lread: 0.66,
             align_ends_type: crate::params::AlignEndsType::default(),
+            flush_right: false,
+            soft_clip_at_reference_ends: true,
+            p_mm_max_read: 1.0,
         };
 
         // Intron from position 2, length 12 (spans positions 2-13 inclusive)
@@ -738,6 +761,9 @@ mod tests {
             align_spliced_mate_map_lmin_over_lmate: 0.66,
             out_filter_score_min_over_lread: 0.66,
             align_ends_type: crate::params::AlignEndsType::default(),
+            flush_right: false,
+            soft_clip_at_reference_ends: true,
+            p_mm_max_read: 1.0,
         };
 
         let motif = scorer.detect_splice_motif(2, 12, &genome);
@@ -782,6 +808,9 @@ mod tests {
             align_spliced_mate_map_lmin_over_lmate: 0.66,
             out_filter_score_min_over_lread: 0.66,
             align_ends_type: crate::params::AlignEndsType::default(),
+            flush_right: false,
+            soft_clip_at_reference_ends: true,
+            p_mm_max_read: 1.0,
         };
 
         let motif = scorer.detect_splice_motif(2, 12, &genome);
@@ -824,6 +853,9 @@ mod tests {
             align_spliced_mate_map_lmin_over_lmate: 0.66,
             out_filter_score_min_over_lread: 0.66,
             align_ends_type: crate::params::AlignEndsType::default(),
+            flush_right: false,
+            soft_clip_at_reference_ends: true,
+            p_mm_max_read: 1.0,
         };
 
         let motif = scorer.detect_splice_motif(2, 12, &genome);
@@ -859,6 +891,9 @@ mod tests {
             align_spliced_mate_map_lmin_over_lmate: 0.66,
             out_filter_score_min_over_lread: 0.66,
             align_ends_type: crate::params::AlignEndsType::default(),
+            flush_right: false,
+            soft_clip_at_reference_ends: true,
+            p_mm_max_read: 1.0,
         };
 
         let (score, gap_type) = scorer.score_gap(0, 5, 0, &genome);
@@ -892,6 +927,9 @@ mod tests {
             align_spliced_mate_map_lmin_over_lmate: 0.66,
             out_filter_score_min_over_lread: 0.66,
             align_ends_type: crate::params::AlignEndsType::default(),
+            flush_right: false,
+            soft_clip_at_reference_ends: true,
+            p_mm_max_read: 1.0,
         };
 
         // Small gap (< align_intron_min) is deletion
@@ -933,6 +971,9 @@ mod tests {
             align_spliced_mate_map_lmin_over_lmate: 0.66,
             out_filter_score_min_over_lread: 0.66,
             align_ends_type: crate::params::AlignEndsType::default(),
+            flush_right: false,
+            soft_clip_at_reference_ends: true,
+            p_mm_max_read: 1.0,
         };
 
         // Gap starting at position 2 (GT), length 26 (>= 21) is splice junction
@@ -972,6 +1013,9 @@ mod tests {
             align_spliced_mate_map_lmin_over_lmate: 0.66,
             out_filter_score_min_over_lread: 0.66,
             align_ends_type: crate::params::AlignEndsType::default(),
+            flush_right: false,
+            soft_clip_at_reference_ends: true,
+            p_mm_max_read: 1.0,
         };
 
         let annotated_score = scorer.score_annotated_junction(0, true);
@@ -1012,6 +1056,9 @@ mod tests {
             align_spliced_mate_map_lmin_over_lmate: 0.66,
             out_filter_score_min_over_lread: 0.66,
             align_ends_type: crate::params::AlignEndsType::default(),
+            flush_right: false,
+            soft_clip_at_reference_ends: true,
+            p_mm_max_read: 1.0,
         };
 
         // CT-AC motif: (1,3,0,1) — reverse complement of GT-AG
@@ -1120,6 +1167,9 @@ mod tests {
             align_spliced_mate_map_lmin_over_lmate: 0.66,
             out_filter_score_min_over_lread: 0.66,
             align_ends_type: crate::params::AlignEndsType::default(),
+            flush_right: false,
+            soft_clip_at_reference_ends: true,
+            p_mm_max_read: 1.0,
         };
 
         // Gap of exactly 589824 starting at position 100 should be splice junction
@@ -1159,6 +1209,24 @@ mod tests {
         assert_eq!(SpliceMotif::GtAt.filter_category(), 3);
     }
 
+    #[test]
+    fn mismatch_nmax_total_takes_the_tightest_of_three_caps() {
+        let mut s = AlignmentScorer::from_params_minimal();
+        s.n_mm_max = 10;
+        s.p_mm_max = 0.3; // fraction of MAPPED length
+        s.p_mm_max_read = 1.0; // fraction of READ length
+
+        // Mapped fraction binds: 0.3 * 20 = 6.
+        assert_eq!(s.mismatch_nmax_total(20, 100), 6);
+        // Absolute cap binds: 0.3 * 100 = 30, capped at 10.
+        assert_eq!(s.mismatch_nmax_total(100, 100), 10);
+
+        // The read-length cap is the term that was previously missing: with a
+        // 100-base alignment it now binds before either of the other two.
+        s.p_mm_max_read = 0.02;
+        assert_eq!(s.mismatch_nmax_total(100, 100), 2);
+    }
+
     #[test]
     fn test_filter_category_from_encoded() {
         assert_eq!(SpliceMotif::filter_category_from_encoded(0), 0); // non-canonical
@@ -1200,6 +1268,9 @@ mod tests {
             align_spliced_mate_map_lmin_over_lmate: 0.66,
             out_filter_score_min_over_lread: 0.66,
             align_ends_type: crate::params::AlignEndsType::default(),
+            flush_right: false,
+            soft_clip_at_reference_ends: true,
+            p_mm_max_read: 1.0,
         };
 
         // Gap of 1001 (> 1000 max) should be deletion, not splice junction
@@ -1252,6 +1323,9 @@ mod tests {
             align_spliced_mate_map_lmin_over_lmate: 0.66,
             out_filter_score_min_over_lread: 0.66,
             align_ends_type: crate::params::AlignEndsType::default(),
+            flush_right: false,
+            soft_clip_at_reference_ends: true,
+            p_mm_max_read: 1.0,
         }
     }
 
diff --git a/src/align/seed.rs b/src/align/seed.rs
index 7c6deb3..572cf5e 100644
--- a/src/align/seed.rs
+++ b/src/align/seed.rs
@@ -48,6 +48,25 @@ impl Seed {
         min_seed_length: usize,
         params: &Parameters,
         debug_name: &str,
+    ) -> Result<Vec<Seed>, Error> {
+        Self::find_seeds_at(read_seq, 0, index, min_seed_length, params, debug_name)
+    }
+
+    /// [`find_seeds`](Self::find_seeds) with the piece's offset inside STAR's
+    /// concatenated `Read1[0]` buffer.
+    ///
+    /// STAR searches one *good piece* at a time (`splitR[0][ip]` = the piece's
+    /// start in the concatenated read, `splitR[1][ip]` = its length). Only the
+    /// `flagDirMap` shortcut in `ReadAlign_mapOneRead.cpp:74` reads the start
+    /// offset, so every other computation here is piece-local and `piece_start`
+    /// is threaded through purely to reproduce that condition.
+    pub fn find_seeds_at(
+        read_seq: &[u8],
+        piece_start: usize,
+        index: &GenomeIndex,
+        min_seed_length: usize,
+        params: &Parameters,
+        debug_name: &str,
     ) -> Result<Vec<Seed>, Error> {
         let mut seeds = Vec::new();
         let read_len = read_seq.len();
@@ -60,13 +79,15 @@ impl Seed {
         //     while(istart*Lstart + Lmapped + seedMapMin < readLen) { ... Lmapped += L; }
 
         // Search L→R (forward direction on read): sparse chain search
-        search_direction_sparse(
+        let flag_dir_map = search_direction_sparse(
             read_seq,
             read_len,
+            piece_start,
             index,
             min_seed_length,
             params,
             false,
+            true,
             debug_name,
             &mut seeds,
         );
@@ -76,15 +97,22 @@ impl Seed {
             return Ok(seeds);
         }
 
-        // Search R→L (reverse direction on read): sparse chain search on RC read
+        // Search R→L (reverse direction on read): sparse chain search on RC read.
+        // `flag_dir_map` is STAR's `flagDirMap` (`ReadAlign_mapOneRead.cpp:50`,
+        // `:61`, `:74`): when the very first L→R search already mapped the piece
+        // to its end, the `istart == 0` chain in the reverse direction would
+        // re-derive the same maximal prefix, so STAR skips it. Chains from
+        // `istart > 0` still run.
         let rc_read = reverse_complement_read(read_seq);
         search_direction_sparse(
             &rc_read,
             read_len,
+            piece_start,
             index,
             min_seed_length,
             params,
             true,
+            flag_dir_map,
             debug_name,
             &mut seeds,
         );
@@ -232,14 +260,19 @@ struct MmpResult {
 fn search_direction_sparse(
     read_seq: &[u8],
     original_read_len: usize,
+    piece_start: usize,
     index: &GenomeIndex,
     min_seed_length: usize,
     params: &Parameters,
     is_rc: bool,
+    run_istart0: bool,
     debug_name: &str,
     seeds: &mut Vec<Seed>,
-) {
+) -> bool {
     let read_len = read_seq.len();
+    // STAR's `flagDirMap`, returned to the caller: stays true unless the first
+    // L→R search maps the piece all the way to its end.
+    let mut flag_dir_map = true;
 
     // STAR (ReadAlign_mapOneRead.cpp lines 41-42):
     //   seedSearchStartLmax = min(P.seedSearchStartLmax, seedSearchStartLmaxOverLread*(Lread-1))
@@ -263,6 +296,12 @@ fn search_direction_sparse(
     let lstart = read_len / nstart; // STAR: Lstart = (splitR[1]-splitR[0]) / Nstart
 
     for istart in 0..nstart {
+        // STAR: `if (flagDirMap || istart>0)` — the reverse direction skips its
+        // istart == 0 chain when the forward direction already mapped the piece
+        // to its end.
+        if istart == 0 && !run_istart0 {
+            continue;
+        }
         let start_pos = (istart * lstart).min(read_len);
         let mut pos = start_pos;
 
@@ -273,12 +312,13 @@ fn search_direction_sparse(
             if pos >= read_len {
                 break;
             }
-            // Stop if remaining bases < seedMapMin (matches STAR's while condition:
-            // istart*Lstart + Lmapped + P.seedMapMin < splitR[1][ip]).
-            // STAR chains continue until only seedMapMin (5) bases remain, NOT
-            // seedSearchStartLmax (50). This allows chains to reach terminal small
-            // exons (e.g. 9M after intron) near the read end.
-            if read_len - pos < min_seed_length {
+            // STAR's while condition, verbatim:
+            //   istart*Lstart + Lmapped + P.seedMapMin < splitR[1][ip]
+            // i.e. keep searching while *more than* seedMapMin bases remain.
+            // rustar-aligner used `remaining < min_seed_length` here, which ran
+            // one extra search when exactly seedMapMin bases were left and
+            // pushed a seed STAR never stores.
+            if read_len - pos <= params.seed_map_min {
                 break;
             }
 
@@ -297,6 +337,18 @@ fn search_direction_sparse(
                 );
             }
 
+            // STAR (`ReadAlign_mapOneRead.cpp:74`): on the very first forward
+            // search of the piece, a match that reaches the piece end means the
+            // reverse direction has nothing new to find from istart == 0.
+            // The comparison is STAR's own — `Shift + L == splitR[1][ip]`, with
+            // `Shift == splitR[0][ip]` at that point — so for a mate that does
+            // not start at offset 0 it is the mate's *global* start that enters
+            // the sum, and the shortcut effectively never fires there.
+            if !is_rc && istart == 0 && pos == start_pos && piece_start + result.advance == read_len
+            {
+                flag_dir_map = false;
+            }
+
             if let Some(mut seed) = result.seed {
                 // Apply seedSearchLmax cap
                 if params.seed_search_lmax > 0 && seed.length > params.seed_search_lmax {
@@ -313,14 +365,16 @@ fn search_direction_sparse(
                 seeds.push(seed);
 
                 if seeds.len() >= params.seed_per_read_nmax {
-                    return;
+                    return flag_dir_map;
                 }
             }
 
             pos += result.advance; // Always advance by MMP length (matches STAR)
-            // Remaining-length check at loop top: stop when < seedMapMin bases remain
+            // Remaining-length check at loop top: stop when <= seedMapMin bases remain
         }
     }
+
+    flag_dir_map
 }
 
 /// Find a seed starting at a specific position in the read.
@@ -505,10 +559,9 @@ fn compare_seq_to_genome(
                 let genome_chunk = &genome_slice_all[genome_start + i..genome_start + simd_end];
                 if let Some(off) = crate::align::simd_scan::find_stop(read_chunk, genome_chunk) {
                     let genome_base = genome_chunk[off];
-                    if genome_base >= 5 {
-                        return (i + off, true);
-                    }
                     let read_base = read_chunk[off];
+                    // A spacer (5) is larger than any read base, so the query
+                    // sorts *before* this suffix. See the scalar branch below.
                     return (i + off, read_base > genome_base);
                 }
                 match_len = simd_end;
@@ -521,19 +574,27 @@ fn compare_seq_to_genome(
         let genome_idx = genome_start + i;
 
         if genome_idx >= index.genome.sequence.len() {
-            // Past end of genome array — treat like padding (STAR: comp_res > 0)
-            return (match_len, true);
+            // Past the end of the genome array. Treated as a spacer, so the
+            // query sorts before this suffix, same as the branch below.
+            return (match_len, false);
         }
 
         let genome_base = index.genome.sequence.base(genome_idx);
-
-        if genome_base >= 5 {
-            // Padding character — STAR returns comp_res > 0 (read > genome)
-            return (match_len, true);
-        }
-
         let read_base = read_seq[read_pos + i];
 
+        // The comparison result is what orders the binary search in
+        // `max_mappable_length`, so it has to agree with the order the suffix
+        // array was built in. STAR compares the raw bytes
+        // (`compareSeqToGenome`, `SuffixArrayFuns.cpp`): read bases are 0-3,
+        // `N` is 4 and the spacer is 5, so a read base is never greater than
+        // either of those and the query always sorts *before* such a suffix.
+        // `comp_res` was hardcoded to `true` for the spacer, the opposite.
+        //
+        // It matters where the sjdb inserts sit in the suffix array. A probe
+        // landing on a spacer inside an insert reported "query is greater",
+        // the search discarded everything below that probe, and the entry
+        // holding the real maximum went with it: on ERR12389696.20597455 the
+        // MMP came back 24 where the genome matches 29.
         if read_base != genome_base {
             return (match_len, read_base > genome_base);
         }
@@ -772,11 +833,23 @@ mod tests {
         Parameters::parse_from(full_args)
     }
 
+    /// Parameters for the toy fixtures below, whose reads are 3-8 bases.
+    ///
+    /// STAR's chain loop is `while (istart*Lstart + Lmapped + P.seedMapMin <
+    /// splitR[1][ip])`, so a piece of at most `seedMapMin` bases is never
+    /// searched at all. At the default `seedMapMin = 5` these reads would
+    /// produce no seeds for that reason alone, which is not what they are
+    /// testing; lower the threshold so the MMP mechanics are what the
+    /// assertions actually exercise.
+    fn toy_params() -> Parameters {
+        params(&["--seedMapMin", "1"])
+    }
+
     #[test]
     fn find_exact_match() {
         let index = make_test_index("ACGTACGT");
         let read = encode_sequence("ACGT");
-        let params = params(&["--runMode", "alignReads"]);
+        let params = toy_params();
 
         let seeds = Seed::find_seeds(&read, &index, 4, &params, "").unwrap();
 
@@ -792,7 +865,7 @@ mod tests {
     fn min_seed_length_filter() {
         let index = make_test_index("AAAAAAAA");
         let read = encode_sequence("AAA");
-        let params = params(&[]);
+        let params = toy_params();
 
         // With min_seed_length=4, should find nothing (read is only 3bp)
         let seeds = Seed::find_seeds(&read, &index, 4, &params, "").unwrap();
@@ -807,7 +880,7 @@ mod tests {
     fn no_match() {
         let index = make_test_index("ACAC");
         let read = encode_sequence("GGGG");
-        let params = params(&[]);
+        let params = toy_params();
 
         let seeds = Seed::find_seeds(&read, &index, 2, &params, "").unwrap();
 
@@ -819,7 +892,7 @@ mod tests {
     fn get_genome_positions() {
         let index = make_test_index("ACGTACGT");
         let read = encode_sequence("ACGT");
-        let params = params(&[]);
+        let params = toy_params();
 
         let seeds = Seed::find_seeds(&read, &index, 4, &params, "").unwrap();
         assert!(!seeds.is_empty());
@@ -838,7 +911,7 @@ mod tests {
     fn test_single_end_mate_id() {
         let index = make_test_index("ACGTACGT");
         let read = encode_sequence("ACGT");
-        let params = params(&[]);
+        let params = toy_params();
 
         let seeds = Seed::find_seeds(&read, &index, 4, &params, "").unwrap();
         assert!(!seeds.is_empty());
@@ -854,7 +927,7 @@ mod tests {
         let index = make_test_index("ACGTACGTTTGGCCAA");
         let mate1 = encode_sequence("ACGT");
         let mate2 = encode_sequence("TTGG");
-        let params = params(&[]);
+        let params = toy_params();
 
         let seeds = Seed::find_paired_seeds(&mate1, &mate2, &index, 4, &params).unwrap();
 
@@ -881,7 +954,7 @@ mod tests {
         let index = make_test_index("ACGTACGT");
         let mate1 = encode_sequence("ACGT");
         let mate2 = encode_sequence("ACGT");
-        let params = params(&[]);
+        let params = toy_params();
 
         let seeds = Seed::find_paired_seeds(&mate1, &mate2, &index, 4, &params).unwrap();
 
@@ -1023,10 +1096,12 @@ mod tests {
         search_direction_sparse(
             &rc_read,
             read.len(),
+            0,
             &index,
             4,
             &params,
             true,
+            true,
             "",
             &mut rc_seeds,
         );
@@ -1088,6 +1163,33 @@ mod tests {
         );
     }
 
+    /// A read base is 0-3, `N` is 4 and the spacer is 5, so a read base is
+    /// never greater than either. When the comparison runs off a chromosome
+    /// into the padding, the query sorts *before* that suffix, and
+    /// `compare_seq_to_genome` must say so: it is this value that orders the
+    /// binary search in `max_mappable_length`, so getting it backwards makes
+    /// the search discard the half of the range holding the real maximum.
+    #[test]
+    fn a_query_running_into_the_padding_sorts_before_that_suffix() {
+        let index = make_test_index("ACGTACGT");
+        let read = encode_sequence("ACGTACGTA");
+
+        // The suffix array entry for forward position 0.
+        let sa_idx = (0..index.suffix_array.len())
+            .find(|&i| {
+                let (pos, rev) = index.suffix_array.decode(index.suffix_array.get(i));
+                pos == 0 && !rev
+            })
+            .expect("forward position 0 must be in the suffix array");
+
+        let (len, comp) = compare_seq_to_genome(&read, 0, &index, sa_idx, 0);
+        assert_eq!(len, 8, "the whole chromosome matches before the padding");
+        assert!(
+            !comp,
+            "at the padding the read base is the smaller one, so the query sorts first"
+        );
+    }
+
     #[test]
     fn test_rc_seed_genome_positions() {
         // Genome: AACCTTGG, read RC = CCAAGGTT
diff --git a/src/align/stitch.rs b/src/align/stitch.rs
index 1c25ecc..b15e586 100644
--- a/src/align/stitch.rs
+++ b/src/align/stitch.rs
@@ -412,8 +412,19 @@ pub struct WindowAlignment {
     /// = length + left_ext_score + right_ext_score (in stitch coords, forward strand).
     /// Computed in stitch_seeds_core after dedup/sort. Default: length as i32.
     pub pre_ext_score: i32,
+    /// Annotated-junction index this seed came from (STAR: `WA_sjA`), or `-1`
+    /// for an ordinary real-genome hit. Set when the SA hit landed in the Gsj
+    /// flanking buffer, which means the seed straddles a known junction.
+    /// Addresses [`SpliceJunctionDb::entry`](crate::junction::SpliceJunctionDb::entry).
+    pub sj_a: i64,
 }
 
+/// One SA hit after Gsj expansion: `(forward_pos, read_offset, length,
+/// sa_pos, sj_a)`. A hit in the real genome yields exactly one of these, with
+/// `sj_a = -1`; a hit straddling the donor/acceptor boundary of the Gsj buffer
+/// yields two, both tagged with the junction they came from.
+type ExpandedHit = (u64, usize, usize, u64, i64);
+
 /// A cluster of seeds mapping to the same genomic region
 #[derive(Debug, Clone)]
 pub struct SeedCluster {
@@ -477,7 +488,7 @@ pub fn cluster_seeds(
     let prepared = index.prepared_junctions.as_slice();
 
     // Expand one raw SA hit into the candidate `(real_fwd_pos, read_offset,
-    // sub_length, sub_sa_pos)` tuples cluster_seeds should consume. Hits
+    // sub_length, sub_sa_pos, sj_a)` tuples cluster_seeds should consume. Hits
     // in the real genome pass through unchanged. Hits in the Gsj
     // flanking buffer are decoded via the sjdb table - yielding one
     // entry for hits confined to a single flank, two entries for hits
@@ -488,45 +499,46 @@ pub fn cluster_seeds(
     // `[Option<T>; 2]` avoids a heap `Vec` allocation on every SA hit expanded
     // in this hot loop (the common `raw_fwd < n_genome_real` case is the
     // overwhelming majority of calls); iterate with `.into_iter().flatten()`.
-    let expand_hit =
-        |sa_pos: u64, strand: bool, length: usize| -> [Option<(u64, usize, usize, u64)>; 2] {
-            let raw_fwd = index.sa_pos_to_forward(sa_pos, strand, length);
-            if raw_fwd < n_genome_real {
-                return [Some((raw_fwd, 0, length, sa_pos)), None];
-            }
-            if sjdb_overhang == 0 || prepared.is_empty() {
-                return [None, None];
-            }
-            let mut decoded = crate::junction::sjdb_insert::decode_gsj_hit(
-                raw_fwd,
-                length,
-                n_genome_real,
-                sjdb_overhang,
-                prepared,
-            );
-            // Reverse-strand hits traverse the donor/acceptor halves in reverse
-            // read order: the leftmost forward bytes (donor flank) align to the
-            // last read bases. Swap the read offsets so each sub-seed's
-            // `read_pos = seed.read_pos + read_offset` lands at the right place
-            // in original-read coords.
-            if strand && decoded.len() == 2 {
-                let acceptor_len = decoded[1].2;
-                decoded[0].1 = acceptor_len;
-                decoded[1].1 = 0;
-            }
-            let mut out = [None, None];
-            for (slot, (real_fwd, read_off, sub_len)) in out.iter_mut().zip(decoded) {
-                let sub_sa_pos = if strand {
-                    n_genome
-                        .saturating_sub(real_fwd)
-                        .saturating_sub(sub_len as u64)
-                } else {
-                    real_fwd
-                };
-                *slot = Some((real_fwd, read_off, sub_len, sub_sa_pos));
-            }
-            out
-        };
+    let expand_hit = |sa_pos: u64, strand: bool, length: usize| -> [Option<ExpandedHit>; 2] {
+        let raw_fwd = index.sa_pos_to_forward(sa_pos, strand, length);
+        if raw_fwd < n_genome_real {
+            // A real-genome hit carries no annotated-junction tag (STAR's
+            // sjA = -1).
+            return [Some((raw_fwd, 0, length, sa_pos, -1)), None];
+        }
+        if sjdb_overhang == 0 || prepared.is_empty() {
+            return [None, None];
+        }
+        let mut decoded = crate::junction::sjdb_insert::decode_gsj_hit(
+            raw_fwd,
+            length,
+            n_genome_real,
+            sjdb_overhang,
+            prepared,
+        );
+        // Reverse-strand hits traverse the donor/acceptor halves in reverse
+        // read order: the leftmost forward bytes (donor flank) align to the
+        // last read bases. Swap the read offsets so each sub-seed's
+        // `read_pos = seed.read_pos + read_offset` lands at the right place
+        // in original-read coords.
+        if strand && decoded.len() == 2 {
+            let acceptor_len = decoded[1].2;
+            decoded[0].1 = acceptor_len;
+            decoded[1].1 = 0;
+        }
+        let mut out = [None, None];
+        for (slot, (real_fwd, read_off, sub_len, sj_a)) in out.iter_mut().zip(decoded) {
+            let sub_sa_pos = if strand {
+                n_genome
+                    .saturating_sub(real_fwd)
+                    .saturating_sub(sub_len as u64)
+            } else {
+                real_fwd
+            };
+            *slot = Some((real_fwd, read_off, sub_len, sub_sa_pos, sj_a as i64));
+        }
+        out
+    };
 
     let anchor_set: Vec<bool> = seeds
         .iter()
@@ -595,7 +607,7 @@ pub fn cluster_seeds(
                 continue;
             }
 
-            for (forward_pos, _read_off, sub_length, _sub_sa_pos) in
+            for (forward_pos, _read_off, sub_length, _sub_sa_pos, _sj_a) in
                 expand_hit(sa_pos, strand, full_length)
                     .into_iter()
                     .flatten()
@@ -766,6 +778,8 @@ pub fn cluster_seeds(
         /// by the donor-side length for the acceptor half of a Gsj
         /// boundary-crossing hit.
         read_pos: usize,
+        /// Annotated-junction index (STAR: `sjA`), or `-1`.
+        sj_a: i64,
     }
 
     let mut win_candidates: Vec<Vec<WinCandidate>> =
@@ -784,7 +798,7 @@ pub fn cluster_seeds(
                 continue;
             }
 
-            for (forward_pos, read_off, sub_length, sub_sa_pos) in
+            for (forward_pos, read_off, sub_length, sub_sa_pos, sj_a) in
                 expand_hit(sa_pos, strand, full_length)
                     .into_iter()
                     .flatten()
@@ -819,6 +833,7 @@ pub fn cluster_seeds(
                     ps_rstart,
                     mate_id: seed.mate_id,
                     read_pos: derived_read_pos,
+                    sj_a,
                 });
             }
         }
@@ -847,7 +862,7 @@ pub fn cluster_seeds(
     // + hash-table rehashes per read (a measured allocator hotspot). Reused
     // storage only; the results are identical.
     let mut by_len: Vec<usize> = Vec::new();
-    let mut diag_ranges: FxHashMap<(i64, u8), Vec<(usize, usize)>> = FxHashMap::default();
+    let mut diag_ranges: FxHashMap<(i64, u8, i64), Vec<(usize, usize)>> = FxHashMap::default();
 
     for win_idx in 0..win_n {
         let candidates = &win_candidates[win_idx];
@@ -860,15 +875,18 @@ pub fn cluster_seeds(
         by_len.extend(0..candidates.len());
         by_len.sort_by(|&a, &b| candidates[b].length.cmp(&candidates[a].length));
 
-        // For each (diagonal, mate_id) pair, track accepted [ps_rstart, ps_rend) ranges.
-        // STAR's assignAlignToWindow checks aFrag==WA[iA][WA_iFrag] before overlap test:
-        // seeds from different fragments are never treated as overlapping duplicates.
+        // For each (diagonal, mate_id, sj_a) triple, track accepted
+        // [ps_rstart, ps_rend) ranges. STAR's assignAlignToWindow checks both
+        // `aFrag == WA[iA][WA_iFrag]` and `sjA == WA[iA][WA_sjA]` before the
+        // overlap test: seeds from different fragments, or derived from
+        // different annotated junctions, are never treated as duplicates even
+        // when they land on the same diagonal.
         diag_ranges.clear();
         for &ci in &by_len {
             let cand = &candidates[ci];
             let diag = cand.forward_pos as i64 - cand.ps_rstart as i64;
             let ps_rend = cand.ps_rstart + cand.length;
-            let key = (diag, cand.mate_id);
+            let key = (diag, cand.mate_id, cand.sj_a);
 
             let blocked = diag_ranges.get(&key).is_some_and(|ranges| {
                 ranges.iter().any(|&(rs, re)| {
@@ -930,6 +948,9 @@ pub fn cluster_seeds(
                 if wa.mate_id != new_mate_id {
                     continue; // STAR: only merge seeds from the same fragment
                 }
+                if wa.sj_a != cand.sj_a {
+                    continue; // STAR: `e.sjA == sjA` guards the overlap test
+                }
                 let wa_ps_rstart = if window.is_reverse {
                     read_len - (wa.length + wa.read_pos)
                 } else {
@@ -968,6 +989,7 @@ pub fn cluster_seeds(
                             is_anchor: is_anchor_seed,
                             mate_id: seed.mate_id,
                             pre_ext_score: length as i32,
+                            sj_a: cand.sj_a,
                         },
                     );
                 }
@@ -1036,6 +1058,7 @@ pub fn cluster_seeds(
                     is_anchor: is_anchor_seed,
                     mate_id: seed.mate_id,
                     pre_ext_score: length as i32,
+                    sj_a: cand.sj_a,
                 },
             );
         }
@@ -1076,6 +1099,11 @@ pub(crate) struct ExonBlock {
     pub(crate) genome_end: u64,   // SA coordinate space (exclusive)
     /// Mate ID: 0=mate1, 1=mate2, 2=SE (STAR: EX_iFrag)
     pub(crate) mate_id: u8,
+    /// Annotated-junction index of the seed this exon was built from
+    /// (STAR: `EX_sjA`), or `-1`. Compared against the incoming seed's `sj_a`
+    /// to recognise that two exons are the two flanks of the same annotated
+    /// junction, which lets the stitcher take that junction verbatim.
+    pub(crate) sj_a: i64,
 }
 
 /// In-progress transcript during recursive search (cheap to clone)
@@ -1092,6 +1120,21 @@ pub(crate) struct WorkingTranscript {
     /// STAR's shiftSJ[isj][0] and shiftSJ[isj][1].
     pub(crate) junction_shifts: Vec<(u32, u32)>,
     pub(crate) n_anchor: u32,
+    /// Score used to *rank* this transcript against its siblings inside the
+    /// recursion: `score` plus the genomic-length penalty.
+    ///
+    /// STAR applies `scoreGenomicLengthLog2scale` in the base case of
+    /// `stitchWindowAligns`, before the dedup and eviction that pick which
+    /// transcripts survive, so a compact alignment can beat a sprawling one of
+    /// equal raw score. rustar-aligner applied it only at finalization, i.e.
+    /// after those decisions had already been made.
+    ///
+    /// Kept separate from `score` rather than folded into it: the PE mate-split
+    /// path builds two per-mate transcripts from one combined `score`, and
+    /// splitting an already-penalised score across mates has no clean meaning.
+    /// `finalize_transcript` still applies the penalty to the final score, so
+    /// nothing is counted twice.
+    pub(crate) rank_score: i32,
     // Tight bounds for extension at finalization
     pub(crate) read_start: usize,
     pub(crate) read_end: usize,
@@ -1111,6 +1154,7 @@ impl WorkingTranscript {
             junction_annotated: Vec::new(),
             junction_shifts: Vec::new(),
             n_anchor: 0,
+            rank_score: 0,
             read_start: 0,
             read_end: 0,
             genome_start: 0,
@@ -1136,6 +1180,68 @@ fn stitch_align_to_transcript(
 ) -> Option<WorkingTranscript> {
     let last_exon = wt.exons.last().unwrap();
 
+    // STAR's "simple stitching if junction belongs to a database"
+    // (`stitchAlignToTranscript.cpp`). When the previous exon and the incoming
+    // piece B carry the *same* annotated-junction tag, sit on the same mate,
+    // and are exactly read-adjacent, the two are the donor and acceptor flanks
+    // of one annotated junction that the Gsj insert split apart. Stitch them
+    // straight onto the annotated boundary: exact donor/acceptor coordinates,
+    // the stored motif and shifts, and the `sjdbScore` bonus, skipping the
+    // general boundary search entirely.
+    //
+    // Without this, the general search flushes a repeat-adjacent junction to
+    // its leftmost position and lands a few bases off the annotated boundary
+    // (e.g. `9M1545N40M` where STAR emits the annotated `11M1545N38M`), and
+    // forfeits the bonus because the shifted coordinates no longer match the
+    // database.
+    //
+    // Receiver coordinates are half-open where STAR's are inclusive, hence
+    // `wa.read_pos == last_exon.read_end` for STAR's `rBstart == rAend + 1`
+    // and `last_exon.genome_end < wa.sa_pos` for `gAend + 1 < gBstart`.
+    if wa.sj_a != -1
+        && last_exon.sj_a == wa.sj_a
+        && last_exon.mate_id == wa.mate_id
+        && wa.read_pos == last_exon.read_end
+        && last_exon.genome_end < wa.sa_pos
+        && let Some(db) = junction_db
+        && let Some(pj) = db.entry(wa.sj_a as usize)
+    {
+        // Too-large repeats around a non-canonical annotated junction: STAR
+        // rejects with -1000006, which is a `None` here — the caller still
+        // takes the EXCLUDE branch, exactly as STAR's recursion does.
+        let exon_a_len = last_exon.read_end - last_exon.read_start;
+        if pj.motif == 0
+            && (wa.length <= pj.shift_right as usize || exon_a_len <= pj.shift_left as usize)
+        {
+            return None;
+        }
+
+        let mut new_wt = wt.clone();
+        new_wt.exons.push(ExonBlock {
+            read_start: wa.read_pos,
+            read_end: wa.read_pos + wa.length,
+            genome_start: wa.sa_pos,
+            genome_end: wa.sa_pos + wa.length as u64,
+            mate_id: wa.mate_id,
+            sj_a: wa.sj_a,
+        });
+        new_wt.n_junction += 1;
+        new_wt
+            .junction_motifs
+            .push(crate::junction::decode_motif(pj.motif));
+        new_wt.junction_annotated.push(true);
+        new_wt
+            .junction_shifts
+            .push((pj.shift_left as u32, pj.shift_right as u32));
+        new_wt.score += wa.length as i32 + scorer.sjdb_score;
+        new_wt.read_end = wa.read_pos + wa.length;
+        new_wt.genome_end = wa.sa_pos + wa.length as u64;
+        if wa.is_anchor {
+            new_wt.n_anchor += 1;
+        }
+        return Some(new_wt);
+    }
+
     // Mate-boundary detection: STAR canonSJ[iex] = -3 (stitchAlignToTranscript.cpp:402)
     // When crossing from mate1 to mate2 (or vice versa), skip junction scoring and
     // check alignMatesGapMax instead.
@@ -1230,6 +1336,7 @@ fn stitch_align_to_transcript(
             genome_start: wa.sa_pos,
             genome_end: wa.sa_pos + wa.length as u64,
             mate_id: wa.mate_id,
+            sj_a: wa.sj_a,
         });
         new_wt.read_end = wa.read_pos + wa.length;
         new_wt.genome_end = wa.sa_pos + wa.length as u64;
@@ -1370,23 +1477,24 @@ fn stitch_align_to_transcript(
         // is motif detection (splice) vs pure positional score (deletion).
         // donor_sa = exclusive end of exon A = STAR's gAend+1. jr_shift = STAR's jR.
         let donor_sa = last_exon.genome_end;
-        let (jr_shift, motif, motif_score, jj_l, jj_r) = scorer.find_best_junction_position(
-            read_seq,
-            last_exon.read_end,
-            donor_sa,
-            read_gap.max(0),
-            genome_gap,
-            &index.genome,
-            cluster.is_reverse,
-            index.genome.n_genome,
-            last_exon.read_end - last_exon.read_start,
-            eff_length,
-        );
+        let (jr_shift, mut motif, motif_score, mut jj_l, mut jj_r) = scorer
+            .find_best_junction_position(
+                read_seq,
+                last_exon.read_end,
+                donor_sa,
+                read_gap.max(0),
+                genome_gap,
+                &index.genome,
+                cluster.is_reverse,
+                index.genome.n_genome,
+                last_exon.read_end - last_exon.read_start,
+                eff_length,
+            );
 
         // Clamp shift: jr_shift = STAR's jR. Lower bound: can't consume entire exon A.
         // Upper bound: scan already limited to < shared+eff_length but clamp for safety.
         let prev_match_len = (last_exon.read_end - last_exon.read_start) as i32;
-        let jr_shift = jr_shift
+        let mut jr_shift = jr_shift
             .max(-prev_match_len)
             .min((eff_length + shared) as i32);
 
@@ -1467,23 +1575,67 @@ fn stitch_align_to_transcript(
         gap_mm += shared_mm;
         d_score += shared_score;
 
+        // Look the junction up in the annotation *before* deciding what kind of
+        // gap this is. STAR resolves `sjdbInd` first and, on a hit, takes the
+        // annotated branch whatever the gap length: `Del >= alignIntronMin`
+        // only decides between junction and deletion for gaps the annotation
+        // does not know (`stitchAlignToTranscript.cpp:198-241`). An annotated
+        // one-base gap is a junction to STAR, not a deletion.
+        //
+        // The lookup also supplies the motif, which is what the mismatch gate
+        // below sees rather than the one scanned from the genome.
+        let junc_donor_sa = (donor_sa as i64 + jr_shift as i64) as u64;
+        let donor_fwd = index.sa_pos_to_forward(junc_donor_sa, cluster.is_reverse, del as usize);
+        let acceptor_fwd = donor_fwd + del as u64 - 1;
+
+        // The metadata-bearing lookup. Junctions inserted by two-pass mode
+        // carry no motif or shift, so they are only visible through the
+        // annotated map below; the two lookups are therefore separate.
+        let sj_entry =
+            junction_db.and_then(|db| db.find(donor_fwd, acceptor_fwd).and_then(|i| db.entry(i)));
+        let is_annotated = sj_entry.is_some()
+            || junction_db.is_some_and(|db| {
+                db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 0)
+                    || db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 1)
+                    || db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 2)
+            });
+        let is_splice = is_annotated || is_splice;
+
         // --- Type-specific scoring and tracking ---
         if is_splice {
-            // Check stitch mismatch limit
+            if let Some(pj) = sj_entry {
+                // The annotated entry's motif wins over the one scanned from
+                // the genome (STAR: `jcan = sjdbMotif[ind]`). This is what the
+                // mismatch gate, `outFilterIntronMotifs` and SJ.out.tab see.
+                motif = crate::junction::decode_motif(pj.motif);
+
+                // Non-canonical annotated junction: snap the boundary onto the
+                // annotated coordinates instead of leaving it wherever the
+                // leftmost-flush scan put it (STAR shifts `jR` by `shiftLeft`).
+                // `jr_shift` is consumed below to move exon A's right edge and
+                // exon B's left edge together, so shifting it here relocates
+                // both coherently and leaves the intron length unchanged.
+                if pj.motif == 0 {
+                    let sl = pj.shift_left as i32;
+                    if (eff_length as i32) <= sl || prev_match_len <= sl {
+                        return None; // STAR -1000006
+                    }
+                    jr_shift += sl;
+                    if (last_exon.read_end as i64 + jr_shift as i64)
+                        >= (eff_read_pos + eff_length) as i64
+                    {
+                        return None; // STAR -1000006
+                    }
+                    jj_l = pj.shift_left as u32;
+                    jj_r = pj.shift_right as u32;
+                }
+            }
+
+            // Check stitch mismatch limit (now against the annotated motif).
             if !scorer.stitch_mismatch_allowed(&motif, gap_mm) {
                 return None;
             }
 
-            let is_annotated = junction_db.is_some_and(|db| {
-                let junc_donor_sa = (donor_sa as i64 + jr_shift as i64) as u64;
-                let donor_fwd =
-                    index.sa_pos_to_forward(junc_donor_sa, cluster.is_reverse, del as usize);
-                let acceptor_fwd = donor_fwd + del as u64 - 1;
-                db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 0)
-                    || db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 1)
-                    || db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 2)
-            });
-
             if is_annotated {
                 d_score += scorer.sjdb_score;
             } else {
@@ -1522,6 +1674,7 @@ fn stitch_align_to_transcript(
             genome_start: b_genome_start,
             genome_end: b_genome_start + b_len as u64,
             mate_id: wa.mate_id,
+            sj_a: wa.sj_a,
         });
     } else {
         // Insertion: read_gap > genome_gap
@@ -1530,14 +1683,15 @@ fn stitch_align_to_transcript(
 
         let mut jr = 0i32; // number of shared bases going to A side
 
+        let genome_offset: u64 = if cluster.is_reverse {
+            index.genome.n_genome
+        } else {
+            0
+        };
+
         if shared > 0 {
             // jR scanning to find optimal insertion placement (STAR lines 265-291)
             let shared_usize = shared as usize;
-            let genome_offset: u64 = if cluster.is_reverse {
-                index.genome.n_genome
-            } else {
-                0
-            };
 
             let mut score1 = 0i32;
             let mut max_score1 = 0i32;
@@ -1566,10 +1720,11 @@ fn stitch_align_to_transcript(
                     }
                 }
 
-                // STAR default (alignInsertionFlush=None): strict > only.
-                // First maximum wins = leftmost insertion in current coordinate space.
-                // For flushRight mode (not yet implemented): Score1 >= maxScore1.
-                if score1 > max_score1 {
+                // `alignInsertionFlush None` (STAR's default): strict `>`, so
+                // the first maximum wins and the insertion sits at its leftmost
+                // position. `Right` accepts ties too, walking the insertion to
+                // the rightmost equally-scoring position.
+                if score1 > max_score1 || (score1 == max_score1 && scorer.flush_right) {
                     max_score1 = score1;
                     jr = jr1 as i32;
                 }
@@ -1606,6 +1761,27 @@ fn stitch_align_to_transcript(
         }
         // shared == 0: simple insertion, no scanning needed, jr stays 0
 
+        // `alignInsertionFlush Right`: walk the insertion further right for as
+        // long as the read keeps matching the genome, so a repeat-ambiguous
+        // insertion lands at the rightmost position it can. Running out of
+        // read on the B side means there is nowhere left to put the insertion
+        // (STAR -1000009).
+        if scorer.flush_right {
+            let jr_limit =
+                (eff_read_pos + eff_length).saturating_sub(last_exon.read_end + ins) as i32;
+            while jr < jr_limit {
+                let r_idx = last_exon.read_end + jr as usize;
+                let g_pos = last_exon.genome_end + jr as u64 + genome_offset;
+                match index.genome.get_base(g_pos) {
+                    Some(gb) if gb < 4 && read_seq.get(r_idx) == Some(&gb) => jr += 1,
+                    _ => break,
+                }
+            }
+            if jr >= jr_limit {
+                return None;
+            }
+        }
+
         let ins_score = scorer.score_ins_open + scorer.score_ins_base * ins as i32;
         d_score += ins_score;
         new_wt.n_gap += 1;
@@ -1631,13 +1807,17 @@ fn stitch_align_to_transcript(
             genome_start: b_genome_start,
             genome_end: eff_genome_pos + eff_length as u64,
             mate_id: wa.mate_id,
+            sj_a: wa.sj_a,
         });
     }
 
-    // Mismatch limit check
+    // Mismatch limit check. STAR's `outFilterMismatchNmaxTotal` is the minimum
+    // of three caps: the absolute one, a fraction of the *mapped* length, and
+    // a fraction of the *read* length. The third was missing here, which made
+    // `--outFilterMismatchNoverReadLmax` unenforceable.
     let total_mm = new_wt.n_mismatch + gap_mm;
     let total_len = new_wt.read_end.max(eff_read_pos + eff_length) - new_wt.read_start;
-    let mm_limit = ((scorer.p_mm_max * total_len as f64) as u32).min(scorer.n_mm_max);
+    let mm_limit = scorer.mismatch_nmax_total(total_len, read_seq.len());
     if total_mm > mm_limit {
         return None;
     }
@@ -1881,6 +2061,50 @@ pub(crate) fn finalize_transcript(
         return None;
     }
 
+    // `--alignSoftClipAtReferenceEnds No`: an alignment may not be
+    // soft-clipped past the end of its chromosome, so discard any that would
+    // be. `Yes`, the default, leaves behaviour unchanged.
+    if !scorer.soft_clip_at_reference_ends {
+        let chr = cluster.chr_idx;
+        let chr_start = index.genome.chr_start[chr];
+        let chr_end = chr_start + index.genome.chr_length[chr];
+        let span_start = wt
+            .genome_start
+            .saturating_sub(left_extend.extend_len as u64);
+        let span_end = wt.genome_end + right_extend.extend_len as u64;
+        let clipped_left = alignment_start - left_extend.extend_len;
+        let clipped_right = read_seq
+            .len()
+            .saturating_sub(alignment_end + right_extend.extend_len);
+        if (clipped_left > 0 && span_start <= chr_start)
+            || (clipped_right > 0 && span_end >= chr_end)
+        {
+            return None;
+        }
+    }
+
+    // Whether a genome gap of `del` bases starting at `donor_end` (SA space) is
+    // a junction rather than a deletion. `--alignIntronMin` decides for gaps the
+    // annotation does not know; an annotated junction is a junction whatever its
+    // length, which is how STAR reads it (`stitchAlignToTranscript.cpp:198-241`)
+    // and the reason yeast's three one-base introns come out as `1N`, not `1D`.
+    let gap_is_junction = |donor_end: u64, del: usize| -> bool {
+        if del >= scorer.align_intron_min as usize && del <= scorer.align_intron_max as usize {
+            return true;
+        }
+        if del == 0 {
+            return false;
+        }
+        let donor_fwd = index.sa_pos_to_forward(donor_end, cluster.is_reverse, del);
+        let acceptor_fwd = donor_fwd + del as u64 - 1;
+        index.junction_db.find(donor_fwd, acceptor_fwd).is_some()
+            || (0..=2).any(|strand| {
+                index
+                    .junction_db
+                    .is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, strand)
+            })
+    };
+
     // STAR finalization check: exon lengths including repeat lengths (shiftSJ)
     // For non-annotated junctions: exon_len >= alignSJoverhangMin + shiftSJ[side]
     // For annotated junctions: exon_len >= alignSJDBoverhangMin
@@ -1896,7 +2120,10 @@ pub(crate) fn finalize_transcript(
             let genome_gap = next_exon.genome_start as i64 - exon.genome_end as i64;
             let read_gap = next_exon.read_start as i64 - exon.read_end as i64;
             let del = genome_gap - read_gap.max(0);
-            if del >= scorer.align_intron_min as i64 && junction_idx < wt.junction_shifts.len() {
+            if del > 0
+                && gap_is_junction(exon.genome_end, del as usize)
+                && junction_idx < wt.junction_shifts.len()
+            {
                 // This is a junction — check exon lengths with repeat
                 let (shift_l, shift_r) = wt.junction_shifts[junction_idx];
                 let is_annotated = wt.junction_annotated[junction_idx];
@@ -1988,9 +2215,7 @@ pub(crate) fn finalize_transcript(
                     append_match(&mut final_cigar, shared);
                 }
                 let del = (genome_gap - read_gap.max(0)) as usize;
-                if del >= scorer.align_intron_min as usize
-                    && del <= scorer.align_intron_max as usize
-                {
+                if gap_is_junction(prev.genome_end, del) {
                     final_cigar.push(Op::new(Kind::Skip, del));
                 } else {
                     final_cigar.push(Op::new(Kind::Deletion, del));
@@ -2338,6 +2563,13 @@ fn stitch_recurse(
                 }
             }
 
+            // STAR applies the genomic-length penalty here, in the base case,
+            // *before* dedup and eviction. Ranking on the raw score instead
+            // lets a sprawling alignment displace a compact one of equal
+            // score, which STAR would never do.
+            wt.rank_score = wt.score
+                + scorer.genomic_length_penalty(wt.genome_end.saturating_sub(wt.genome_start));
+
             // Dedup via blocks_overlap: drop if subset of existing higher-score transcript.
             // Use same_structure guard: only dedup transcripts with same number of exon
             // blocks. A non-spliced path should never be killed by a spliced one here
@@ -2361,11 +2593,11 @@ fn stitch_recurse(
 
                 // Only dedup transcripts with same number of exon blocks (junctions).
                 let same_structure = wt.exons.len() == existing.exons.len();
-                if same_structure && overlap >= wt_len && existing.score >= wt.score {
+                if same_structure && overlap >= wt_len && existing.rank_score >= wt.rank_score {
                     dominated = true;
                     break;
                 }
-                if same_structure && overlap >= ex_len && wt.score >= existing.score {
+                if same_structure && overlap >= ex_len && wt.rank_score >= existing.rank_score {
                     remove_indices.push(idx);
                 }
             }
@@ -2380,8 +2612,8 @@ fn stitch_recurse(
                 } else if let Some(worst_idx) = transcripts
                     .iter()
                     .enumerate()
-                    .min_by_key(|(_, t)| t.score)
-                    .filter(|(_, t)| t.score < wt.score)
+                    .min_by_key(|(_, t)| t.rank_score)
+                    .filter(|(_, t)| t.rank_score < wt.rank_score)
                     .map(|(i, _)| i)
                 {
                     // STAR-faithful eviction: keep the N best transcripts.
@@ -2407,6 +2639,7 @@ fn stitch_recurse(
             genome_start: wa.sa_pos,
             genome_end: wa.sa_pos + wa.length as u64,
             mate_id: wa.mate_id,
+            sj_a: wa.sj_a,
         });
         new_wt.score = wa.length as i32;
         new_wt.read_start = wa.read_pos;
@@ -2629,6 +2862,9 @@ pub(crate) fn split_combined_wt(
             junction_annotated: m1_ja,
             junction_shifts: m1_js,
             n_anchor: 0,
+            // Per-mate transcripts are never ranked against each other inside
+            // the recursion; `finalize_transcript` applies the penalty.
+            rank_score: m1_score,
             read_start: m1_read_start,
             read_end: m1_read_end,
             genome_start: m1_genome_start,
@@ -2644,6 +2880,9 @@ pub(crate) fn split_combined_wt(
             junction_annotated: m2_ja,
             junction_shifts: m2_js,
             n_anchor: 0,
+            // Per-mate transcripts are never ranked against each other inside
+            // the recursion; `finalize_transcript` applies the penalty.
+            rank_score: m2_score,
             read_start: m2_read_start,
             read_end: m2_read_end,
             genome_start: m2_genome_start,
@@ -2784,6 +3023,74 @@ pub(crate) fn stitch_seeds_with_jdb_debug(
     transcripts
 }
 
+/// Collapse window alignments that are redundant copies of one another before
+/// the recursive stitcher runs.
+///
+/// For each diagonal (`genome_pos - positive_strand_read_start`), overlapping
+/// entries are merged into intervals and only the longest entry of each merged
+/// interval is kept. Without this the stitcher's include/exclude recursion
+/// explodes on the many redundant seeds that cover the same diagonal.
+///
+/// The grouping key is `(diagonal, mate_id, sjA)`, which is what STAR's
+/// `assignAlignToWindow` compares before treating two aligns as duplicates: it
+/// guards the overlap test with both `aFrag == WA[iA][WA_iFrag]` and
+/// `sjA == WA[iA][WA_sjA]`. Entries split out of *different* annotated
+/// junctions are therefore never duplicates of each other, even when they sit
+/// on the same diagonal and cover the same read bases.
+fn dedup_wa_by_diagonal(wa_entries: &mut Vec<WindowAlignment>, read_len: usize, is_rev: bool) {
+    use rustc_hash::{FxHashMap, FxHashSet};
+
+    type DiagKey = (i64, u8, i64);
+    type DiagSeeds = Vec<(usize, usize, usize)>;
+    let mut diag_seeds: FxHashMap<DiagKey, DiagSeeds> = FxHashMap::default();
+    for (idx, wa) in wa_entries.iter().enumerate() {
+        let ps = if is_rev {
+            read_len - (wa.length + wa.read_pos)
+        } else {
+            wa.read_pos
+        };
+        let diag = wa.genome_pos as i64 - ps as i64;
+        diag_seeds
+            .entry((diag, wa.mate_id, wa.sj_a))
+            .or_default()
+            .push((ps, ps + wa.length, idx));
+    }
+
+    let mut keep_indices = FxHashSet::default();
+    for (_key, mut seeds) in diag_seeds {
+        seeds.sort_unstable();
+        let mut merged_end = seeds[0].1;
+        let mut best_idx = seeds[0].2;
+        let mut best_len = seeds[0].1 - seeds[0].0;
+
+        for &(s, e, idx) in &seeds[1..] {
+            if s <= merged_end {
+                // Overlapping: extend the interval, track the longest entry in it.
+                merged_end = merged_end.max(e);
+                let len = e - s;
+                if len > best_len {
+                    best_len = len;
+                    best_idx = idx;
+                }
+            } else {
+                // Disjoint: commit the previous interval's best and start a new one.
+                keep_indices.insert(best_idx);
+                merged_end = e;
+                best_len = e - s;
+                best_idx = idx;
+            }
+        }
+        keep_indices.insert(best_idx);
+    }
+
+    let mut idx = 0usize;
+    wa_entries.retain(|_| {
+        let keep = keep_indices.contains(&idx);
+        idx += 1;
+        keep
+    });
+}
+
 /// Shared core: preprocessing + recursive stitcher, returns working transcripts + context.
 #[allow(clippy::too_many_arguments)]
 pub(crate) fn stitch_seeds_core(
@@ -2805,72 +3112,7 @@ pub(crate) fn stitch_seeds_core(
     // of repetitive seeds, matching STAR's WA_Anchor=2 "last anchor" logic.
     let mut wa_entries: Vec<WindowAlignment> = cluster.alignments.clone();
 
-    // Diagonal dedup: for each diagonal (genome_pos - ps_rstart), merge overlapping
-    // seeds into intervals, keeping only the longest seed per merged interval.
-    // This prevents combinatorial explosion in the recursive stitcher when many
-    // redundant seeds cover the same diagonal region.
-    // Uses positive-strand coordinates consistent with cluster_seeds overlap detection.
-    {
-        use rustc_hash::{FxHashMap, FxHashSet};
-        let read_len = read_seq.len();
-        let is_rev = cluster.is_reverse;
-        // For each (diagonal, mate_id) pair, find the longest seed per merged interval.
-        // STAR's assignAlignToWindow checks aFrag==WA[iA][WA_iFrag] before overlap test:
-        // seeds from different fragments are never treated as duplicates.
-        type DiagMateKey = (i64, u8);
-        type DiagSeeds = Vec<(usize, usize, usize)>;
-        let mut diag_seeds: FxHashMap<DiagMateKey, DiagSeeds> = FxHashMap::default();
-        for (idx, wa) in wa_entries.iter().enumerate() {
-            let ps = if is_rev {
-                read_len - (wa.length + wa.read_pos)
-            } else {
-                wa.read_pos
-            };
-            let diag = wa.genome_pos as i64 - ps as i64;
-            diag_seeds
-                .entry((diag, wa.mate_id))
-                .or_default()
-                .push((ps, ps + wa.length, idx));
-        }
-
-        let mut keep_indices = FxHashSet::default();
-        for (_diag, mut seeds) in diag_seeds {
-            // Sort by start position
-            seeds.sort_unstable();
-            // Merge intervals, keeping the index of the longest seed in each merged group
-            let mut merged_end = seeds[0].1;
-            let mut best_idx = seeds[0].2;
-            let mut best_len = seeds[0].1 - seeds[0].0;
-
-            for &(s, e, idx) in &seeds[1..] {
-                if s <= merged_end {
-                    // Overlapping — extend and track longest
-                    merged_end = merged_end.max(e);
-                    let len = e - s;
-                    if len > best_len {
-                        best_len = len;
-                        best_idx = idx;
-                    }
-                } else {
-                    // New interval — commit previous best
-                    keep_indices.insert(best_idx);
-                    merged_end = e;
-                    best_len = e - s;
-                    best_idx = idx;
-                }
-            }
-            // Commit last group
-            keep_indices.insert(best_idx);
-        }
-
-        // Retain only the kept indices
-        let mut idx = 0usize;
-        wa_entries.retain(|_| {
-            let keep = keep_indices.contains(&idx);
-            idx += 1;
-            keep
-        });
-    }
+    dedup_wa_by_diagonal(&mut wa_entries, read_seq.len(), cluster.is_reverse);
 
     // STAR-faithful coordinate conversion for stitching:
     // STAR stores WA_gStart in FORWARD genome coordinates (converting RC positions via
@@ -2946,8 +3188,15 @@ pub(crate) fn stitch_seeds_core(
         );
         for (i, wa) in wa_entries.iter().enumerate().take(30) {
             eprintln!(
-                "  wa[{}]: read_pos={}, sa_pos={}, genome_pos={}, length={}, anchor={}, mate={}",
-                i, wa.read_pos, wa.sa_pos, wa.genome_pos, wa.length, wa.is_anchor, wa.mate_id
+                "  wa[{}]: read_pos={}, sa_pos={}, genome_pos={}, length={}, anchor={}, mate={}, sjA={}",
+                i,
+                wa.read_pos,
+                wa.sa_pos,
+                wa.genome_pos,
+                wa.length,
+                wa.is_anchor,
+                wa.mate_id,
+                wa.sj_a
             );
         }
         if wa_entries.len() > 30 {
@@ -3141,6 +3390,65 @@ mod tests {
     use crate::index::sa_index::SaIndex;
     use crate::index::suffix_array::SuffixArray;
 
+    fn wa(read_pos: usize, length: usize, genome_pos: u64, sj_a: i64) -> WindowAlignment {
+        WindowAlignment {
+            seed_idx: 0,
+            read_pos,
+            length,
+            genome_pos,
+            sa_pos: genome_pos,
+            n_rep: 1,
+            is_anchor: true,
+            mate_id: 2,
+            pre_ext_score: length as i32,
+            sj_a,
+        }
+    }
+
+    /// A read whose last exon is a few bases long produces two window
+    /// alignments covering the same donor-side bases: one tagged with the
+    /// annotated junction the read arrives on, one tagged with the junction it
+    /// leaves by. STAR keeps both, because `assignAlignToWindow` only treats
+    /// two aligns as duplicates when their `sjA` match as well as their
+    /// fragment and diagonal.
+    ///
+    /// Keying the dedup on the diagonal alone dropped the departing copy, which
+    /// left the micro-exon entry without a partner carrying the same `sjA`, so
+    /// the terminal exon came out soft-clipped instead of spliced.
+    #[test]
+    fn same_diagonal_entries_from_different_junctions_both_survive() {
+        let mut entries = vec![
+            wa(61, 36, 14_434_354, 348),
+            wa(61, 36, 14_434_354, 349),
+            wa(97, 3, 14_436_873, 349),
+        ];
+        dedup_wa_by_diagonal(&mut entries, 100, false);
+
+        assert_eq!(entries.len(), 3, "no entry may be dropped: {entries:?}");
+        let tags: Vec<i64> = entries.iter().map(|e| e.sj_a).collect();
+        assert!(tags.contains(&348) && tags.contains(&349));
+        assert!(
+            entries.iter().any(|e| e.length == 3 && e.sj_a == 349),
+            "the micro-exon entry must keep a partner with the same sjA"
+        );
+    }
+
+    /// The dedup still collapses what it is there to collapse: redundant
+    /// overlapping copies on one diagonal that carry the same `sjA`.
+    #[test]
+    fn same_diagonal_overlapping_entries_collapse_to_the_longest() {
+        let mut entries = vec![
+            wa(0, 20, 1000, -1),
+            wa(5, 40, 1005, -1),
+            wa(80, 10, 1080, -1),
+        ];
+        dedup_wa_by_diagonal(&mut entries, 100, false);
+
+        assert_eq!(entries.len(), 2);
+        assert_eq!(entries[0].length, 40, "longest of the merged interval wins");
+        assert_eq!(entries[1].read_pos, 80, "the disjoint interval survives");
+    }
+
     fn make_simple_index() -> GenomeIndex {
         // Simple genome: ACGTACGTNN (10 bases)
         let seq = vec![0, 1, 2, 3, 0, 1, 2, 3, 4, 4];
@@ -3245,6 +3553,7 @@ mod tests {
                 is_anchor: true,
                 mate_id: 2,
                 pre_ext_score: 5,
+                sj_a: -1,
             },
             WindowAlignment {
                 seed_idx: 1,
@@ -3256,6 +3565,7 @@ mod tests {
                 is_anchor: true,
                 mate_id: 2,
                 pre_ext_score: 5,
+                sj_a: -1,
             },
         ];
 
@@ -3544,6 +3854,9 @@ mod tests {
             align_spliced_mate_map_lmin_over_lmate: 0.66,
             out_filter_score_min_over_lread: 0.66,
             align_ends_type: crate::params::AlignEndsType::default(),
+            flush_right: false,
+            soft_clip_at_reference_ends: true,
+            p_mm_max_read: 1.0,
         };
 
         // Left overhang (prev.length) = 3, below min of 5
@@ -3588,6 +3901,9 @@ mod tests {
             align_spliced_mate_map_lmin_over_lmate: 0.66,
             out_filter_score_min_over_lread: 0.66,
             align_ends_type: crate::params::AlignEndsType::default(),
+            flush_right: false,
+            soft_clip_at_reference_ends: true,
+            p_mm_max_read: 1.0,
         };
 
         // Both overhangs >= 5
@@ -3698,4 +4014,347 @@ mod tests {
         let replacement_correct = baseline + scorer.sjdb_score;
         assert_ne!(additive_buggy, replacement_correct);
     }
+
+    // -----------------------------------------------------------------
+    // sjAB fast path (STAR "simple stitching if junction belongs to a
+    // database")
+    // -----------------------------------------------------------------
+
+    fn sjab_cluster() -> SeedCluster {
+        SeedCluster {
+            alignments: Vec::new(),
+            chr_idx: 0,
+            genome_start: 0,
+            genome_end: 400,
+            is_reverse: false,
+            anchor_idx: 0,
+            anchor_bin: 0,
+        }
+    }
+
+    /// A working transcript holding one exon that came from junction `sj_a`,
+    /// covering read `[0, exon_len)` at genome `[0, exon_len)`.
+    fn sjab_wt(exon_len: usize, sj_a: i64) -> WorkingTranscript {
+        let mut wt = WorkingTranscript::new();
+        wt.exons.push(ExonBlock {
+            read_start: 0,
+            read_end: exon_len,
+            genome_start: 0,
+            genome_end: exon_len as u64,
+            mate_id: 2,
+            sj_a,
+        });
+        wt.score = exon_len as i32;
+        wt.read_end = exon_len;
+        wt.genome_end = exon_len as u64;
+        wt
+    }
+
+    /// The acceptor-side seed of the same junction: read-adjacent to the exon
+    /// above, but far away on the genome.
+    fn sjab_wa(read_pos: usize, length: usize, sa_pos: u64, sj_a: i64) -> WindowAlignment {
+        WindowAlignment {
+            seed_idx: 1,
+            read_pos,
+            length,
+            genome_pos: sa_pos,
+            sa_pos,
+            n_rep: 1,
+            is_anchor: false,
+            mate_id: 2,
+            pre_ext_score: length as i32,
+            sj_a,
+        }
+    }
+
+    /// A one-entry table whose stored coordinates are deliberately nowhere near
+    /// the seeds being stitched. The `sjA` path addresses entries by *index*,
+    /// trusting the tag the Gsj hit carried, so it still fires; the general
+    /// path looks entries up by *coordinate* and finds nothing. That is what
+    /// makes the two paths distinguishable here.
+    fn sjab_db(motif: u8, shift_left: u8, shift_right: u8) -> crate::junction::SpliceJunctionDb {
+        crate::junction::SpliceJunctionDb::from_prepared(vec![crate::junction::PreparedJunction {
+            chr_idx: 0,
+            start_pos: 50_000,
+            end_pos: 50_180,
+            motif,
+            shift_left,
+            shift_right,
+            strand: 1,
+        }])
+    }
+
+    #[test]
+    fn sjab_shortcut_uses_annotated_junction() {
+        use crate::align::score::{AlignmentScorer, SpliceMotif};
+
+        let index = make_index_with_seq(&[0, 1, 2, 3, 0, 1, 2, 3, 0, 1]);
+        let scorer = AlignmentScorer::from_params_minimal();
+        let db = sjab_db(3, 0, 0); // GC/AG, no micro-repeat
+        let cluster = sjab_cluster();
+
+        let wt = sjab_wt(20, 0);
+        let wa = sjab_wa(20, 30, 200, 0);
+        let out = stitch_align_to_transcript(
+            &wt,
+            &wa,
+            &[0u8; 50],
+            &index,
+            &scorer,
+            &cluster,
+            Some(&db),
+            0,
+            "t",
+        )
+        .expect("annotated shortcut should stitch");
+
+        // Exon B is taken verbatim at the annotated boundary — no shifting.
+        assert_eq!(out.exons.len(), 2);
+        assert_eq!(out.exons[1].read_start, 20);
+        assert_eq!(out.exons[1].genome_start, 200);
+        assert_eq!(out.exons[1].sj_a, 0);
+
+        // The junction is recorded as annotated, with the *stored* motif
+        // rather than one re-derived from the genome.
+        assert_eq!(out.n_junction, 1);
+        assert_eq!(out.junction_annotated, vec![true]);
+        assert_eq!(out.junction_motifs, vec![SpliceMotif::GcAg]);
+
+        // Score is match-per-base plus the annotated bonus, and specifically
+        // not the GC/AG motif penalty the general path would have applied.
+        assert_eq!(out.score, 20 + 30 + scorer.sjdb_score);
+        assert_ne!(out.score, 20 + 30 + scorer.score_gap_gcag);
+    }
+
+    #[test]
+    fn sjab_shortcut_requires_matching_tags_and_adjacency() {
+        use crate::align::score::AlignmentScorer;
+
+        let index = make_index_with_seq(&[0, 1, 2, 3, 0, 1, 2, 3, 0, 1]);
+        let scorer = AlignmentScorer::from_params_minimal();
+        let db = sjab_db(3, 0, 0);
+        let cluster = sjab_cluster();
+        let read = [0u8; 50];
+
+        // The fast path's fingerprint is that it takes the *stored* motif and
+        // leaves exon B exactly where the annotation puts it. The general path
+        // can also mark this junction annotated (it does its own coordinate
+        // lookup), so `junction_annotated` alone would not tell them apart.
+        let signature = |wt: &WorkingTranscript, wa: &WindowAlignment| {
+            stitch_align_to_transcript(wt, wa, &read, &index, &scorer, &cluster, Some(&db), 0, "t")
+                .map(|t| (t.junction_motifs[0], t.exons[1].genome_start))
+        };
+        let fast_path = Some((crate::align::score::SpliceMotif::GcAg, 200u64));
+
+        // Baseline: the shortcut fires.
+        assert_eq!(
+            signature(&sjab_wt(20, 0), &sjab_wa(20, 30, 200, 0)),
+            fast_path
+        );
+
+        // Untagged seed: no shortcut, so the stored GC/AG motif is not adopted.
+        assert_ne!(
+            signature(&sjab_wt(20, 0), &sjab_wa(20, 30, 200, -1)),
+            fast_path
+        );
+        // Tags disagree: the two pieces are flanks of *different* junctions.
+        assert_ne!(
+            signature(&sjab_wt(20, 0), &sjab_wa(20, 30, 200, 1)),
+            fast_path
+        );
+        // Not read-adjacent: a gap in the read is not what the Gsj split makes.
+        assert_ne!(
+            signature(&sjab_wt(20, 0), &sjab_wa(21, 30, 200, 0)),
+            fast_path
+        );
+    }
+
+    #[test]
+    fn noncanonical_annotated_junction_snaps_to_annotated_coords() {
+        use crate::align::score::AlignmentScorer;
+
+        let index = make_index_with_seq(&[0, 1, 2, 3, 0, 1, 2, 3, 0, 1]);
+        let scorer = AlignmentScorer::from_params_minimal();
+        let cluster = sjab_cluster();
+        let read = [0u8; 50];
+        // No `sj_a` tag, so this exercises the general coordinate-lookup path.
+        let wt = sjab_wt(20, -1);
+        let wa = sjab_wa(20, 30, 200, -1);
+
+        let run = |db: Option<&crate::junction::SpliceJunctionDb>| {
+            stitch_align_to_transcript(&wt, &wa, &read, &index, &scorer, &cluster, db, 0, "t")
+        };
+
+        // Calibrate: where does the unannotated boundary search land?
+        let base = run(None).expect("baseline stitch should succeed");
+        let donor_end = base.exons[0].genome_end;
+        let acceptor_start = base.exons[1].genome_start;
+        let intron_len = acceptor_start - donor_end;
+
+        // Annotate exactly that junction as non-canonical with a 2-base
+        // micro-repeat to its left.
+        let shift_left = 2u8;
+        let annotated = crate::junction::SpliceJunctionDb::from_prepared(vec![
+            crate::junction::PreparedJunction {
+                chr_idx: 0,
+                start_pos: donor_end,
+                end_pos: acceptor_start - 1,
+                motif: 0,
+                shift_left,
+                shift_right: 0,
+                strand: 0,
+            },
+        ]);
+        let snapped = run(Some(&annotated)).expect("annotated stitch should succeed");
+
+        // Both boundaries move right by shift_left; the intron keeps its length.
+        assert_eq!(snapped.exons[0].genome_end, donor_end + shift_left as u64);
+        assert_eq!(
+            snapped.exons[1].genome_start,
+            acceptor_start + shift_left as u64
+        );
+        assert_eq!(
+            snapped.exons[1].genome_start - snapped.exons[0].genome_end,
+            intron_len,
+            "snapping must relocate the junction, not resize the intron"
+        );
+
+        // It is recorded as annotated, with the stored shifts.
+        assert_eq!(snapped.junction_annotated, vec![true]);
+        assert_eq!(snapped.junction_shifts, vec![(shift_left as u32, 0)]);
+
+        // A canonical annotated junction at the same coordinates, carrying the
+        // same micro-repeat, is *not* snapped: STAR only shifts non-canonical
+        // entries. Canonical entries store `start_pos + shift_left`, so the
+        // pre-shift positions are backed off to land on the same stored key.
+        let canonical = crate::junction::SpliceJunctionDb::from_prepared(vec![
+            crate::junction::PreparedJunction {
+                chr_idx: 0,
+                start_pos: donor_end - shift_left as u64,
+                end_pos: acceptor_start - 1 - shift_left as u64,
+                motif: 1,
+                shift_left,
+                shift_right: 0,
+                strand: 1,
+            },
+        ]);
+        let unsnapped = run(Some(&canonical)).expect("canonical annotated stitch should succeed");
+        assert_eq!(unsnapped.exons[0].genome_end, donor_end);
+        assert_eq!(unsnapped.exons[1].genome_start, acceptor_start);
+        assert_eq!(unsnapped.junction_annotated, vec![true]);
+    }
+
+    /// A gap shorter than `--alignIntronMin` is a deletion unless the
+    /// annotation knows it, in which case it is a junction however short.
+    ///
+    /// STAR resolves `sjdbInd` before testing `Del >= alignIntronMin`
+    /// (`stitchAlignToTranscript.cpp:198-241`), so the length test only decides
+    /// for gaps the annotation has never heard of. Yeast R64-1-1 has three
+    /// one-base introns; before this, reads crossing them scored as deletions
+    /// and came out `1D` where STAR writes `1N`.
+    #[test]
+    fn a_short_annotated_gap_is_a_junction_not_a_deletion() {
+        use crate::align::score::AlignmentScorer;
+
+        let index = make_index_with_seq(&[0, 1, 2, 3, 0, 1, 2, 3, 0, 1]);
+        let scorer = AlignmentScorer::from_params_minimal();
+        assert!(
+            scorer.align_intron_min > 5,
+            "the fixture's 5-base gap must be below the intron threshold"
+        );
+        let cluster = sjab_cluster();
+        let read = [0u8; 50];
+
+        // Exon A covers read[0,20) at genome[0,20); the next seed sits five
+        // bases further along the genome, so the gap is a 5-base deletion.
+        let wt = sjab_wt(20, -1);
+        let wa = sjab_wa(20, 20, 25, -1);
+
+        let plain =
+            stitch_align_to_transcript(&wt, &wa, &read, &index, &scorer, &cluster, None, 0, "t")
+                .expect("unannotated stitch should succeed");
+        assert_eq!(plain.n_gap, 1, "unannotated: a deletion");
+        assert_eq!(plain.n_junction, 0);
+
+        // The same gap, annotated.
+        let donor_end = plain.exons[0].genome_end;
+        let acceptor_start = plain.exons[1].genome_start;
+        let db = crate::junction::SpliceJunctionDb::from_prepared(vec![
+            crate::junction::PreparedJunction {
+                chr_idx: 0,
+                start_pos: donor_end,
+                end_pos: acceptor_start - 1,
+                motif: 1,
+                shift_left: 0,
+                shift_right: 0,
+                strand: 1,
+            },
+        ]);
+        let annotated = stitch_align_to_transcript(
+            &wt,
+            &wa,
+            &read,
+            &index,
+            &scorer,
+            &cluster,
+            Some(&db),
+            0,
+            "t",
+        )
+        .expect("annotated stitch should succeed");
+        assert_eq!(annotated.n_junction, 1, "annotated: a junction");
+        assert_eq!(annotated.n_gap, 0);
+        assert_eq!(annotated.junction_annotated, vec![true]);
+        assert!(
+            annotated.score > plain.score,
+            "the annotated form earns sjdbScore rather than a deletion penalty: \
+             {} vs {}",
+            annotated.score,
+            plain.score
+        );
+    }
+
+    #[test]
+    fn too_large_repeat_around_annotated_junction_rejects() {
+        use crate::align::score::AlignmentScorer;
+
+        let index = make_index_with_seq(&[0, 1, 2, 3, 0, 1, 2, 3, 0, 1]);
+        let scorer = AlignmentScorer::from_params_minimal();
+        let cluster = sjab_cluster();
+        let read = [0u8; 50];
+
+        // Non-canonical annotated junction whose micro-repeats swallow the
+        // pieces being stitched: STAR rejects with -1000006, we return None.
+        let db = sjab_db(0, 40, 40);
+        let out = stitch_align_to_transcript(
+            &sjab_wt(20, 0),
+            &sjab_wa(20, 30, 200, 0),
+            &read,
+            &index,
+            &scorer,
+            &cluster,
+            Some(&db),
+            0,
+            "t",
+        );
+        assert!(
+            out.is_none(),
+            "a seed shorter than shiftRight must reject the annotated shortcut"
+        );
+
+        // Same junction with small repeats: accepted.
+        let db_ok = sjab_db(0, 1, 1);
+        let out_ok = stitch_align_to_transcript(
+            &sjab_wt(20, 0),
+            &sjab_wa(20, 30, 200, 0),
+            &read,
+            &index,
+            &scorer,
+            &cluster,
+            Some(&db_ok),
+            0,
+            "t",
+        );
+        assert!(out_ok.is_some());
+    }
 }
diff --git a/src/chimeric/detect.rs b/src/chimeric/detect.rs
index 4997049..33d564c 100644
--- a/src/chimeric/detect.rs
+++ b/src/chimeric/detect.rs
@@ -1112,6 +1112,7 @@ mod tests {
                 is_anchor: true,
                 mate_id: 2,
                 pre_ext_score: (genome_end - genome_start) as i32,
+                sj_a: -1,
             }],
             chr_idx,
             genome_start,
diff --git a/src/index/io.rs b/src/index/io.rs
index ce6de64..2385f39 100644
--- a/src/index/io.rs
+++ b/src/index/io.rs
@@ -69,7 +69,7 @@ impl GenomeIndex {
         //     db would be empty), losing all `sjdbScore` bonuses and annotated
         //     junction recognition. Keyed on the stored (post-sjdbPrepare) donor/
         //     acceptor coordinates, matching what the stitch scan produces.
-        let junction_db = if let Some(ref gtf_path) = params.sjdb_gtf_file {
+        let mut junction_db = if let Some(ref gtf_path) = params.sjdb_gtf_file {
             SpliceJunctionDb::from_gtf_configured(
                 gtf_path,
                 &genome,
@@ -78,20 +78,27 @@ impl GenomeIndex {
                 &params.sjdb_gtf_tag_exon_parent_transcript,
             )?
         } else if !prepared_junctions.is_empty() {
-            let raw: Vec<(usize, u64, u64, u8)> = prepared_junctions
-                .iter()
-                .map(|j| (j.chr_idx, j.stored_start(), j.stored_end(), j.strand))
-                .collect();
             log::info!(
                 "No GTF at align time; loaded {} annotated junctions from index sjdbInfo.txt",
-                raw.len()
+                prepared_junctions.len()
             );
-            SpliceJunctionDb::from_raw_junctions(&raw)
+            SpliceJunctionDb::from_prepared(prepared_junctions.clone())
         } else {
             log::info!("No GTF file provided, all junctions will be novel");
             SpliceJunctionDb::empty()
         };
 
+        // `sj_a` tags come out of `decode_gsj_hit`, which indexes the junction
+        // array stored in the index. Whenever that array exists it must also be
+        // what `find()` searches, otherwise a tag would address a different
+        // junction than the one the SA hit came from — silently wrong CIGARs
+        // rather than a missed optimisation. An align-time GTF may rebuild the
+        // annotated lookup map, but never the table.
+        if !prepared_junctions.is_empty() {
+            junction_db.set_table(prepared_junctions.clone());
+        }
+        let junction_db = junction_db;
+
         log::info!(
             "Junction database loaded: {} annotated junctions",
             junction_db.len()
diff --git a/src/index/mod.rs b/src/index/mod.rs
index 88ee62d..bfccbf9 100644
--- a/src/index/mod.rs
+++ b/src/index/mod.rs
@@ -327,8 +327,6 @@ impl GenomeIndex {
         }
 
         let (junction_db, prepared_junctions) = if !raw.is_empty() {
-            let jdb = SpliceJunctionDb::from_raw_junctions(&raw);
-
             let prepared: Vec<PreparedJunction> = raw
                 .iter()
                 .map(|&(chr_idx, intron_start, intron_end, strand)| {
@@ -343,6 +341,10 @@ impl GenomeIndex {
                 })
                 .collect();
             let prepared = sjdb_insert::sort_and_dedup(prepared);
+            // Build the annotated db from the prepared array itself, so the
+            // `sjA` table and the Gsj buffer are the same list in the same
+            // order by construction.
+            let jdb = SpliceJunctionDb::from_prepared(prepared.clone());
 
             let gsj =
                 sjdb_insert::build_gsj(&prepared, &genome, n_genome_real, params.sjdb_overhang)?;
diff --git a/src/junction/gtf.rs b/src/junction/gtf.rs
index 411b4f5..4411341 100644
--- a/src/junction/gtf.rs
+++ b/src/junction/gtf.rs
@@ -204,9 +204,16 @@ pub fn extract_junctions_configured(
             let intron_start_local_1b = exon1.end + 1;
             let intron_end_local_1b = exon2.start.saturating_sub(1);
 
-            if intron_end_local_1b <= intron_start_local_1b {
+            // STAR skips a pair only when the exons touch or overlap
+            // (`GTF_transcriptGeneSJ.cpp:123`, `exonLoci[iex][exS] <=
+            // exonLoci[iex-1][exE] + 1`). A one-base intron — `exon2.start ==
+            // exon1.end + 2` — is a junction it keeps, and yeast has three of
+            // them. Requiring `end > start` here dropped exactly those, so the
+            // aligner never saw them as annotated and emitted `1D` where STAR
+            // emits `1N`.
+            if intron_end_local_1b < intron_start_local_1b {
                 log::warn!(
-                    "Invalid junction coordinates: {intron_start_local_1b}-{intron_end_local_1b} (possibly overlapping exons)"
+                    "Invalid junction coordinates: {intron_start_local_1b}-{intron_end_local_1b} (overlapping exons)"
                 );
                 continue;
             }
@@ -360,6 +367,54 @@ mod tests {
         assert_eq!(strand, 1);
     }
 
+    /// A one-base intron is a junction, and touching exons are not.
+    ///
+    /// STAR skips a pair only when `exon2.start <= exon1.end + 1`
+    /// (`GTF_transcriptGeneSJ.cpp:123`), so `exon2.start == exon1.end + 2`
+    /// yields a one-base intron it keeps. Yeast R64-1-1 has three; requiring
+    /// `end > start` here dropped all three from the index, and reads crossing
+    /// them then came out with `1D` where STAR writes `1N`.
+    #[test]
+    fn a_one_base_intron_is_a_junction_and_touching_exons_are_not() {
+        let genome = Genome {
+            transform_blocks: None,
+            sequence: vec![0; 1000].into(),
+            n_genome: 1000,
+            n_genome_real: 1000,
+            n_chr_real: 1,
+            chr_start: vec![0, 1000],
+            chr_length: vec![1000],
+            chr_name: vec!["chr1".to_string()],
+        };
+        let exon = |start: u64, end: u64, transcript: &str| GtfRecord {
+            seqname: "chr1".to_string(),
+            feature: "exon".to_string(),
+            start,
+            end,
+            strand: '+',
+            attributes: vec![
+                ("gene_id".to_string(), "G1".to_string()),
+                ("transcript_id".to_string(), transcript.to_string()),
+            ]
+            .into_iter()
+            .collect(),
+        };
+
+        // T1: exons 100-200 and 202-300 — a one-base intron at 201.
+        let junctions =
+            extract_junctions_from_exons(vec![exon(100, 200, "T1"), exon(202, 300, "T1")], &genome)
+                .unwrap();
+        assert_eq!(junctions.len(), 1, "the one-base intron must survive");
+        let (_, start, end, _) = junctions[0];
+        assert_eq!((start, end), (200, 200), "0-based, inclusive, length one");
+
+        // T2: exons 100-200 and 201-300 — adjacent, no intron between them.
+        let touching =
+            extract_junctions_from_exons(vec![exon(100, 200, "T2"), exon(201, 300, "T2")], &genome)
+                .unwrap();
+        assert!(touching.is_empty(), "touching exons are not a junction");
+    }
+
     #[test]
     fn test_extract_junctions_multiple_transcripts() {
         let genome = Genome {
diff --git a/src/junction/mod.rs b/src/junction/mod.rs
index 13776ce..2761e80 100644
--- a/src/junction/mod.rs
+++ b/src/junction/mod.rs
@@ -11,7 +11,8 @@ mod sj_output;
 pub mod sjdb_insert;
 
 pub use sj_output::SpliceJunctionStats;
-pub(crate) use sj_output::{SjKey, encode_motif};
+pub(crate) use sj_output::{SjKey, decode_motif, encode_motif};
+pub use sjdb_insert::PreparedJunction;
 
 use crate::params::Parameters;
 
@@ -55,6 +56,15 @@ pub struct NovelJunctionKey {
 pub struct SpliceJunctionDb {
     /// Map: (chr_idx, intron_start, intron_end, strand) → annotated
     junctions: HashMap<JunctionKey, JunctionInfo>,
+    /// STAR's `mapGen.sjdb*` arrays: the annotated junctions in the exact
+    /// order they occupy the Gsj buffer, sorted by `(stored_start,
+    /// stored_end)`. Index `i` here is STAR's `sjA` / `sjdbInd`, i.e. the
+    /// value [`sjdb_insert::decode_gsj_hit`] derives from a Gsj SA hit.
+    ///
+    /// Empty when the database was built without motif/shift metadata; in
+    /// that case [`find`](Self::find) always returns `None` and the
+    /// annotated-junction fast paths simply do not engage.
+    table: Vec<PreparedJunction>,
 }
 
 impl SpliceJunctionDb {
@@ -62,7 +72,102 @@ impl SpliceJunctionDb {
     pub fn empty() -> Self {
         Self {
             junctions: HashMap::new(),
+            table: Vec::new(),
+        }
+    }
+
+    /// Build from an already sorted-and-deduplicated [`PreparedJunction`]
+    /// list — the form `sjdbInfo.txt` is read back into and the form
+    /// `genomeGenerate` produces before building the Gsj buffer.
+    ///
+    /// The HashMap is keyed on the *stored* (post-`sjdbPrepare`) donor and
+    /// acceptor coordinates, which is what the stitch-time scan produces.
+    pub fn from_prepared(prepared: Vec<PreparedJunction>) -> Self {
+        let mut junctions = HashMap::with_capacity(prepared.len());
+        for j in &prepared {
+            junctions.insert(
+                JunctionKey {
+                    chr_idx: j.chr_idx,
+                    intron_start: j.stored_start(),
+                    intron_end: j.stored_end(),
+                    strand: j.strand,
+                },
+                JunctionInfo { annotated: true },
+            );
+        }
+        Self {
+            junctions,
+            table: prepared,
+        }
+    }
+
+    /// Replace the `sjA`-addressable table without touching the annotated
+    /// lookup map.
+    ///
+    /// Needed because `sj_a` tags come from [`sjdb_insert::decode_gsj_hit`],
+    /// which indexes the junction list stored in the index (`sjdbInfo.txt`).
+    /// When a GTF is *also* supplied at align time the map is rebuilt from
+    /// that GTF, but the table must keep addressing the index's array or the
+    /// tags would point at the wrong junctions.
+    pub fn set_table(&mut self, prepared: Vec<PreparedJunction>) {
+        self.table = prepared;
+    }
+
+    /// STAR's `binarySearch2` (`ReadAlign_stitchAlignToTranscript.cpp` →
+    /// `binarySearch2.h`): the index of the junction whose stored donor and
+    /// acceptor are exactly `x` and `y`, or `None`.
+    ///
+    /// Coordinates are genome-absolute, so no chromosome index is needed:
+    /// STAR's `sjdbStart` / `sjdbEnd` are already unique across contigs.
+    pub fn find(&self, x: u64, y: u64) -> Option<usize> {
+        let n = self.table.len();
+        if n == 0 || x > self.table[n - 1].stored_start() || x < self.table[0].stored_start() {
+            return None;
         }
+        let (mut i1, mut i2) = (0usize, n - 1);
+        while i2 > i1 + 1 {
+            let i3 = usize::midpoint(i1, i2);
+            if self.table[i3].stored_start() > x {
+                i2 = i3;
+            } else {
+                i1 = i3;
+            }
+        }
+        let i3 = if x == self.table[i1].stored_start() {
+            i1
+        } else if x == self.table[i2].stored_start() {
+            i2
+        } else {
+            return None;
+        };
+        // Scan the run of equal `stored_start` values (backward then forward)
+        // for a matching `stored_end`.
+        for jj in (0..=i3).rev() {
+            if x != self.table[jj].stored_start() {
+                break;
+            } else if y == self.table[jj].stored_end() {
+                return Some(jj);
+            }
+        }
+        for jj in i3..n {
+            if x != self.table[jj].stored_start() {
+                return None;
+            } else if y == self.table[jj].stored_end() {
+                return Some(jj);
+            }
+        }
+        None
+    }
+
+    /// The junction at `sjA` index `i`, or `None` when the table is absent
+    /// or the index is out of range.
+    pub fn entry(&self, i: usize) -> Option<&PreparedJunction> {
+        self.table.get(i)
+    }
+
+    /// Number of entries in the `sjA`-addressable table.
+    pub fn table_len(&self) -> usize {
+        self.table.len()
     }
 
     /// Build junction database from GTF file with configurable GTF attribute names.
@@ -81,7 +186,36 @@ impl SpliceJunctionDb {
         let raw = gtf::extract_junctions_configured(exons, genome, transcript_tag)?;
         log::info!("Extracted {} annotated junctions from GTF", raw.len());
 
-        Ok(Self::from_raw_junctions(&raw))
+        // Keep the historical (raw-coordinate) annotated lookup map, but also
+        // derive the motif/shift/strand table so the annotated-junction fast
+        // paths have something to address. Building the table here makes the
+        // align-time GTF path carry the same metadata the index-loaded path
+        // already had.
+        let mut db = Self::from_raw_junctions(&raw);
+        db.set_table(Self::prepare_table(&raw, genome));
+        Ok(db)
+    }
+
+    /// Run every raw `(chr_idx, intron_start, intron_end, strand)` through
+    /// `sjdbPrepare`'s motif detection and micro-repeat shift computation,
+    /// then apply STAR's post-dedup sort. Produces exactly the array
+    /// `genomeGenerate` writes to `sjdbInfo.txt`.
+    fn prepare_table(raw: &[(usize, u64, u64, u8)], genome: &Genome) -> Vec<PreparedJunction> {
+        let n_genome_real = genome.n_genome_real;
+        let prepared: Vec<PreparedJunction> = raw
+            .iter()
+            .map(|&(chr_idx, intron_start, intron_end, strand)| {
+                sjdb_insert::prepare_junction(
+                    chr_idx,
+                    intron_start,
+                    intron_end,
+                    strand,
+                    genome,
+                    n_genome_real,
+                )
+            })
+            .collect();
+        sjdb_insert::sort_and_dedup(prepared)
     }
 
     /// Build junction database from GTF file (default STAR attribute names).
@@ -94,6 +228,10 @@ impl SpliceJunctionDb {
     /// the `genomeGenerate` path so it can share the parsed GTF with
     /// `TranscriptomeIndex` and the `sjdb_insert` pipeline without
     /// re-parsing the file.
+    ///
+    /// The resulting database has **no** `sjA` table; callers that need
+    /// [`find`](Self::find) must follow up with [`set_table`](Self::set_table)
+    /// or use [`from_prepared`](Self::from_prepared) instead.
     pub fn from_raw_junctions(raw: &[(usize, u64, u64, u8)]) -> Self {
         let mut junctions = HashMap::with_capacity(raw.len());
         for &(chr_idx, intron_start, intron_end, strand) in raw {
@@ -105,7 +243,10 @@ impl SpliceJunctionDb {
             };
             junctions.insert(key, JunctionInfo { annotated: true });
         }
-        Self { junctions }
+        Self {
+            junctions,
+            table: Vec::new(),
+        }
     }
 
     /// Check if a junction is annotated in the GTF.
@@ -408,6 +549,104 @@ mod tests {
         assert_eq!(novel_junctions[0].0.intron_start, 300);
     }
 
+    /// Build a canonical (motif != 0) prepared junction whose stored
+    /// coordinates are exactly `(start, end)`.
+    fn pj(start: u64, end: u64, strand: u8) -> PreparedJunction {
+        PreparedJunction {
+            chr_idx: 0,
+            start_pos: start,
+            end_pos: end,
+            motif: 1,
+            shift_left: 0,
+            shift_right: 0,
+            strand,
+        }
+    }
+
+    #[test]
+    fn find_round_trips_every_prepared_junction() {
+        // Deliberately unsorted input: `sort_and_dedup` establishes the
+        // (stored_start, stored_end) order `find`'s binary search needs.
+        let prepared = sjdb_insert::sort_and_dedup(vec![
+            pj(900, 1000, 1),
+            pj(100, 200, 1),
+            pj(500, 600, 2),
+            pj(300, 400, 1),
+        ]);
+        let db = SpliceJunctionDb::from_prepared(prepared.clone());
+        assert_eq!(db.table_len(), 4);
+
+        // The invariant the sjAB fast path and the annotated snap both rely
+        // on: index i in the table is addressable by its own stored coords.
+        for (i, j) in prepared.iter().enumerate() {
+            assert_eq!(
+                db.find(j.stored_start(), j.stored_end()),
+                Some(i),
+                "junction {i} did not round-trip"
+            );
+            assert_eq!(db.entry(i).unwrap().stored_start(), j.stored_start());
+        }
+
+        // Misses in every direction: below the first, above the last, a
+        // start that exists with a wrong end, and an end that exists with a
+        // wrong start.
+        assert_eq!(db.find(50, 200), None);
+        assert_eq!(db.find(2000, 3000), None);
+        assert_eq!(db.find(100, 201), None);
+        assert_eq!(db.find(101, 200), None);
+    }
+
+    #[test]
+    fn find_disambiguates_a_run_of_equal_starts() {
+        // Several junctions sharing a donor: the backward/forward scan around
+        // the binary-search landing point must pick the right acceptor.
+        let prepared = sjdb_insert::sort_and_dedup(vec![
+            pj(100, 700, 1),
+            pj(100, 200, 1),
+            pj(100, 500, 1),
+            pj(100, 300, 1),
+            pj(900, 1000, 1),
+        ]);
+        let db = SpliceJunctionDb::from_prepared(prepared.clone());
+        for (i, j) in prepared.iter().enumerate() {
+            assert_eq!(db.find(j.stored_start(), j.stored_end()), Some(i));
+        }
+        assert_eq!(db.find(100, 400), None);
+    }
+
+    #[test]
+    fn find_is_inert_without_a_table() {
+        // `from_raw_junctions` carries no motif/shift metadata, so the
+        // annotated fast paths must simply not engage rather than misfire.
+        let db = SpliceJunctionDb::from_raw_junctions(&[(0, 100, 200, 1)]);
+        assert!(db.is_annotated(0, 100, 200, 1));
+        assert_eq!(db.table_len(), 0);
+        assert_eq!(db.find(100, 200), None);
+        assert!(db.entry(0).is_none());
+    }
+
+    #[test]
+    fn from_prepared_keys_the_map_on_stored_coordinates() {
+        // Non-canonical junction: stored coords are the shifted ones, which is
+        // what the stitch-time scan produces.
+        let noncan = PreparedJunction {
+            chr_idx: 0,
+            start_pos: 100,
+            end_pos: 200,
+            motif: 0,
+            shift_left: 3,
+            shift_right: 0,
+            strand: 0,
+        };
+        assert_eq!(noncan.stored_start(), 100);
+        assert_eq!(noncan.original_start(), 103);
+
+        let db = SpliceJunctionDb::from_prepared(vec![noncan]);
+        assert!(db.is_annotated(0, 100, 200, 0));
+        assert!(!db.is_annotated(0, 103, 203, 0));
+        assert_eq!(db.find(100, 200), Some(0));
+    }
+
     #[test]
     fn test_db_keyed_in_genome_absolute_zero_based_multi_chr() {
         use crate::junction::gtf::{GtfRecord, extract_junctions_configured};
diff --git a/src/junction/sj_output.rs b/src/junction/sj_output.rs
index 9fa41dd..65bc577 100644
--- a/src/junction/sj_output.rs
+++ b/src/junction/sj_output.rs
@@ -401,6 +401,21 @@ pub(crate) fn encode_motif(motif: SpliceMotif) -> u8 {
     }
 }
 
+/// Inverse of [`encode_motif`]: turn a stored `sjdbMotif` code back into a
+/// [`SpliceMotif`]. Out-of-range codes map to `NonCanonical`, matching
+/// [`SpliceMotif::filter_category_from_encoded`].
+pub(crate) fn decode_motif(encoded: u8) -> SpliceMotif {
+    match encoded {
+        1 => SpliceMotif::GtAg,
+        2 => SpliceMotif::CtAc,
+        3 => SpliceMotif::GcAg,
+        4 => SpliceMotif::CtGc,
+        5 => SpliceMotif::AtAc,
+        6 => SpliceMotif::GtAt,
+        _ => SpliceMotif::NonCanonical,
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use crate::params::Parameters;
diff --git a/src/junction/sjdb_insert.rs b/src/junction/sjdb_insert.rs
index e358622..5b04b17 100644
--- a/src/junction/sjdb_insert.rs
+++ b/src/junction/sjdb_insert.rs
@@ -140,7 +140,8 @@ pub fn read_sjdb_info_tab(path: &Path, genome: &Genome) -> Result<SjdbInfoTab, E
 /// should use the original `(fwd_pos, length)` directly and not invoke
 /// this function.
 ///
-/// Returns a vector of `(real_fwd_pos, read_offset, length)` tuples:
+/// Returns a vector of `(real_fwd_pos, read_offset, length, junction_idx)`
+/// tuples:
 /// - **One** entry when the hit lies entirely within a single junction's
 ///   donor flank or acceptor flank — the read↔genome correspondence is a
 ///   single contiguous run with `read_offset = 0` and `length` unchanged.
@@ -151,6 +152,11 @@ pub fn read_sjdb_info_tab(path: &Path, genome: &Genome) -> Result<SjdbInfoTab, E
 ///   read window but at non-adjacent genome positions — the downstream
 ///   stitch DP can then chain them via its splice branch.
 ///
+/// `junction_idx` is the entry's index in `junctions` — STAR's `sjA`, the tag
+/// its stitcher carries on every align derived from the sjdb region. Both
+/// halves of a boundary-crossing hit share the same index, since they are the
+/// two flanks of one junction.
+///
 /// Returns an empty vector when the hit's Gsj junction index is out of
 /// bounds, when the hit's slot offset lands on the trailing spacer byte,
 /// or when the hit would extend past the spacer.
@@ -160,7 +166,7 @@ pub fn decode_gsj_hit(
     n_genome_real: u64,
     sjdb_overhang: u32,
     junctions: &[PreparedJunction],
-) -> Vec<(u64, usize, usize)> {
+) -> Vec<(u64, usize, usize, usize)> {
     if fwd_pos < n_genome_real {
         return Vec::new();
     }
@@ -200,8 +206,8 @@ pub fn decode_gsj_hit(
         let donor_real = donor_genome_start + slot_offset;
         let acceptor_real = acceptor_genome_start;
         vec![
-            (donor_real, 0, donor_len),
-            (acceptor_real, donor_len, acceptor_len),
+            (donor_real, 0, donor_len, junction_idx),
+            (acceptor_real, donor_len, acceptor_len, junction_idx),
         ]
     }
 }
@@ -1052,9 +1058,9 @@ mod tests {
         let out = decode_gsj_hit(hit_pos, 12, n_genome_real, overhang, &junctions);
         assert_eq!(out.len(), 2);
         // Donor: real_fwd = 990 + 6 = 996, read_offset = 0, sub_len = 4.
-        assert_eq!(out[0], (996, 0, 4));
+        assert_eq!(out[0], (996, 0, 4, 0));
         // Acceptor: real_fwd = 2001, read_offset = 4, sub_len = 8.
-        assert_eq!(out[1], (2001, 4, 8));
+        assert_eq!(out[1], (2001, 4, 8, 0));
     }
 
     #[test]
@@ -1086,9 +1092,9 @@ mod tests {
         let out = decode_gsj_hit(hit_pos, 10, n_genome_real, overhang, &junctions);
         assert_eq!(out.len(), 2);
         // original_start = 1003 → donor flank [993..1003). Hit at slot 5 = 993+5 = 998.
-        assert_eq!(out[0], (998, 0, 5));
+        assert_eq!(out[0], (998, 0, 5, 0));
         // original_end = 2003 → acceptor flank starts at 2004.
-        assert_eq!(out[1], (2004, 5, 5));
+        assert_eq!(out[1], (2004, 5, 5, 0));
     }
 
     #[test]
diff --git a/src/params/mod.rs b/src/params/mod.rs
index 0536a85..7b38133 100644
--- a/src/params/mod.rs
+++ b/src/params/mod.rs
@@ -55,6 +55,8 @@ impl std::str::FromStr for RunMode {
     }
 }
 
+// ---------------------------------------------------------------------------
+
 impl std::fmt::Display for RunMode {
     fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
         match self {
@@ -713,6 +715,12 @@ pub struct Parameters {
     #[arg(long = "outFilterMismatchNoverLmax", default_value_t = 0.3)]
     pub out_filter_mismatch_nover_lmax: f64,
 
+    /// Max mismatches per pair as a fraction of *read* length (both mates
+    /// summed), as opposed to `outFilterMismatchNoverLmax`, which is a
+    /// fraction of the *mapped* length.
+    #[arg(long = "outFilterMismatchNoverReadLmax", default_value_t = 1.0)]
+    pub out_filter_mismatch_nover_read_lmax: f64,
+
     /// Min alignment score (absolute)
     #[arg(long = "outFilterScoreMin", default_value_t = 0)]
     pub out_filter_score_min: i32,
@@ -818,6 +826,22 @@ pub struct Parameters {
           default_values_t = vec![0, -1, 0, 0], allow_hyphen_values = true)]
     pub align_sj_stitch_mismatch_nmax: Vec<i32>,
 
+    /// Allow protrusion of one mate's alignment past the other's start:
+    /// `ConcordantPair` or `<Nbases> ConcordantPair`.
+    #[arg(long = "alignEndsProtrude", num_args = 1..=2,
+          default_values_t = vec!["0".to_string(), "ConcordantPair".to_string()])]
+    pub align_ends_protrude: Vec<String>,
+
+    /// Allow the soft-clipping of alignments past the ends of chromosomes
+    /// (`Yes`, default) or forbid it (`No`).
+    #[arg(long = "alignSoftClipAtReferenceEnds", default_value = "Yes")]
+    pub align_soft_clip_at_reference_ends: String,
+
+    /// Which end of a repeat-ambiguous insertion the insertion is flushed to:
+    /// `None` (default, no flushing) or `Right`.
+    #[arg(long = "alignInsertionFlush", default_value = "None")]
+    pub align_insertion_flush: String,
+
     /// Splice junction penalty (canonical)
     #[arg(long = "scoreGap", default_value_t = 0)]
     pub score_gap: i32,
@@ -887,6 +911,14 @@ pub struct Parameters {
     #[arg(long = "seedPerWindowNmax", default_value_t = 50)]
     pub seed_per_window_nmax: usize,
 
+    /// Max number of one-seed loci per window.
+    #[arg(long = "seedNoneLociPerWindow", default_value_t = 10)]
+    pub seed_none_loci_per_window: usize,
+
+    /// Min length of the seed sequences split by the seed search.
+    #[arg(long = "seedSplitMin", default_value_t = 12)]
+    pub seed_split_min: usize,
+
     /// Max distance between seed search start positions (defines Nstart = readLen/seedSearchStartLmax + 1)
     #[arg(long = "seedSearchStartLmax", default_value_t = 50)]
     pub seed_search_start_lmax: usize,
@@ -915,6 +947,12 @@ pub struct Parameters {
     #[arg(long = "alignTranscriptsPerWindowNmax", default_value_t = 100)]
     pub align_transcripts_per_window_nmax: usize,
 
+    /// Max number of different alignments per read to consider. STAR stops
+    /// collecting once the running total plus one window's worth would exceed
+    /// this.
+    #[arg(long = "alignTranscriptsPerReadNmax", default_value_t = 10000)]
+    pub align_transcripts_per_read_nmax: usize,
+
     // ── Splice junction database ────────────────────────────────────────
     /// GTF file for splice junction annotations
     #[arg(long = "sjdbGTFfile")]
@@ -1551,6 +1589,45 @@ impl Parameters {
             ));
         }
 
+        // Validate --alignSoftClipAtReferenceEnds.
+        if !matches!(
+            params.align_soft_clip_at_reference_ends.as_str(),
+            "Yes" | "No"
+        ) {
+            return Err(command.error(
+                ErrorKind::InvalidValue,
+                format!(
+                    "unknown --alignSoftClipAtReferenceEnds '{}'; expected Yes or No",
+                    params.align_soft_clip_at_reference_ends
+                ),
+            ));
+        }
+        // Validate --alignInsertionFlush.
+        if !matches!(params.align_insertion_flush.as_str(), "None" | "Right") {
+            return Err(command.error(
+                ErrorKind::InvalidValue,
+                format!(
+                    "unknown --alignInsertionFlush '{}'; expected None or Right",
+                    params.align_insertion_flush
+                ),
+            ));
+        }
+        // Validate --alignEndsProtrude: an optional base count followed by
+        // the protrusion type.
+        {
+            let p = &params.align_ends_protrude;
+            let kind = p.last().map_or("ConcordantPair", String::as_str);
+            let n_ok = p.len() < 2 || p[0].parse::<u32>().is_ok();
+            if !n_ok || !matches!(kind, "ConcordantPair") {
+                return Err(command.error(
+                    ErrorKind::InvalidValue,
+                    format!(
+                        "invalid --alignEndsProtrude '{}'; expected [<Nbases>] ConcordantPair",
+                        p.join(" ")
+                    ),
+                ));
+            }
+        }
         // ── STARsolo validation ─────────────────────────────────────────
         if params.run_mode == RunMode::AlignReads && params.solo_enabled() {
             // CB_UMI_Complex needs one CB position + whitelist per segment.
@@ -1881,6 +1958,53 @@ mod tests {
         Parameters::try_parse_from(&full)
     }
 
+    #[test]
+    fn new_align_knobs_parse_and_reject() {
+        let p = try_parse(&["--readFilesIn", "r.fq"]).unwrap();
+        assert_eq!(p.align_ends_type, "Local");
+        assert_eq!(p.align_soft_clip_at_reference_ends, "Yes");
+        assert_eq!(p.align_insertion_flush, "None");
+        assert_eq!(p.align_transcripts_per_read_nmax, 10000);
+        assert!((p.out_filter_mismatch_nover_read_lmax - 1.0).abs() < f64::EPSILON);
+        assert_eq!(p.seed_none_loci_per_window, 10);
+        assert_eq!(p.seed_split_min, 12);
+
+        assert!(try_parse(&["--readFilesIn", "r.fq", "--alignEndsType", "EndToEnd"]).is_ok());
+        assert!(try_parse(&["--readFilesIn", "r.fq", "--alignInsertionFlush", "Right"]).is_ok());
+        assert!(
+            try_parse(&[
+                "--readFilesIn",
+                "r.fq",
+                "--alignSoftClipAtReferenceEnds",
+                "No"
+            ])
+            .is_ok()
+        );
+        assert!(
+            try_parse(&[
+                "--readFilesIn",
+                "r.fq",
+                "--alignEndsProtrude",
+                "5",
+                "ConcordantPair"
+            ])
+            .is_ok()
+        );
+
+        // Off-menu values are rejected loudly rather than silently ignored.
+        assert!(try_parse(&["--readFilesIn", "r.fq", "--alignEndsType", "Global"]).is_err());
+        assert!(try_parse(&["--readFilesIn", "r.fq", "--alignInsertionFlush", "Left"]).is_err());
+        assert!(
+            try_parse(&[
+                "--readFilesIn",
+                "r.fq",
+                "--alignSoftClipAtReferenceEnds",
+                "Maybe"
+            ])
+            .is_err()
+        );
+    }
+
     #[test]
     fn defaults() {
         let p = try_parse(&["--readFilesIn", "reads.fq"]).unwrap();
diff --git a/test/run_tests.sh b/test/run_tests.sh
index fbfa7ba..5b75599 100755
--- a/test/run_tests.sh
+++ b/test/run_tests.sh
@@ -44,6 +44,15 @@ TEST_CASES=(
     "yeast_1k_paired:yeast:ERR12389696_sub_1_1000.fastq.gz,ERR12389696_sub_2_1000.fastq.gz:paired:--outSAMtype SAM"
     "yeast_1k_twopass:yeast:ERR12389696_sub_1_1000.fastq.gz:single:--twopassMode Basic --outSAMtype SAM"
     "yeast_10k:yeast:ERR12389696_sub_1_10k.fastq.gz:single:--outSAMtype SAM"
+    # Annotated-junction tier. The tiers above use real RNA-seq reads, of which
+    # only a small fraction cross an intron -- yeast has a few hundred -- so the
+    # annotated-junction code paths (the sjA shortcut, annotated-boundary
+    # snapping) are barely exercised and a divergence in them is easy to miss.
+    # These reads are drawn from spliced transcripts and every one of them
+    # crosses an annotated junction. Generate with:
+    #   test/simulate_junction_reads.py <genome.fa> <annotation.gtf> \
+    #       test/data/small/yeast/reads/sim_junction_1k.fastq 1000 100
+    "yeast_junction_gtf:yeast_gtf:sim_junction_1k.fastq:single:--outSAMtype SAM"
 )
 
 # ==============================================================================
@@ -350,6 +359,7 @@ run_test_case() {
     # Parse test specification
     IFS=':' read -r name dataset reads mode extra_args <<< "$test_spec"
 
+
     local test_dir="$RESULTS_DIR/${TIMESTAMP}_${name}"
     local star_dir="$test_dir/star"
     local rustar_aligner_dir="$test_dir/rustar-aligner-aligner"
diff --git a/test/simulate_junction_reads.py b/test/simulate_junction_reads.py
new file mode 100644
index 0000000..8e0b21b
--- /dev/null
+++ b/test/simulate_junction_reads.py
@@ -0,0 +1,156 @@
+#!/usr/bin/env python3
+"""Simulate junction-spanning reads from a genome + GTF.
+
+The existing test tiers all use an unannotated index, which makes the
+annotated-junction code paths (`sjA` shortcut, annotated-boundary snapping)
+invisible: they simply never run. This generates reads that are guaranteed to
+cross annotated exon-exon junctions, so those paths are exercised and any
+divergence from STAR shows up as a junction count or CIGAR difference.
+
+Reads are drawn from spliced transcript sequences, centred on junctions, with a
+fixed seed so the FASTQ is reproducible byte-for-byte.
+
+Usage:
+    simulate_junction_reads.py <genome.fa> <annotation.gtf> <out.fastq> [n_reads] [read_len]
+"""
+
+import gzip
+import random
+import sys
+from collections import defaultdict
+
+
+def read_fasta(path):
+    """Return {contig_name: sequence}. Names are truncated at first whitespace."""
+    opener = gzip.open if path.endswith(".gz") else open
+    seqs, name, chunks = {}, None, []
+    with opener(path, "rt") as fh:
+        for line in fh:
+            if line.startswith(">"):
+                if name is not None:
+                    seqs[name] = "".join(chunks)
+                name = line[1:].split()[0]
+                chunks = []
+            else:
+                chunks.append(line.strip())
+    if name is not None:
+        seqs[name] = "".join(chunks)
+    return seqs
+
+
+def read_exons(path):
+    """Return {transcript_id: (contig, strand, [(start, end), ...])}, 1-based inclusive."""
+    opener = gzip.open if path.endswith(".gz") else open
+    exons = defaultdict(list)
+    meta = {}
+    with opener(path, "rt") as fh:
+        for line in fh:
+            if line.startswith("#"):
+                continue
+            f = line.rstrip("\n").split("\t")
+            if len(f) < 9 or f[2] != "exon":
+                continue
+            attrs = f[8]
+            key = 'transcript_id "'
+            i = attrs.find(key)
+            if i < 0:
+                continue
+            j = attrs.find('"', i + len(key))
+            tid = attrs[i + len(key) : j]
+            exons[tid].append((int(f[3]), int(f[4])))
+            meta[tid] = (f[0], f[6])
+    out = {}
+    for tid, blocks in exons.items():
+        blocks.sort()
+        contig, strand = meta[tid]
+        out[tid] = (contig, strand, blocks)
+    return out
+
+
+COMPLEMENT = str.maketrans("ACGTNacgtn", "TGCANtgcan")
+
+
+def revcomp(s):
+    return s.translate(COMPLEMENT)[::-1]
+
+
+def main():
+    if len(sys.argv) < 4:
+        sys.exit(__doc__)
+    genome_path, gtf_path, out_path = sys.argv[1:4]
+    n_reads = int(sys.argv[4]) if len(sys.argv) > 4 else 1000
+    read_len = int(sys.argv[5]) if len(sys.argv) > 5 else 100
+
+    genome = read_fasta(genome_path)
+    transcripts = read_exons(gtf_path)
+
+    # Collect every (transcript, junction) pair with enough flank on both sides
+    # to host a read that actually crosses the junction.
+    candidates = []
+    for tid, (contig, strand, blocks) in transcripts.items():
+        if contig not in genome or len(blocks) < 2:
+            continue
+        # Offset of each junction within the spliced transcript.
+        offset = 0
+        for k in range(len(blocks) - 1):
+            offset += blocks[k][1] - blocks[k][0] + 1
+            candidates.append((tid, offset))
+
+    if not candidates:
+        sys.exit("no multi-exon transcripts found; nothing to simulate")
+
+    rng = random.Random(20260728)
+    rng.shuffle(candidates)
+
+    spliced_cache = {}
+
+    def spliced(tid):
+        if tid not in spliced_cache:
+            contig, strand, blocks = transcripts[tid]
+            seq = "".join(genome[contig][s - 1 : e] for s, e in blocks)
+            spliced_cache[tid] = (seq, strand)
+        return spliced_cache[tid]
+
+    # Yeast has only a few hundred introns, so one read per junction is not
+    # enough to be interesting. Cycle the junction list, varying the overhang
+    # each pass, until the requested count is reached or a pass places nothing.
+    written = 0
+    with open(out_path, "w") as out:
+        while written < n_reads:
+            placed_this_pass = 0
+            for tid, junction_offset in candidates:
+                if written >= n_reads:
+                    break
+                seq, _strand = spliced(tid)
+                if len(seq) < read_len:
+                    continue
+                # Put the junction in the middle third of the read so both
+                # overhangs comfortably clear alignSJDBoverhangMin.
+                overhang = rng.randint(read_len // 3, 2 * read_len // 3)
+                start = junction_offset - overhang
+                if start < 0 or start + read_len > len(seq):
+                    continue
+                read = seq[start : start + read_len].upper()
+                if "N" in read:
+                    continue
+                # Half the reads on the opposite strand, so both orientations
+                # of the annotated-junction path get exercised.
+                if written % 2 == 1:
+                    read = revcomp(read)
+                out.write(f"@sim_{written}_{tid}_j{junction_offset}\n")
+                out.write(read + "\n+\n" + "I" * read_len + "\n")
+                written += 1
+                placed_this_pass += 1
+            if placed_this_pass == 0:
+                break
+
+    print(f"wrote {written} junction-spanning reads to {out_path}", file=sys.stderr)
+    if written < n_reads:
+        print(
+            f"note: only {written} of {n_reads} requested reads were placeable",
+            file=sys.stderr,
+        )
+
+
+if __name__ == "__main__":
+    main()