From bb21247def7ce4476b74c97fb60c6c2766609794 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 27 Aug 2026 21:12:40 +0200 Subject: [PATCH] fix(filter): take outFilterMismatchNoverLmax over the mapped length STAR divides the mismatch count by the transcript's mapped length (trBest->rLength in ReadAlign_mappedFilter.cpp). This divided by the read length, which is the denominator of a different parameter, outFilterMismatchNoverReadLmax. The two agree on a fully aligned read and diverge on a soft-clipped one, in the permissive direction: a 50-base alignment of a 100-base read carrying 16 mismatches read as 0.16 rather than as 0.32, so it passed the 0.3 default that STAR uses to reject it. The filter was loosest exactly where the alignment was weakest. Measured on the nf-core/rnaseq test data (50 000 pairs), the change is output-neutral: uniquely mapped 41665 and multi-mapped 961 either way, because the absolute --outFilterMismatchNmax cap of 10 dominates at default settings. It bites when that cap is raised, which is what the unit tests cover. The ratio is now a named function rather than an expression inside the retain closure, so it can be tested directly. Closes #238. --- CHANGELOG.md | 7 +++ src/align/read_align.rs | 107 ++++++++++++++++++++++++++++++++++++++-- 2 files changed, 111 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 14b4346..d090367 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -42,6 +42,13 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Bug fixes +- `--outFilterMismatchNoverLmax` is taken over the *mapped* length, as STAR + takes it (`trBest->rLength` in `ReadAlign_mappedFilter.cpp`), rather than + over the read length. The two agree on a fully aligned read and diverge on a + soft-clipped one, in the permissive direction: a 50-base alignment carrying + 16 mismatches read as 0.16 rather than 0.32, so it passed the 0.3 default + that STAR applies to reject it. Output-neutral on the nf-core/rnaseq test + data, where the absolute `--outFilterMismatchNmax` cap dominates. Closes #238. - 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..13079ab 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -132,6 +132,30 @@ pub enum PairedAlignmentResult { }, } +/// Bases of the read that `t` actually aligns, summed over its exons. +fn mapped_read_length(t: &Transcript) -> usize { + t.exons + .iter() + .map(|e| e.read_end.saturating_sub(e.read_start)) + .sum() +} + +/// Mismatches per aligned base, the ratio `--outFilterMismatchNoverLmax` caps. +/// +/// STAR divides by the transcript's mapped length (`trBest->rLength` in +/// `ReadAlign_mappedFilter.cpp`), not by the read length. The distinction only +/// shows up on a soft-clipped alignment, and it matters in the permissive +/// direction: dividing by the read length made the filter looser the more of +/// the read had been clipped away, so a 50-base alignment carrying 16 +/// mismatches read as 0.16 rather than as 0.32 (#238). +/// +/// `--outFilterMismatchNoverReadLmax` is the parameter that *is* taken over the +/// read length, and it keeps its own denominator. +fn mismatch_rate_over_mapped(t: &Transcript) -> f64 { + let mapped = mapped_read_length(t).max(1); + f64::from(t.n_mismatch) / (mapped as f64) +} + /// Align a read to the genome. /// /// # Algorithm @@ -437,8 +461,13 @@ pub fn align_read( return false; } - // Relative mismatch count (mismatches / read_length) - let mismatch_rate = t.n_mismatch as f64 / read_length; + // Relative mismatch count. STAR divides by the *mapped* length + // (`trBest->rLength` in `ReadAlign_mappedFilter.cpp`), not by the read + // length: a soft-clipped alignment is judged on the part that actually + // aligned. Dividing by the read length made the filter more permissive + // the more of the read was clipped away (#238). + let mapped_length = mapped_read_length(t); + let mismatch_rate = mismatch_rate_over_mapped(t); if mismatch_rate > params.out_filter_mismatch_nover_lmax { *filter_reasons.entry("mismatch_rate").or_insert(0) += 1; log::debug!( @@ -447,7 +476,7 @@ pub fn align_read( mismatch_rate * 100.0, params.out_filter_mismatch_nover_lmax * 100.0, t.n_mismatch, - read_length, + mapped_length, t.score ); return false; @@ -1550,6 +1579,78 @@ pub(crate) fn pe_junctions_consistent(left: &Transcript, right: &Transcript) -> #[cfg(test)] mod tests { use super::*; + + /// STAR's denominator is the mapped length, so a soft-clipped alignment is + /// judged on what aligned. The read-length denominator was looser exactly + /// where the alignment was weakest (#238). + #[test] + fn mismatch_rate_is_taken_over_the_mapped_length() { + // 100-base read, 50 bases aligned, 16 mismatches. + let mut t = make_transcript_for_rate(0, 1_000, 1_050, 0, 50); + t.n_mismatch = 16; + + assert_eq!(mapped_read_length(&t), 50); + assert!((mismatch_rate_over_mapped(&t) - 0.32).abs() < 1e-9); + + // Over the full read the same alignment would read 0.16, which is + // under the 0.3 default and so would have been kept. + let read_length = 100.0; + assert!((f64::from(t.n_mismatch) / read_length - 0.16).abs() < 1e-9); + } + + /// Two exons: the mapped length is their sum, not the genomic span. + #[test] + fn mapped_length_sums_exons_and_ignores_the_intron() { + let mut t = make_transcript_for_rate(0, 1_000, 5_060, 0, 30); + t.exons.push(crate::align::transcript::Exon { + genome_start: 5_000, + genome_end: 5_060, + read_start: 30, + read_end: 90, + i_frag: 0, + }); + t.n_mismatch = 9; + assert_eq!(mapped_read_length(&t), 90); + assert!((mismatch_rate_over_mapped(&t) - 0.1).abs() < 1e-9); + } + + /// A transcript with no exons cannot divide by zero. + #[test] + fn mismatch_rate_of_an_empty_transcript_is_finite() { + let mut t = make_transcript_for_rate(0, 0, 0, 0, 0); + t.exons.clear(); + t.n_mismatch = 3; + assert!(mismatch_rate_over_mapped(&t).is_finite()); + } + + fn make_transcript_for_rate( + chr_idx: usize, + genome_start: u64, + genome_end: u64, + read_start: usize, + read_end: usize, + ) -> Transcript { + Transcript { + chr_idx, + genome_start, + genome_end, + is_reverse: false, + exons: vec![crate::align::transcript::Exon { + genome_start, + genome_end, + read_start, + read_end, + i_frag: 0, + }], + cigar: Vec::new(), + score: 0, + n_mismatch: 0, + n_gap: 0, + n_junction: 0, + junction_motifs: Vec::new(), + junction_annotated: Vec::new(), + } + } use crate::genome::Genome; use crate::index::packed_array::PackedArray; use crate::index::sa_index::SaIndex;