diff --git a/Cargo.lock b/Cargo.lock index da20215e..f84de282 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -385,9 +385,9 @@ dependencies = [ [[package]] name = "crossbeam-epoch" -version = "0.9.18" +version = "0.9.20" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5b82ac4a3c2ca9c3460964f020e1402edd5753411d7737aa39c3714ad1b5420e" +checksum = "2d6914041f254d6e9176c01941b21115dcfb7089e55135a35411081bd106ef3f" dependencies = [ "crossbeam-utils", ] @@ -1037,9 +1037,9 @@ dependencies = [ [[package]] name = "libz-sys" -version = "1.1.25" +version = "1.1.31" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d52f4c29e2a68ac30c9087e1b772dc9f44a2b66ed44edf2266cf2be9b03dafc1" +checksum = "ff20308ecabd635a2dc1dd5834a86c462cca20efba1d5d3ddab1f8aea3455cb2" dependencies = [ "cc", "cmake", diff --git a/examples/flagstat_counter_probe.rs b/examples/flagstat_counter_probe.rs new file mode 100644 index 00000000..2c485849 --- /dev/null +++ b/examples/flagstat_counter_probe.rs @@ -0,0 +1,40 @@ +//! Emit raw production flagstat counters for independent differential validation. +//! Each row calls the library implementation on a fresh counter. + +use rustqc::rna::rseqc::flagstat; + +fn main() -> std::io::Result<()> { + use std::io::Write; + + let mut out = std::io::BufWriter::new(std::io::stdout().lock()); + writeln!(out, "flags\ttid\tmtid\tmapq\ttotal\tprimary\tsecondary\tsupplementary\tduplicates\tprimary_duplicates\tmapped\tprimary_mapped\tpaired\tread1\tread2\tproperly_paired\tboth_mapped\tsingletons\tmate_diff_chr\tmate_diff_chr_mapq5")?; + for flags in 0..4096u16 { + for (tid, mtid, mapq) in [(0, 0, 60), (0, 1, 60), (1, 0, 4)] { + let mut c = flagstat::FlagstatCounts::default(); + c.add(flags, tid, mtid, mapq); + write!(out, "{flags}\t{tid}\t{mtid}\t{mapq}")?; + for value in [ + c.total, + c.primary, + c.secondary, + c.supplementary, + c.duplicates, + c.primary_duplicates, + c.mapped, + c.primary_mapped, + c.paired, + c.read1, + c.read2, + c.properly_paired, + c.both_mapped, + c.singletons, + c.mate_diff_chr, + c.mate_diff_chr_mapq5, + ] { + write!(out, "\t{},{}", value[0], value[1])?; + } + writeln!(out)?; + } + } + out.flush() +} diff --git a/src/rna/rseqc/accumulators.rs b/src/rna/rseqc/accumulators.rs index b91a409e..fdbfbe85 100644 --- a/src/rna/rseqc/accumulators.rs +++ b/src/rna/rseqc/accumulators.rs @@ -12,6 +12,7 @@ use indexmap::IndexMap; use rust_htslib::bam; use super::bam_stat::{BamStatResult, GcDepthBin}; +use super::flagstat::FlagstatCounts; /// Default GC-depth bin size in base pairs (matches upstream samtools default). const GCD_BIN_SIZE: u64 = 20_000; @@ -157,30 +158,20 @@ pub struct BamStatAccum { // --- samtools flagstat additional fields --- /// Secondary alignments (0x100) — counted independently of QC/dup. pub secondary: u64, - /// Supplementary alignments (0x800) — counted independently of QC/dup. + /// Supplementary alignments (0x800), excluding secondary records, regardless of QC/dup. pub supplementary: u64, - /// All mapped records (not 0x4), regardless of QC/dup. - pub mapped: u64, /// Paired reads (0x1), regardless of QC/dup. pub paired_flagstat: u64, - /// Read1 in pair (0x40), regardless of QC/dup — for flagstat. - pub read1_flagstat: u64, - /// Read2 in pair (0x80), regardless of QC/dup — for flagstat. - pub read2_flagstat: u64, /// First fragments for samtools stats: primary reads that are not "last fragments". pub first_fragments: u64, /// Last fragments for samtools stats: primary reads with 0x80 flag. pub last_fragments: u64, /// Properly paired reads (0x1 + 0x2), regardless of QC/dup. pub properly_paired: u64, - /// Both mates mapped (paired + both !unmapped). - pub both_mapped: u64, - /// Singletons (paired, this mapped, mate unmapped). - pub singletons: u64, - /// Paired, both mapped, different reference. - pub mate_diff_chr: u64, /// Paired, both mapped, different reference, MAPQ >= 5. pub mate_diff_chr_mapq5: u64, + /// The `samtools flagstat` counters, split by QC-pass/fail. + pub flagstat: FlagstatCounts, // --- samtools idxstats additional fields --- /// Per-reference (tid) mapped and unmapped counts. @@ -319,17 +310,12 @@ impl Default for BamStatAccum { proper_pair_diff_chrom: 0, secondary: 0, supplementary: 0, - mapped: 0, paired_flagstat: 0, - read1_flagstat: 0, - read2_flagstat: 0, first_fragments: 0, last_fragments: 0, properly_paired: 0, - both_mapped: 0, - singletons: 0, - mate_diff_chr: 0, mate_diff_chr_mapq5: 0, + flagstat: FlagstatCounts::default(), chrom_counts: HashMap::new(), unplaced_unmapped: 0, total_len: 0, @@ -409,15 +395,13 @@ impl BamStatAccum { // ================================================================= // samtools flagstat counters (count ALL records, no early returns) // ================================================================= + self.flagstat.add(flags, tid, record.mtid(), mapq); + // samtools gives SECONDARY priority for records carrying both flags. if is_secondary { self.secondary += 1; - } - if is_supplementary { + } else if is_supplementary { self.supplementary += 1; } - if is_mapped { - self.mapped += 1; - } // samtools stats: "1st fragments" / "last fragments" count primary reads only // For paired reads: read2 flag -> last, everything else -> 1st // For SE reads (no PAIRED flag): all counted as 1st fragments @@ -432,27 +416,12 @@ impl BamStatAccum { // (secondary/supplementary are excluded from paired/read1/read2/properly-paired counts) if is_paired && is_primary { self.paired_flagstat += 1; - if flags & BAM_FREAD1 != 0 { - self.read1_flagstat += 1; - } - if flags & BAM_FREAD2 != 0 { - self.read2_flagstat += 1; - } if flags & BAM_FPROPER_PAIR != 0 { self.properly_paired += 1; } let mate_unmapped = flags & BAM_FMUNMAP != 0; - if is_mapped && !mate_unmapped { - self.both_mapped += 1; - if tid != record.mtid() { - self.mate_diff_chr += 1; - if mapq >= 5 { - self.mate_diff_chr_mapq5 += 1; - } - } - } - if is_mapped && mate_unmapped { - self.singletons += 1; + if is_mapped && !mate_unmapped && tid != record.mtid() && mapq >= 5 { + self.mate_diff_chr_mapq5 += 1; } } @@ -1193,17 +1162,12 @@ impl BamStatAccum { // samtools flagstat fields self.secondary += other.secondary; self.supplementary += other.supplementary; - self.mapped += other.mapped; self.paired_flagstat += other.paired_flagstat; - self.read1_flagstat += other.read1_flagstat; - self.read2_flagstat += other.read2_flagstat; self.first_fragments += other.first_fragments; self.last_fragments += other.last_fragments; self.properly_paired += other.properly_paired; - self.both_mapped += other.both_mapped; - self.singletons += other.singletons; - self.mate_diff_chr += other.mate_diff_chr; self.mate_diff_chr_mapq5 += other.mate_diff_chr_mapq5; + self.flagstat.merge(&other.flagstat); // samtools idxstats fields for (tid, (m, u)) in other.chrom_counts { @@ -2336,9 +2300,9 @@ fn merge_vec_arrays(target: &mut Vec<[u64; N]>, source: Vec<[u64 if source.len() > target.len() { target.resize(source.len(), [0u64; N]); } - for (i, arr) in source.into_iter().enumerate() { - for j in 0..N { - target[i][j] += arr[j]; + for (target_row, source_row) in target.iter_mut().zip(source) { + for (target_count, source_count) in target_row.iter_mut().zip(source_row) { + *target_count += source_count; } } } @@ -2372,17 +2336,12 @@ impl BamStatAccum { // samtools flagstat fields secondary: self.secondary, supplementary: self.supplementary, - mapped: self.mapped, paired_flagstat: self.paired_flagstat, - read1_flagstat: self.read1_flagstat, - read2_flagstat: self.read2_flagstat, first_fragments: self.first_fragments, last_fragments: self.last_fragments, properly_paired: self.properly_paired, - both_mapped: self.both_mapped, - singletons: self.singletons, - mate_diff_chr: self.mate_diff_chr, mate_diff_chr_mapq5: self.mate_diff_chr_mapq5, + flagstat: self.flagstat, // samtools idxstats fields chrom_counts: self.chrom_counts, unplaced_unmapped: self.unplaced_unmapped, diff --git a/src/rna/rseqc/bam_stat.rs b/src/rna/rseqc/bam_stat.rs index 2cd06acb..255e6a24 100644 --- a/src/rna/rseqc/bam_stat.rs +++ b/src/rna/rseqc/bam_stat.rs @@ -10,6 +10,8 @@ use std::collections::HashMap; use std::io::Write; use std::path::Path; +use super::flagstat::FlagstatCounts; + // ============================================================================ // Result struct // ============================================================================ @@ -57,30 +59,20 @@ pub struct BamStatResult { // --- samtools flagstat fields --- /// Secondary alignments (0x100). pub secondary: u64, - /// Supplementary alignments (0x800). + /// Supplementary alignments (0x800) that are not also secondary (0x100). pub supplementary: u64, - /// All mapped records (not 0x4). - pub mapped: u64, /// Paired reads (0x1). pub paired_flagstat: u64, - /// Read1 in pair (0x40). - pub read1_flagstat: u64, - /// Read2 in pair (0x80). - pub read2_flagstat: u64, /// First fragments for samtools stats (read1 or unpaired). pub first_fragments: u64, /// Last fragments for samtools stats (read2). pub last_fragments: u64, /// Properly paired reads (0x1 + 0x2). pub properly_paired: u64, - /// Both mates mapped. - pub both_mapped: u64, - /// Singletons (this mapped, mate unmapped). - pub singletons: u64, - /// Paired, both mapped, different reference. - pub mate_diff_chr: u64, /// Paired, both mapped, different reference, MAPQ >= 5. pub mate_diff_chr_mapq5: u64, + /// The `samtools flagstat` counters, split by QC-pass/fail. + pub flagstat: FlagstatCounts, // --- samtools idxstats fields --- /// Per-reference (tid) mapped and unmapped counts. @@ -209,17 +201,12 @@ impl Default for BamStatResult { proper_pair_diff_chrom: 0, secondary: 0, supplementary: 0, - mapped: 0, paired_flagstat: 0, - read1_flagstat: 0, - read2_flagstat: 0, first_fragments: 0, last_fragments: 0, properly_paired: 0, - both_mapped: 0, - singletons: 0, - mate_diff_chr: 0, mate_diff_chr_mapq5: 0, + flagstat: FlagstatCounts::default(), chrom_counts: HashMap::new(), unplaced_unmapped: 0, total_len: 0, diff --git a/src/rna/rseqc/flagstat.rs b/src/rna/rseqc/flagstat.rs index cb1e370d..18d850b9 100644 --- a/src/rna/rseqc/flagstat.rs +++ b/src/rna/rseqc/flagstat.rs @@ -9,33 +9,167 @@ use anyhow::{Context, Result}; use log::debug; use super::bam_stat::BamStatResult; +use crate::rna::bam_flags::*; // ============================================================================ -// Output formatting +// Counting // ============================================================================ -/// Format a flagstat count line with the standard ` + 0 ` format. -fn fmt_line(count: u64, description: &str) -> String { - format!("{count} + 0 {description}") +/// samtools flagstat counters. Each is `[QC-passed, QC-failed]`, split on +/// 0x200 like samtools' `bam_flagstat_t`. +#[derive(Debug, Clone, Default)] +pub struct FlagstatCounts { + /// All records. + pub total: [u64; 2], + /// Neither secondary nor supplementary. + pub primary: [u64; 2], + /// Secondary (0x100). + pub secondary: [u64; 2], + /// Supplementary (0x800), excluding secondary records. + pub supplementary: [u64; 2], + /// Duplicates (0x400). + pub duplicates: [u64; 2], + /// Primary duplicates. + pub primary_duplicates: [u64; 2], + /// Mapped (not 0x4). + pub mapped: [u64; 2], + /// Primary mapped. + pub primary_mapped: [u64; 2], + /// Primary and paired (0x1). + pub paired: [u64; 2], + /// Primary, paired, read1 (0x40). + pub read1: [u64; 2], + /// Primary, paired, read2 (0x80). + pub read2: [u64; 2], + /// Primary, paired, proper pair (0x2) and mapped. + pub properly_paired: [u64; 2], + /// Primary, paired, this and mate mapped. + pub both_mapped: [u64; 2], + /// Primary, paired, mapped with mate unmapped. + pub singletons: [u64; 2], + /// Both mapped, mate on a different reference. + pub mate_diff_chr: [u64; 2], + /// As `mate_diff_chr`, with MAPQ >= 5. + pub mate_diff_chr_mapq5: [u64; 2], } -/// Format a flagstat percentage as `(XX.XX% : N/A)`. -/// -/// samtools uses `:` separated `passed : failed` percentages. Since we do not -/// separate QC-pass/fail, the failed side is always `N/A`. -fn fmt_pct(numerator: u64, denominator: u64) -> String { - if denominator == 0 { - "(N/A : N/A)".to_string() - } else { - let pct = 100.0 * numerator as f64 / denominator as f64; - format!("({pct:.2}% : N/A)") +impl FlagstatCounts { + /// Count one record, following samtools' `flagstat_loop` (bam_stat.c). + pub fn add(&mut self, flags: u16, tid: i32, mtid: i32, mapq: u8) { + let w = usize::from(flags & BAM_FQCFAIL != 0); + let is_secondary = flags & BAM_FSECONDARY != 0; + let is_supplementary = flags & BAM_FSUPPLEMENTARY != 0; + let is_mapped = flags & BAM_FUNMAP == 0; + let is_dup = flags & BAM_FDUP != 0; + + self.total[w] += 1; + if is_secondary { + self.secondary[w] += 1; + } else if is_supplementary { + self.supplementary[w] += 1; + } + if is_mapped { + self.mapped[w] += 1; + } + if is_dup { + self.duplicates[w] += 1; + } + if is_secondary || is_supplementary { + return; + } + + self.primary[w] += 1; + if is_mapped { + self.primary_mapped[w] += 1; + } + if is_dup { + self.primary_duplicates[w] += 1; + } + if flags & BAM_FPAIRED != 0 { + self.paired[w] += 1; + if flags & BAM_FREAD1 != 0 { + self.read1[w] += 1; + } + if flags & BAM_FREAD2 != 0 { + self.read2[w] += 1; + } + if is_mapped && flags & BAM_FPROPER_PAIR != 0 { + self.properly_paired[w] += 1; + } + let mate_mapped = flags & BAM_FMUNMAP == 0; + if is_mapped && mate_mapped { + self.both_mapped[w] += 1; + if tid != mtid { + self.mate_diff_chr[w] += 1; + if mapq >= 5 { + self.mate_diff_chr_mapq5[w] += 1; + } + } + } + if is_mapped && !mate_mapped { + self.singletons[w] += 1; + } + } + } + + /// Add another set of counters into this one. + pub fn merge(&mut self, other: &FlagstatCounts) { + let pairs = [ + (&mut self.total, other.total), + (&mut self.primary, other.primary), + (&mut self.secondary, other.secondary), + (&mut self.supplementary, other.supplementary), + (&mut self.duplicates, other.duplicates), + (&mut self.primary_duplicates, other.primary_duplicates), + (&mut self.mapped, other.mapped), + (&mut self.primary_mapped, other.primary_mapped), + (&mut self.paired, other.paired), + (&mut self.read1, other.read1), + (&mut self.read2, other.read2), + (&mut self.properly_paired, other.properly_paired), + (&mut self.both_mapped, other.both_mapped), + (&mut self.singletons, other.singletons), + (&mut self.mate_diff_chr, other.mate_diff_chr), + (&mut self.mate_diff_chr_mapq5, other.mate_diff_chr_mapq5), + ]; + for (mine, theirs) in pairs { + mine[0] += theirs[0]; + mine[1] += theirs[1]; + } } } +// ============================================================================ +// Output formatting +// ============================================================================ + +/// Format a flagstat count line as ` + `. +fn fmt_line(count: [u64; 2], description: &str) -> String { + format!("{} + {} {description}", count[0], count[1]) +} + +/// Format a flagstat percentage as `(XX.XX% : YY.YY%)`, one per QC column, +/// `N/A` where the denominator is zero. +fn fmt_pct(numerator: [u64; 2], denominator: [u64; 2]) -> String { + let pct = |n: u64, d: u64| { + if d == 0 { + "N/A".to_string() + } else { + // samtools bam_stat.c divides in float before multiplying by 100.0. + format!("{:.2}%", (n as f32 / d as f32) as f64 * 100.0) + } + }; + format!( + "({} : {})", + pct(numerator[0], denominator[0]), + pct(numerator[1], denominator[1]) + ) +} + /// Write samtools flagstat-compatible output. /// /// The output format matches `samtools flagstat` exactly, with one line per -/// statistic. Each line has the form ` + 0 `. +/// statistic. Each line has the form ` + `. /// /// # Arguments /// * `result` - The computed BAM statistics @@ -46,95 +180,87 @@ pub fn write_flagstat(result: &BamStatResult, output_path: &Path) -> Result<()> let mut out = std::fs::File::create(output_path) .with_context(|| format!("Failed to create flagstat file: {}", output_path.display()))?; - // Total = primary + secondary + supplementary (all QC states) - let total = result.total_records; - - // Primary count for percentage denominators - let primary = result.primary_count; + let fs = &result.flagstat; // Line 1: total (QC-passed + QC-failed) writeln!( out, "{}", - fmt_line(total, "in total (QC-passed reads + QC-failed reads)") + fmt_line(fs.total, "in total (QC-passed reads + QC-failed reads)") )?; // Line 2: primary - writeln!(out, "{}", fmt_line(primary, "primary"))?; + writeln!(out, "{}", fmt_line(fs.primary, "primary"))?; // Line 3: secondary - writeln!(out, "{}", fmt_line(result.secondary, "secondary"))?; + writeln!(out, "{}", fmt_line(fs.secondary, "secondary"))?; // Line 4: supplementary - writeln!(out, "{}", fmt_line(result.supplementary, "supplementary"))?; + writeln!(out, "{}", fmt_line(fs.supplementary, "supplementary"))?; // Line 5: duplicates - writeln!(out, "{}", fmt_line(result.duplicates, "duplicates"))?; + writeln!(out, "{}", fmt_line(fs.duplicates, "duplicates"))?; // Line 6: primary duplicates writeln!( out, "{}", - fmt_line(result.primary_duplicates, "primary duplicates") + fmt_line(fs.primary_duplicates, "primary duplicates") )?; // Line 7: mapped writeln!( out, "{} {}", - fmt_line(result.mapped, "mapped"), - fmt_pct(result.mapped, total) + fmt_line(fs.mapped, "mapped"), + fmt_pct(fs.mapped, fs.total) )?; // Line 8: primary mapped writeln!( out, "{} {}", - fmt_line(result.primary_mapped, "primary mapped"), - fmt_pct(result.primary_mapped, primary) + fmt_line(fs.primary_mapped, "primary mapped"), + fmt_pct(fs.primary_mapped, fs.primary) )?; // Line 9: paired in sequencing - writeln!( - out, - "{}", - fmt_line(result.paired_flagstat, "paired in sequencing") - )?; + writeln!(out, "{}", fmt_line(fs.paired, "paired in sequencing"))?; // Line 10: read1 - writeln!(out, "{}", fmt_line(result.read1_flagstat, "read1"))?; + writeln!(out, "{}", fmt_line(fs.read1, "read1"))?; // Line 11: read2 - writeln!(out, "{}", fmt_line(result.read2_flagstat, "read2"))?; + writeln!(out, "{}", fmt_line(fs.read2, "read2"))?; // Line 12: properly paired (percentage of paired reads) writeln!( out, "{} {}", - fmt_line(result.properly_paired, "properly paired"), - fmt_pct(result.properly_paired, result.paired_flagstat) + fmt_line(fs.properly_paired, "properly paired"), + fmt_pct(fs.properly_paired, fs.paired) )?; // Line 13: with itself and mate mapped writeln!( out, "{}", - fmt_line(result.both_mapped, "with itself and mate mapped") + fmt_line(fs.both_mapped, "with itself and mate mapped") )?; // Line 14: singletons (percentage of paired reads) writeln!( out, "{} {}", - fmt_line(result.singletons, "singletons"), - fmt_pct(result.singletons, result.paired_flagstat) + fmt_line(fs.singletons, "singletons"), + fmt_pct(fs.singletons, fs.paired) )?; // Line 15: with mate mapped to a different chr writeln!( out, "{}", - fmt_line(result.mate_diff_chr, "with mate mapped to a different chr") + fmt_line(fs.mate_diff_chr, "with mate mapped to a different chr") )?; // Line 16: with mate mapped to a different chr (mapQ>=5) @@ -142,7 +268,7 @@ pub fn write_flagstat(result: &BamStatResult, output_path: &Path) -> Result<()> out, "{}", fmt_line( - result.mate_diff_chr_mapq5, + fs.mate_diff_chr_mapq5, "with mate mapped to a different chr (mapQ>=5)" ) )?; @@ -195,4 +321,310 @@ mod tests { "Missing mapped line in flagstat output" ); } + + // ------------------------------------------------------------------------ + // Oracle: samtools 1.22.1 `flagstat_loop` and `out_fmt_default` + // (bam_stat.c:48-79 and 124-143), transcribed. + // ------------------------------------------------------------------------ + + /// samtools' `bam_flagstat_t`; index 0 is QC-passed, 1 is QC-failed. + #[derive(Debug, Default)] + struct SamtoolsFlagstat { + n_reads: [u64; 2], + n_mapped: [u64; 2], + n_pair_all: [u64; 2], + n_pair_map: [u64; 2], + n_pair_good: [u64; 2], + n_sgltn: [u64; 2], + n_read1: [u64; 2], + n_read2: [u64; 2], + n_dup: [u64; 2], + n_diffchr: [u64; 2], + n_diffhigh: [u64; 2], + n_secondary: [u64; 2], + n_supp: [u64; 2], + n_primary: [u64; 2], + n_pmapped: [u64; 2], + n_pdup: [u64; 2], + } + + impl SamtoolsFlagstat { + fn add(&mut self, flag: u16, tid: i32, mtid: i32, qual: u8) { + let w = usize::from(flag & BAM_FQCFAIL != 0); + self.n_reads[w] += 1; + if flag & BAM_FSECONDARY != 0 { + self.n_secondary[w] += 1; + } else if flag & BAM_FSUPPLEMENTARY != 0 { + self.n_supp[w] += 1; + } else { + self.n_primary[w] += 1; + if flag & BAM_FPAIRED != 0 { + self.n_pair_all[w] += 1; + if flag & BAM_FPROPER_PAIR != 0 && flag & BAM_FUNMAP == 0 { + self.n_pair_good[w] += 1; + } + if flag & BAM_FREAD1 != 0 { + self.n_read1[w] += 1; + } + if flag & BAM_FREAD2 != 0 { + self.n_read2[w] += 1; + } + if flag & BAM_FMUNMAP != 0 && flag & BAM_FUNMAP == 0 { + self.n_sgltn[w] += 1; + } + if flag & BAM_FUNMAP == 0 && flag & BAM_FMUNMAP == 0 { + self.n_pair_map[w] += 1; + if mtid != tid { + self.n_diffchr[w] += 1; + if qual >= 5 { + self.n_diffhigh[w] += 1; + } + } + } + } + if flag & BAM_FUNMAP == 0 { + self.n_pmapped[w] += 1; + } + if flag & BAM_FDUP != 0 { + self.n_pdup[w] += 1; + } + } + if flag & BAM_FUNMAP == 0 { + self.n_mapped[w] += 1; + } + if flag & BAM_FDUP != 0 { + self.n_dup[w] += 1; + } + } + + fn render(&self) -> String { + // samtools computes `(float)n / total * 100.0`. + fn percent(n: u64, total: u64) -> String { + if total == 0 { + "N/A".to_string() + } else { + format!("{:.2}%", (n as f32 / total as f32) as f64 * 100.0) + } + } + let line = |c: [u64; 2], d: &str| format!("{} + {} {d}\n", c[0], c[1]); + let pct = |c: [u64; 2], t: [u64; 2], d: &str| { + format!( + "{} + {} {d} ({} : {})\n", + c[0], + c[1], + percent(c[0], t[0]), + percent(c[1], t[1]) + ) + }; + [ + line(self.n_reads, "in total (QC-passed reads + QC-failed reads)"), + line(self.n_primary, "primary"), + line(self.n_secondary, "secondary"), + line(self.n_supp, "supplementary"), + line(self.n_dup, "duplicates"), + line(self.n_pdup, "primary duplicates"), + pct(self.n_mapped, self.n_reads, "mapped"), + pct(self.n_pmapped, self.n_primary, "primary mapped"), + line(self.n_pair_all, "paired in sequencing"), + line(self.n_read1, "read1"), + line(self.n_read2, "read2"), + pct(self.n_pair_good, self.n_pair_all, "properly paired"), + line(self.n_pair_map, "with itself and mate mapped"), + pct(self.n_sgltn, self.n_pair_all, "singletons"), + line(self.n_diffchr, "with mate mapped to a different chr"), + line( + self.n_diffhigh, + "with mate mapped to a different chr (mapQ>=5)", + ), + ] + .concat() + } + } + + /// (tid, mtid, mapq) placements: mate on the same contig, mate on another + /// contig, mate on another contig with MAPQ below 5. + const PLACEMENTS: [(i32, i32, u8); 3] = [(0, 0, 60), (0, 1, 60), (1, 0, 3)]; + + fn record(flag: u16, (tid, mtid, mapq): (i32, i32, u8)) -> bam::Record { + use rust_htslib::bam::record::{Cigar, CigarString}; + let mut rec = bam::Record::new(); + rec.set( + b"r", + Some(&CigarString(vec![Cigar::Match(4)])), + b"ACGT", + &[30; 4], + ); + rec.set_flags(flag); + rec.set_tid(tid); + rec.set_pos(9); + rec.set_mtid(mtid); + rec.set_mpos(9); + rec.set_mapq(mapq); + rec + } + + fn render_flagstat(accum: BamStatAccum, name: &str) -> String { + let path = std::env::temp_dir().join(format!( + "rustqc_test_flagstat_{name}_{}.txt", + std::process::id() + )); + write_flagstat(&accum.into_result(), &path).expect("Failed to write flagstat"); + let contents = std::fs::read_to_string(&path).unwrap(); + let _ = std::fs::remove_file(&path); + contents + } + + /// Each of the 4,096 flags alone, in every placement, against the oracle. + #[test] + fn test_flagstat_matches_samtools_flagstat_loop_per_flag() { + let mut mismatches = Vec::new(); + for placement in PLACEMENTS { + for flag in 0..4096u16 { + let mut accum = BamStatAccum::default(); + accum.process_read(&record(flag, placement), 30); + let mut oracle = SamtoolsFlagstat::default(); + oracle.add(flag, placement.0, placement.1, placement.2); + + let got = render_flagstat(accum, "per_flag"); + let want = oracle.render(); + if got != want { + let diff: Vec = want + .lines() + .zip(got.lines()) + .filter(|(w, g)| w != g) + .map(|(w, g)| format!("want `{w}`, got `{g}`")) + .collect(); + mismatches.push(format!("{flag:#06x} {placement:?}: {}", diff.join("; "))); + } + } + } + assert!( + mismatches.is_empty(), + "{} of {} (flag, placement) cases differ from samtools' flagstat_loop; first:\n{}", + mismatches.len(), + 4096 * PLACEMENTS.len(), + mismatches[..mismatches.len().min(8)].join("\n") + ); + } + + /// All flags and placements in one file, split across two merged + /// accumulators as the per-contig parallel path does. + #[test] + fn test_flagstat_matches_samtools_flagstat_loop_merged() { + let mut parts = [BamStatAccum::default(), BamStatAccum::default()]; + let mut oracle = SamtoolsFlagstat::default(); + for placement in PLACEMENTS { + for flag in 0..4096u16 { + parts[placement.0 as usize].process_read(&record(flag, placement), 30); + oracle.add(flag, placement.0, placement.1, placement.2); + } + } + let [mut accum, other] = parts; + accum.merge(other); + assert_eq!(render_flagstat(accum, "merged"), oracle.render()); + } + + #[test] + fn test_flagstat_dual_flags_with_qcfail_and_duplicate() { + let mut accum = BamStatAccum::default(); + let mut direct = FlagstatCounts::default(); + let flags = [ + 0, + BAM_FDUP, + BAM_FSECONDARY, + BAM_FSUPPLEMENTARY, + BAM_FSECONDARY | BAM_FSUPPLEMENTARY, + BAM_FSECONDARY | BAM_FSUPPLEMENTARY | BAM_FDUP, + ]; + for qc in [0, BAM_FQCFAIL] { + for flag in flags { + let flag = flag | qc; + accum.process_read(&record(flag, PLACEMENTS[0]), 30); + direct.add(flag, 0, 0, 60); + } + } + for fs in [&accum.flagstat, &direct] { + assert_eq!(fs.total, [6, 6]); + assert_eq!(fs.primary, [2, 2]); + assert_eq!(fs.secondary, [3, 3]); + assert_eq!(fs.supplementary, [1, 1]); + assert_eq!(fs.duplicates, [2, 2]); + assert_eq!(fs.primary_duplicates, [1, 1]); + assert_eq!(fs.mapped, [6, 6]); + assert_eq!(fs.primary_mapped, [2, 2]); + for w in 0..2 { + assert_eq!( + fs.primary[w] + fs.secondary[w] + fs.supplementary[w], + fs.total[w] + ); + } + } + let result = accum.into_result(); + assert_eq!(result.total_records, 12); + assert_eq!(result.primary_count, 4); + assert_eq!(result.secondary, 6); + assert_eq!(result.supplementary, 2); + assert_eq!( + result.primary_count + result.secondary + result.supplementary, + result.total_records + ); + } + + #[test] + fn test_flagstat_partitioned_counts_match_serial() { + for partitions in [1, 2, 7, 31] { + let mut parts = vec![FlagstatCounts::default(); partitions]; + let mut oracle = SamtoolsFlagstat::default(); + for (i, (flag, placement)) in (0..4096u16) + .flat_map(|flag| PLACEMENTS.map(|placement| (flag, placement))) + .enumerate() + { + parts[i % partitions].add(flag, placement.0, placement.1, placement.2); + oracle.add(flag, placement.0, placement.1, placement.2); + } + let mut merged = FlagstatCounts::default(); + for part in parts.into_iter().rev() { + merged.merge(&part); + } + for w in 0..2 { + assert_eq!( + merged.primary[w] + merged.secondary[w] + merged.supplementary[w], + merged.total[w] + ); + } + let mut accum = BamStatAccum::default(); + accum.flagstat.merge(&merged); + assert_eq!( + render_flagstat(accum, "partitioned"), + oracle.render(), + "{partitions} partitions" + ); + } + } + + #[test] + fn test_flagstat_large_counter_percentages_match_samtools_c() { + // samtools 1.22.1 and 1.24 bam_stat.c percent(): + // sprintf(buffer, "%.2f%%", (float)n / total * 100.0). + // Expected strings were checked by compiling that C function unchanged. + let cases = [ + (0, 0, "N/A"), + (0, 1_000_000_000, "0.00%"), + (1, 32, "3.12%"), + (1, 6, "16.67%"), + (16_777_217, 16_777_218, "100.00%"), + (500_049_999, 1_000_000_000, "50.01%"), + (999_950_000, 1_000_000_000, "99.99%"), + (499_050_000, 1_000_000_000, "49.90%"), + (499_150_000, 1_000_000_000, "49.92%"), + (i64::MAX as u64 - 1, i64::MAX as u64, "100.00%"), + ]; + for (n, d, expected) in cases { + assert_eq!( + fmt_pct([n, n], [d, d]), + format!("({expected} : {expected})"), + "{n}/{d}" + ); + } + } } diff --git a/src/rna/rseqc/stats.rs b/src/rna/rseqc/stats.rs index 20abf8c0..5d81f9e5 100644 --- a/src/rna/rseqc/stats.rs +++ b/src/rna/rseqc/stats.rs @@ -219,11 +219,8 @@ pub fn write_stats(result: &BamStatResult, output_path: &Path) -> Result<()> { // Quality and length stats sn_no_comment(&mut out, "reads QC failed:", filtered)?; - sn_no_comment( - &mut out, - "non-primary alignments:", - result.secondary + result.supplementary, - )?; + // samtools stats reports nreads_secondary separately from supplementary. + sn_no_comment(&mut out, "non-primary alignments:", result.secondary)?; sn_no_comment(&mut out, "supplementary alignments:", result.supplementary)?; sn( &mut out, @@ -525,13 +522,13 @@ fn write_gcc( for cycle in 0..max_cycles { let mut combined = [0u64; 6]; if cycle < fbc.len() { - for j in 0..6 { - combined[j] += fbc[cycle][j]; + for (count, &value) in combined.iter_mut().zip(&fbc[cycle]) { + *count += value; } } if cycle < lbc.len() { - for j in 0..6 { - combined[j] += lbc[cycle][j]; + for (count, &value) in combined.iter_mut().zip(&lbc[cycle]) { + *count += value; } } // Denominator is A+C+G+T only (first 4 elements), matching upstream diff --git a/tests/expected/samtools/README.md b/tests/expected/samtools/README.md new file mode 100644 index 00000000..81d0cfb9 --- /dev/null +++ b/tests/expected/samtools/README.md @@ -0,0 +1,10 @@ +# samtools reference outputs + +`all_flags.flagstat` is `samtools flagstat` 1.22.1 output (1.24 prints the same) for the records `all_flags_sam_lines()` in `tests/integration_test.rs` writes. To regenerate: + +```bash +{ printf '@HD\tVN:1.6\tSO:coordinate\n@SQ\tSN:chr1\tLN:1000\n@SQ\tSN:chr2\tLN:1000\n' + awk 'BEGIN{for(f=0;f<4096;f++){printf "f%da\t%d\tchr1\t10\t60\t4M\t=\t10\t0\tACGT\tIIII\n",f,f; printf "f%db\t%d\tchr1\t10\t60\t4M\tchr2\t10\t0\tACGT\tIIII\n",f,f}; for(f=0;f<4096;f++) printf "f%dc\t%d\tchr2\t10\t3\t4M\tchr1\t10\t0\tACGT\tIIII\n",f,f}' +} > all_flags.sam +samtools flagstat all_flags.sam > tests/expected/samtools/all_flags.flagstat +``` diff --git a/tests/expected/samtools/all_flags.flagstat b/tests/expected/samtools/all_flags.flagstat new file mode 100644 index 00000000..7ddfb9b7 --- /dev/null +++ b/tests/expected/samtools/all_flags.flagstat @@ -0,0 +1,16 @@ +6144 + 6144 in total (QC-passed reads + QC-failed reads) +1536 + 1536 primary +3072 + 3072 secondary +1536 + 1536 supplementary +3072 + 3072 duplicates +768 + 768 primary duplicates +3072 + 3072 mapped (50.00% : 50.00%) +768 + 768 primary mapped (50.00% : 50.00%) +768 + 768 paired in sequencing +384 + 384 read1 +384 + 384 read2 +192 + 192 properly paired (25.00% : 25.00%) +192 + 192 with itself and mate mapped +192 + 192 singletons (25.00% : 25.00%) +128 + 128 with mate mapped to a different chr +64 + 64 with mate mapped to a different chr (mapQ>=5) diff --git a/tests/integration_test.rs b/tests/integration_test.rs index c7ea3acd..212bce7d 100644 --- a/tests/integration_test.rs +++ b/tests/integration_test.rs @@ -12,26 +12,9 @@ use std::time::{SystemTime, UNIX_EPOCH}; use rust_htslib::bam; use rust_htslib::bam::header::HeaderRecord; -/// Helper: get the path to the rustqc binary. -/// -/// Uses the `CARGO_BIN_EXE_rustqc` env var when available (set by cargo test), -/// otherwise falls back to building and using `target/debug/rustqc`. +/// Use the binary Cargo built for this integration target and profile. fn rustqc_binary() -> String { - // cargo test sets this for integration tests automatically - if let Ok(path) = std::env::var("CARGO_BIN_EXE_rustqc") { - return path; - } - - // Fallback: build and use debug binary - let status = Command::new("cargo") - .args(["build"]) - .status() - .expect("Failed to run cargo build"); - assert!(status.success(), "cargo build failed"); - - let path = Path::new("target/debug/rustqc"); - assert!(path.exists(), "Binary not found at {:?}", path); - path.to_str().unwrap().to_string() + env!("CARGO_BIN_EXE_rustqc").to_string() } /// Helper: run rustqc rna on test data and return the output directory @@ -1026,3 +1009,259 @@ fn test_dup_check_parallel_uses_global_duplicate_state() { let _ = fs::remove_dir_all(root); } + +// =================================================================== +// samtools flagstat parity +// =================================================================== + +/// The two plain mapped reads every flagstat case sits between, so that some +/// read lands in a gene and `rustqc rna` produces its outputs. +const FLAGSTAT_BACKGROUND: [&str; 2] = [ + "r0\t0\tchr1\t1001\t60\t50M\t*\t0\t0\tACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTAC\t??????????????????????????????????????????????????", + "r1\t0\tchr1\t1201\t60\t50M\t*\t0\t0\tACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTAC\t??????????????????????????????????????????????????", +]; + +/// Run `rustqc rna` on the background reads plus one record with `flag` (and +/// mate field `rnext`) and return its flagstat and stats outputs. +fn samtools_for_one_record(label: &str, flag: u16, rnext: &str) -> (String, String) { + let root = unique_test_dir(label); + let bam_path = root.join("s.bam"); + let gtf_path = root.join("a.gtf"); + let outdir = root.join("out"); + fs::create_dir_all(&outdir).unwrap(); + fs::write( + >f_path, + "chr1\tt\tgene\t900\t3000\t.\t+\t.\tgene_id \"g1\"; gene_name \"g1\";\n\ + chr1\tt\texon\t900\t3000\t.\t+\t.\tgene_id \"g1\"; transcript_id \"t1\"; exon_number \"1\";\n", + ) + .unwrap(); + let pnext = if rnext == "*" { 0 } else { 1001 }; + let case = format!( + "r2\t{flag}\tchr1\t1001\t60\t50M\t{rnext}\t{pnext}\t0\tACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTAC\t??????????????????????????????????????????????????" + ); + write_bam_fixture( + &bam_path, + &[("chr1", 100_000), ("chr2", 100_000)], + &[FLAGSTAT_BACKGROUND[0], &case, FLAGSTAT_BACKGROUND[1]], + ); + + let output = run_rustqc_custom(&outdir, &bam_path, >f_path, 1, true, &[]); + assert!( + output.status.success(), + "rustqc failed on flag {flag:#06x}:\n{}", + String::from_utf8_lossy(&output.stderr) + ); + let flagstat = fs::read_to_string(outdir.join("samtools/s.flagstat")).unwrap(); + let stats = fs::read_to_string(outdir.join("samtools/s.stats")).unwrap(); + let _ = fs::remove_dir_all(root); + (flagstat, stats) +} + +fn flagstat_for_one_record(label: &str, flag: u16, rnext: &str) -> String { + samtools_for_one_record(label, flag, rnext).0 +} + +/// A QC-failed read (0x200) belongs in samtools' second column. +/// Expected text is `samtools flagstat` 1.22.1 output for this input. +#[test] +fn test_flagstat_qc_failed_column() { + let expected = "\ +2 + 1 in total (QC-passed reads + QC-failed reads) +2 + 1 primary +0 + 0 secondary +0 + 0 supplementary +0 + 0 duplicates +0 + 0 primary duplicates +2 + 1 mapped (100.00% : 100.00%) +2 + 1 primary mapped (100.00% : 100.00%) +0 + 0 paired in sequencing +0 + 0 read1 +0 + 0 read2 +0 + 0 properly paired (N/A : N/A) +0 + 0 with itself and mate mapped +0 + 0 singletons (N/A : N/A) +0 + 0 with mate mapped to a different chr +0 + 0 with mate mapped to a different chr (mapQ>=5) +"; + assert_eq!(flagstat_for_one_record("fs-qcfail", 0x200, "*"), expected); +} + +/// "properly paired" needs the read itself mapped (samtools bam_stat.c:61), +/// so PAIRED|PROPER_PAIR|UNMAP|MUNMAP|READ1 (0x4f) does not count. +/// Expected text is `samtools flagstat` 1.22.1 output for this input. +#[test] +fn test_flagstat_properly_paired_requires_mapped_read() { + let expected = "\ +3 + 0 in total (QC-passed reads + QC-failed reads) +3 + 0 primary +0 + 0 secondary +0 + 0 supplementary +0 + 0 duplicates +0 + 0 primary duplicates +2 + 0 mapped (66.67% : N/A) +2 + 0 primary mapped (66.67% : N/A) +1 + 0 paired in sequencing +1 + 0 read1 +0 + 0 read2 +0 + 0 properly paired (0.00% : N/A) +0 + 0 with itself and mate mapped +0 + 0 singletons (0.00% : N/A) +0 + 0 with mate mapped to a different chr +0 + 0 with mate mapped to a different chr (mapQ>=5) +"; + assert_eq!( + flagstat_for_one_record("fs-proper-unmapped", 0x4f, "="), + expected + ); +} + +/// A QC-failed duplicate (0x600) is a duplicate in the QC-failed column, so +/// "primary duplicates" never exceeds "duplicates". +/// Expected text is `samtools flagstat` 1.22.1 output for this input. +#[test] +fn test_flagstat_qc_failed_duplicate_counts_as_duplicate() { + let expected = "\ +2 + 1 in total (QC-passed reads + QC-failed reads) +2 + 1 primary +0 + 0 secondary +0 + 0 supplementary +0 + 1 duplicates +0 + 1 primary duplicates +2 + 1 mapped (100.00% : 100.00%) +2 + 1 primary mapped (100.00% : 100.00%) +0 + 0 paired in sequencing +0 + 0 read1 +0 + 0 read2 +0 + 0 properly paired (N/A : N/A) +0 + 0 with itself and mate mapped +0 + 0 singletons (N/A : N/A) +0 + 0 with mate mapped to a different chr +0 + 0 with mate mapped to a different chr (mapQ>=5) +"; + assert_eq!( + flagstat_for_one_record("fs-qcfail-dup", 0x600, "*"), + expected + ); +} + +/// SECONDARY takes priority over SUPPLEMENTARY, matching samtools 1.22.1. +#[test] +fn test_flagstat_secondary_supplementary_overlap_counts_secondary_only() { + let expected = "\ +3 + 0 in total (QC-passed reads + QC-failed reads) +2 + 0 primary +1 + 0 secondary +0 + 0 supplementary +0 + 0 duplicates +0 + 0 primary duplicates +3 + 0 mapped (100.00% : N/A) +2 + 0 primary mapped (100.00% : N/A) +0 + 0 paired in sequencing +0 + 0 read1 +0 + 0 read2 +0 + 0 properly paired (N/A : N/A) +0 + 0 with itself and mate mapped +0 + 0 singletons (N/A : N/A) +0 + 0 with mate mapped to a different chr +0 + 0 with mate mapped to a different chr (mapQ>=5) +"; + assert_eq!(flagstat_for_one_record("fs-sec-supp", 0x900, "*"), expected); +} + +/// SAM lines for every 12-bit flag in three placements: mate on the same +/// contig, mate on the other contig (MAPQ 60), and mate on the other contig +/// with MAPQ 3. `tests/expected/samtools/all_flags.flagstat` is samtools +/// 1.22.1's output for exactly these records (see the README there). +fn all_flags_sam_lines() -> Vec { + let mut chr1 = Vec::new(); + let mut chr2 = Vec::new(); + for flag in 0..4096u16 { + chr1.push(format!( + "f{flag}a\t{flag}\tchr1\t10\t60\t4M\t=\t10\t0\tACGT\tIIII" + )); + chr1.push(format!( + "f{flag}b\t{flag}\tchr1\t10\t60\t4M\tchr2\t10\t0\tACGT\tIIII" + )); + chr2.push(format!( + "f{flag}c\t{flag}\tchr2\t10\t3\t4M\tchr1\t10\t0\tACGT\tIIII" + )); + } + chr1.extend(chr2); + chr1 +} + +/// Every flag, against the unmodified samtools 1.22.1 golden output, on +/// both the serial and per-contig parallel paths. +#[test] +fn test_flagstat_matches_samtools_on_every_flag() { + let expected = fs::read_to_string("tests/expected/samtools/all_flags.flagstat").unwrap(); + + let lines = all_flags_sam_lines(); + let lines: Vec<&str> = lines.iter().map(String::as_str).collect(); + for threads in [1, 2] { + let root = unique_test_dir(&format!("fs-all-flags-t{threads}")); + let bam_path = root.join("all_flags.bam"); + let gtf_path = root.join("all_flags.gtf"); + let outdir = root.join("out"); + fs::create_dir_all(&outdir).unwrap(); + write_test_gtf(>f_path, &["chr1", "chr2"]); + write_bam_fixture(&bam_path, &[("chr1", 1000), ("chr2", 1000)], &lines); + + let output = run_rustqc_custom(&outdir, &bam_path, >f_path, threads, true, &[]); + assert!( + output.status.success(), + "rustqc failed with --threads {threads}:\n{}", + String::from_utf8_lossy(&output.stderr) + ); + let flagstat = fs::read_to_string(outdir.join("samtools/all_flags.flagstat")).unwrap(); + assert_eq!(flagstat, expected, "flagstat with --threads {threads}"); + let _ = fs::remove_dir_all(root); + } +} + +/// A dual-flagged QC-failed duplicate is secondary in the failed column, +/// stays a duplicate, and contributes to neither supplementary nor primary. +#[test] +fn test_flagstat_secondary_supplementary_qcfail_duplicate() { + let (flagstat, stats) = samtools_for_one_record("fs-sec-supp-qcfail-dup", 0xf00, "*"); + let expected = "\ +2 + 1 in total (QC-passed reads + QC-failed reads) +2 + 0 primary +0 + 1 secondary +0 + 0 supplementary +0 + 1 duplicates +0 + 0 primary duplicates +2 + 1 mapped (100.00% : 100.00%) +2 + 0 primary mapped (100.00% : N/A) +0 + 0 paired in sequencing +0 + 0 read1 +0 + 0 read2 +0 + 0 properly paired (N/A : N/A) +0 + 0 with itself and mate mapped +0 + 0 singletons (N/A : N/A) +0 + 0 with mate mapped to a different chr +0 + 0 with mate mapped to a different chr (mapQ>=5) +"; + assert_eq!(flagstat, expected); + assert!( + stats.contains("SN\tnon-primary alignments:\t1\n"), + "{stats}" + ); + assert!( + stats.contains("SN\tsupplementary alignments:\t0\n"), + "{stats}" + ); +} + +#[test] +fn test_stats_supplementary_is_not_non_primary() { + let (_, stats) = samtools_for_one_record("stats-supp-only", 0x800, "*"); + assert!( + stats.contains("SN\tnon-primary alignments:\t0\n"), + "{stats}" + ); + assert!( + stats.contains("SN\tsupplementary alignments:\t1\n"), + "{stats}" + ); +}