From f119afe23b984c3a1d8115072bb58fff4565e700 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 13 Aug 2026 09:21:04 +0200 Subject: [PATCH] fix: count secondary+supplementary records as secondary only MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit samtools gives the SECONDARY (0x100) flag priority over SUPPLEMENTARY (0x800): a record carrying both is counted as secondary and never reaches the supplementary counter (bam_stat.c flagstat_loop and stats.c both return early on secondary). RustQC tested the two bits independently, so dual-flagged records were counted twice — inflating the flagstat supplementary total and breaking the invariant `primary + secondary + supplementary == total`. The same conflation affected the samtools stats SN section: "non-primary alignments" is `nreads_secondary` in samtools, not secondary+supplementary. Verified against samtools 1.24 on a synthetic BAM with 4 primary pairs, 2 secondary-only, 3 supplementary-only and 5 dual-flagged records: 18 total / 8 primary / 7 secondary / 3 supplementary, non-primary 7. Closes #125 Co-Authored-By: Claude Opus 5 (1M context) --- src/rna/rseqc/accumulators.rs | 13 +++++-- src/rna/rseqc/bam_stat.rs | 4 ++- src/rna/rseqc/flagstat.rs | 64 +++++++++++++++++++++++++++++++++++ src/rna/rseqc/stats.rs | 8 ++--- 4 files changed, 80 insertions(+), 9 deletions(-) diff --git a/src/rna/rseqc/accumulators.rs b/src/rna/rseqc/accumulators.rs index b91a409e..49f84b19 100644 --- a/src/rna/rseqc/accumulators.rs +++ b/src/rna/rseqc/accumulators.rs @@ -157,7 +157,9 @@ 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) that are *not* also secondary — + /// counted independently of QC/dup. Matches samtools, which gives the + /// SECONDARY flag priority over SUPPLEMENTARY. pub supplementary: u64, /// All mapped records (not 0x4), regardless of QC/dup. pub mapped: u64, @@ -409,10 +411,15 @@ impl BamStatAccum { // ================================================================= // samtools flagstat counters (count ALL records, no early returns) // ================================================================= + // samtools gives the SECONDARY flag priority: a record carrying both + // 0x100 and 0x800 is counted as secondary only, never as supplementary + // (bam_stat.c flagstat_loop / stats.c both return early on secondary). + // Counting the two bits independently would inflate the supplementary + // total and break the `primary + secondary + supplementary == total` + // invariant. if is_secondary { self.secondary += 1; - } - if is_supplementary { + } else if is_supplementary { self.supplementary += 1; } if is_mapped { diff --git a/src/rna/rseqc/bam_stat.rs b/src/rna/rseqc/bam_stat.rs index 2cd06acb..5e788a00 100644 --- a/src/rna/rseqc/bam_stat.rs +++ b/src/rna/rseqc/bam_stat.rs @@ -57,7 +57,9 @@ pub struct BamStatResult { // --- samtools flagstat fields --- /// Secondary alignments (0x100). pub secondary: u64, - /// Supplementary alignments (0x800). + /// Supplementary alignments (0x800) that are not also secondary (0x100). + /// samtools gives the SECONDARY flag priority, so dual-flagged records + /// count only towards `secondary`. pub supplementary: u64, /// All mapped records (not 0x4). pub mapped: u64, diff --git a/src/rna/rseqc/flagstat.rs b/src/rna/rseqc/flagstat.rs index cb1e370d..2fdc8ad3 100644 --- a/src/rna/rseqc/flagstat.rs +++ b/src/rna/rseqc/flagstat.rs @@ -195,4 +195,68 @@ mod tests { "Missing mapped line in flagstat output" ); } + + /// A record carrying both SECONDARY (0x100) and SUPPLEMENTARY (0x800) must + /// be counted as secondary only, matching `samtools flagstat`. Counting it + /// in both totals inflates the supplementary count and breaks the + /// `primary + secondary + supplementary == total` invariant. + #[test] + fn test_dual_flagged_reads_count_as_secondary_only() { + // 4 primary pairs (8 records), 2 secondary-only, 3 supplementary-only, + // 5 secondary+supplementary. Reference values come from + // `samtools flagstat` 1.24 on the same records: + // 18 in total / 8 primary / 7 secondary / 3 supplementary + let mut sam = String::from("@HD\tVN:1.6\tSO:coordinate\n@SQ\tSN:chr1\tLN:20000\n"); + // Records are emitted in ascending coordinate order (the accumulators + // assume coordinate-sorted input). + let mut pos = 100u64; + let push = |sam: &mut String, name: String, flag: u16, pos: &mut u64| { + sam.push_str(&format!( + "{name}\t{flag}\tchr1\t{pos}\t30\t10M\t*\t0\t0\tACGTACGTAC\tIIIIIIIIII\n" + )); + *pos += 20; + }; + for i in 0..4 { + push(&mut sam, format!("p{i}"), 99, &mut pos); + push(&mut sam, format!("p{i}"), 147, &mut pos); + } + for i in 0..2 { + push(&mut sam, format!("s{i}"), 256, &mut pos); + } + for i in 0..3 { + push(&mut sam, format!("u{i}"), 2048, &mut pos); + } + for i in 0..5 { + push(&mut sam, format!("d{i}"), 2304, &mut pos); + } + + let tmp_path = std::env::temp_dir().join("rustqc_test_flagstat_dual_flags.sam"); + std::fs::write(&tmp_path, sam).expect("Failed to write test SAM"); + + let mut reader = bam::Reader::from_path(&tmp_path).expect("Failed to open test SAM"); + let mut accum = BamStatAccum::default(); + let mut record = bam::Record::new(); + while let Some(res) = reader.read(&mut record) { + res.expect("Error reading SAM record"); + accum.process_read(&record, 30); + } + let _ = std::fs::remove_file(&tmp_path); + + let result = accum.into_result(); + assert_eq!(result.total_records, 18, "total records"); + assert_eq!(result.primary_count, 8, "primary count"); + assert_eq!( + result.secondary, 7, + "secondary must include dual-flagged records" + ); + assert_eq!( + result.supplementary, 3, + "supplementary must exclude dual-flagged records" + ); + assert_eq!( + result.primary_count + result.secondary + result.supplementary, + result.total_records, + "primary + secondary + supplementary must equal total" + ); + } } diff --git a/src/rna/rseqc/stats.rs b/src/rna/rseqc/stats.rs index 20abf8c0..dfa9a215 100644 --- a/src/rna/rseqc/stats.rs +++ b/src/rna/rseqc/stats.rs @@ -219,11 +219,9 @@ 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` here: stats.c returns early on + // the SECONDARY flag, so supplementary records never reach this counter. + sn_no_comment(&mut out, "non-primary alignments:", result.secondary)?; sn_no_comment(&mut out, "supplementary alignments:", result.supplementary)?; sn( &mut out,