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..0f5e362 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.** 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 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. + +**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..85c5f10 100644 --- a/src/chimeric/detect.rs +++ b/src/chimeric/detect.rs @@ -662,10 +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 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; @@ -713,6 +716,9 @@ fn place( jr += 1; } + // At least 0: the mate-spacer skip above can leave `jr_best` at `jr_max`, no further. + let new_l1 = jr_max - jr_best; + // Shift the junction. let jr_best = jr_best as usize; let (j0, j1); @@ -731,15 +737,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 +1477,46 @@ 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. STAR keeps such + /// a candidate (the overhang filters drop it later), and so does `place`. + #[test] + fn place_keeps_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"]); + 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 /// 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.