feat(genome): GTF-derived splice-junction insertion into Genome + SA (Phase 17.X)

[ewels]

Summary

Ports STAR's sjdbPrepare.cpp + sjdbBuildIndex.cpp to ruSTAR. At genomeGenerate, GTF-derived splice junctions are now baked into the Genome binary and indexed by the suffix array, matching STAR's layout byte-for-byte.

Stacked on scverse#7 (TranscriptomeSAM). Merge scverse#7 first, then rebase this branch onto main.

Validation — 15/15 files match STAR at genomeGenerate

test/diff_star_transcriptome.sh against STAR 2.7.11b on the synthetic 2-chr / 4-transcript / 3-junction test case:

✓ chrName.txt
✓ chrLength.txt
✓ chrStart.txt
✓ chrNameLength.txt
✓ transcriptInfo.tab
✓ exonInfo.tab
✓ geneInfo.tab
✓ exonGeTrInfo.tab
✓ sjdbList.fromGTF.out.tab
✓ sjdbInfo.txt
✓ sjdbList.out.tab
✓ Genome
✓ genomeParameters.txt
~ SA size-only match (51929 bytes) — content differs (expected, different SA algorithm)
~ SAindex size-only match (6027 bytes) — content differs (expected, different SA algorithm)
RESULT: 15/15 files identical (or size-matched) to STAR's output

SA/SAindex content divergence is expected: STAR's bucketed parallel SA construction breaks ties differently than ruSTAR's comparator-based sort. Size parity proves we indexed the same genome+Gsj length. Byte-for-byte SA parity would require reimplementing STAR's SA algorithm — out of scope for this PR.

Round-trip verified: ruSTAR can load the extended index it built and align reads end-to-end (load path correctly reads n_genome=524585 from genomeFileSizes, confirming Gsj bytes are included alongside the real genome bytes).

What changed

src/junction/sjdb_insert.rs (new module, 19 unit tests)

• compute_shifts — STAR's sjdbShiftLeft / sjdbShiftRight (repeat length across intron boundary, capped at 255, stops at N).
• PreparedJunction + prepare_junction — combines motif detection, shifts, strand derivation. stored_start/end (for sjdbInfo.txt) vs original_start/end (for Gsj extraction) follow STAR's convention where canonical motifs store the pre-shift ORIGINAL coord and non-canonical ones store the SHIFTED coord.
• sort_and_dedup — STAR's post-sort + cross-strand dedup (sjdbPrepare.cpp:141-192).
• build_gsj — concatenated flanking sequences, 2*overhang+1 bytes per junction (donor + acceptor + spacer=5).
• write_sjdb_info_tab + write_sjdb_list_out_tab — byte-exact file writers (header <n>\t<overhang>\n + six tab-separated fields per row; chr-local 1-based coords + strand char).

src/genome/mod.rs — Genome::append_sjdb()

Extends the forward buffer with Gsj, rebuilds the RC over the extended range. chr_start[n_chr_real] stays pinned at the pre-sjdb forward total (matches STAR's chrStart.txt — verified against the 10k yeast docker output). n_genome grows to include Gsj.

src/index/mod.rs — GenomeIndex::build orchestration

1. Parse GTF once (shared between junction_db, transcriptome, and sjdb_insert — was parsed twice before).
2. Convert chr-local 1-based GTF intron coords to 0-based absolute genome offsets for prepare_junction.
3. Sort + dedup.
4. Build Gsj and call Genome::append_sjdb() BEFORE the SA is built (STAR indexes Gsj alongside the real genome in a single SA, sjdbBuildIndex.cpp:293).
5. Stash sorted prepared junctions in a new prepared_junctions field so write() can emit sjdbInfo.txt + sjdbList.out.tab.

src/index/io.rs — load path

Reads n_genome from genomeParameters.txt's genomeFileSizes line. The chr_start boundary is only the pre-sjdb forward total, so with sjdb-extended indices the two values diverge and chr_start cannot be used as the total genome size.

Other

• src/junction/mod.rs — SpliceJunctionDb::from_raw_junctions helper (shared by from_gtf and the build path, avoiding the GTF-parse-twice problem).
• test/diff_star_transcriptome.sh — extended from 10 files to 15 (added Genome, sjdbInfo.txt, sjdbList.out.tab, plus size-only check for SA/SAindex).

Pre-existing bugs flagged (separate PRs)

Documented during the transcriptome-BAM validation in scverse#7 — not introduced or resolved here:

1. QUAL encoding in transcriptome BAM — ruSTAR double-encodes Phred+33.
2. Default SAM attributes — ruSTAR emits AS nM in addition to STAR's default NH HI.
3. FASTQ parser — fails on some specific reads in the 50k yeast set.

Test plan

• cargo test — 359 unit + 4 integration + 1 doc test, all passing
• cargo clippy --all-targets — 0 new warnings (27 pre-existing)
• cargo fmt --check — clean
• test/diff_star_transcriptome.sh — 15/15 files identical/size-matched against STAR 2.7.11b
• Round-trip: ruSTAR reads its own sjdb-extended index and aligns reads end-to-end

🤖 Generated with Claude Code

[ewels]

Moved to upstream scverse#8.