Skip to content

Commit 0ab80c5

Browse files
pinin4fjordsclaude
andcommitted
fix(stats): route reads exceeding --outFilterMultimapNmax to too_many_loci
The existing AlignmentStats::record_alignment routing arm for n_alignments > max_multimaps was unreachable on the PE path - the too-many-loci filter returned n_for_mapq = 0, so the caller hit record_alignment(0, ...) and incremented the unmapped bucket. Log.final.out always reported 0 in the too-many-loci field and the per-bucket sum fell 3-of-50000 short of input on the yeast PE test profile (STAR records exactly 3 there on the same data). Return joint_pairs.len() as n_for_mapq from the PE filter site and route through record_alignment with that count, mirroring the SE caller pattern so the _ => arm fires. The read still isn't emitted in the BAM with outSAMmultNmax=-1, matching STAR; only the accounting changes. Fixes #53 Co-Authored-By: Claude <noreply@anthropic.com>
1 parent 70be24d commit 0ab80c5

3 files changed

Lines changed: 57 additions & 3 deletions

File tree

src/align/read_align.rs

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1084,10 +1084,11 @@ pub fn align_paired_read(
10841084

10851085
// Step 3: TooManyLoci check (post-dedup, matching STAR's ordering: multMapSelect → dedup → TooManyLoci).
10861086
if joint_pairs.len() > params.out_filter_multimap_nmax as usize {
1087+
let n_loci = joint_pairs.len();
10871088
return Ok((
10881089
Vec::new(),
10891090
pe_chimeric,
1090-
0,
1091+
n_loci,
10911092
Some(UnmappedReason::TooManyLoci),
10921093
));
10931094
}

src/lib.rs

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1454,8 +1454,8 @@ fn align_reads_paired_end<W: AlignmentWriter + ?Sized>(
14541454
.collect();
14551455

14561456
if results.is_empty() {
1457-
// Both mates unmapped
1458-
stats.record_alignment(0, max_multimaps);
1457+
let n_for_stats = if n_for_mapq > 0 { n_for_mapq } else { 0 };
1458+
stats.record_alignment(n_for_stats, max_multimaps);
14591459
stats.record_unmapped_reason(
14601460
unmapped_reason.unwrap_or(crate::stats::UnmappedReason::Other),
14611461
);

src/stats.rs

Lines changed: 53 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -827,6 +827,59 @@ mod tests {
827827
assert_eq!(stats.unmapped_mismatches.load(Ordering::Relaxed), 1);
828828
}
829829

830+
#[test]
831+
fn test_too_many_loci_routed_and_totals_conserved() {
832+
let stats = AlignmentStats::new();
833+
let max_multimaps = 3;
834+
835+
stats.record_alignment(5, max_multimaps);
836+
stats.record_unmapped_reason(UnmappedReason::TooManyLoci);
837+
838+
assert_eq!(stats.total_reads.load(Ordering::Relaxed), 1);
839+
assert_eq!(stats.too_many_loci.load(Ordering::Relaxed), 1);
840+
assert_eq!(stats.uniquely_mapped.load(Ordering::Relaxed), 0);
841+
assert_eq!(stats.multi_mapped.load(Ordering::Relaxed), 0);
842+
assert_eq!(stats.unmapped.load(Ordering::Relaxed), 0);
843+
assert_eq!(stats.unmapped_other.load(Ordering::Relaxed), 0);
844+
assert_eq!(stats.unmapped_short.load(Ordering::Relaxed), 0);
845+
assert_eq!(stats.unmapped_mismatches.load(Ordering::Relaxed), 0);
846+
}
847+
848+
#[test]
849+
fn test_bucket_sum_conserved_with_too_many_loci() {
850+
let stats = AlignmentStats::new();
851+
let max_multimaps = 10;
852+
853+
for _ in 0..3 {
854+
stats.record_alignment(1, max_multimaps);
855+
}
856+
for _ in 0..2 {
857+
stats.record_alignment(5, max_multimaps);
858+
}
859+
for _ in 0..4 {
860+
stats.record_alignment(0, max_multimaps);
861+
stats.record_unmapped_reason(UnmappedReason::TooShort);
862+
}
863+
for _ in 0..2 {
864+
stats.record_alignment(15, max_multimaps);
865+
stats.record_unmapped_reason(UnmappedReason::TooManyLoci);
866+
}
867+
868+
let total = stats.total_reads.load(Ordering::Relaxed);
869+
let unique = stats.uniquely_mapped.load(Ordering::Relaxed);
870+
let multi = stats.multi_mapped.load(Ordering::Relaxed);
871+
let unmapped_short = stats.unmapped_short.load(Ordering::Relaxed);
872+
let unmapped_other = stats.unmapped_other.load(Ordering::Relaxed);
873+
let unmapped_mismatches = stats.unmapped_mismatches.load(Ordering::Relaxed);
874+
let too_many_loci = stats.too_many_loci.load(Ordering::Relaxed);
875+
876+
assert_eq!(total, 11);
877+
assert_eq!(too_many_loci, 2);
878+
let bucket_sum =
879+
unique + multi + unmapped_short + unmapped_other + unmapped_mismatches + too_many_loci;
880+
assert_eq!(bucket_sum, total);
881+
}
882+
830883
#[test]
831884
fn test_splice_motif_aggregation() {
832885
use crate::align::score::SpliceMotif;

0 commit comments

Comments
 (0)