diff --git a/CHANGELOG.md b/CHANGELOG.md index 14b4346b..7c20e6cb 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/Cargo.lock b/Cargo.lock index 9f77c86d..c98c649e 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -174,9 +174,9 @@ checksum = "1e748733b7cbc798e1434b6ac524f0c1ff2ab456fe201501e6497c8417a4fc33" [[package]] name = "caps-sa" -version = "0.6.1" +version = "0.7.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "be62f8675f16876a7f29cfd8b26cc18dcc52eca05a8428d530cc7036b2269f30" +checksum = "e46cbd8870dc17f488813e4c5499ff1bdc2feedb9a650ca2b597933943101394" dependencies = [ "rayon", "tempfile", diff --git a/Cargo.toml b/Cargo.toml index 8a4638f9..831e694d 100644 --- a/Cargo.toml +++ b/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 diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 66f40888..8f154482 100644 --- a/DIVERGENCE.md +++ b/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 `### ` 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 `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. diff --git a/src/genome/mod.rs b/src/genome/mod.rs index 174ba42c..134a229c 100644 --- a/src/genome/mod.rs +++ b/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: `### \n` where commandLineFull is - // " -- -- ...". 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::>() - .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: `### \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::>() + .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))?; diff --git a/src/index/mod.rs b/src/index/mod.rs index 88ee62d0..89fa21b1 100644 --- a/src/index/mod.rs +++ b/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"); diff --git a/src/index/sa_build.rs b/src/index/sa_build.rs index a23a53d5..113e72c9 100644 --- a/src/index/sa_build.rs +++ b/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; diff --git a/src/index/sa_index.rs b/src/index/sa_index.rs index 9761b105..5f15c3af 100644 --- a/src/index/sa_index.rs +++ b/src/index/sa_index.rs @@ -139,6 +139,7 @@ impl SaIndex { } let num_indices = Self::calculate_num_indices(nbases) as usize; let absent_mask: u64 = 1u64 << (gstrand_bit + 2); + let n_mark_mask: u64 = 1u64 << (gstrand_bit + 1); // `isaStep` from STAR: `nSA / 4^nbases`. With 5.9 B SA entries // and `nbases = 14`, this is ~22 — i.e. on average we expect @@ -161,6 +162,20 @@ impl SaIndex { .map(|_| AtomicU64::new(u64::MAX)) .collect(); + // STAR's `SAiMarkNbit` marks: when a suffix's k-mer prefix hits an + // N at level `iL4`, STAR ORs `SAiMarkNmaskC` onto the *last written* + // entry `SAi[start[iL1] + ind0[iL1]]` for every level `iL1 >= iL4` + // (`genomeSAindex.cpp:140-145`). The target slots are collected here + // as a shared bitset and OR-ed into `firsts` after the gap-fill. + // Genome-only indexes come out byte-identical to STAR's; on sjdb + // builds a handful of mark bits can land one slot away from STAR's + // because STAR places them against the pre-insertion base genome — + // see DIVERGENCE.md §3.1a. + let marks: Vec = (0..num_indices.div_ceil(64)) + .into_par_iter() + .map(|_| AtomicU64::new(0)) + .collect(); + let sa_mask: u64 = if sa_word_length == 64 { u64::MAX } else { @@ -306,11 +321,75 @@ impl SaIndex { (hi, next_kmer, next_il4) }; + // Absolute-index variant of `calc_kmer` for the chunk- + // boundary backward scan below: reads one SA entry at an + // absolute SA index via its own small pread (page-cached, + // and only a handful of entries are ever visited). + let calc_kmer_abs = |idx: usize| -> std::io::Result<(u64, i32)> { + let bit = idx as u64 * sa_word_length as u64; + let byte = bit / 8; + let shift = (bit % 8) as u32; + let mut b = [0u8; 16]; + let _ = read_at(sa_file, &mut b, byte)?; + let word = u64::from_le_bytes(b[0..8].try_into().unwrap()); + let packed = (word >> shift) & sa_mask; + let pos = packed & gstrand_mask; + let is_reverse = (packed >> gstrand_bit) != 0; + let genome_pos = if is_reverse { + pos as usize + n_genome + } else { + pos as usize + }; + let mut kmer: u64 = 0; + for ii in 0..nbases as usize { + if genome_pos + ii >= genome_seq.len() { + return Ok((kmer << (2 * (nbases as usize - ii)), ii as i32)); + } + let g = genome_seq[genome_pos + ii]; + if g >= 4 { + return Ok((kmer << (2 * (nbases as usize - ii)), ii as i32)); + } + kmer = (kmer << 2) | (g as u64); + } + Ok((kmer, -1)) + }; + // Per-chunk last-written kmer index at each level. - // `None` means "nothing written in this chunk yet"; - // the first iteration always writes (matches STAR's - // `isa == 0` special-case). + // `None` means "nothing written yet anywhere before this + // point" (only possible at the very start of the SA; + // matches STAR's `isa == 0` special-case). + // + // For chunks after the first, seed `ind0_local` with the + // state STAR's serial scan would have on entering this + // chunk: for each level `iL`, the level-`iL` prefix of the + // most recent suffix before `chunk_start` whose first + // `iL+1` bases contain no N. Suffix prefixes are sorted and + // level validity is prefix-monotone, so a short backward + // walk resolves all levels (typically 1-2 entries). Without + // this seed, the N-mark targets at chunk boundaries would + // be unknown and cross-chunk `SAiMarkNbit` marks would be + // dropped. let mut ind0_local: [Option; 32] = [None; 32]; + { + let mut lo = 0usize; // levels < lo are resolved + let mut j = chunk_start; + while lo < nbases as usize && j > 0 { + j -= 1; + let (kmer, il4) = calc_kmer_abs(j)?; + let valid = if il4 < 0 { + nbases as usize + } else { + il4 as usize + }; + if valid > lo { + for (il, slot) in ind0_local.iter_mut().enumerate().take(valid).skip(lo) + { + *slot = Some(kmer >> (2 * (nbases as usize - 1 - il))); + } + lo = valid; + } + } + } let mut i: usize = 0; let (mut ind_full, mut il4) = calc_kmer(i); @@ -318,15 +397,25 @@ impl SaIndex { let sa_idx = (chunk_start + i) as u64; for il in 0..nbases as usize { if il as i32 == il4 { - // N at level `il`. STAR sets the N flag - // on `ind0[il1]` for `il1 >= il4`; we - // don't track the N flag in this version - // (our `hierarchical_lookup` doesn't - // consult it), so just break out of the - // level loop. This means our SAindex - // file is not byte-identical to STAR's - // in the N-bit positions — documented - // limitation. + // N at level `il`: STAR ORs `SAiMarkNmaskC` + // onto the last-written entry at every level + // `il1 >= il4` (`genomeSAindex.cpp:140-145`). + // Record the target slots; the bit is OR-ed + // into the packed values after the gap-fill. + // `None` (nothing written yet at that level + // anywhere in the SA) is skipped — STAR never + // hits that case on a real genome because the + // lexicographically smallest suffixes are + // N-free. + for (il1, slot0) in + ind0_local.iter().enumerate().take(nbases as usize).skip(il) + { + if let Some(prev) = slot0 { + let slot = (genome_sa_index_start[il1] + prev) as usize; + marks[slot / 64] + .fetch_or(1u64 << (slot % 64), Ordering::Relaxed); + } + } break; } let ind_pref = ind_full >> (2 * (nbases as usize - 1 - il)); @@ -385,9 +474,24 @@ impl SaIndex { } } + // Apply the collected `SAiMarkNbit` marks. Marked slots are always + // present slots (they were the last-written `ind0` entry when the + // mark was recorded), so this ORs the N bit onto a first-occurrence + // `sa_idx` value — exactly STAR's `SAi[..] | SAiMarkNmaskC`. + for (w, mword) in marks.iter().enumerate() { + let mut m = mword.load(Ordering::Relaxed); + while m != 0 { + let slot = w * 64 + m.trailing_zeros() as usize; + firsts[slot].fetch_or(n_mark_mask, Ordering::Relaxed); + m &= m - 1; + } + } + drop(marks); + // Final sequential pack into the output `PackedArray`. // Every slot in `firsts` is now valid (either the - // first-occurrence `sa_idx` or `next | absent_mask`). + // first-occurrence `sa_idx`, possibly with the N-mark bit, + // or `next | absent_mask`). let mut data = PackedArray::new(sai_word_length, num_indices); for (i, slot) in firsts.iter().enumerate() { data.write(i, slot.load(Ordering::Relaxed)); @@ -466,21 +570,25 @@ impl SaIndex { "Building SA index: nbases={nbases}, num_indices={num_indices}, word_length={word_length}" ); - // Initialize packed array with "absent" markers - let mut data = PackedArray::new(word_length, num_indices as usize); - let absent_marker = (1u64 << (gstrand_bit + 2)) | ((1u64 << gstrand_bit) - 1); - - for i in 0..num_indices as usize { - data.write(i, absent_marker); - } + let absent_mask: u64 = 1u64 << (gstrand_bit + 2); + let n_mark_mask: u64 = 1u64 << (gstrand_bit + 1); + + // STAR-exact serial build (`genomeSAindex.cpp` reference + // algorithm): record the first-occurrence `sa_idx` of every + // present k-mer, collect `SAiMarkNbit` targets when a suffix's + // prefix hits an N, then gap-fill absent slots with + // `next_present_sa_idx | absent_mask` and tail-fill with + // `n_sa | absent_mask` — the same encoding `build_parallel` + // emits, byte-identical to STAR's SAindex. + let mut firsts: Vec = vec![u64::MAX; num_indices as usize]; + let mut marks: Vec = vec![0u64; (num_indices as usize).div_ceil(64)]; + let mut ind0: [Option; 32] = [None; 32]; - // Iterate through SA and record first occurrence of each k-mer. // Inner k-loop maintains `kmer_idx` **incrementally** — one - // base read per k iteration (vs. the original `O(k²)` read - // pattern that re-scanned the prefix for every k). When an N - // is encountered at position `genome_pos + (k - 1)` we - // `break` rather than `continue`: every longer k-mer at this - // same `genome_pos` necessarily includes that N too. + // base read per k iteration. When an N (or the genome end) is + // hit at level `k - 1`, STAR ORs the N-mark onto the last + // written entry of every level `>= k - 1` and stops: every + // longer k-mer at this `genome_pos` includes that N too. for sa_idx in 0..sa.len() { let sa_entry = sa.get(sa_idx); let (pos, is_reverse) = sa.decode(sa_entry); @@ -492,23 +600,64 @@ impl SaIndex { let mut kmer_idx: u64 = 0; for k in 1..=nbases { - if genome_pos + (k as usize) > genome.sequence.len() { + let il = (k - 1) as usize; + let past_end = genome_pos + (k as usize) > genome.sequence.len(); + if past_end || genome.sequence.base(genome_pos + il) >= 4 { + for (il1, prev0) in ind0.iter().enumerate().take(nbases as usize).skip(il) { + if let Some(prev) = prev0 { + let slot = (genome_sa_index_start[il1] + prev) as usize; + marks[slot / 64] |= 1u64 << (slot % 64); + } + } break; } - let next_base = genome.sequence.base(genome_pos + (k - 1) as usize); - if next_base >= 4 { - break; + kmer_idx = (kmer_idx << 2) | (genome.sequence.base(genome_pos + il) as u64); + + let is_new = match ind0[il] { + None => true, + Some(prev) => kmer_idx > prev, + }; + if is_new { + let sai_pos = (genome_sa_index_start[il] + kmer_idx) as usize; + if firsts[sai_pos] == u64::MAX { + firsts[sai_pos] = sa_idx as u64; + } + ind0[il] = Some(kmer_idx); } - kmer_idx = (kmer_idx << 2) | (next_base as u64); + } + } - let sai_pos = genome_sa_index_start[(k - 1) as usize] + kmer_idx; - let current_entry = data.read(sai_pos as usize); - let is_absent = (current_entry >> (gstrand_bit + 2)) & 1 != 0; - if is_absent { - data.write(sai_pos as usize, sa_idx as u64); + // Backward gap-fill per level, then apply N-marks, then pack. + for (il, &level_start_raw) in genome_sa_index_start + .iter() + .enumerate() + .take(nbases as usize) + { + let level_start = level_start_raw as usize; + let level_size = 4u64.pow(il as u32 + 1) as usize; + let mut next_present: u64 = sa.len() as u64; + for off in (0..level_size).rev() { + let slot = level_start + off; + if firsts[slot] == u64::MAX { + firsts[slot] = next_present | absent_mask; + } else { + next_present = firsts[slot]; } } } + for (w, &mword) in marks.iter().enumerate() { + let mut m = mword; + while m != 0 { + let slot = w * 64 + m.trailing_zeros() as usize; + firsts[slot] |= n_mark_mask; + m &= m - 1; + } + } + + let mut data = PackedArray::new(word_length, num_indices as usize); + for (i, &v) in firsts.iter().enumerate() { + data.write(i, v); + } Ok(SaIndex { nbases, diff --git a/src/io/log.rs b/src/io/log.rs index 5757687d..4cd4df61 100644 --- a/src/io/log.rs +++ b/src/io/log.rs @@ -135,6 +135,76 @@ pub fn write_log_out( Ok(()) } +/// Write a STAR-shaped `Log.out` for `genomeGenerate`. +/// +/// STAR writes its run log to `Log.out` during +/// `genomeGenerate` and copies it into the genome directory at the end +/// (`genomeGenerate.cpp`), so a STAR-built index directory always +/// contains a `Log.out`. The content is a free-form run log (timestamps, +/// disk-space notes) that can never byte-match across runs; this mirrors +/// the structure — version header, command-line/parameter sections, and +/// the final `DONE: Genome generation, EXITING` line — so tooling that +/// expects the file finds a familiar shape. +pub fn write_genome_generate_log( + path: &Path, + params: &Parameters, + time_start: chrono::DateTime, + time_finish: chrono::DateTime, +) -> std::io::Result<()> { + let file = std::fs::File::create(path)?; + let mut out = BufWriter::new(file); + + let short_fmt = "%b %e %H:%M:%S"; // "Feb 10 17:11:26" + + writeln!(out, "STAR version={}", env!("CARGO_PKG_VERSION"))?; + writeln!( + out, + "STAR compilation time,server,dir={} :", + time_start.format("%Y-%m-%dT%H:%M:%S%:z") + )?; + writeln!(out, "STAR git: ")?; + + let cmd = params.command_line.as_deref().unwrap_or(""); + let pairs = cli_params(cmd); + + writeln!(out, "##### Command Line:")?; + writeln!(out, "{cmd}")?; + + writeln!(out, "###### All USER parameters from Command Line:")?; + for (k, v) in &pairs { + writeln!(out, "{k:<30}{v} ~RE-DEFINED")?; + } + writeln!(out, "##### Finished reading parameters from all sources")?; + writeln!(out)?; + writeln!( + out, + "##### Final user re-defined parameters-----------------:" + )?; + for (k, v) in &pairs { + writeln!(out, "{k:<34}{v}")?; + } + writeln!(out)?; + writeln!(out, "-------------------------------")?; + writeln!(out, "##### Final effective command line:")?; + writeln!(out, "{cmd}")?; + writeln!(out, "----------------------------------------")?; + writeln!(out)?; + + writeln!( + out, + "{} ... starting to generate Genome files", + time_start.format(short_fmt) + )?; + writeln!( + out, + "{} ..... finished successfully", + time_finish.format(short_fmt) + )?; + writeln!(out, "DONE: Genome generation, EXITING")?; + + Ok(()) +} + /// Write STAR-compatible `Log.progress.out`. /// /// STAR updates this file periodically during alignment; for short runs (and diff --git a/src/junction/gtf.rs b/src/junction/gtf.rs index 411b4f5b..9bd8703a 100644 --- a/src/junction/gtf.rs +++ b/src/junction/gtf.rs @@ -201,16 +201,28 @@ pub fn extract_junctions_configured( let exon1 = &exons[i]; let exon2 = &exons[i + 1]; - let intron_start_local_1b = exon1.end + 1; - let intron_end_local_1b = exon2.start.saturating_sub(1); - - if intron_end_local_1b <= intron_start_local_1b { - log::warn!( - "Invalid junction coordinates: {intron_start_local_1b}-{intron_end_local_1b} (possibly overlapping exons)" - ); + // STAR (`GTF_transcriptGeneSJ.cpp:123-134`): touching exons + // (`exS <= exE+1`) silently produce no junction; overlapping + // exons (`exS <= exE`) additionally warn. Anything else is a + // junction — including 1-base introns (`exS == exE+2`), which + // STAR keeps. + if exon2.start <= exon1.end + 1 { + if exon2.start <= exon1.end { + log::warn!( + "Overlapping exons in GTF: {}:{}-{} and {}-{}", + exon1.seqname, + exon1.start, + exon1.end, + exon2.start, + exon2.end + ); + } continue; } + let intron_start_local_1b = exon1.end + 1; + let intron_end_local_1b = exon2.start - 1; + let intron_start = chr_off + intron_start_local_1b - 1; let intron_end = chr_off + intron_end_local_1b - 1; diff --git a/src/junction/sjdb_insert.rs b/src/junction/sjdb_insert.rs index e358622f..dbfb3fb6 100644 --- a/src/junction/sjdb_insert.rs +++ b/src/junction/sjdb_insert.rs @@ -113,6 +113,7 @@ pub fn read_sjdb_info_tab(path: &Path, genome: &Genome) -> Result db_strand, + _ => 0, + }, } } -/// Sort a prepared junction list into STAR's post-dedup order and apply -/// the cross-strand deduplication that STAR does after its second sort -/// (`sjdbPrepare.cpp:141-192`). +/// Sort a prepared junction list into STAR's post-dedup order, +/// replicating both dedup passes of `sjdbPrepare.cpp`: /// -/// STAR's first-pass (intra-strand) dedup collapses duplicate sjdb -/// entries from the same source at the same `(start, end, strand)`. -/// rustar-aligner's `SpliceJunctionDb` already deduplicates on that key at the -/// HashMap level, so those first-pass branches never trigger here; the -/// second-pass cross-strand collision dedup does. +/// **Pass 1** (`sjdbPrepare.cpp:75-123`): sort by left-shifted +/// coordinates, partitioned by the raw source strand (`'+'`, `'-'`, +/// `'.'` sort as separate blocks). Junctions whose shifted coordinates +/// coincide within a strand block are alternative representations of +/// the same splice event inside a repeat; STAR keeps one — preferring +/// canonical motifs, then the smallest left shift. (STAR also compares +/// source priority here; every rustar source currently shares one +/// priority, so those branches are omitted.) /// -/// Dedup rules when two surviving junctions share `(stored_start, -/// stored_end)` but have different strand assignments: +/// **Pass 2** (`sjdbPrepare.cpp:125-191`): re-sort survivors by their +/// stored (motif-restored) coordinates and collapse entries that share +/// `(stored_start, stored_end)` across strand blocks: /// -/// - Undefined strand vs defined strand → keep the defined-strand one. +/// - Defined-strand entry beats a `'.'`-source entry. /// - Both non-canonical → collapse to a single entry with strand = 0 /// (undefined). /// - One canonical + one not → keep the canonical one. /// - Both canonical but on correct vs wrong strand relative to motif — /// keep the one whose strand matches `2 - motif % 2`. pub fn sort_and_dedup(mut junctions: Vec) -> Vec { + // Pass 1: shifted coords, strand-partitioned ('+' block, then '-', + // then '.', mirroring STAR's `shift1` of 0 / nGenomeReal / 2n). + let strand_block = |j: &PreparedJunction| match j.src_strand { + 1 => 0u8, + 2 => 1u8, + _ => 2u8, + }; + junctions.sort_by(|a, b| { + strand_block(a) + .cmp(&strand_block(b)) + .then_with(|| a.start_pos.cmp(&b.start_pos)) + .then_with(|| a.end_pos.cmp(&b.end_pos)) + }); + + let mut pass1: Vec = Vec::with_capacity(junctions.len()); + for j in junctions { + match pass1.last_mut() { + Some(last) + if strand_block(last) == strand_block(&j) + && last.start_pos == j.start_pos + && last.end_pos == j.end_pos => + { + // sjdbPrepare.cpp:116-121 (equal priority): the new + // junction wins if it is canonical and the old one is + // not, or if both have the same canonicality and the + // new one has the smaller left shift. + if (j.motif > 0 && last.motif == 0) + || ((j.motif > 0) == (last.motif > 0) && j.shift_left < last.shift_left) + { + *last = j; + } + } + _ => pass1.push(j), + } + } + + // Pass 2: stored coords, cross-strand collapse. + let mut junctions = pass1; junctions.sort_by(|a, b| { a.stored_start() .cmp(&b.stored_start()) @@ -397,31 +450,27 @@ pub fn sort_and_dedup(mut junctions: Vec) -> Vec Option { - // Strand 0 = undefined. - if old.strand > 0 && new.strand == 0 { - return None; // keep old + // sjdbPrepare.cpp:154-159 — STAR compares the RECORDED strand of the + // old junction (`mapGen.sjdbStrand`, our derived `strand`) against the + // RAW source strand of the new one (`sjdbLoci.str`, our `src_strand`). + if old.strand > 0 && new.src_strand == 0 { + return None; // new junction strand is not defined — keep old } - if old.strand == 0 && new.strand > 0 { - return Some(new.clone()); // replace + if old.strand == 0 && new.src_strand != 0 { + return Some(new.clone()); // old junction strand is not defined — replace } - // Both non-canonical → collapse to undefined strand on the old one. + // Both non-canonical → collapse to undefined strand on the old one + // (sjdbPrepare.cpp:160-163). if old.motif == 0 && new.motif == 0 { let mut merged = old.clone(); merged.strand = 0; return Some(merged); } - // One canonical, one not: prefer canonical. - if old.motif > 0 && new.motif == 0 { - return None; - } - if old.motif == 0 && new.motif > 0 { - return Some(new.clone()); - } - // Both canonical with defined strands. Keep the one on the correct - // strand for its motif (2 - motif % 2). If the old one is on the - // correct strand, skip the new one; otherwise replace. - let old_expected = 2 - (old.motif % 2); - if old.strand == old_expected { + // sjdbPrepare.cpp:164-170: keep the old junction when it is canonical + // and the new one is not, OR when the old junction sits on the correct + // strand for its motif (`motif % 2 == 2 - strand`); otherwise the new + // junction is on the correct strand and replaces it. + if (old.motif > 0 && new.motif == 0) || (old.motif % 2 == 2 - old.strand) { None } else { Some(new.clone()) @@ -735,6 +784,7 @@ mod tests { shift_left, shift_right: 0, strand, + src_strand: strand, } } @@ -917,6 +967,7 @@ mod tests { shift_left: 0, shift_right: 1, strand: 1, + src_strand: 1, }, // Non-canonical: stored = shifted (139_187..139_217). PreparedJunction { @@ -927,6 +978,7 @@ mod tests { shift_left: 0, shift_right: 0, strand: 1, + src_strand: 1, }, ]; write_sjdb_info_tab(tmp.path(), &junctions, 99).unwrap(); @@ -951,6 +1003,7 @@ mod tests { shift_left: 0, shift_right: 0, strand: 1, + src_strand: 1, }; // Non-canonical with shift_left=3 — STAR writes // `stored + shift_left + 1`, which is `original + 1`. @@ -962,6 +1015,7 @@ mod tests { shift_left: 3, shift_right: 0, strand: 0, + src_strand: 0, }; write_sjdb_list_out_tab(tmp.path(), &[canon, noncan], &genome).unwrap(); let bytes = std::fs::read(tmp.path()).unwrap(); @@ -992,6 +1046,7 @@ mod tests { shift_left: 0, shift_right: 0, strand: 2, + src_strand: 2, }; write_sjdb_list_out_tab(tmp.path(), &[pj_b], &genome).unwrap(); let bytes = std::fs::read(tmp.path()).unwrap(); @@ -1001,19 +1056,39 @@ mod tests { #[test] fn dedup_prefers_strand_matching_motif() { - // motif=1 (GT/AG +) stored on wrong strand (2) vs a competing - // motif=2 (CT/AC -) on its correct strand. Same stored coords. - // STAR keeps the one whose strand matches `2 - motif%2`. - // For motif=1: expected strand = 2 - 1%2 = 1. - // For motif=2: expected strand = 2 - 2%2 = 2. - let old_wrong = pj(0, 100, 200, 1, 0, 2); // motif 1 wants strand 1, has 2 - let new_right = pj(0, 100, 200, 2, 0, 2); // motif 2 wants strand 2, has 2 + // Same stored coords reached from different source-strand blocks + // (so pass 1 keeps both and the cross-strand pass 2 decides): + // motif=2 (CT/AC, wants strand 2) annotated on '+' vs the same + // motif annotated on '-'. STAR keeps the one whose strand matches + // `2 - motif%2` (sjdbPrepare.cpp:164-170). + let old_wrong = pj(0, 100, 200, 2, 0, 1); // motif 2 wants strand 2, has 1 + let new_right = pj(0, 100, 200, 2, 0, 2); // motif 2 on its correct strand let out = sort_and_dedup(vec![old_wrong, new_right.clone()]); assert_eq!(out.len(), 1); - // `old_wrong` is on wrong strand for its motif — STAR replaces. + // `old_wrong` is on the wrong strand for its motif — STAR replaces. assert_eq!(out[0], new_right); } + #[test] + fn dedup_pass1_merges_shifted_duplicates_same_strand() { + // Two same-strand junctions whose LEFT-SHIFTED coords coincide are + // repeat-shifted copies of one splice event. STAR's first dedup + // pass keeps one — canonical first, then smallest shift_left + // (sjdbPrepare.cpp:107-122). + let shifted_more = pj(0, 100, 200, 1, 5, 1); + let shifted_less = pj(0, 100, 200, 1, 2, 1); + let out = sort_and_dedup(vec![shifted_more, shifted_less.clone()]); + assert_eq!(out.len(), 1); + assert_eq!(out[0], shifted_less); + + // Canonical beats non-canonical regardless of shift. + let noncan = pj(0, 300, 400, 0, 1, 1); + let canon = pj(0, 300, 400, 3, 4, 1); + let out = sort_and_dedup(vec![noncan, canon.clone()]); + assert_eq!(out.len(), 1); + assert_eq!(out[0], canon); + } + #[test] fn decode_gsj_hit_outside_buffer_returns_empty() { let junctions = vec![pj(0, 1000, 2000, 1, 0, 1)]; @@ -1133,6 +1208,7 @@ mod tests { shift_left: 0, shift_right: 1, strand: 1, + src_strand: 1, }, PreparedJunction { chr_idx: 0, @@ -1142,6 +1218,7 @@ mod tests { shift_left: 3, shift_right: 0, strand: 0, + src_strand: 0, }, ]; let tmp = tempfile::NamedTempFile::new().unwrap(); diff --git a/src/lib.rs b/src/lib.rs index 5086fcbb..7b60e9b9 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -108,6 +108,8 @@ fn genome_generate(params: &Parameters) -> anyhow::Result<()> { ); } + let time_start = chrono::Local::now(); + info!("Building genome index (streaming SA + on-the-fly SAindex)..."); // Streaming path: opens SA file early, packs each caps-sa emit // directly to disk + into the SAindex builder, never holding the @@ -137,6 +139,29 @@ fn genome_generate(params: &Parameters) -> anyhow::Result<()> { GenomeIndex::generate_streaming(&orig_params)?; } + // STAR writes `Log.out` to `` during genomeGenerate + // and copies it into the genome directory at the end, so a STAR-built + // index directory always contains one; mirror that. The genomeDir copy + // is a second independent write, not `fs::copy` from the prefix file: + // concurrent genomeGenerate processes sharing a working directory (the + // integration-test harness does this) race on `Log.out`, and on + // Windows `CopyFileEx` opens its source without write sharing, turning + // that race into a sharing-violation error (os error 32). Two plain + // creates use share-all flags and cannot collide. + let time_finish = chrono::Local::now(); + crate::io::log::write_genome_generate_log( + ¶ms.output_path("Log.out"), + params, + time_start, + time_finish, + )?; + crate::io::log::write_genome_generate_log( + ¶ms.genome_dir.join("Log.out"), + params, + time_start, + time_finish, + )?; + info!("Genome generation complete!"); Ok(()) } diff --git a/src/quant/transcriptome.rs b/src/quant/transcriptome.rs index 2feb761e..8a7658f4 100644 --- a/src/quant/transcriptome.rs +++ b/src/quant/transcriptome.rs @@ -572,10 +572,9 @@ impl TranscriptomeIndex { )); } } - // STAR sorts by sjStart only (funCompareUint2 on the first uint64). - // Keep it stable so gene list order across duplicates matches - // transcript-insertion order. - junctions.sort_by_key(|&(s, _, _, _, _)| s); + // STAR sorts by (sjStart, sjEnd) — `funCompareUint2` compares TWO + // uint64s (`GTF_transcriptGeneSJ.cpp:140`). + junctions.sort_by_key(|&(s, e, _, _, _)| (s, e)); let strand_char = |s: u8| match s { 1 => '+', @@ -583,38 +582,42 @@ impl TranscriptomeIndex { _ => '.', }; - // Dedup pass: merge genes across identical (chr, start, end, strand). + // Collapse pass (`GTF_transcriptGeneSJ.cpp:145-158`): a new output + // row starts whenever (start, end, strand) differs from the + // PREVIOUS sorted entry; otherwise the gene joins the current + // row's gene set. STAR stores genes in a `std::set`, so the + // comma-joined list is ascending and duplicate-free. let mut i = 0; while i < junctions.len() { let (sj_start, sj_end, chr_idx, strand, gene1) = junctions[i]; + let mut genes: std::collections::BTreeSet = std::collections::BTreeSet::new(); + genes.insert(gene1); + let mut j = i + 1; + while j < junctions.len() { + let (s2, e2, _c2, st2, g2) = junctions[j]; + if s2 == sj_start && e2 == sj_end && st2 == strand { + genes.insert(g2); + j += 1; + } else { + break; + } + } + let chr_offset = genome.chr_start[chr_idx]; let start_1based = sj_start + 1 - chr_offset; let end_1based = (sj_end + 1) - chr_offset; write!( out, - "{}\t{}\t{}\t{}\t{}", + "{}\t{}\t{}\t{}", genome.chr_name[chr_idx], start_1based, end_1based, strand_char(strand), - gene1 ) .map_err(|e| Error::io(e, &path))?; - - // Append genes from subsequent entries with the same key. - let mut j = i + 1; - let mut seen: std::collections::BTreeSet = std::collections::BTreeSet::new(); - seen.insert(gene1); - while j < junctions.len() { - let (s2, e2, c2, st2, g2) = junctions[j]; - if s2 == sj_start && e2 == sj_end && c2 == chr_idx && st2 == strand { - if seen.insert(g2) { - write!(out, ",{g2}").map_err(|e| Error::io(e, &path))?; - } - j += 1; - } else { - break; - } + for (k, g) in genes.iter().enumerate() { + let sep = if k == 0 { '\t' } else { ',' }; + write!(out, "{sep}{g}").map_err(|e| Error::io(e, &path))?; } writeln!(out).map_err(|e| Error::io(e, &path))?; i = j; diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index 86691d31..dc94067f 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -102,7 +102,12 @@ fn build_index(fasta: &Path, genome_dir: &Path, sa_nbases: &str, gtf: Option<&Pa .arg("--genomeFastaFiles") .arg(fasta) .arg("--genomeSAindexNbases") - .arg(sa_nbases); + .arg(sa_nbases) + // Per-test prefix: genomeGenerate writes `Log.out`, and the + // default `./` prefix would make concurrently running test + // processes share one file in the crate directory. + .arg("--outFileNamePrefix") + .arg(genome_dir.join("run_")); if let Some(g) = gtf { cmd.arg("--sjdbGTFfile") .arg(g) diff --git a/tests/phase9_threading.rs b/tests/phase9_threading.rs index 3386a144..718a8540 100644 --- a/tests/phase9_threading.rs +++ b/tests/phase9_threading.rs @@ -76,6 +76,8 @@ fn test_single_thread_alignment() { .arg(&fasta_path) .arg("--genomeSAindexNbases") .arg("5") + .arg("--outFileNamePrefix") + .arg(genome_dir.join("run_")) .assert() .success(); @@ -135,6 +137,8 @@ fn test_multi_thread_alignment() { .arg(&fasta_path) .arg("--genomeSAindexNbases") .arg("5") + .arg("--outFileNamePrefix") + .arg(genome_dir.join("run_")) .assert() .success(); @@ -191,6 +195,8 @@ fn test_thread_count_consistency() { .arg(&fasta_path) .arg("--genomeSAindexNbases") .arg("5") + .arg("--outFileNamePrefix") + .arg(genome_dir.join("run_")) .assert() .success(); diff --git a/tests/transcriptome_sam.rs b/tests/transcriptome_sam.rs index 17bbfaec..54a4c42d 100644 --- a/tests/transcriptome_sam.rs +++ b/tests/transcriptome_sam.rs @@ -166,6 +166,8 @@ fn transcriptome_sam_end_to_end_smoke_test() { gtf_path.to_str().unwrap(), "--genomeSAindexNbases", "5", + "--outFileNamePrefix", + genome_dir.join("run_").to_str().unwrap(), ]) .assert() .success();