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

CHANGELOG.md

@@ -13,6 +13,14 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 ### Other changes
 
+- Production genome indexing now enables caps-sa 0.7's bounded geometric LCP
+  memoization through its stable policy API. On the complete ruSTAR-shaped
+  GRCh38 plus GENCODE v50 fixture (6.56 billion symbols, 6.18 billion retained
+  suffixes, 1.40 million segments, 32 physical cores), the final caps-sa 0.7
+  implementation built the SA in 172.953 s versus 267.592 s for its original
+  0.7 baseline: 35.4% faster, with peak RSS reduced from 10,512,408 to
+  9,169,892 KiB. The complete output hash was unchanged.
+
 - `cluster_seeds` reuses its window-bin map across reads on a thread instead
   of rebuilding it per read. Merging two windows re-keys every bin in the
   merged span, so the per-read pre-sizing was only a floor and the map

Cargo.toml

@@ -46,7 +46,7 @@ chrono = "0.4"
 tempfile = "3"
 bitflags = { version = "2.12.1", features = ["std"] }
 shlex = "2.0.1"
-caps-sa = "0.6"
+caps-sa = "0.7"
 # mimalloc as the global allocator. Two reasons:
 #   1. **Memory return**: glibc malloc creates one arena per worker
 #      thread (rayon spawns ~num_cpus workers + sub-threads) and

DIVERGENCE.md

@@ -89,16 +89,42 @@ On the 10k yeast PE benchmark, 4 reads differ in alignment score (AS) because ST
 
 **What STAR does.** At `genomeGenerate`, STAR writes a `### <commandLineFull>` header line reproducing the full command line that built the index.
 
-**What rustar-aligner does.** rustar-aligner emits a fixed skeleton containing the parameters it knows at invocation (`--runMode`, `--runThreadN`, `--genomeDir`, `--genomeFastaFiles`, `--genomeSAindexNbases`, `--sjdbGTFfile`, `--sjdbOverhang`). The remaining value lines of `genomeParameters.txt` match STAR's `genomeParametersWrite.cpp` order and tab/space formatting.
+**What rustar-aligner does.** rustar-aligner echoes its own actual command line after `### ` (falling back to a parameter skeleton for API callers constructed without one). Every value line of `genomeParameters.txt` matches STAR's `genomeParametersWrite.cpp` order and tab/space formatting, including the *effective* `sjdbOverhang` (0 when the index has no sjdb, mirroring `mapGen.sjdbOverhang`).
 
-**Why.** The header is informational; reproducing an arbitrary STAR invocation's exact argv byte-for-byte serves no functional purpose and the index loads identically either way.
+**Why.** The header is informational; the binary path and argument spacing can never byte-match an arbitrary STAR invocation, and the index loads identically either way.
 
-**Impact.** The `###` header line will not byte-match an arbitrary STAR run. No effect on alignment, index loading, or any downstream tool.
+**Impact.** The `###` header line will not byte-match a STAR run (different `argv[0]` and spacing). Every other line matches byte-for-byte. No effect on alignment, index loading, or any downstream tool.
 
 **Source.** `src/genome/mod.rs` (`genomeParameters.txt` writer).
 
 ---
 
+### 3.1a `SAindex` N-mark bits adjacent to junction-flank k-mers
+
+**What STAR does.** With `--sjdbGTFfile`, STAR builds the base-genome `SAindex` first (`genomeSAindex.cpp`) and then *patches* it while inserting junction-flank suffixes (`sjdbBuildIndex.cpp`). `SAiMarkNbit` marks — "suffixes for this k-mer slot may border an N" — are placed against the **base-genome** k-mer landscape: a mark lands on the last k-mer that was present *before* the junction flanks were inserted, marks are silently dropped when an inserted flank suffix takes over a slot's first-occurrence value (`sjdbBuildIndex.cpp:228-231` overwrites the packed value, flags included), and flank suffixes that touch the inter-junction spacer get marks via a separate T-fill backward-scan rule (`sjdbBuildIndex.cpp:262-284`).
+
+**What rustar-aligner does.** rustar-aligner builds the final genome+flank text in one pass and replicates `genomeSAindex.cpp`'s serial mark semantics over that final text: the mark lands on the last k-mer *present in the final index* before the N-run.
+
+**Why.** On GRCh38 + GENCODE v49 this changes a handful of bits (2 slots out of 357,913,940 on the measured build): exactly the slots where a k-mer became present only via a junction flank. STAR's placement there is an artifact of its incremental patch, not a semantic choice; reproducing it would mean simulating the two-phase build. Both placements are valid conservative markers — the bit only widens seed-search bounds near Ns.
+
+**Impact.** ≤ a few bytes of the ~1.5 GB `SAindex` differ on sjdb builds (indexes built *without* a GTF are byte-identical). STAR loads either file and produces identical alignments (verified on 100k read pairs). No effect on any coordinate, count, or emitted record.
+
+**Source.** `src/index/sa_index.rs` (`build_parallel`, `build`).
+
+---
+
+### 3.1b `Log.out` in the genome directory
+
+**What STAR does.** `genomeGenerate` writes its free-form run log to `<outFileNamePrefix>Log.out` and copies it into the genome directory, so a STAR-built index always contains a `Log.out`.
+
+**What rustar-aligner does.** The same — a STAR-shaped `Log.out` (version header, command-line/parameter sections, phase timestamps, `DONE: Genome generation, EXITING`) is written to the output prefix and copied into the genome directory.
+
+**Impact.** The file's *content* is a run log (timestamps, host-specific paths) and can never byte-match across runs or tools; only its presence and shape are mirrored. Nothing loads it at align time.
+
+**Source.** `src/io/log.rs` (`write_genome_generate_log`), `src/lib.rs` (`genome_generate`).
+
+---
+
 ### 3.2 `CellReads.stats` row order
 
 **What STAR does.** `--soloCellReadStats CB` emits its rows by iterating a libc++ `std::unordered_map`, so the order is a hash-table walk rather than a sort. At the map sizes this produces, libc++ chains new entries at the head of their bucket and walks buckets in order, which comes out as the reverse of each barcode's first appearance in read order.

src/genome/mod.rs

@@ -343,7 +343,18 @@ impl Genome {
     /// - `chrStart.txt` — chromosome start positions + final n_genome entry
     /// - `chrNameLength.txt` — tab-separated name + length
     /// - `genomeParameters.txt` — key-value pairs of genome generation parameters
-    pub fn write_index_files(&self, dir: &Path, params: &Parameters) -> Result<(), Error> {
+    ///
+    /// `effective_sjdb_overhang` is the overhang actually baked into the
+    /// genome (STAR's `mapGen.sjdbOverhang`): `params.sjdb_overhang` when
+    /// sjdb junctions were inserted, `0` when the index has no sjdb —
+    /// STAR writes the effective value, not the parameter, into
+    /// `genomeParameters.txt`, and its loader trusts it at align time.
+    pub fn write_index_files(
+        &self,
+        dir: &Path,
+        params: &Parameters,
+        effective_sjdb_overhang: u32,
+    ) -> Result<(), Error> {
         use std::fs;
         use std::io::Write;
 
@@ -393,7 +404,7 @@ impl Genome {
         // trailing whitespace on vector values). STAR's loader reads these
         // keys via `<<` streaming; the leading `###` comment lines are
         // skipped.
-        self.write_genome_parameters_txt(dir, params)?;
+        self.write_genome_parameters_txt(dir, params, effective_sjdb_overhang)?;
 
         // --genomeTransformType Haploid: the block map for reverse conversion.
         if let Some(blocks) = &self.transform_blocks {
@@ -405,40 +416,42 @@ impl Genome {
         Ok(())
     }
 
-    fn write_genome_parameters_txt(&self, dir: &Path, params: &Parameters) -> Result<(), Error> {
+    fn write_genome_parameters_txt(
+        &self,
+        dir: &Path,
+        params: &Parameters,
+        effective_sjdb_overhang: u32,
+    ) -> Result<(), Error> {
         use std::fs;
         use std::io::Write;
 
         let path = dir.join("genomeParameters.txt");
         let mut f = fs::File::create(&path).map_err(|e| Error::io(e, &path))?;
 
-        // STAR writes: `### <commandLineFull>\n` where commandLineFull is
-        // "<argv[0]>   --<name1> <val1>   --<name2> <val2> ...".  We emit
-        // the same skeleton using our known-at-invocation parameters.
-        // Not exposed for retrospective exact-byte match against an arbitrary
-        // STAR run's commandLineFull — see `DIVERGENCE.md` (§3.1) for the short
-        // list of parameters we echo.
-        let fasta_list = params
-            .genome_fasta_files
-            .iter()
-            .map(|p| p.display().to_string())
-            .collect::<Vec<_>>()
-            .join(" ");
-        let gtf = params
-            .sjdb_gtf_file
-            .as_ref()
-            .map_or_else(|| "-".to_string(), |p| p.display().to_string());
-        writeln!(
-            f,
-            "### STAR   --runMode genomeGenerate      --runThreadN {thr}   --genomeDir {dir}   --genomeFastaFiles {fa}      --genomeSAindexNbases {sai}   --sjdbGTFfile {gtf}   --sjdbOverhang {ov}",
-            thr = params.run_thread_n,
-            dir = dir.display(),
-            fa = fasta_list,
-            sai = params.genome_sa_index_nbases,
-            gtf = gtf,
-            ov = params.sjdb_overhang,
-        )
-        .map_err(|e| Error::io(e, &path))?;
+        // STAR writes: `### <commandLineFull>\n` — an echo of the actual
+        // invocation. Emit the real command line (STAR's loader skips
+        // `###` comment lines, and byte-matching an arbitrary STAR run's
+        // argv is impossible anyway — see `DIVERGENCE.md` §3.1); fall
+        // back to a parameter skeleton for callers constructed without
+        // a command line.
+        if let Some(cmd) = params.command_line.as_deref() {
+            writeln!(f, "### {cmd}").map_err(|e| Error::io(e, &path))?;
+        } else {
+            let fasta_list = params
+                .genome_fasta_files
+                .iter()
+                .map(|p| p.display().to_string())
+                .collect::<Vec<_>>()
+                .join(" ");
+            writeln!(
+                f,
+                "### STAR   --runMode genomeGenerate      --runThreadN {thr}   --genomeDir {dir}   --genomeFastaFiles {fa}",
+                thr = params.run_thread_n,
+                dir = dir.display(),
+                fa = fasta_list,
+            )
+            .map_err(|e| Error::io(e, &path))?;
+        }
 
         // GstrandBit: floor(log2(nGenome + limitSjdbInsertNsj*sjdbLength))+1,
         // clamped at a minimum of 32. STAR's default limitSjdbInsertNsj is
@@ -469,7 +482,7 @@ impl Genome {
         writeln!(f, "genomeTransformType\tNone").map_err(|e| Error::io(e, &path))?;
         writeln!(f, "genomeTransformVCF\t-").map_err(|e| Error::io(e, &path))?;
 
-        writeln!(f, "sjdbOverhang\t{}", params.sjdb_overhang).map_err(|e| Error::io(e, &path))?;
+        writeln!(f, "sjdbOverhang\t{effective_sjdb_overhang}").map_err(|e| Error::io(e, &path))?;
 
         // sjdbFileChrStartEnd: empty vector → `-` plus STAR's trailing space.
         writeln!(f, "sjdbFileChrStartEnd\t- ").map_err(|e| Error::io(e, &path))?;

src/index/mod.rs

@@ -103,14 +103,11 @@ impl GenomeIndex {
     ///    `gsj` to the genome.
     /// 2. Write genome files (`Genome`, `chrInfo`, `genomeParameters.txt`,
     ///    etc.) immediately — no dependency on the SA.
-    /// 3. Open `genome_dir/SA` through a [`PackedStreamWriter`]; build
-    ///    a [`SaIndexBuilder`][sa_index::SaIndexBuilder]. The caps-sa
-    ///    emit callback feeds each entry to **both**, so the SA file
-    ///    grows as construction progresses and the SAindex is built
-    ///    on the fly. Total peak RSS during this phase ≈
-    ///    `genome.sequence` + caps-sa scratch + `SaIndex.data` —
-    ///    ~5 GB on the human genome vs the ~47 GB the in-memory
-    ///    path peaked at.
+    /// 3. Open `genome_dir/SA` through a [`PackedStreamWriter`]; the
+    ///    caps-sa emit callback bit-packs each entry straight into
+    ///    the SA file, so the SA file grows as construction
+    ///    progresses and the ~25 GB SA `PackedArray` never has to be
+    ///    materialised in RAM.
     /// 4. Finalise the SA writer (flush partial-byte + padding).
     /// 5. **Build the SAindex in parallel** from the on-disk SA via
     ///    mmap + [`SaIndex::build_parallel`]. caps-sa's phase-4
@@ -139,7 +136,12 @@ impl GenomeIndex {
         }
 
         log::info!("Writing genome files to {}...", dir.display());
-        genome.write_index_files(dir, params)?;
+        let effective_sjdb_overhang = if prepared_junctions.is_empty() {
+            0
+        } else {
+            params.sjdb_overhang
+        };
+        genome.write_index_files(dir, params, effective_sjdb_overhang)?;
 
         let gstrand_bit = SuffixArray::calculate_gstrand_bit(genome.n_genome);
         let gstrand_mask = (1u64 << gstrand_bit) - 1;
@@ -402,7 +404,8 @@ impl GenomeIndex {
     /// Write index files to directory.
     pub fn write(&self, dir: &Path, params: &Parameters) -> Result<(), Error> {
         // Write genome files
-        self.genome.write_index_files(dir, params)?;
+        self.genome
+            .write_index_files(dir, params, self.sjdb_overhang)?;
 
         // Write SA file
         let sa_path = dir.join("SA");

src/index/sa_build.rs

@@ -633,7 +633,13 @@ fn dispatch_caps_sa_segmented(
     );
 
     if use_ext_mem(n) {
-        let opts = caps_sa_ext_mem_opts(temp_dir);
+        // The production segmented genome-plus-junction layout contains many
+        // long shared contexts across phase-4 partition merges. Enable caps-sa's
+        // bounded, partition-local geometric LCP memoization with its measured
+        // defaults. The policy remains explicit here: caps-sa itself defaults
+        // to the direct kernel for generic inputs.
+        let opts = caps_sa_ext_mem_opts(temp_dir)
+            .lcp_memoization(caps_sa::LcpMemoizationPolicy::geometric());
         // Predicate accepts ACGT only (rejects N at 4, spacer at 5).
         // Borrows `original` via `&[u8]` — `Send + Sync` is satisfied.
         let original_ref: &[u8] = original;