feat(params): accept --limitGenomeGenerateRAM for STAR CLI compatibility#36
Closed
pinin4fjords wants to merge 4 commits into
Closed
feat(params): accept --limitGenomeGenerateRAM for STAR CLI compatibility#36pinin4fjords wants to merge 4 commits into
--limitGenomeGenerateRAM for STAR CLI compatibility#36pinin4fjords wants to merge 4 commits into
Conversation
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 scverse#25
Contributor
Author
|
I'm unsure to what extent this makes sense really. I get Claude's rationale that it would be good not to break any CLI options though, even if they do nothing. |
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>
Contributor
Author
|
Added a runtime warning when the flag is set to a non-default value. Verified on macOS/aarch64: Default and no-flag invocations stay silent (no spurious warning when the user didn't explicitly ask for a cap). A user passing a real RAM cap now sees that the cap isn't applied rather than having to read the PR body to find out. Pushed as a follow-up commit ( |
Collaborator
|
I merged this, and extended it to take bytes and human readable values, like |
Psy-Fer
added a commit
that referenced
this pull request
May 29, 2026
* 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>
Collaborator
|
merged in #77 |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
3 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Accepts the STAR-compatible
--limitGenomeGenerateRAMflag at the CLI parser. The value is parsed and stored but not currently enforced — rustar's memory management is independent. This unblocks every wrapper that passes the flag derived from job resources (e.g. nf-core/rnaseq'sSTAR_GENOMEGENERATENextflow module).Before
After
Test plan
cargo buildcargo clippy --all-targets -- -D warningscargo fmt --checkcargo test --libA future PR can wire the value into the suffix-array build path if real RAM capping is desired.
Fixes #25