Window bins capped by the index's genomeChrBinNbits; chimeric place() cannot hang

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

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

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;

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, &params(&[])).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]

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");

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(&params).genome_chr_bin_nbits;
     params.redefine_window_params(index.genome.n_genome);
 
     // Use the transcriptome index loaded alongside the genome (populated

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"];

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.