fix(sjdb): key `SpliceJunctionDb` in genome-absolute 0-based coordinates by pinin4fjords · Pull Request #45 · scverse/rustar-aligner

pinin4fjords · 2026-05-12T17:38:48Z

Summary

SpliceJunctionDb was keyed in chromosome-local 1-based coordinates (extracted straight from GTF column 5) but consulted in two different coordinate spaces at runtime:

Call site	Convention	On chr 0 (chr_start=0)	On chr N+
stitch-time (`src/align/stitch.rs:1305-1314`)	genome-absolute 0-based	off by 1 -> miss	off by `chr_start[N]` -> miss
stats-time (`src/lib.rs:1860-1894`)	genome-absolute 1-based	matches -> hit	off by `chr_start[N]` -> miss

Neither matched the DB. Smoking gun: on the same WT_REP2 BAM, SJ.out.tab had 2 rows annotated=1 (stats-time accidentally hits chr-0) while Log.final.out reported Number of splices: Annotated (sjdb) | 0 (stitch-time always misses). The stitch-time miss is the load-bearing one - it drops sjdb_score and pulls in the stricter align_sj_overhang_min gate, costing ~50 % of GT/AG splices on the nf-core/rnaseq test profile.

Fix

Normalise the DB to genome-absolute 0-based at construction (matching the existing convention used by prepared_junctions in src/index/mod.rs and by SpliceJunctionStats). Single source of truth contained in the GTF extraction step.

extract_junctions_configured in src/junction/gtf.rs now adds chr_start[chr_idx] and subtracts 1 before returning each (intron_start, intron_end). SpliceJunctionDb::from_raw_junctions consumes those tuples verbatim.
Stats-time call sites at src/lib.rs:1860 and src/lib.rs:796 now record (genome_pos, genome_pos + intron_len - 1) instead of (genome_pos + 1, genome_pos + intron_len). genome_pos is the 0-based first intronic base, which is what detect_splice_motif already expects.
Stitch-time call site at src/align/stitch.rs is unchanged - it already passes genome-absolute 0-based.
SJ.out.tab writer in src/junction/sj_output.rs now adds + 1 after subtracting chr_start[chr_idx] so the chr-local 1-based output is unchanged.
prepared_junctions construction in src/index/mod.rs is simplified: the chr-local-to-absolute conversion that used to live there now happens upstream in the GTF extractor.

Test plan

New unit test test_db_keyed_in_genome_absolute_zero_based_multi_chr builds a SpliceJunctionDb on a 2-chromosome toy genome (chr_start[1] = 1000) and asserts is_annotated succeeds for the chr-1 junction at its genome-absolute 0-based key. The pre-fix chr-local key and the pre-fix stitch-time off-by-one key must both miss
Existing tests that encoded the chr-local 1-based assumption have been updated to genome-absolute 0-based equivalents (those were testing the bug)
cargo build
cargo clippy --lib -- -D warnings
cargo fmt --check
cargo test (384 lib tests + integration suites pass)

After this fix, Number of splices: Annotated (sjdb) in Log.final.out matches the count of annotated=1 rows in SJ.out.tab on the same BAM, and the per-sample Annotated count is in the same order as STAR's on equivalent inputs.

Fixes #27

The DB was keyed in chromosome-local 1-based coords (straight from the GTF) but consulted in genome-absolute 0-based at stitch time, so every sjdb lookup during alignment missed. Annotated splice events never got the sjdb_score bonus and the stricter overhang gate fired in their place, dropping ~50 % of GT/AG splices on the test profile. The stats-time site queried in genome-absolute 1-based-equivalent, so it accidentally matched on chr 0 (chr_start=0) but missed on every other chromosome -- producing inconsistent answers between SJ.out.tab and Log.final.out on the same BAM. Normalise the DB to genome-absolute 0-based at construction, matching the convention prepared_junctions and SpliceJunctionStats already use. Update the stats-time call site to query in the new space. Stitch-time needs no change. Fixes scverse#27 Co-Authored-By: Claude <noreply@anthropic.com>

The previous commit normalised the DB to genome-absolute 0-based and updated the stats-time query. Stitch-time was left unchanged on the assumption that it already passed genome-absolute 0-based — half right. donor_fwd was correct (first intron base, 0-based) but acceptor_fwd = donor_fwd + del landed on the first base AFTER the intron, while the DB keys store the last intron base. Lookup missed every annotated junction. Subtract 1 from the acceptor to land on the last intron base. After this, Log.final.out Annotated (sjdb) goes from 0 to a non-zero count that's consistent with the annotated=1 row count in SJ.out.tab on the same BAM. Co-Authored-By: Claude <noreply@anthropic.com>

pinin4fjords · 2026-05-12T18:01:26Z

Verified end-to-end on macOS/aarch64 against the rebuilt fix branch (PE yeast + --twopassMode Basic --sjdbGTFfile genes.gtf):

=== rustar Log.final.out splice section ===
Number of splices: Total            | 366
Number of splices: Annotated (sjdb) | 95         <- was 0 pre-fix
Number of splices: GT/AG            | 266
Number of splices: GC/AG            | 2
Number of splices: AT/AC            | 5
Number of splices: Non-canonical    | 93

=== rustar SJ.out.tab annotated col distribution ===
    14  annotated=0
     2  annotated=1

The two counters are now consistent — the pre-fix Log.final.out = 0 vs SJ.out.tab = 2 annotated smoking gun is gone. Both fixes in this branch are load-bearing: the coord-space normalisation, plus the acceptor_fwd - 1 follow-up on the stitch-time lookup (without the latter, donor_fwd + del lands one base past the last intron base and every annotated lookup still misses).

Remaining gap (separate from this PR, tracked as #47): total splice count is still 366 vs STAR's ~720 — pass-1 isn't using sjdb-derived junctions as alignment candidates, only at scoring time. That's a deeper change.

Pass-1 alignment was silently dropping SA hits that fall in the appended splice-junction flanking buffer (Gsj). Those hits encode candidate splice events for annotated junctions: when a read seed straddles the donor/acceptor boundary of a Gsj slot, the matching genome bytes live in two non-adjacent real-genome positions. STAR decodes these via g1 >= P.nGenomeReal and feeds the (donor, acceptor) pair into the seed pipeline; rustar-aligner was discarding them at position_to_chr (which returns None for positions past the last real chromosome). Plumb n_genome_real onto Genome (pinned at the pre-sjdb forward total) and reload prepared junctions + sjdbOverhang from sjdbInfo.txt on the index-load path. cluster_seeds now expands each SA hit through decode_gsj_hit: real-genome hits pass through unchanged, single-flank Gsj hits are skipped (the equivalent real-genome SA entry already exists in the same range), and boundary-crossing Gsj hits split into two virtual seeds with the donor-side read run and the acceptor-side read run pointing at their respective real-genome positions. The existing stitch DP then chains them via its splice branch. Verified on the yeast SRR6357072 reproducer from the issue: Number of splices: Total jumps from 366 to 631 (target ~720) and GT/AG from 266 to 531 with non-canonical unchanged at 93. The remaining gap is gated on PR #45's annotated-junction lookup landing. Fixes #47 Co-Authored-By: Claude <noreply@anthropic.com>

@rg

* feat(params): accept --limitGenomeGenerateRAM for STAR CLI compatibility STAR exposes --limitGenomeGenerateRAM as a memory cap for genomeGenerate. Pipelines that wrap STAR commonly derive a value from their task resources and pass it through. Currently rustar rejects the flag at the CLI parser, breaking those pipelines. Accept the flag with STAR's default of 31 GiB. The value is parsed but not enforced — rustar's memory management is independent. Fixes #25 * fix(chimeric): create output parent directory before writing The chimeric output writer constructs its path as <outFileNamePrefix>/Chimeric.out.junction, treating the prefix as a directory regardless of whether it ends in `/`. In two-pass mode the chim writer fires before any other output creates that directory, so the file open fails with "No such file or directory" and the entire run aborts. Call create_dir_all on the parent of the chim output path before opening the file. Without two-pass mode the bug was masked because another output writer happened to create the dir first. Fixes #35 * fix(bam): subtract Phred+33 offset from QUAL before writing to BAM The BAM binary QUAL field stores raw Phred values (0-93) per SAM spec §4.2.3. rustar was passing FASTQ ASCII bytes (Phred+33 encoded, e.g. 'B' = ASCII 66 = Phred 33) straight to noodles' QualityScores::from, which expects raw Phred values. Every QUAL byte in a rustar BAM was therefore +33 above what the spec mandates; samtools view rendered the SAM text with another +33 on top, producing biologically impossible quality strings like 'cccgg...'. Subtract 33 from each FASTQ ASCII byte at every BAM-write call site before constructing the QualityScores. Use saturating_sub so a malformed FASTQ byte < 33 clamps to 0 rather than underflowing. samtools stats average_quality on the nf-core/rnaseq test profile now reads 35.3 (matching STAR) rather than the previous 68.3. Fixes #34 * fix(transcriptome-bam): write per-record RG:Z: tag matching @rg header When --outSAMattrRGline is supplied, rustar writes the @rg header to both the genome and transcriptome BAMs, but only stamps the per-record RG:Z: tag on the genome BAM. The transcriptome BAM was emitting 0 RG:Z: records vs STAR's 1-per-record output. Add the RG tag stamp to the transcriptome record builder (paired-end and single-end paths) so every record carries RG:Z:<id> matching the @rg header, byte-symmetric with STAR. Fixes #32 * feat(bam): populate @pg header with PN/VN/CL fields Per SAM spec §1.3, the @pg line conventionally carries PN (program name), VN (version), and CL (command line) alongside ID. rustar was emitting only ID:rustar-aligner, leaving downstream provenance tools (MultiQC's program-version table, dx-toolkit lineage tracking) with a blank entry. Expand the header writer to emit: @pg\tID:rustar-aligner\tPN:rustar-aligner\tVN:<cargo pkg version>\tCL:<args> The full command line is captured in main() before clap parses it, then threaded into Parameters via a new (skip) field so it reaches the SAM header builder. Version comes from CARGO_PKG_VERSION at compile time. This matches STAR's @pg format and gives downstream tools the provenance they need. Fixes #33 (the @pg header gap; AS divergence is a separate item). * chore(params): trim verbose doc comment * chore(chimeric): remove narrative test comment * chore(bam): remove narrative inline comment * chore(transcriptome-bam): remove narrative comments * chore(bam): remove narrative comments * fix(transcriptome-bam): populate mate fields on paired-end records build_transcriptome_records_pe was stamping SEGMENTED + FIRST/LAST_SEGMENT but leaving PROPERLY_SEGMENTED, MATE_REVERSE_COMPLEMENTED, RNEXT, PNEXT, and TLEN unset on every transcriptome record. Salmon's alignment-mode fragment-length inference reads those mate fields; with them missing it falls back to its 250+/-25 prior, distorting paired-end TPMs. Walk the projected (rec1, rec2) pairs after the existing flag-stamping loops and mirror the mate-bookkeeping the genome-space build_paired_mate_record already does: PROPERLY_SEGMENTED + MATE_REVERSE_COMPLEMENTED + mate ref id + mate position + signed TLEN (leftmost +, rightmost -). Fixes #22 * fix(bam): emit SAM-spec NM:i: tag distinct from STAR-internal nM:i: --outSAMattributes NM was being routed to the existing nM writer (substitutions only), so requests for the SAM-spec edit-distance tag silently produced wrong values. Downstream tools that read NM:i: (samtools stats, picard, MultiQC) saw nothing. Treat NM and nM as distinct attribute tokens. When NM is requested, compute SAM-spec edit distance per the SAM v1 spec section 1.4: substitutions + inserted bases + deleted bases (excluding intron N skips). Keep nM emission unchanged when explicitly requested. Fixes #29 * feat(bam): emit XS:A: strand tag when --outSAMstrandField intronMotif The flag was accepted at the CLI parser but never produced any XS:A: tags on the output BAM. Downstream tools that need strand inference (StringTie, Cufflinks, rseqc infer_experiment.py) saw nothing. Mirror STAR's two coupling paths: --outSAMattributes XS auto-enables --outSAMstrandField intronMotif, and vice versa. For each output record with at least one canonical intron motif, map the motif to + or - per STAR's convention; non-canonical or mixed motifs omit the tag. Fixes #30 * feat(output): emit Log.out, Log.progress.out, _STARpass1/SJ.out.tab STAR writes three log files alongside Log.final.out and keeps two-pass intermediates in <prefix>_STARpass1/. rustar was only writing Log.final.out and emitting the pass-1 SJ tab as <prefix>SJ.pass1.out.tab at the top level. Add minimal Log.out (parameters dump + per-phase timestamps) and Log.progress.out (timestamp + mapping speed) writers next to the existing Log.final.out writer. Move the pass-1 SJ tab into <prefix>_STARpass1/SJ.out.tab and mkdir the parent first. The Log.out content is intentionally a stub matching the file's existence rather than STAR's full verbosity; that's a follow-up. Fixes #28 Co-Authored-By: Claude <noreply@anthropic.com> * fix(sjdb): key SpliceJunctionDb in genome-absolute 0-based coordinates The DB was keyed in chromosome-local 1-based coords (straight from the GTF) but consulted in genome-absolute 0-based at stitch time, so every sjdb lookup during alignment missed. Annotated splice events never got the sjdb_score bonus and the stricter overhang gate fired in their place, dropping ~50 % of GT/AG splices on the test profile. The stats-time site queried in genome-absolute 1-based-equivalent, so it accidentally matched on chr 0 (chr_start=0) but missed on every other chromosome -- producing inconsistent answers between SJ.out.tab and Log.final.out on the same BAM. Normalise the DB to genome-absolute 0-based at construction, matching the convention prepared_junctions and SpliceJunctionStats already use. Update the stats-time call site to query in the new space. Stitch-time needs no change. Fixes #27 Co-Authored-By: Claude <noreply@anthropic.com> * revert(stats): drop Log.out / Log.progress.out stub writers The previous commit added writers that mimicked STAR's section-header structure with placeholder content (a Debug-format params dump + three timestamps for Log.out, a header + final-timestamp line for Log.progress.out). That closes the file-existence gap but makes the output look like STAR-equivalent verbose logging when it isn't — downstream tools parsing for memory usage, per-chunk progress, or warnings would silently get nothing. Drop the stub writers. The SJ.pass1.out.tab → <prefix>_STARpass1/SJ.out.tab placement fix stays — that's a real path-parity change with no fakery. Real STAR-equivalent Log.out / Log.progress.out content (per-phase progress, warnings, memory, periodic updates during long runs) is a separate engineering effort, not a stub. Co-Authored-By: Claude <noreply@anthropic.com> * chore(params): warn at runtime when --limitGenomeGenerateRAM is set Accepting the flag silently could mislead a user who passes a non-default cap into thinking rustar enforces it. Emit a warning whenever the value differs from STAR's default so the user knows the cap isn't applied. Co-Authored-By: Claude <noreply@anthropic.com> * fix(sjdb): off-by-one on stitch-time acceptor coordinate The previous commit normalised the DB to genome-absolute 0-based and updated the stats-time query. Stitch-time was left unchanged on the assumption that it already passed genome-absolute 0-based — half right. donor_fwd was correct (first intron base, 0-based) but acceptor_fwd = donor_fwd + del landed on the first base AFTER the intron, while the DB keys store the last intron base. Lookup missed every annotated junction. Subtract 1 from the acceptor to land on the last intron base. After this, Log.final.out Annotated (sjdb) goes from 0 to a non-zero count that's consistent with the annotated=1 row count in SJ.out.tab on the same BAM. Co-Authored-By: Claude <noreply@anthropic.com> * fix: decode Gsj-region SA hits into split splice seeds Pass-1 alignment was silently dropping SA hits that fall in the appended splice-junction flanking buffer (Gsj). Those hits encode candidate splice events for annotated junctions: when a read seed straddles the donor/acceptor boundary of a Gsj slot, the matching genome bytes live in two non-adjacent real-genome positions. STAR decodes these via g1 >= P.nGenomeReal and feeds the (donor, acceptor) pair into the seed pipeline; rustar-aligner was discarding them at position_to_chr (which returns None for positions past the last real chromosome). Plumb n_genome_real onto Genome (pinned at the pre-sjdb forward total) and reload prepared junctions + sjdbOverhang from sjdbInfo.txt on the index-load path. cluster_seeds now expands each SA hit through decode_gsj_hit: real-genome hits pass through unchanged, single-flank Gsj hits are skipped (the equivalent real-genome SA entry already exists in the same range), and boundary-crossing Gsj hits split into two virtual seeds with the donor-side read run and the acceptor-side read run pointing at their respective real-genome positions. The existing stitch DP then chains them via its splice branch. Verified on the yeast SRR6357072 reproducer from the issue: Number of splices: Total jumps from 366 to 631 (target ~720) and GT/AG from 266 to 531 with non-canonical unchanged at 93. The remaining gap is gated on PR #45's annotated-junction lookup landing. Fixes #47 Co-Authored-By: Claude <noreply@anthropic.com> * fix(stitch): use sjdb_score in place of motif_score on annotated junctions STAR's Transcript_alignScore.cpp adds either sjdbScore (annotated) or motifScore (unannotated) to the alignment score, not both. rustar was adding both on annotated junctions, over-crediting the score by motif_score per junction. Currently invisible because #27 keeps is_annotated false in practice; becomes the dominant remaining AS divergence on annotated splices as soon as #27 lands. Fixes #50 Co-Authored-By: Claude <noreply@anthropic.com> * fix(score): score_annotated_junction helper uses replacement not addition The helper at src/align/score.rs:137-143 had the same additive bug as the stitch-time call site this PR fixes — adding sjdb_score to the motif score rather than replacing it. Currently only called from its own tests, but a landmine if any other caller picks it up. Switch to replacement semantics matching the stitch path. Update the tests that locked in the old (wrong) additive result. Co-Authored-By: Claude <noreply@anthropic.com> * fix(stats): route reads exceeding --outFilterMultimapNmax to too_many_loci The existing AlignmentStats::record_alignment routing arm for n_alignments > max_multimaps was unreachable on the PE path - the too-many-loci filter returned n_for_mapq = 0, so the caller hit record_alignment(0, ...) and incremented the unmapped bucket. Log.final.out always reported 0 in the too-many-loci field and the per-bucket sum fell 3-of-50000 short of input on the yeast PE test profile (STAR records exactly 3 there on the same data). Return joint_pairs.len() as n_for_mapq from the PE filter site and route through record_alignment with that count, mirroring the SE caller pattern so the _ => arm fires. The read still isn't emitted in the BAM with outSAMmultNmax=-1, matching STAR; only the accounting changes. Fixes #53 Co-Authored-By: Claude <noreply@anthropic.com> * refactor(stats): drop redundant n_for_mapq conditional in PE caller n_for_mapq is usize, so `if n > 0 { n } else { 0 }` is identical to passing n_for_mapq directly to record_alignment. Match the SE caller shape: the variable already carries the right count from the PE TooManyLoci return site, and record_alignment routes via max_multimaps. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * fix test * upstream changes * chore: switch sam attributes to bitflags * Extend PR #36 for human readable ram entries, like 64G * add missing n_genome_real for #51 * Fix CigarOp to cigar::op::new * fmt fix * Update run_tests.sh --------- Co-authored-by: Jonathan Manning <jonathan.manning@seqera.io> Co-authored-by: Claude <noreply@anthropic.com> Co-authored-by: Philipp A. <flying-sheep@web.de>

Psy-Fer · 2026-05-29T12:15:58Z

merged in #77

pinin4fjords and others added 2 commits May 12, 2026 18:38

Merge branch 'main' into pr/pinin4fjords/45

0ab761f

Psy-Fer mentioned this pull request May 29, 2026

Integration/pr batch 1 #77

Merged

7 tasks

Psy-Fer closed this May 29, 2026

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

fix(sjdb): key `SpliceJunctionDb` in genome-absolute 0-based coordinates#45

fix(sjdb): key `SpliceJunctionDb` in genome-absolute 0-based coordinates#45
pinin4fjords wants to merge 3 commits into
scverse:mainfrom
pinin4fjords:fix/sjdb-coord-space

pinin4fjords commented May 12, 2026

Uh oh!

pinin4fjords commented May 12, 2026 •

edited

Loading

Uh oh!

Psy-Fer commented May 29, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

Conversation

pinin4fjords commented May 12, 2026

Summary

Fix

Test plan

Uh oh!

pinin4fjords commented May 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Psy-Fer commented May 29, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

pinin4fjords commented May 12, 2026 •

edited

Loading