genomeGenerate: libsais builder, --genomeType SuperTranscriptome, --genomeTransformOutput SAM

[BenjaminDEMAILLE]

genomeGenerate: libsais suffix-array builder, selected by --limitGenomeGenerateRAM, plus four genome-index parameters.

What changed

sa_build::build_libsais builds the suffix array with libsais instead of caps-sa, with automatic alphabet-width selection (u8/u16 for few segments, i32 large-alphabet for many-segment genomes).

Selection is driven by --limitGenomeGenerateRAM, not an environment variable. An earlier revision of this PR gated libsais behind RUSTAR_USE_LIBSAIS=1, which reads as debug scaffolding and is not how STAR exposes a build choice. libsais is now used when its estimated peak (17 × N bytes) fits the RAM limit, and caps-sa's external-memory path is used otherwise. That turns --limitGenomeGenerateRAM — previously warned-and-ignored — into a real flag.

--genomeTransformOutput SAM reports alignments in the original genome's coordinates. --genomeTransformType Haploid bakes a VCF's variants into the genome so reads carrying those alleles align without mismatches, but every reported coordinate then refers to a genome nobody else has. Each alignment is mapped back through the conversion blocks written at build time: an indel baked into the sequence reappears as an I/D CIGAR operation at the original position, and junction motifs are reclassified against the original genome rather than the transformed one.

Transcripts are converted immediately after alignment, before anything reads them, so records, junction counts and statistics share one coordinate space. Converting at the writer instead would leave SJ.out.tab and the counters describing a genome the SAM file does not. The SAM header comes from the original genome as well — original coordinates under transformed reference lengths would produce a file no downstream tool could read correctly.

A paired result converts as a unit: if one mate fails to convert the pair is dropped whole rather than reported as a half-mapped read that never existed. Insert size is recomputed from the converted mates.

Combinations whose other outputs would remain in transformed space are refused: --quantMode, --soloType, --outWigType. --genomeTransformOutput SJ and Quant are unimplemented and rejected.

--genomeType SuperTranscriptome condenses the genome to the union of its annotated exons: overlapping exons merge, the merged intervals concatenate, and overlapping transcripts group into superTranscripts, one per condensed chromosome (st0, st1, ...). The index then covers only exonic sequence, and introns cannot be crossed because they are not present. Requires --sjdbGTFfile; Transcriptome is still refused.

The condensed genome comes back with its annotation already addressed to it, so the transcriptome tables, the junction extractor and the suffix array consume it unchanged. Four files are written alongside the index, as STAR writes them: superTranscriptSequences.fasta, transcriptSequences.fasta, superTranscriptSJcollapsed.tsv, and fullGenome/conversionToFullGenome.tsv.

D20, deliberately not reproduced. STAR condenses the sequence before it fills the reverse-complement half of its genome buffer, so every minus-strand exon reads the allocation's spacer bytes and minus-strand superTranscripts come out as uninitialised memory. Here that half is already filled, so a minus-strand exon yields its actual reverse complement. Reproducing STAR would mean writing code whose correct behaviour is to emit garbage and a test asserting it. Recorded in docs-old/dev/divergences.md; the test asserts the hand-derived reverse complement and that no 0x00 byte reaches the output.

One behaviour worth stating because it looks like a gap: a transcript's own intron leaves no junction behind. With nothing else annotated between two exons they end up contiguous in the condensed genome, so there is no gap left to record. A junction survives only when another transcript's exon separates them. Both cases have a test.

Also adds --genomeSuffixLengthMax and --sjdbInsertSave to the parameter surface, accepted and documented inert; off-menu values are rejected rather than silently ignored.

Why

Index construction was the slowest step with no output risk attached: a suffix array is unique for a given text, so a faster builder cannot change results. build_libsais sorts the same per-segment sentinel-transformed text the existing caps-sa sentinel arm uses, so its packed output is byte-for-byte identical to build(), and hence to STAR 2.7.11b.

The RAM gate exists because the two builders trade differently. libsais builds in memory; caps-sa's default path spills to disk. On a 48 MB / 24-chromosome genome:

Builder	Construction time	Peak RSS	SA output
caps-sa	12.8 s	~585 MB	396,000,003 bytes
libsais	4.2 s	~1.49 GB	byte-identical

~3.0× faster at ~2.5× the peak RAM. Extrapolated to human (2N ≈ 6.2 B), libsais would need ~50 GB while caps-sa external-memory stays at a few GB — hence the gate rather than a flat default.

Verification

Byte-identity checked at four levels, all green:

• unit — build_libsais == caps-sa build() on single-chromosome, multi-chromosome-with-N, and repetitive fixtures
• large-alphabet — the i32 arm forced on small fixtures
• streaming — build_libsais driven exactly as GenomeIndex::build does, matching the caps-sa streaming byte stream
• real genome — the 48 MB genomeGenerate run above: identical SA and SAindex

The back-transform is verified end to end: a genome with a 5-base deletion at position 5001, one read at 6000. Without the flag the reported POS is 5996; with it, 6001 — the original coordinate — and the CIGAR is unchanged, since the read spans no variant. Unit tests cover an exon straddling a block boundary coming back out as 10M5D10M, a seamless split collapsing again, a motif read from the original genome, a junction canonical only in the transformed genome turning non-canonical, and an alignment outside the block map being dropped rather than guessed at.

libsais is pinned with default-features = false (no OpenMP); libsais-sys builds its C core via cc. CI green on all 5 platforms including windows-x86_64. On-disk index format unchanged.

SuperTranscriptome is verified end to end: the 20 kb test genome with its two annotated 50 bp exons condenses to a single 100 bp st0; a read straddling the 200 bp intron aligns 50M at position 26, and the conversion map names the two source intervals at 10000 and 10250. Unit tests cover exon merging, the condensed annotation, disjoint transcripts becoming separate superTranscripts, junction collapsing, and D20.

585 lib + 23 integration tests green, clippy -D warnings clean, cargo fmt --check clean, MSRV 1.89.

Note: docs-old/dev/divergences.md is created here with the D20 entry and also on #152 with the solo entries. Whichever merges second merges the two entry lists.

Bookkeeping note: the four flags above are listed in NOT_YET_ACCEPTED on #147. Whichever merges second moves them out of that list in the same commit.

🤖 Generated with Claude Code