Skip to content

perf(build): enable geometric LCP memoization for genome indexing - #228

Open
rob-p wants to merge 9 commits into
mainfrom
codex/geometric-memoization
Open

perf(build): enable geometric LCP memoization for genome indexing#228
rob-p wants to merge 9 commits into
mainfrom
codex/geometric-memoization

Conversation

@rob-p

@rob-p rob-p commented Aug 13, 2026 •
edited
Loading

Copy link
Copy Markdown
Contributor

Summary

This PR enables caps-sa's bounded geometric LCP memoization for ruSTAR's
production segmented external-memory genome-index construction.

  • depend on the published caps-sa = "0.7" crates.io release
  • select the measured defaults through the final stable API,
    LcpMemoizationPolicy::geometric()
  • leave caps-sa's direct kernel as the default for generic callers
  • leave ruSTAR's small in-memory and sentinel-transform paths unchanged

Why this is ruSTAR-specific

The production genome-plus-junction layout contains many long shared contexts
across the phase-4 partition cascades. Each cascade owns an independent bounded
memoization table, activates it lazily after learning enough long exact LCP
intervals, and compares a short prefix normally before consulting the table.
There is no shared-table write contention.

Full-genome benchmark

Complete ruSTAR-shaped GENCODE Human v50 fixture:

  • GRCh38 primary assembly plus comprehensive annotation
  • sjdbOverhang=100
  • actual ruSTAR GTF parsing, junction preparation/deduplication, and
    splice-junction flank generation
  • complete forward/reverse-complement layout
  • STAR segment boundaries and boundary ordering
  • ACGT-start filtering
  • 32 pinned physical cores and u64 indices
  • 6,557,611,930 input symbols, 6,176,694,310 retained suffixes, and about
    1.40 million segments

Memoization policy A/B

The isolated policy comparison that motivated this ruSTAR change was:

Path Build time User CPU Peak RSS
Direct 291.710 s 8,029.81 s 10,672,468 KiB
Geometric memoization 267.128 s 7,607.24 s 9,696,720 KiB

That is an 8.43% wall-time improvement and a 5.26% reduction in user CPU.
Both paths emitted the same count and streaming hash
(e81c8f9881e322148741a23c92ae2000). An earlier independent pair measured
295.048 s versus 270.063 s (8.47%).

Final caps-sa 0.7 release

The pre-release caps-sa audit then optimized the surrounding ruSTAR-shaped
construction as well: direct position/LCP array flow, task-local ping-pong
sorting, and a compact SegmentedText boundary directory for the 1.40 million
segments.

caps-sa 0.7 candidate Build time User CPU Peak RSS Phase 1 Phase 4
Original 0.7 baseline 267.592 s 7,622.88 s 10,512,408 KiB 106.847 s 155.794 s
Published 0.7.0 172.953 s 4,731.05 s 9,169,892 KiB 49.029 s 118.743 s

The final published implementation is 35.4% faster than that original 0.7
baseline and uses 12.8% less peak RSS. It emitted the same full output hash.
On the focused chr21 backbone plus every annotation-derived flank, all
359,616,038 output positions also matched the stored baseline byte-for-byte.

Dependency and API

caps-sa 0.7.0 is published on crates.io. Cargo.lock resolves the registry
package with checksum
e46cbd8870dc17f488813e4c5499ff1bdc2feedb9a650ca2b597933943101394.
The earlier temporary Git SHA pin has been removed.

ruSTAR uses the final concise API:

caps_sa_ext_mem_opts(temp_dir)
    .lcp_memoization(caps_sa::LcpMemoizationPolicy::geometric())

The policy remains explicit in ruSTAR; caps-sa continues to default to the
direct kernel for general inputs.

Validation

Validated locally against the published crates.io package recorded in
Cargo.lock:

  • cargo check --all-targets --locked
  • cargo test --release --locked: 590 library tests plus all integration and
    doc-test targets passed
  • all five segmented-vs-sentinel STAR-ordering oracle tests passed
  • cargo clippy --all-targets --locked -- -D warnings
  • cargo fmt --check

The release gate is now satisfied: caps-sa 0.7.0 is published, caps-sa main
and v0.7.0 resolve to the audited release commit, and this PR consumes the
registry release rather than a Git dependency.

Full-genome benchmark vs upstream STAR (GRCh38 + sjdb)

Measured on the final branch state (commits ea1a359..814d7ce included):
GRCh38 primary assembly plus GENCODE v49 basic annotation,
--sjdbOverhang 100, --runThreadN 32 pinned to 32 physical cores
(taskset -c 0-31, 2× EPYC 9575F), x86-64-v4 build, inputs prewarmed into
the page cache, runs sequential on an otherwise idle machine,
/usr/bin/time -v:

genomeGenerate Wall time User CPU Peak RSS
ruSTAR (this PR) 4:52.3 4,735 s 17.7 GB
STAR 2.7.11b 11:55.5 8,266 s 36.5 GB

ruSTAR builds the complete splice-aware human index 2.45× faster than
upstream STAR in less than half the memory.

STAR index compatibility: byte-identity evidence

We compared the two GRCh38+sjdb genome directories file by file and fixed
every divergence this PR could close:

  • SAindex N-mark flag (ea1a359): ruSTAR never set STAR's SAiMarkNbit,
    leaving 58,474 byte diffs. Both SAindex builders now replicate
    genomeSAindex.cpp's mark semantics; the parallel builder seeds per-chunk
    ind0 state with a backward scan so cross-chunk marks land on STAR's
    slots, and the in-memory builder adopts STAR's absent-slot encoding.
  • Junction set (d12bcff): ruSTAR dropped 1-base introns that STAR keeps
    and lacked STAR's first dedup pass over left-shifted coordinates, so
    repeat-shifted duplicates survived. sjdbList.fromGTF.out.tab also used a
    start-only sort that broke its adjacency merge, leaving 143,860 excess
    lines. All three paths now match GTF_transcriptGeneSJ.cpp and
    sjdbPrepare.cpp branch for branch.
  • Metadata (814d7ce): genomeParameters.txt now records the effective
    sjdbOverhang (0 for a no-sjdb index, where the old value of 100 could
    mislead STAR's loader) and echoes the real command line, and
    genomeGenerate writes a Log.out into the genome directory as STAR
    does.

Byte-level comparison of the two GRCh38+sjdb genome directories after these
fixes:

File Result
Genome, SA (25.1 GB), chr*, sjdbInfo.txt, sjdbList.out.tab, sjdbList.fromGTF.out.tab, transcriptInfo.tab, exonInfo.tab, exonGeTrInfo.tab, geneInfo.tab byte-identical
SAindex (1.5 GB) 4 bytes differ (2 N-mark bits, see below)
genomeParameters.txt identical except the ### argv comment line
Log.out present in both; free-form run log

Indexes built without a GTF are byte-identical across every file, SAindex
included; we verified this on GRCh38 and on yeast.

The 4 remaining SAindex bytes are an artifact of STAR's two-phase build:
it patches the base-genome SAindex during junction insertion
(sjdbBuildIndex.cpp:228-231, 262-284) and never updates N-marks for
k-mers that only became present via junction flanks. ruSTAR's single pass
over the final text marks the current slot instead. Both placements are
valid conservative markers, and DIVERGENCE.md §3.1a documents the
difference.

Functional cross-check

Upstream STAR 2.7.11b aligned 100k SRR21186103 paired-end reads against the
ruSTAR-built index and against its own. Aligned.out.sam (209,729 records)
and SJ.out.tab are byte-identical, excluding only the path-bearing
@PG/@CO header lines; both runs report 89.50% uniquely mapped reads and
106,768 splices. An index generated with ruSTAR is a drop-in replacement
for a STAR-generated one.

🤖 Generated with Claude Code

https://claude.ai/code/session_01HS4Mudeqj3Y1jk5p3fhUWt

@rob-p rob-p changed the title Enable geometric LCP memoization for genome indexing perf(build): enable geometric LCP memoization for genome indexing Aug 13, 2026
@rob-p
rob-p marked this pull request as ready for review August 13, 2026 14:42
@rob-p rob-p mentioned this pull request Aug 13, 2026
Open
4 tasks
BenjaminDEMAILLE and others added 5 commits August 27, 2026 14:01
@BenjaminDEMAILLE
Merge main into codex/geometric-memoization
fcc1f2a
@rob-p @claude
fix(index): set STAR's SAiMarkNbit N-flag when building SAindex
ea1a359 
STAR ORs SAiMarkNmaskC onto the last-written SAindex entry at every
level >= iL4 whenever a suffix's k-mer prefix hits an N
(genomeSAindex.cpp:140-145). rustar never set the bit, leaving 58,474
byte diffs against STAR's SAindex on GRCh38.

Implement the mark semantics in both builders:

- build_parallel: collect mark targets in a shared bitset and OR them
  into the packed values after the gap-fill. Chunk workers seed their
  per-level ind0 state with a short backward scan over the preceding
  SA entries so cross-chunk marks land on the same slots as STAR's
  serial scan.
- build (in-memory serial): follow the reference algorithm directly,
  and adopt STAR's absent-slot encoding (next-present | absent mask,
  tail-filled with nSA | absent mask) that build_parallel already used.

Genome-only GRCh38 and yeast indexes are now byte-identical to
STAR 2.7.11b's, SAindex included. On sjdb builds a handful of mark
bits still differ by one slot because STAR places them against the
pre-insertion base genome (see DIVERGENCE.md 3.1a).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HS4Mudeqj3Y1jk5p3fhUWt
@rob-p @claude
fix(sjdb): match STAR's junction set and fromGTF list exactly
d12bcff 
Three divergences from STAR's GTF-junction pipeline left the inserted
junction set 4 entries larger than STAR's on GRCh38 + GENCODE v49
(469,461 vs 469,457) and sjdbList.fromGTF.out.tab 143,860 lines longer:

- junction/gtf.rs: keep 1-base introns. STAR only skips touching
  (exS <= exE+1, silent) and overlapping (exS <= exE, warning) exon
  pairs (GTF_transcriptGeneSJ.cpp:123-134); rustar also dropped
  length-1 introns, losing 2 junctions STAR keeps.
- sjdb_insert.rs: implement STAR's first dedup pass over left-shifted
  coordinates, partitioned by the raw source strand
  (sjdbPrepare.cpp:75-123). Repeat-shifted duplicates of one splice
  event now collapse (canonical first, then smallest left shift);
  rustar previously kept 6 such duplicates. PreparedJunction gains a
  src_strand field because both of STAR's passes branch on the
  annotation strand, not the motif-derived one, and the cross-strand
  pass now mirrors sjdbPrepare.cpp:154-170 branch-for-branch.
- transcriptome.rs: sort sjdbList.fromGTF.out.tab by (start, end) —
  funCompareUint2 compares two uint64s, not one — collapse rows the
  way GTF_transcriptGeneSJ.cpp:145-158 does, and emit gene lists as
  ascending sets.

Genome, SA, sjdbInfo.txt, sjdbList.out.tab, and
sjdbList.fromGTF.out.tab are now byte-identical to STAR 2.7.11b's on
the full GRCh38 + GENCODE v49 build.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HS4Mudeqj3Y1jk5p3fhUWt
@rob-p @claude
fix(genomeGenerate): STAR-faithful genomeParameters.txt and Log.out
814d7ce 
- genomeParameters.txt: write the EFFECTIVE sjdbOverhang (STAR's
  mapGen.sjdbOverhang — 0 when the index has no sjdb) instead of the
  parameter value. A no-GTF index previously claimed overhang 100,
  which STAR's loader would trust at align time. The ### header now
  echoes the actual command line like STAR does, falling back to the
  old skeleton for API callers constructed without one.
- Log.out: genomeGenerate now writes a STAR-shaped run log to the
  output prefix and copies it into the genome directory, matching
  STAR's manifest (genomeGenerate.cpp copies its log there).
- index/mod.rs: thread the effective overhang into write_index_files;
  drop a stale doc reference to the removed on-the-fly SAindexBuilder
  flow.
- DIVERGENCE.md: update 3.1 for the argv echo, add 3.1a documenting
  the remaining sjdb-build SAindex N-mark bits (STAR places them
  against the pre-insertion base genome; 2 slots out of 357,913,940
  on GRCh38 + GENCODE v49) and 3.1b for Log.out content.

With this and the two preceding commits, a full GRCh38 + GENCODE v49
sjdb index differs from STAR 2.7.11b's only in the ### argv line,
Log.out content, and those 4 SAindex bytes; STAR aligns 100k read
pairs byte-identically against either index.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HS4Mudeqj3Y1jk5p3fhUWt
@rob-p @claude
fix(windows): avoid Log.out sharing violation in genomeGenerate
272615c 
The genomeDir Log.out was produced with fs::copy from <prefix>Log.out.
Integration tests run genomeGenerate as concurrent processes in one
working directory, so several of them raced on ./Log.out; on Windows,
CopyFileEx opens its source without write sharing, and a concurrent
File::create from another process fails with a sharing violation
(os error 32). Three sjdb tests failed this way on windows-x86_64 CI.

Write the genomeDir Log.out as a second independent create instead of
copying — std's File::create uses share-all flags on Windows, so plain
concurrent creates cannot collide. Also give every genomeGenerate
invocation in the integration tests a per-test --outFileNamePrefix so
test processes stop sharing ./Log.out at all (and stop littering the
crate directory).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HS4Mudeqj3Y1jk5p3fhUWt
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

No reviews

Assignees

No one assigned

Labels

None yet

Projects

None yet

Milestone

No milestone

Development

Successfully merging this pull request may close these issues.

2 participants

@rob-p @BenjaminDEMAILLE