From f4f4bdcba39cec64514876d9b2825596bd46357b Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 27 Aug 2026 19:21:43 +0200 Subject: [PATCH] fix(stats): paired-end unmapped reads split between too short and other STAR decides the bucket in ReadAlign_mappedFilter.cpp: nW == 0, meaning no good window at all, is unmappedOther; only a read that had a window whose best transcript failed the score or length thresholds is unmappedShort. The single-end path here already made that distinction. The paired-end path returned TooShort for every unmapped pair, so the "other" bucket in Log.final.out was permanently zero and "too short" carried both populations. Measured on the nf-core/rnaseq test data the issue reports (50 000 pairs, S. cerevisiae chrI + GFP transgene, --outFilterMultimapNmax 20): before too short 7374 other 0 after too short 3778 other 3596 STAR too short 3766 other 3609 The totals were already right (7374 against STAR's 7375), which is why this showed up as mis-categorisation rather than as a counting bug. tests/unmapped_buckets.rs plants all three populations: pairs that map, pairs that seed but cannot clear the score threshold, and pairs with no genomic seed at all. Checked against the control: without this change the paired-end test fails and the single-end one still passes, which is the shape the bug had. Closes #48. --- CHANGELOG.md | 8 ++ src/align/read_align.rs | 15 ++- tests/unmapped_buckets.rs | 223 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 245 insertions(+), 1 deletion(-) create mode 100644 tests/unmapped_buckets.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index 14b4346..7d0a24e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -42,6 +42,14 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Bug fixes +- `Log.final.out` splits unmapped paired-end reads between `too short` and + `other` again. STAR calls a read `other` when no good window was found at + all and `too short` only when a window existed whose best transcript failed + the score or length thresholds (`ReadAlign_mappedFilter.cpp`); the + paired-end path reported every unmapped pair as `too short`, so the `other` + bucket was permanently zero. On the nf-core/rnaseq test data (50 000 pairs, + yeast chrI + GFP) this moves rustar from `too short 7374 / other 0` to + `too short 3778 / other 3596`, against STAR's `3766 / 3609`. Closes #48. - Read names are cut at `--readNameSeparator` (default `/`), as STAR does. A read named `foo/1` was previously emitted as `foo/1` where STAR emits `foo`. diff --git a/src/align/read_align.rs b/src/align/read_align.rs index ec66a16..ac262fa 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -1306,7 +1306,20 @@ pub fn align_paired_read( )) } } - (None, None) => Ok((Vec::new(), pe_chimeric, 0, Some(UnmappedReason::TooShort))), + // STAR splits this case in two (`ReadAlign_mappedFilter.cpp`): a read + // with no good window at all is `unmappedOther` ("other"), and only a + // read that *had* a window whose best transcript failed the score or + // length thresholds is `unmappedShort` ("too short"). Reporting every + // unmapped pair as too short empties the "other" bucket, which is what + // #48 measured on paired-end data. + (None, None) => { + let reason = if clusters.is_empty() { + UnmappedReason::Other + } else { + UnmappedReason::TooShort + }; + Ok((Vec::new(), pe_chimeric, 0, Some(reason))) + } } } diff --git a/tests/unmapped_buckets.rs b/tests/unmapped_buckets.rs new file mode 100644 index 0000000..617219c --- /dev/null +++ b/tests/unmapped_buckets.rs @@ -0,0 +1,223 @@ +//! `Log.final.out` has to tell "too short" and "other" apart. +//! +//! STAR splits them in `ReadAlign_mappedFilter.cpp`: a read with no good +//! window at all is `unmappedOther`, and only a read that *had* a window whose +//! best transcript failed the score or length thresholds is `unmappedShort`. +//! Folding both into "too short" leaves the "other" bucket permanently zero, +//! which is what issue #48 measured on paired-end data. + +use assert_cmd::cargo::cargo_bin_cmd; +use std::fs; +use std::io::Write; +use std::path::{Path, PathBuf}; +use tempfile::TempDir; + +const READ_LEN: usize = 100; + +fn lcg_seq(seed: u32, length: usize) -> Vec { + let bases: [u8; 4] = *b"ACGT"; + let mut state = seed; + let mut seq = Vec::with_capacity(length); + for _ in 0..length { + state = state.wrapping_mul(1_103_515_245).wrapping_add(12345); + seq.push(bases[((state >> 16) & 3) as usize]); + } + seq +} + +fn write_fasta(dir: &Path, genome: &[u8]) -> PathBuf { + let path = dir.join("genome.fa"); + let mut f = fs::File::create(&path).unwrap(); + writeln!(f, ">chr1").unwrap(); + f.write_all(genome).unwrap(); + writeln!(f).unwrap(); + path +} + +fn build_index(fasta: &Path, genome_dir: &Path) { + fs::create_dir_all(genome_dir).unwrap(); + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "genomeGenerate", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--genomeFastaFiles", + fasta.to_str().unwrap(), + "--genomeSAindexNbases", + "7", + ]) + .assert() + .success(); +} + +fn parse_final_log(path: &Path) -> Vec<(String, u64)> { + fs::read_to_string(path) + .unwrap() + .lines() + .filter_map(|l| { + let (label, value) = l.split_once('|')?; + let value = value.trim(); + if value.ends_with('%') { + return None; + } + Some((label.trim().to_string(), value.parse::().ok()?)) + }) + .collect() +} + +fn get(rows: &[(String, u64)], label: &str) -> u64 { + rows.iter() + .find(|(l, _)| l == label) + .map(|(_, v)| *v) + .unwrap_or_else(|| panic!("{label} missing from Log.final.out")) +} + +/// Paired reads of three kinds: mapping cleanly, matching the genome over a +/// stretch too short to pass the score filter, and matching nothing at all. +fn write_reads(dir: &Path, genome: &[u8]) -> (PathBuf, PathBuf) { + let p1 = dir.join("r1.fq"); + let p2 = dir.join("r2.fq"); + let mut f1 = fs::File::create(&p1).unwrap(); + let mut f2 = fs::File::create(&p2).unwrap(); + + let mut emit = |name: &str, s1: &[u8], s2: &[u8]| { + writeln!(f1, "@{name}").unwrap(); + f1.write_all(s1).unwrap(); + writeln!(f1, "\n+\n{}", "I".repeat(s1.len())).unwrap(); + writeln!(f2, "@{name}").unwrap(); + f2.write_all(s2).unwrap(); + writeln!(f2, "\n+\n{}", "I".repeat(s2.len())).unwrap(); + }; + + // 20 pairs that map: mate2 is the reverse complement of a downstream slice. + let rc = |s: &[u8]| -> Vec { + s.iter() + .rev() + .map(|&b| match b { + b'A' => b'T', + b'C' => b'G', + b'G' => b'C', + _ => b'A', + }) + .collect() + }; + for i in 0..20usize { + let start = 2_000 + i * 300; + let m1 = &genome[start..start + READ_LEN]; + let m2 = rc(&genome[start + 150..start + 150 + READ_LEN]); + emit(&format!("ok{i}"), m1, &m2); + } + + // 20 pairs that seed but cannot clear the score threshold: 30 genomic + // bases then 70 bases belonging to no chromosome. A window exists, so + // STAR calls these "too short", not "other". + for i in 0..20usize { + let start = 9_000 + i * 137; + let mut m1 = genome[start..start + 30].to_vec(); + m1.extend_from_slice(&lcg_seq(4_242 + i as u32, READ_LEN - 30)); + let mut m2 = genome[start + 200..start + 230].to_vec(); + m2.extend_from_slice(&lcg_seq(9_191 + i as u32, READ_LEN - 30)); + emit(&format!("short{i}"), &m1, &m2); + } + + // 20 pairs with no genomic seed at all, which is the "other" bucket. + for i in 0..20usize { + let m1 = lcg_seq(700_000 + i as u32, READ_LEN); + let m2 = lcg_seq(900_000 + i as u32, READ_LEN); + emit(&format!("other{i}"), &m1, &m2); + } + + (p1, p2) +} + +#[test] +fn paired_end_unmapped_reads_split_between_too_short_and_other() { + let dir = TempDir::new().unwrap(); + let root = dir.path(); + let genome = lcg_seq(88888, 20_000); + let genome_dir = root.join("genome"); + build_index(&write_fasta(root, &genome), &genome_dir); + let (r1, r2) = write_reads(root, &genome); + + let prefix = format!("{}/pe_", root.display()); + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + r1.to_str().unwrap(), + r2.to_str().unwrap(), + "--outFileNamePrefix", + &prefix, + ]) + .assert() + .success(); + + let rows = parse_final_log(Path::new(&format!("{prefix}Log.final.out"))); + let input = get(&rows, "Number of input reads"); + let unique = get(&rows, "Uniquely mapped reads number"); + let too_short = get(&rows, "Number of reads unmapped: too short"); + let other = get(&rows, "Number of reads unmapped: other"); + + assert_eq!(input, 60, "fixture size"); + assert!(unique > 0, "the mapping pairs must map: {rows:?}"); + assert!( + other > 0, + "reads with no genomic seed belong in `other`, which was permanently \ + zero before this split: too_short={too_short}, other={other}" + ); + assert!( + too_short > 0, + "reads that seed but fail the score threshold belong in `too short`: \ + too_short={too_short}, other={other}" + ); + // Every unmapped read lands in exactly one bucket. + let multi = get(&rows, "Number of reads mapped to multiple loci"); + let too_many = get(&rows, "Number of reads mapped to too many loci"); + let mismatches = get(&rows, "Number of reads unmapped: too many mismatches"); + assert_eq!( + unique + multi + too_many + mismatches + too_short + other, + input, + "the buckets have to add up to the input: {rows:?}" + ); +} + +#[test] +fn single_end_unmapped_reads_split_between_too_short_and_other() { + // The single-end path already distinguished the two; this keeps it that + // way, and makes the pair of tests say which path is being checked. + let dir = TempDir::new().unwrap(); + let root = dir.path(); + let genome = lcg_seq(88888, 20_000); + let genome_dir = root.join("genome"); + build_index(&write_fasta(root, &genome), &genome_dir); + let (r1, _r2) = write_reads(root, &genome); + + let prefix = format!("{}/se_", root.display()); + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + r1.to_str().unwrap(), + "--outFileNamePrefix", + &prefix, + ]) + .assert() + .success(); + + let rows = parse_final_log(Path::new(&format!("{prefix}Log.final.out"))); + assert!( + get(&rows, "Number of reads unmapped: other") > 0, + "single-end `other` bucket: {rows:?}" + ); + assert!( + get(&rows, "Number of reads unmapped: too short") > 0, + "single-end `too short` bucket: {rows:?}" + ); +}