From ee13ae72af54ebe18bf4cbced7e06e6dae979533 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 04:22:35 -0700 Subject: [PATCH 1/3] Stabilize summary output and sampling logs --- modkit-core/src/modbam_util/subcommands.rs | 18 +- modkit-core/src/summarize.rs | 5 +- modkit-core/src/writers.rs | 385 +++++++++++++++++---- 3 files changed, 345 insertions(+), 63 deletions(-) diff --git a/modkit-core/src/modbam_util/subcommands.rs b/modkit-core/src/modbam_util/subcommands.rs index 509b270b..7d13ad20 100644 --- a/modkit-core/src/modbam_util/subcommands.rs +++ b/modkit-core/src/modbam_util/subcommands.rs @@ -1771,6 +1771,10 @@ pub struct ModSummarize { interval_size: u32, } +fn sampling_fraction_log_message(fraction: f64) -> String { + format!("sampling {}% of reads", fraction * 100f64) +} + impl ModSummarize { pub fn run(&self) -> anyhow::Result<()> { let _handle = init_logging(self.log_filepath.as_ref()); @@ -1807,7 +1811,9 @@ impl ModSummarize { (None, Some(num_reads)) => { info!("sampling {num_reads} reads from BAM") } - (Some(pct), None) => info!("sampling {pct}% of reads"), + (Some(fraction), None) => { + info!("{}", sampling_fraction_log_message(fraction)) + } (Some(_), Some(_)) => unreachable!(), }); @@ -2473,3 +2479,13 @@ impl ModMode { } } } + +#[cfg(test)] +mod tests { + use super::sampling_fraction_log_message; + + #[test] + fn summary_sampling_fraction_log_uses_percentage_units() { + assert_eq!(sampling_fraction_log_message(0.1), "sampling 10% of reads"); + } +} diff --git a/modkit-core/src/summarize.rs b/modkit-core/src/summarize.rs index 04e90d7b..dd97f5b4 100644 --- a/modkit-core/src/summarize.rs +++ b/modkit-core/src/summarize.rs @@ -48,8 +48,9 @@ pub struct ModSummary<'a> { impl<'a> ModSummary<'a> { pub(crate) fn mod_bases(&self) -> String { - self.mod_call_counts - .keys() + [DnaBase::A, DnaBase::C, DnaBase::G, DnaBase::T] + .into_iter() + .filter(|base| self.mod_call_counts.contains_key(base)) .map(|d| d.char().to_string()) .collect::>() .join(",") diff --git a/modkit-core/src/writers.rs b/modkit-core/src/writers.rs index e1b16e93..e449ea00 100644 --- a/modkit-core/src/writers.rs +++ b/modkit-core/src/writers.rs @@ -567,6 +567,38 @@ pub struct TableWriter { writer: BufWriter, } +/// Stable writer-local ordering: letter codes lexicographically, followed by +/// ChEBI identifiers numerically. +fn compare_summary_mod_codes( + left: &ModCodeRepr, + right: &ModCodeRepr, +) -> Ordering { + match (left, right) { + (ModCodeRepr::Code(left), ModCodeRepr::Code(right)) => left.cmp(right), + (ModCodeRepr::Code(_), ModCodeRepr::ChEbi(_)) => Ordering::Less, + (ModCodeRepr::ChEbi(_), ModCodeRepr::Code(_)) => Ordering::Greater, + (ModCodeRepr::ChEbi(left), ModCodeRepr::ChEbi(right)) => { + left.cmp(right) + } + } +} + +fn compare_summary_base_states( + left: &BaseState, + right: &BaseState, +) -> Ordering { + match (left, right) { + (BaseState::Canonical(left), BaseState::Canonical(right)) => { + left.cmp(right) + } + (BaseState::Canonical(_), BaseState::Modified(_)) => Ordering::Less, + (BaseState::Modified(_), BaseState::Canonical(_)) => Ordering::Greater, + (BaseState::Modified(left), BaseState::Modified(right)) => { + compare_summary_mod_codes(left, right) + } + } +} + impl TableWriter { pub fn new() -> Self { let out = BufWriter::new(std::io::stdout()); @@ -608,7 +640,14 @@ impl<'a, W: Write> OutWriter> for TableWriter { { metadata_table.add_row(row![ format!("modification_codes_for_{base}"), - format!("{}", codes.iter().copied().sorted().join(",")) + format!( + "{}", + codes + .iter() + .copied() + .sorted_by(compare_summary_mod_codes) + .join(",") + ) ]); } @@ -625,23 +664,14 @@ impl<'a, W: Write> OutWriter> for TableWriter { "all_frac", ]); - let iter = item + for (canonical_base, mut mod_codes) in item .per_base_mod_codes .into_iter() .sorted_by(|(a, _), (b, _)| a.cmp(b)) - .map(|(primary_base, mod_codes)| { - let pass_counts = item.mod_call_counts.get(&primary_base); - let filtered_counts = - item.filtered_mod_call_counts.get(&primary_base); - (primary_base, pass_counts, filtered_counts, mod_codes) - }); - for ( - canonical_base, - pass_mod_to_counts, - filtered_counts, - mut mod_codes, - ) in iter { + let pass_mod_to_counts = item.mod_call_counts.get(&canonical_base); + let filtered_counts = + item.filtered_mod_call_counts.get(&canonical_base); let total_pass_calls = pass_mod_to_counts .map(|counts| counts.values().sum::()) .unwrap_or(0); @@ -650,62 +680,51 @@ impl<'a, W: Write> OutWriter> for TableWriter { .unwrap_or(0); let total_calls = total_filtered_calls + total_pass_calls; - let mut seen_canonical = false; if let Some(pass_counts) = pass_mod_to_counts { - for (base_state, pass_counts) in - pass_counts.iter().sorted_by(|(a, _), (b, _)| a.cmp(b)) + for base_state in pass_counts.keys() { + if let BaseState::Modified(mod_code) = base_state { + mod_codes.insert(*mod_code); + } + } + } + + let mut add_row = |base_state: BaseState, label: String| { + if let Some(pass_count) = pass_mod_to_counts + .and_then(|counts| counts.get(&base_state)) { - let label = match base_state { - BaseState::Canonical(_) => { - seen_canonical = true; - format!("-") // could be a const.. - } - BaseState::Modified(repr) => { - mod_codes.remove(repr); - format!("{repr}") - } - }; - let filtered = *item - .filtered_mod_call_counts - .get(&canonical_base) - .and_then(|filtered_counts| { - filtered_counts.get(&base_state) - }) - .unwrap_or(&0); - let all_counts = *pass_counts + filtered; + let filtered = filtered_counts + .and_then(|counts| counts.get(&base_state)) + .copied() + .unwrap_or(0); + let all_counts = *pass_count + filtered; let all_frac = all_counts as f32 / total_calls as f32; let pass_frac = - *pass_counts as f32 / total_pass_calls as f32; + *pass_count as f32 / total_pass_calls as f32; report_table.add_row(row![ canonical_base.char(), label, - pass_counts, + pass_count, pass_frac, all_counts, all_frac, ]); + } else { + report_table.add_row(row![ + canonical_base.char(), + label, + 0u64, + 0f32, + 0u64, + 0f32 + ]); } - } + }; - if !seen_canonical { - report_table.add_row(row![ - canonical_base.char(), - format!("-"), - 0u64, - 0f32, - 0u64, - 0f32 - ]); - } - for mod_code in mod_codes { - report_table.add_row(row![ - canonical_base.char(), - format!("{mod_code}"), - 0u64, - 0f32, - 0u64, - 0f32 - ]); + add_row(BaseState::Canonical(canonical_base), "-".to_string()); + for mod_code in + mod_codes.into_iter().sorted_by(compare_summary_mod_codes) + { + add_row(BaseState::Modified(mod_code), mod_code.to_string()); } } let mut report_emitted = report_table.print(&mut self.writer)?; @@ -816,21 +835,33 @@ impl<'a, W: Write> OutWriter> for TsvWriter { let mut report = String::new(); let mod_called_bases = item.mod_bases(); report.push_str(&format!("mod_bases\t{}\n", mod_called_bases)); - for (dna_base, read_count) in item.reads_with_mod_calls { + for (dna_base, read_count) in item + .reads_with_mod_calls + .into_iter() + .sorted_by(|(a, _), (b, _)| a.cmp(b)) + { report.push_str(&format!( "count_reads_{}\t{}\n", dna_base.char(), read_count )); } - for (canonical_base, mod_counts) in item.mod_call_counts { + for (canonical_base, mod_counts) in item + .mod_call_counts + .into_iter() + .sorted_by(|(a, _), (b, _)| a.cmp(b)) + { let total_calls = mod_counts.values().sum::() as f64; let total_filtered_calls = item .filtered_mod_call_counts .get(&canonical_base) .map(|filtered_counts| filtered_counts.values().sum::()) .unwrap_or(0); - for (base_state, counts) in mod_counts { + for (base_state, counts) in + mod_counts.into_iter().sorted_by(|(left, _), (right, _)| { + compare_summary_base_states(left, right) + }) + { let label = match base_state { BaseState::Canonical(_) => format!("unmodified"), BaseState::Modified(repr) => format!("modified_{repr}"), @@ -1473,3 +1504,237 @@ where Ok(total_rows as u64) } } + +#[cfg(test)] +mod tests { + use std::collections::{HashMap, HashSet}; + use std::io::{BufWriter, Write}; + + use super::{OutWriter, TableWriter, TsvWriter}; + use crate::mod_base_code::{BaseState, DnaBase, ModCodeRepr}; + use crate::summarize::ModSummary; + + const EXPECTED_TSV: &str = concat!( + "mod_bases\tA,C,G,T\n", + "count_reads_A\t1\n", + "count_reads_C\t2\n", + "count_reads_G\t3\n", + "count_reads_T\t4\n", + "A_pass_calls_unmodified\t8\n", + "A_pass_frac_unmodified\t1\n", + "A_fail_calls_unmodified\t1\n", + "A_total_mod_calls\t8\n", + "A_total_fail_mod_calls\t1\n", + "C_pass_calls_unmodified\t6\n", + "C_pass_frac_unmodified\t0.75\n", + "C_fail_calls_unmodified\t1\n", + "C_pass_calls_modified_m\t2\n", + "C_pass_frac_modified_m\t0.25\n", + "C_fail_calls_modified_m\t1\n", + "C_total_mod_calls\t8\n", + "C_total_fail_mod_calls\t2\n", + "G_pass_calls_unmodified\t5\n", + "G_pass_frac_unmodified\t1\n", + "G_fail_calls_unmodified\t0\n", + "G_total_mod_calls\t5\n", + "G_total_fail_mod_calls\t0\n", + "T_pass_calls_unmodified\t4\n", + "T_pass_frac_unmodified\t1\n", + "T_fail_calls_unmodified\t1\n", + "T_total_mod_calls\t4\n", + "T_total_fail_mod_calls\t1\n", + "total_reads_used\t10\n", + ); + + const EXPECTED_TABLE: &str = concat!( + "# bases A,C,G,T \n", + "# total_reads_used 10 \n", + "# count_reads_A 1 \n", + "# count_reads_C 2 \n", + "# count_reads_G 3 \n", + "# count_reads_T 4 \n", + "# pass_threshold_A 0.1 \n", + "# pass_threshold_C 0.2 \n", + "# pass_threshold_G 0.3 \n", + "# pass_threshold_T 0.4 \n", + "# modification_codes_for_A a \n", + "# modification_codes_for_C f,h,m \n", + " base code pass_count pass_frac all_count all_frac \n", + " A - 8 1 9 1 \n", + " A a 0 0 0 0 \n", + " C - 6 0.75 7 0.7 \n", + " C f 0 0 0 0 \n", + " C h 0 0 0 0 \n", + " C m 2 0.25 3 0.3 \n", + ); + + const EXPECTED_MIXED_CODE_TABLE: &str = concat!( + "# bases C \n", + "# total_reads_used 1 \n", + "# count_reads_C 1 \n", + "# pass_threshold_C 0.5 \n", + "# modification_codes_for_C m,123 \n", + " base code pass_count pass_frac all_count all_frac \n", + " C - 1 0.5 1 0.5 \n", + " C m 0 0 0 0 \n", + " C 123 1 0.5 1 0.5 \n", + ); + + fn summary_with_insertion_order(reverse: bool) -> ModSummary<'static> { + let base_order = if reverse { + [DnaBase::T, DnaBase::G, DnaBase::C, DnaBase::A] + } else { + [DnaBase::A, DnaBase::C, DnaBase::G, DnaBase::T] + }; + + let mut reads_with_mod_calls = HashMap::new(); + let mut mod_call_counts = HashMap::new(); + let mut filtered_mod_call_counts = HashMap::new(); + let mut per_base_thresholds = HashMap::new(); + let mut per_base_mod_codes = HashMap::new(); + for base in base_order { + reads_with_mod_calls.insert( + base, + match base { + DnaBase::A => 1, + DnaBase::C => 2, + DnaBase::G => 3, + DnaBase::T => 4, + }, + ); + per_base_thresholds.insert( + base, + match base { + DnaBase::A => 0.1, + DnaBase::C => 0.2, + DnaBase::G => 0.3, + DnaBase::T => 0.4, + }, + ); + + let canonical = BaseState::Canonical(base); + let mut pass = HashMap::new(); + let mut filtered = HashMap::new(); + if base == DnaBase::C && reverse { + pass.insert(BaseState::Modified('m'.into()), 2); + pass.insert(canonical, 6); + filtered.insert(BaseState::Modified('m'.into()), 1); + filtered.insert(canonical, 1); + } else { + pass.insert( + canonical, + match base { + DnaBase::A => 8, + DnaBase::C => 6, + DnaBase::G => 5, + DnaBase::T => 4, + }, + ); + if base == DnaBase::C { + pass.insert(BaseState::Modified('m'.into()), 2); + } + filtered.insert( + canonical, + match base { + DnaBase::A | DnaBase::C | DnaBase::T => 1, + DnaBase::G => 0, + }, + ); + if base == DnaBase::C { + filtered.insert(BaseState::Modified('m'.into()), 1); + } + } + mod_call_counts.insert(base, pass); + filtered_mod_call_counts.insert(base, filtered); + + if matches!(base, DnaBase::A | DnaBase::C) { + let codes: Vec = match (base, reverse) { + (DnaBase::A, false) => vec!['a'.into()], + (DnaBase::A, true) => vec!['a'.into()], + (DnaBase::C, false) => { + vec!['f'.into(), 'h'.into(), 'm'.into()] + } + (DnaBase::C, true) => { + vec!['m'.into(), 'h'.into(), 'f'.into()] + } + _ => unreachable!(), + }; + per_base_mod_codes.insert(base, HashSet::from_iter(codes)); + } + } + + ModSummary::new( + reads_with_mod_calls, + mod_call_counts, + filtered_mod_call_counts, + 10, + per_base_thresholds, + None, + per_base_mod_codes, + ) + } + + fn render_tsv(summary: ModSummary<'static>) -> String { + let mut writer = TsvWriter { writer: Vec::new() }; + OutWriter::write(&mut writer, summary).unwrap(); + String::from_utf8(writer.writer).unwrap() + } + + fn render_table(summary: ModSummary<'static>) -> String { + let mut writer = + TableWriter { writer: BufWriter::new(Vec::::new()) }; + OutWriter::write(&mut writer, summary).unwrap(); + writer.writer.flush().unwrap(); + String::from_utf8(writer.writer.into_inner().unwrap()).unwrap() + } + + fn mixed_code_summary(reverse: bool) -> ModSummary<'static> { + let canonical = BaseState::Canonical(DnaBase::C); + let chebi = BaseState::Modified(ModCodeRepr::ChEbi(123)); + let pass_counts = if reverse { + HashMap::from([(chebi, 1), (canonical, 1)]) + } else { + HashMap::from([(canonical, 1), (chebi, 1)]) + }; + let mod_codes = if reverse { + HashSet::from([ModCodeRepr::ChEbi(123), 'm'.into()]) + } else { + HashSet::from(['m'.into(), ModCodeRepr::ChEbi(123)]) + }; + ModSummary::new( + HashMap::from([(DnaBase::C, 1)]), + HashMap::from([(DnaBase::C, pass_counts)]), + HashMap::new(), + 1, + HashMap::from([(DnaBase::C, 0.5)]), + None, + HashMap::from([(DnaBase::C, mod_codes)]), + ) + } + + #[test] + fn summary_tsv_is_deterministic_across_insertion_orders() { + let forward = render_tsv(summary_with_insertion_order(false)); + let reverse = render_tsv(summary_with_insertion_order(true)); + + assert_eq!(forward, reverse); + assert_eq!(forward, EXPECTED_TSV); + } + + #[test] + fn summary_table_is_deterministic_across_insertion_orders() { + let forward = render_table(summary_with_insertion_order(false)); + let reverse = render_table(summary_with_insertion_order(true)); + + assert_eq!(forward, reverse); + assert_eq!(forward, EXPECTED_TABLE); + } + + #[test] + fn summary_table_uses_one_code_and_chebi_order() { + let forward = render_table(mixed_code_summary(false)); + let reverse = render_table(mixed_code_summary(true)); + assert_eq!(forward, reverse); + assert_eq!(forward, EXPECTED_MIXED_CODE_TABLE); + } +} From ba73cfc6130095014a103072a6b4f54004db9876 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 04:43:39 -0700 Subject: [PATCH 2/3] Report filtered-only summary states --- modkit-core/src/summarize.rs | 6 +- modkit-core/src/writers.rs | 309 ++++++++++++++++++++++++----------- 2 files changed, 215 insertions(+), 100 deletions(-) diff --git a/modkit-core/src/summarize.rs b/modkit-core/src/summarize.rs index dd97f5b4..c167e54c 100644 --- a/modkit-core/src/summarize.rs +++ b/modkit-core/src/summarize.rs @@ -50,7 +50,11 @@ impl<'a> ModSummary<'a> { pub(crate) fn mod_bases(&self) -> String { [DnaBase::A, DnaBase::C, DnaBase::G, DnaBase::T] .into_iter() - .filter(|base| self.mod_call_counts.contains_key(base)) + .filter(|base| { + self.mod_call_counts.contains_key(base) + || self.filtered_mod_call_counts.contains_key(base) + || self.per_base_mod_codes.contains_key(base) + }) .map(|d| d.char().to_string()) .collect::>() .join(",") diff --git a/modkit-core/src/writers.rs b/modkit-core/src/writers.rs index e449ea00..7d682cf8 100644 --- a/modkit-core/src/writers.rs +++ b/modkit-core/src/writers.rs @@ -583,19 +583,47 @@ fn compare_summary_mod_codes( } } -fn compare_summary_base_states( - left: &BaseState, - right: &BaseState, -) -> Ordering { - match (left, right) { - (BaseState::Canonical(left), BaseState::Canonical(right)) => { - left.cmp(right) - } - (BaseState::Canonical(_), BaseState::Modified(_)) => Ordering::Less, - (BaseState::Modified(_), BaseState::Canonical(_)) => Ordering::Greater, - (BaseState::Modified(left), BaseState::Modified(right)) => { - compare_summary_mod_codes(left, right) - } +fn summary_output_bases(item: &ModSummary<'_>) -> Vec { + [DnaBase::A, DnaBase::C, DnaBase::G, DnaBase::T] + .into_iter() + .filter(|base| { + item.mod_call_counts.contains_key(base) + || item.filtered_mod_call_counts.contains_key(base) + || item.per_base_mod_codes.contains_key(base) + }) + .collect() +} + +fn summary_mod_codes(item: &ModSummary<'_>, base: DnaBase) -> Vec { + let mut codes = item + .per_base_mod_codes + .get(&base) + .into_iter() + .flatten() + .copied() + .collect::>(); + for counts in [ + item.mod_call_counts.get(&base), + item.filtered_mod_call_counts.get(&base), + ] + .into_iter() + .flatten() + { + codes.extend(counts.keys().filter_map(|state| match state { + BaseState::Canonical(_) => None, + BaseState::Modified(code) => Some(*code), + })); + } + codes.sort_by(compare_summary_mod_codes); + codes.dedup(); + codes +} + +fn summary_fraction(numerator: u64, denominator: u64) -> f32 { + if denominator == 0 { + 0.0 + } else { + numerator as f32 / denominator as f32 } } @@ -608,6 +636,7 @@ impl TableWriter { impl<'a, W: Write> OutWriter> for TableWriter { fn write(&mut self, item: ModSummary<'a>) -> AnyhowResult { + let output_bases = summary_output_bases(&item); let mut metadata_table = Table::new(); let metadata_format = FormatBuilder::new().padding(1, 1).left_border('#').build(); @@ -635,20 +664,14 @@ impl<'a, W: Write> OutWriter> for TableWriter { if let Some(region) = item.region { metadata_table.add_row(row!["region", region.to_string()]); } - for (base, codes) in - item.per_base_mod_codes.iter().sorted_by(|(a, _), (b, _)| a.cmp(b)) - { - metadata_table.add_row(row![ - format!("modification_codes_for_{base}"), - format!( - "{}", - codes - .iter() - .copied() - .sorted_by(compare_summary_mod_codes) - .join(",") - ) - ]); + for base in output_bases.iter().copied() { + let codes = summary_mod_codes(&item, base); + if !codes.is_empty() { + metadata_table.add_row(row![ + format!("modification_codes_for_{base}"), + codes.into_iter().join(",") + ]); + } } let emitted = metadata_table.print(&mut self.writer)?; @@ -664,11 +687,8 @@ impl<'a, W: Write> OutWriter> for TableWriter { "all_frac", ]); - for (canonical_base, mut mod_codes) in item - .per_base_mod_codes - .into_iter() - .sorted_by(|(a, _), (b, _)| a.cmp(b)) - { + for canonical_base in output_bases { + let mod_codes = summary_mod_codes(&item, canonical_base); let pass_mod_to_counts = item.mod_call_counts.get(&canonical_base); let filtered_counts = item.filtered_mod_call_counts.get(&canonical_base); @@ -680,50 +700,30 @@ impl<'a, W: Write> OutWriter> for TableWriter { .unwrap_or(0); let total_calls = total_filtered_calls + total_pass_calls; - if let Some(pass_counts) = pass_mod_to_counts { - for base_state in pass_counts.keys() { - if let BaseState::Modified(mod_code) = base_state { - mod_codes.insert(*mod_code); - } - } - } - let mut add_row = |base_state: BaseState, label: String| { - if let Some(pass_count) = pass_mod_to_counts + let pass_count = pass_mod_to_counts .and_then(|counts| counts.get(&base_state)) - { - let filtered = filtered_counts - .and_then(|counts| counts.get(&base_state)) - .copied() - .unwrap_or(0); - let all_counts = *pass_count + filtered; - let all_frac = all_counts as f32 / total_calls as f32; - let pass_frac = - *pass_count as f32 / total_pass_calls as f32; - report_table.add_row(row![ - canonical_base.char(), - label, - pass_count, - pass_frac, - all_counts, - all_frac, - ]); - } else { - report_table.add_row(row![ - canonical_base.char(), - label, - 0u64, - 0f32, - 0u64, - 0f32 - ]); - } + .copied() + .unwrap_or(0); + let filtered = filtered_counts + .and_then(|counts| counts.get(&base_state)) + .copied() + .unwrap_or(0); + let all_counts = pass_count.saturating_add(filtered); + let all_frac = summary_fraction(all_counts, total_calls); + let pass_frac = summary_fraction(pass_count, total_pass_calls); + report_table.add_row(row![ + canonical_base.char(), + label, + pass_count, + pass_frac, + all_counts, + all_frac, + ]); }; add_row(BaseState::Canonical(canonical_base), "-".to_string()); - for mod_code in - mod_codes.into_iter().sorted_by(compare_summary_mod_codes) - { + for mod_code in mod_codes { add_row(BaseState::Modified(mod_code), mod_code.to_string()); } } @@ -833,11 +833,12 @@ impl<'a, W: Write> OutWriter> for TsvWriter { format will require the --tsv option" ); let mut report = String::new(); + let output_bases = summary_output_bases(&item); let mod_called_bases = item.mod_bases(); report.push_str(&format!("mod_bases\t{}\n", mod_called_bases)); for (dna_base, read_count) in item .reads_with_mod_calls - .into_iter() + .iter() .sorted_by(|(a, _), (b, _)| a.cmp(b)) { report.push_str(&format!( @@ -846,33 +847,45 @@ impl<'a, W: Write> OutWriter> for TsvWriter { read_count )); } - for (canonical_base, mod_counts) in item - .mod_call_counts - .into_iter() - .sorted_by(|(a, _), (b, _)| a.cmp(b)) - { - let total_calls = mod_counts.values().sum::() as f64; - let total_filtered_calls = item - .filtered_mod_call_counts - .get(&canonical_base) - .map(|filtered_counts| filtered_counts.values().sum::()) + for canonical_base in output_bases { + let pass_counts = item.mod_call_counts.get(&canonical_base); + let filtered_counts = + item.filtered_mod_call_counts.get(&canonical_base); + let total_calls = pass_counts + .map(|counts| counts.values().sum::()) .unwrap_or(0); - for (base_state, counts) in - mod_counts.into_iter().sorted_by(|(left, _), (right, _)| { - compare_summary_base_states(left, right) - }) - { - let label = match base_state { - BaseState::Canonical(_) => format!("unmodified"), - BaseState::Modified(repr) => format!("modified_{repr}"), + let total_filtered_calls = filtered_counts + .map(|counts| counts.values().sum::()) + .unwrap_or(0); + let mut states = vec![( + BaseState::Canonical(canonical_base), + "unmodified".to_string(), + )]; + states.extend( + summary_mod_codes(&item, canonical_base).into_iter().map( + |mod_code| { + ( + BaseState::Modified(mod_code), + format!("modified_{mod_code}"), + ) + }, + ), + ); + + for (base_state, label) in states { + let counts = pass_counts + .and_then(|counts| counts.get(&base_state)) + .copied() + .unwrap_or(0); + let filtered = filtered_counts + .and_then(|counts| counts.get(&base_state)) + .copied() + .unwrap_or(0); + let pass_fraction = if total_calls == 0 { + 0.0 + } else { + counts as f64 / total_calls as f64 }; - let filtered = *item - .filtered_mod_call_counts - .get(&canonical_base) - .and_then(|filtered_counts| { - filtered_counts.get(&base_state) - }) - .unwrap_or(&0); report.push_str(&format!( "{}_pass_calls_{}\t{}\n", canonical_base.char(), @@ -883,7 +896,7 @@ impl<'a, W: Write> OutWriter> for TsvWriter { "{}_pass_frac_{}\t{}\n", canonical_base.char(), label, - counts as f64 / total_calls + pass_fraction )); report.push_str(&format!( "{}_fail_calls_{}\t{}\n", @@ -895,7 +908,7 @@ impl<'a, W: Write> OutWriter> for TsvWriter { report.push_str(&format!( "{}_total_mod_calls\t{}\n", canonical_base.char(), - total_calls as u64 + total_calls )); report.push_str(&format!( "{}_total_fail_mod_calls\t{}\n", @@ -1523,11 +1536,20 @@ mod tests { "A_pass_calls_unmodified\t8\n", "A_pass_frac_unmodified\t1\n", "A_fail_calls_unmodified\t1\n", + "A_pass_calls_modified_a\t0\n", + "A_pass_frac_modified_a\t0\n", + "A_fail_calls_modified_a\t0\n", "A_total_mod_calls\t8\n", "A_total_fail_mod_calls\t1\n", "C_pass_calls_unmodified\t6\n", "C_pass_frac_unmodified\t0.75\n", "C_fail_calls_unmodified\t1\n", + "C_pass_calls_modified_f\t0\n", + "C_pass_frac_modified_f\t0\n", + "C_fail_calls_modified_f\t0\n", + "C_pass_calls_modified_h\t0\n", + "C_pass_frac_modified_h\t0\n", + "C_fail_calls_modified_h\t0\n", "C_pass_calls_modified_m\t2\n", "C_pass_frac_modified_m\t0.25\n", "C_fail_calls_modified_m\t1\n", @@ -1566,6 +1588,8 @@ mod tests { " C f 0 0 0 0 \n", " C h 0 0 0 0 \n", " C m 2 0.25 3 0.3 \n", + " G - 5 1 5 1 \n", + " T - 4 1 5 1 \n", ); const EXPECTED_MIXED_CODE_TABLE: &str = concat!( @@ -1580,6 +1604,52 @@ mod tests { " C 123 1 0.5 1 0.5 \n", ); + const EXPECTED_MIXED_CODE_TSV: &str = concat!( + "mod_bases\tC\n", + "count_reads_C\t1\n", + "C_pass_calls_unmodified\t1\n", + "C_pass_frac_unmodified\t0.5\n", + "C_fail_calls_unmodified\t0\n", + "C_pass_calls_modified_m\t0\n", + "C_pass_frac_modified_m\t0\n", + "C_fail_calls_modified_m\t0\n", + "C_pass_calls_modified_123\t1\n", + "C_pass_frac_modified_123\t0.5\n", + "C_fail_calls_modified_123\t0\n", + "C_total_mod_calls\t2\n", + "C_total_fail_mod_calls\t0\n", + "total_reads_used\t1\n", + ); + + const EXPECTED_FILTERED_ONLY_TSV: &str = concat!( + "mod_bases\tC\n", + "count_reads_C\t1\n", + "C_pass_calls_unmodified\t0\n", + "C_pass_frac_unmodified\t0\n", + "C_fail_calls_unmodified\t1\n", + "C_pass_calls_modified_m\t0\n", + "C_pass_frac_modified_m\t0\n", + "C_fail_calls_modified_m\t3\n", + "C_pass_calls_modified_123\t0\n", + "C_pass_frac_modified_123\t0\n", + "C_fail_calls_modified_123\t2\n", + "C_total_mod_calls\t0\n", + "C_total_fail_mod_calls\t6\n", + "total_reads_used\t1\n", + ); + + const EXPECTED_FILTERED_ONLY_TABLE: &str = concat!( + "# bases C \n", + "# total_reads_used 1 \n", + "# count_reads_C 1 \n", + "# pass_threshold_C 0.5 \n", + "# modification_codes_for_C m,123 \n", + " base code pass_count pass_frac all_count all_frac \n", + " C - 0 0 1 0.16666667 \n", + " C m 0 0 3 0.5 \n", + " C 123 0 0 2 0.33333334 \n", + ); + fn summary_with_insertion_order(reverse: bool) -> ModSummary<'static> { let base_order = if reverse { [DnaBase::T, DnaBase::G, DnaBase::C, DnaBase::A] @@ -1712,6 +1782,27 @@ mod tests { ) } + fn filtered_only_summary() -> ModSummary<'static> { + let canonical = BaseState::Canonical(DnaBase::C); + let methyl = BaseState::Modified('m'.into()); + let chebi = BaseState::Modified(ModCodeRepr::ChEbi(123)); + ModSummary::new( + HashMap::from([(DnaBase::C, 1)]), + HashMap::new(), + HashMap::from([( + DnaBase::C, + HashMap::from([(canonical, 1), (methyl, 3), (chebi, 2)]), + )]), + 1, + HashMap::from([(DnaBase::C, 0.5)]), + None, + HashMap::from([( + DnaBase::C, + HashSet::from(['m'.into(), ModCodeRepr::ChEbi(123)]), + )]), + ) + } + #[test] fn summary_tsv_is_deterministic_across_insertion_orders() { let forward = render_tsv(summary_with_insertion_order(false)); @@ -1737,4 +1828,24 @@ mod tests { assert_eq!(forward, reverse); assert_eq!(forward, EXPECTED_MIXED_CODE_TABLE); } + + #[test] + fn summary_tsv_uses_one_code_and_chebi_order() { + let forward = render_tsv(mixed_code_summary(false)); + let reverse = render_tsv(mixed_code_summary(true)); + assert_eq!(forward, reverse); + assert_eq!(forward, EXPECTED_MIXED_CODE_TSV); + } + + #[test] + fn summary_outputs_filtered_only_states_exactly() { + assert_eq!( + render_tsv(filtered_only_summary()), + EXPECTED_FILTERED_ONLY_TSV + ); + assert_eq!( + render_table(filtered_only_summary()), + EXPECTED_FILTERED_ONLY_TABLE + ); + } } From 413cffd1bc1451901a5d58c3b63acfafb3210b5b Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Tue, 4 Aug 2026 22:55:24 -0700 Subject: [PATCH 3/3] Exclude unmapped calls from CpG summaries --- modkit-core/src/modbam_util/subcommands.rs | 6 +- modkit/tests/test_summary.rs | 75 +++++++++++++++++++++- 2 files changed, 79 insertions(+), 2 deletions(-) diff --git a/modkit-core/src/modbam_util/subcommands.rs b/modkit-core/src/modbam_util/subcommands.rs index 7d13ad20..fa6ee56f 100644 --- a/modkit-core/src/modbam_util/subcommands.rs +++ b/modkit-core/src/modbam_util/subcommands.rs @@ -1938,7 +1938,11 @@ impl ModSummarize { multi_progress.clone(), )?; - if self.mapped_only || self.matched_only || self.motif.is_some() { + if self.mapped_only + || self.matched_only + || self.motif.is_some() + || self.cpg + { qual_hist } else if let Some(nr) = self.num_reads { if qual_hist.ok_records < nr { diff --git a/modkit/tests/test_summary.rs b/modkit/tests/test_summary.rs index bcfae7d1..559fe59d 100644 --- a/modkit/tests/test_summary.rs +++ b/modkit/tests/test_summary.rs @@ -5,11 +5,84 @@ use crate::common::{ use anyhow::Context; use mod_kit::mod_bam::{CollapseMethod, EdgeFilter}; use mod_kit::mod_base_code::{BaseState, DnaBase}; -use std::collections::HashSet; +use rust_htslib::bam::{self, record::Aux, Read}; +use std::collections::{HashMap, HashSet}; use std::path::Path; +use std::process::Command; +use tempfile::tempdir; mod common; +fn run_cpg_summary(bam_fp: &Path) -> HashMap { + let output = Command::new(env!("CARGO_BIN_EXE_modkit")) + .args([ + "summary", + bam_fp.to_str().unwrap(), + "--reference", + "../tests/resources/CGI_ladder_3.6kb_ref.fa", + "--cpg", + "--filter-threshold", + "0", + "--tsv", + "--suppress-progress", + "--threads", + "1", + "--io-threads", + "1", + ]) + .output() + .unwrap(); + assert!( + output.status.success(), + "summary failed: {}", + String::from_utf8_lossy(&output.stderr) + ); + String::from_utf8(output.stdout) + .unwrap() + .lines() + .map(|line| { + let (key, value) = line.split_once('\t').unwrap(); + (key.to_string(), value.to_string()) + }) + .collect() +} + +#[test] +fn test_summary_cpg_excludes_unmapped_calls() { + let mapped_bam = + Path::new("../tests/resources/bc_anchored_10_reads.sorted.bam"); + let temp_dir = tempdir().unwrap(); + let mixed_bam = temp_dir.path().join("mapped-and-unmapped.bam"); + let mut reader = bam::Reader::from_path(mapped_bam).unwrap(); + let header = bam::Header::from_template(reader.header()); + let mut writer = + bam::Writer::from_path(&mixed_bam, &header, bam::Format::Bam).unwrap(); + for record in reader.records() { + writer.write(&record.unwrap()).unwrap(); + } + + let mut unmapped = bam::Record::new(); + unmapped.set(b"unmapped", None, b"CCCC", &[255; 4]); + unmapped.set_unmapped(); + unmapped.set_tid(-1); + unmapped.set_pos(-1); + unmapped.set_mtid(-1); + unmapped.set_mpos(-1); + unmapped.push_aux(b"MM", Aux::String("C+m?,0,0,0,0;")).unwrap(); + unmapped.push_aux(b"ML", Aux::ArrayU8((&[255; 4][..]).into())).unwrap(); + unmapped.push_aux(b"MN", Aux::I32(4)).unwrap(); + writer.write(&unmapped).unwrap(); + drop(writer); + bam::index::build(&mixed_bam, None, bam::index::Type::Bai, 1).unwrap(); + + let mixed_summary = run_cpg_summary(&mixed_bam); + let mapped_summary = run_cpg_summary(mapped_bam); + assert_eq!(mixed_summary["total_reads_used"], "10"); + assert_eq!(mixed_summary["count_reads_C"], "10"); + assert_eq!(mixed_summary["C_total_mod_calls"], "77"); + assert_eq!(mixed_summary, mapped_summary); +} + /// tests that the summary from a BAM is the same if run with the BAI /// or without (just taking the first N reads). In this case we use all /// of the reads in the BAM