From 94b5a370218d866af13464a7ae1a6f23f0a4bbad Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 3 Nov 2026 00:26:32 +0100 Subject: [PATCH 1/2] Window bins capped by the index; chimeric place() cannot hang redefine_window_params skipped the genomeChrBinNbits cap when both alignIntronMax and alignMatesGapMax were 0, and read the command-line value rather than the index's. Genomes of many small contigs lost reads and reported repeats as unique. The cap now always applies with the index value, and the bin size uses STAR's integer divisions. place() cast a negative block length to usize, which stalled the run on one read. Such a candidate is now rejected. STAR stops with a bus error on the same read (DIVERGENCE.md 1.9, needs maintainer sign-off). Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 12 +++++ DIVERGENCE.md | 12 +++++ src/align/read_align.rs | 5 ++- src/chimeric/detect.rs | 55 +++++++++++++++++++++-- src/index/mod.rs | 2 +- src/lib.rs | 1 + src/params/mod.rs | 90 ++++++++++++++++++++++--------------- tests/alignment_features.rs | 82 ++++++++++++++++++++++++++++++++- 8 files changed, 217 insertions(+), 42 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1ab8f0e..a7c91d8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,18 @@ Sections commonly used: Features, Bug fixes, Other changes. ## [Unreleased] +### Bug fixes + +- **Window bins are capped by the index's `genomeChrBinNbits`, as STAR.** The cap + was skipped when `--alignIntronMax` and `--alignMatesGapMax` were both 0, and + read the command-line value rather than the index's, so genomes of many small + contigs lost reads as "too short" and reported repeats as unique. The bin size + now also uses STAR's integer divisions (`--alignIntronMax 1450` gives 8, not 9). +- **Chimeric detection no longer hangs on one read.** `place()` turned a negative + block length into a huge one when the trailing segment ended before the leading + segment's last block (seen with `--chimJunctionOverhangMin 1`). Such a candidate + is now rejected. STAR crashes on the same read (`DIVERGENCE.md` ยง1.9). + ### Other changes - Third-party notices completed. `LICENSES/` gains the licences of code diff --git a/DIVERGENCE.md b/DIVERGENCE.md index f0890fb..5c7a49a 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -171,6 +171,18 @@ A rustar-aligner extension with no STAR equivalent, off by default (`No`). With **Source.** `src/solo/count.rs` (`write_cellranger_summary`), `src/solo/gene.rs` (`Region`). +### 1.9 A chimeric candidate with no junction to place is rejected + +**Needs maintainer sign-off.** + +**What STAR does.** Observed output, source not read: STAR 2.7.11b stops with a bus error (exit 138) on a read pair with `--chimMultimapNmax 20 --chimJunctionOverhangMin 1`, where rustar-aligner on `main` stalled in `place()` because the trailing segment's new length was not positive. + +**What rustar-aligner does.** It rejects that candidate, reports no chimera for it, and carries on with exit 0. + +**Why.** A run should not die on one read, and there is no STAR output to match. Stopping with an error instead would be a one-line change. + +**Source.** `src/chimeric/detect.rs` (`place`). + --- ## 2. Cases where rustar-aligner outperforms STAR diff --git a/src/align/read_align.rs b/src/align/read_align.rs index 8fc7feb..740ad92 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -1429,8 +1429,9 @@ fn try_pair_transcripts( // Genomic span check, only when alignMatesGapMax is set. With the default 0 STAR // has no span limit beyond the window itself (stitchAlignToTranscript tests // alignMatesGapMax only when it is > 0), so a pair whose span exceeds - // winBinNbits*winAnchorDistNbins (~589 kb) is still valid if both mates sit in - // one window. + // 2^winBinNbits * winAnchorDistNbins (589,824 unless the index's + // genomeChrBinNbits caps winBinNbits) is still valid if both mates sit in one + // window. let span = right.genome_end - left.genome_start; if params.align_mates_gap_max > 0 && span > params.align_mates_gap_max as u64 { return None; diff --git a/src/chimeric/detect.rs b/src/chimeric/detect.rs index 9ba58a8..af414a6 100644 --- a/src/chimeric/detect.rs +++ b/src/chimeric/detect.rs @@ -666,6 +666,11 @@ fn place( let m = ro_start1 + l1; if m > ro_start0 { m - ro_start0 - 1 } else { 0 } }; + // The trailing segment ends before the leading one's last block: no junction + // to scan for, and a length below would go negative. + if jr_max == 0 { + return None; + } let (mut motif, mut jr_best) = (0i32, 0i64); let (mut j_score, mut j_score_best) = (0i32, -999_999i32); let mut jr = 0i64; @@ -713,6 +718,13 @@ fn place( jr += 1; } + // The mate-spacer skip above can leave `jr_best` at `jr_max`, with nothing left + // of the trailing segment. + let new_l1 = jr_max - jr_best; + if new_l1 <= 0 { + return None; + } + // Shift the junction. let jr_best = jr_best as usize; let (j0, j1); @@ -731,15 +743,14 @@ fn place( } { let x = &mut t1.exons[e1]; - let new_l = (ro_start1 + x.l as i64 - ro_start0 - jr_best as i64 - 1) as usize; if t1.str_ == 0 { let d = (ro_start0 + jr_best as i64 + 1 - ro_start1) as usize; x.r += d; x.g += d as u64; - x.l = new_l; + x.l = new_l1 as usize; j1 = x.g.wrapping_sub(1); } else { - x.l = new_l; + x.l = new_l1 as usize; j1 = x.g + x.l as u64; } } @@ -1472,6 +1483,44 @@ mod tests { assert_eq!(capped.tr[0].len(), 2); } + /// The trailing segment's first block ends before the leading segment's last + /// block starts. The junction shift used to make its length negative, cast to + /// a huge `usize`, and the loops over that length never ended. + #[test] + fn place_rejects_a_trailing_block_ending_before_the_last_leading_block() { + let gn = genome(); + let fwd = slice(&gn, 100, 100); + let read = ChimRead::new(&fwd, [100, 0], "r", [&fwd, &[]]); + let t0 = wintr(vec![exon(100, 0, 25, 0), exon(200, 25, 35, 0)], vec![1], 10); + let t1 = wintr( + vec![ + exon(CHR_PAD + 100, 15, 5, 0), + exon(CHR_PAD + 200, 25, 75, 0), + ], + vec![1], + 10, + ); + assert!(place(t0, t1, 0, &read, &gn, ¶ms(&[])).is_none()); + } + + /// The scan skips offset 30 (mate 1's length) and so can settle on + /// `jr_best == jr_max`, leaving a trailing block of length 0. + #[test] + fn place_rejects_a_junction_leaving_no_trailing_bases() { + let gn = genome(); + let mut fwd = slice(&gn, 100, 61); + fwd[30] = 5; + let read = ChimRead::new(&fwd, [30, 30], "r", [&fwd[..30], &fwd[31..]]); + // The last scanned base (offset 31) must favour the leading segment. + let g1 = (CHR_PAD + 100..CHR_PAD + 200) + .find(|&g| gn.get_base(g + 11).unwrap() != fwd[31]) + .unwrap(); + let t0 = wintr(vec![exon(100, 0, 40, 0)], vec![], 10); + let t1 = wintr(vec![exon(g1, 20, 12, 0)], vec![], 10); + let prm = params(&["--chimScoreJunctionNonGTAG", "0"]); + assert!(place(t0, t1, 0, &read, &gn, &prm).is_none()); + } + /// A single-end read, half chr0 and half chr1: the old detector pins the /// better half and finds the other as its partner. #[test] diff --git a/src/index/mod.rs b/src/index/mod.rs index badff9d..9e65b77 100644 --- a/src/index/mod.rs +++ b/src/index/mod.rs @@ -508,7 +508,7 @@ impl GenomeIndex { /// `genomeSAindexNbases`, `genomeChrBinNbits` and `genomeSAsparseD` are taken /// from its `genomeParameters.txt` (`Genome_genomeLoad.cpp:104-107`), and /// `genomeFastaFiles` keeps its mapping-time default `-`. -fn loaded_genome_params(params: &Parameters) -> Parameters { +pub(crate) fn loaded_genome_params(params: &Parameters) -> Parameters { let mut p = params.clone(); p.genome_fasta_files = vec![std::path::PathBuf::from("-")]; let path = params.genome_dir.join("genomeParameters.txt"); diff --git a/src/lib.rs b/src/lib.rs index 1d0c993..ef48647 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -411,6 +411,7 @@ fn align_reads(params: &Parameters) -> anyhow::Result<()> { // Redefine window parameters based on genome size (STAR's Genome_genomeLoad.cpp) let mut params = params.clone(); + params.genome_chr_bin_nbits = crate::index::loaded_genome_params(¶ms).genome_chr_bin_nbits; params.redefine_window_params(index.genome.n_genome); // Use the transcriptome index loaded alongside the genome (populated diff --git a/src/params/mod.rs b/src/params/mod.rs index 47e4ecf..f8600dc 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1713,59 +1713,51 @@ impl Parameters { } } - /// Compute the default window distance: 2^winBinNbits * winAnchorDistNbins. - /// Used for max_cluster_dist and as default alignIntronMax (when 0). + /// The window distance: 2^winBinNbits * winAnchorDistNbins. pub fn win_bin_window_dist(&self) -> u64 { (1u64 << self.win_bin_nbits) * self.win_anchor_dist_nbins as u64 } - /// Redefine window parameters based on genome size and intron/gap limits. - /// Ports STAR's Genome_genomeLoad.cpp logic that recomputes winBinNbits, - /// winFlankNbins, and winAnchorDistNbins after loading the genome. + /// Redefine window parameters based on genome size and intron/gap limits, + /// as STAR reports doing after loading the genome (its `Log.out`). /// - /// IMPORTANT: winBinNbits is only redefined when alignIntronMax > 0 OR - /// alignMatesGapMax > 0. When both are 0, winBinNbits stays at its default (16). + /// `genome_chr_bin_nbits` must already hold the index's value. Windows are not + /// widened beyond it, including with alignIntronMax = alignMatesGapMax = 0: bins + /// wider than the padding between contigs would mix contigs in one window. pub fn redefine_window_params(&mut self, n_genome: u64) { let intron_max = self.align_intron_max as u64; let gap_max = self.align_mates_gap_max as u64; + let redefine = intron_max > 0 || gap_max > 0; - if intron_max == 0 && gap_max == 0 { - // STAR: no redefinition when both are 0. Log effective max intron. - let max_intron = (1u64 << self.win_bin_nbits) * self.win_anchor_dist_nbins as u64; - log::info!( - "alignIntronMax=alignMatesGapMax=0, max intron ~= (2^winBinNbits)*winAnchorDistNbins={max_intron}" + if redefine { + // STAR: max(max(4, alignIntronMax), alignMatesGapMax==0 ? 1000 : alignMatesGapMax) + let max_span = std::cmp::max( + std::cmp::max(4u64, intron_max), + if gap_max == 0 { 1000 } else { gap_max }, ); - return; - } - - // STAR: max(max(4, alignIntronMax), alignMatesGapMax==0 ? 1000 : alignMatesGapMax) - let max_span = std::cmp::max( - std::cmp::max(4u64, intron_max), - if gap_max == 0 { 1000 } else { gap_max }, - ); - // winBinNbits = floor(log2(max_span / 4) + 0.5) - self.win_bin_nbits = ((max_span as f64 / 4.0).log2() + 0.5).floor() as u32; + // Integer divisions: --alignIntronMax 1450 gives winBinNbits=8 in STAR's + // Log.out, where 1450/4 = 362.5 would round to 9. + self.win_bin_nbits = (((max_span / 4) as f64).log2() + 0.5).floor() as u32; + let genome_based = (((n_genome / 40000 + 1) as f64).log2() + 0.5).floor() as u32; + self.win_bin_nbits = self.win_bin_nbits.max(genome_based); + } - // max with genome-based value: floor(log2(nGenome/40000 + 1) + 0.5) - let genome_based = ((n_genome as f64 / 40000.0 + 1.0).log2() + 0.5).floor() as u32; - self.win_bin_nbits = self.win_bin_nbits.max(genome_based); + self.win_bin_nbits = self.win_bin_nbits.min(self.genome_chr_bin_nbits); - // Cap at genomeChrBinNbits - if self.win_bin_nbits > self.genome_chr_bin_nbits { - self.win_bin_nbits = self.genome_chr_bin_nbits; + if redefine { + let max_gap = std::cmp::max(intron_max, gap_max); + self.win_flank_nbins = (max_gap / (1u64 << self.win_bin_nbits) + 1) as u32; + self.win_anchor_dist_nbins = 2 * self.win_flank_nbins; } - // Redefine winFlankNbins and winAnchorDistNbins - let max_gap = std::cmp::max(intron_max, gap_max); - self.win_flank_nbins = (max_gap / (1u64 << self.win_bin_nbits) + 1) as u32; - self.win_anchor_dist_nbins = 2 * self.win_flank_nbins; - log::info!( - "Redefined window params: winBinNbits={}, winFlankNbins={}, winAnchorDistNbins={}", + "Window params: winBinNbits={}, winFlankNbins={}, winAnchorDistNbins={}, \ + max window span (2^winBinNbits)*winAnchorDistNbins={}", self.win_bin_nbits, self.win_flank_nbins, - self.win_anchor_dist_nbins + self.win_anchor_dist_nbins, + self.win_bin_window_dist() ); } @@ -3475,6 +3467,34 @@ mod tests { ); } + /// Expected values are those in STAR 2.7.11b's `Log.out` for the same options. + #[test] + fn window_params_follow_star_log() { + let run = |chr_bin: u32, n_genome: u64, extra: &[&str]| { + let mut args = vec!["rustar-aligner", "--readFilesIn", "r.fq"]; + args.extend_from_slice(extra); + let mut p = Parameters::parse_from(args); + p.genome_chr_bin_nbits = chr_bin; + p.redefine_window_params(n_genome); + (p.win_bin_nbits, p.win_flank_nbins, p.win_anchor_dist_nbins) + }; + let n = 1_619_968; + // No intron or gap limit: only the cap applies. + assert_eq!(run(11, n, &[]), (11, 4, 9)); + assert_eq!(run(18, n, &[]), (16, 4, 9)); + assert_eq!(run(11, n, &["--winBinNbits", "5"]), (5, 4, 9)); + // Capped before the flanks are derived. + assert_eq!(run(11, n, &["--alignIntronMax", "100000"]), (11, 49, 98)); + // 1450 / 4 is 362, not 362.5. + assert_eq!(run(11, n, &["--alignIntronMax", "1450"]), (8, 6, 12)); + assert_eq!(run(11, n, &["--alignMatesGapMax", "1450"]), (8, 6, 12)); + // 14450688 / 40000 is 361, not 361.27. + assert_eq!( + run(12, 14_450_688, &["--alignIntronMax", "1000"]), + (8, 4, 8) + ); + } + #[test] fn chim_out_type_follows_star() { let base = ["--readFilesIn", "r.fq", "--chimSegmentMin", "12"]; diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index 7e71c04..0311e6c 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -3154,8 +3154,88 @@ fn align_pair( .collect() } +/// Contigs far smaller than the default 2^16 window bin: an index built with +/// `--genomeChrBinNbits 11` must cap winBinNbits at 11 at mapping time, although +/// the mapping command line does not repeat that option. Otherwise contigs share +/// a bin, reads are lost as "too short" and repeats are reported as unique. +/// Expected numbers are STAR 2.7.11b's on this genome and these reads. +#[test] +fn test_window_bins_capped_by_the_index_chr_bins() { + let tmpdir = TempDir::new().unwrap(); + let sequence = xorshift_seq(0x9E37_79B9_7F4A_7C15, 700_000); + let repeat = &sequence[..200]; + let mut fasta = String::new(); + let mut reads = String::new(); + let mut add_read = |name: String, seq: &[u8]| { + let seq = String::from_utf8_lossy(seq); + reads += &format!("@{name}\n{seq}\n+\n{}\n", "I".repeat(seq.len())); + }; + let mut start = 200; + for i in 0..200usize { + let len = 2000 + (i * 13) % 1500; + let mut contig = sequence[start..start + len].to_vec(); + start += len; + if i < 4 { + contig[500..700].copy_from_slice(repeat); + } + fasta += &format!(">c{i}\n{}\n", String::from_utf8_lossy(&contig)); + for st in [0, 300, 900, len - 100] { + let read = &contig[st..st + 100]; + if i % 2 == 0 { + add_read(format!("r{i}_{st}"), read); + } else { + add_read(format!("r{i}_{st}"), &rc(read)); + } + } + } + add_read("rep".to_string(), &repeat[20..120]); + let fa_path = tmpdir.path().join("contigs.fa"); + fs::write(&fa_path, fasta).unwrap(); + let fq_path = tmpdir.path().join("reads.fq"); + fs::write(&fq_path, reads).unwrap(); + let genome_dir = tmpdir.path().join("idx"); + fs::create_dir_all(&genome_dir).unwrap(); + cargo_bin_cmd!("rustar-aligner") + .arg("--runMode") + .arg("genomeGenerate") + .arg("--genomeDir") + .arg(&genome_dir) + .arg("--genomeFastaFiles") + .arg(&fa_path) + .args(["--genomeSAindexNbases", "8", "--genomeChrBinNbits", "11"]) + .arg("--outFileNamePrefix") + .arg(genome_dir.join("run_")) + .assert() + .success(); + let prefix = tmpdir.path().join("out_"); + cargo_bin_cmd!("rustar-aligner") + .arg("--genomeDir") + .arg(&genome_dir) + .arg("--readFilesIn") + .arg(&fq_path) + .arg("--outFileNamePrefix") + .arg(&prefix) + .assert() + .success(); + + let log = fs::read_to_string(tmpdir.path().join("out_Log.final.out")).unwrap(); + let value = |label: &str| { + log.lines() + .find(|l| l.contains(label)) + .and_then(|l| l.split('|').nth(1)) + .map(|v| v.trim().to_string()) + .unwrap() + }; + assert_eq!(value("Uniquely mapped reads number"), "800"); + assert_eq!(value("Number of reads mapped to multiple loci"), "1"); + assert_eq!(value("Number of reads unmapped: too short"), "0"); + let sam = fs::read_to_string(tmpdir.path().join("out_Aligned.out.sam")).unwrap(); + let rep = sam.lines().filter(|l| l.starts_with("rep\t")).count(); + assert_eq!(rep, 4); +} + /// A pair whose mates are both spliced across 420 kb introns spans 840 kb, -/// more than winBinNbits*winAnchorDistNbins (589,824). With alignIntronMax and +/// more than 2^winBinNbits * winAnchorDistNbins (589,824). With alignIntronMax and /// alignMatesGapMax at 0, STAR puts the whole pair in one window and reports it; /// rustar-aligner used to reject any pair spanning more than that distance. /// Expected alignments are STAR 2.7.11b's on this genome. From 74d8a023f5589cb33885f25ff8af553b374c4dc4 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 3 Nov 2026 09:04:10 +0100 Subject: [PATCH 2/2] Chimeric place(): reject only a negative trailing length, as STAR crashes The guard also rejected a trailing block of length exactly 0, which STAR keeps (its jRmax clamps to 0 and the new length is 0, not negative). That only mattered with --chimJunctionOverhangMin 0 or on the multimap path with min_score <= 0. Only the case whose length would wrap is rejected now, and DIVERGENCE.md 1.9 cites STAR's source for it. Co-Authored-By: Claude Opus 5.5 --- DIVERGENCE.md | 4 ++-- src/chimeric/detect.rs | 26 +++++++++++--------------- 2 files changed, 13 insertions(+), 17 deletions(-) diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 5c7a49a..0f5e362 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -175,9 +175,9 @@ A rustar-aligner extension with no STAR equivalent, off by default (`No`). With **Needs maintainer sign-off.** -**What STAR does.** Observed output, source not read: STAR 2.7.11b stops with a bus error (exit 138) on a read pair with `--chimMultimapNmax 20 --chimJunctionOverhangMin 1`, where rustar-aligner on `main` stalled in `place()` because the trailing segment's new length was not positive. +**What STAR does.** When the trailing segment's first block ends before the leading segment's last block starts, `jRmax` is clamped to 0 and the new block length `roStart1+L-roStart0-jRbest-1` (`ReadAlign_chimericDetectionOld.cpp:154-155`, `:253-256`; the same formula in `ChimericAlign_chimericStitching.cpp`) is negative in a `uint`, so it wraps. STAR 2.7.11b then stops with a bus error (exit 138) on a read pair with `--chimMultimapNmax 20 --chimJunctionOverhangMin 1`, where rustar-aligner on `main` stalled in `place()` on the same wrapped length. A block that ends exactly where the other starts gets length 0, which STAR keeps. -**What rustar-aligner does.** It rejects that candidate, reports no chimera for it, and carries on with exit 0. +**What rustar-aligner does.** It rejects only the negative case, reports no chimera for it, and carries on with exit 0. The length-0 case is kept, as in STAR, and the overhang filters decide it. **Why.** A run should not die on one read, and there is no STAR output to match. Stopping with an error instead would be a one-line change. diff --git a/src/chimeric/detect.rs b/src/chimeric/detect.rs index af414a6..85c5f10 100644 --- a/src/chimeric/detect.rs +++ b/src/chimeric/detect.rs @@ -662,15 +662,13 @@ fn place( } }; - let jr_max = { - let m = ro_start1 + l1; - if m > ro_start0 { m - ro_start0 - 1 } else { 0 } - }; - // The trailing segment ends before the leading one's last block: no junction - // to scan for, and a length below would go negative. - if jr_max == 0 { + // The trailing segment ends before the leading one's last block starts: its + // new length below would be negative, which crashes STAR (DIVERGENCE.md 1.9). + // Ending exactly there leaves length 0, which STAR keeps, and so do we. + if ro_start1 + l1 <= ro_start0 { return None; } + let jr_max = ro_start1 + l1 - ro_start0 - 1; let (mut motif, mut jr_best) = (0i32, 0i64); let (mut j_score, mut j_score_best) = (0i32, -999_999i32); let mut jr = 0i64; @@ -718,12 +716,8 @@ fn place( jr += 1; } - // The mate-spacer skip above can leave `jr_best` at `jr_max`, with nothing left - // of the trailing segment. + // At least 0: the mate-spacer skip above can leave `jr_best` at `jr_max`, no further. let new_l1 = jr_max - jr_best; - if new_l1 <= 0 { - return None; - } // Shift the junction. let jr_best = jr_best as usize; @@ -1504,9 +1498,10 @@ mod tests { } /// The scan skips offset 30 (mate 1's length) and so can settle on - /// `jr_best == jr_max`, leaving a trailing block of length 0. + /// `jr_best == jr_max`, leaving a trailing block of length 0. STAR keeps such + /// a candidate (the overhang filters drop it later), and so does `place`. #[test] - fn place_rejects_a_junction_leaving_no_trailing_bases() { + fn place_keeps_a_junction_leaving_no_trailing_bases() { let gn = genome(); let mut fwd = slice(&gn, 100, 61); fwd[30] = 5; @@ -1518,7 +1513,8 @@ mod tests { let t0 = wintr(vec![exon(100, 0, 40, 0)], vec![], 10); let t1 = wintr(vec![exon(g1, 20, 12, 0)], vec![], 10); let prm = params(&["--chimScoreJunctionNonGTAG", "0"]); - assert!(place(t0, t1, 0, &read, &gn, &prm).is_none()); + let placed = place(t0, t1, 0, &read, &gn, &prm).expect("kept, as in STAR"); + assert_eq!(placed.t1.exons[0].l, 0); } /// A single-end read, half chr0 and half chr1: the old detector pins the