From 697de7be719148b94a63108299d3cc3f4bb38283 Mon Sep 17 00:00:00 2001 From: ljwharbers Date: Wed, 2 Sep 2026 11:56:21 +0200 Subject: [PATCH 1/2] Keep records with conflicting mod probabilities; fix pileup on PacBio tags PacBio Jasmine (>= 26.1.3) calls 5mC and 5hmC with independent models and writes them as separate `C+m?` and `C+h?`/`G-h?` sub-tags, so at a CpG the probabilities can sum to well over 1.0. modkit rejected the entire record (`conflict-explicit-prob-greater-than-one`) in adjust-mods, update-tags, call-mods, extract, summary, sample-probs and the older pileup paths, which removed 32-65% of the reads of affected samples, mostly silently. - mod_bam: positions whose probabilities from separate sub-tags sum to more than one are now dropped instead of the whole record. When the caller has a collapse method (`--ignore h`, `--convert h m`) the position is resolved with it first (h removed / summed and saturated at 1.0), so the 5mC call is kept. The number of dropped positions is counted and reported at the end of pileup, adjust-mods, update-tags, extract and friends; check-tags still reports the records as invalid. - pileup: the same positions are skipped in the htslib-based general worker and in the optimized `BaseModsAdapter` (previously accepted unchecked). - pileup: a record without MM/ML tags no longer aborts the whole interval (the aggregator only logged the failure at DEBUG and dropped all rows of the interval). Failed intervals are now reported as errors. - pileup: `BaseModsAdapter` accepts opposite-strand sub-tags (`T-a.`), fixing threshold estimation on PacBio HiFi BAMs, and `pileup` falls back to the general workers when such calls are present instead of failing every record in the optimized workers. - tests: synthetic PacBio-style BAMs (tests/make_pacbio_style_tags.py) with unit and integration tests; CHANGELOG and troubleshooting docs updated. Co-Authored-By: Claude Fable 5.1 --- CHANGELOG.md | 6 + book/src/troubleshooting.md | 27 ++ modkit-core/src/adjust.rs | 3 +- modkit-core/src/extract/subcommand.rs | 2 + modkit-core/src/extract/util.rs | 9 +- modkit-core/src/mod_bam.rs | 417 ++++++++++++++++-- modkit-core/src/modbam_util/check_tags.rs | 18 +- modkit-core/src/modbam_util/subcommands.rs | 1 + modkit-core/src/pileup/base_mods_adapter.rs | 113 ++++- modkit-core/src/pileup/pileup_processor.rs | 38 +- modkit-core/src/pileup/subcommand.rs | 79 +++- modkit-core/src/read_cache.rs | 4 +- modkit-core/src/read_ids_to_base_mod_probs.rs | 12 +- modkit-core/src/util.rs | 7 + modkit/tests/test_pacbio_style_tags.rs | 222 ++++++++++ tests/make_pacbio_style_tags.py | 79 ++++ tests/resources/pacbio_style_tags.bam | Bin 0 -> 5345 bytes tests/resources/pacbio_style_tags.bam.bai | Bin 0 -> 360 bytes .../pacbio_style_tags_no_untagged.bam | Bin 0 -> 5030 bytes .../pacbio_style_tags_no_untagged.bam.bai | Bin 0 -> 360 bytes 20 files changed, 966 insertions(+), 71 deletions(-) create mode 100644 modkit/tests/test_pacbio_style_tags.rs create mode 100644 tests/make_pacbio_style_tags.py create mode 100644 tests/resources/pacbio_style_tags.bam create mode 100644 tests/resources/pacbio_style_tags.bam.bai create mode 100644 tests/resources/pacbio_style_tags_no_untagged.bam create mode 100644 tests/resources/pacbio_style_tags_no_untagged.bam.bai diff --git a/CHANGELOG.md b/CHANGELOG.md index d74fa17f..afda1224 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,12 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Unreleased] +### Fixes + - [all] Records where the base modification probabilities from separate MM sub-tags sum to more than one at a position (e.g. PacBio Jasmine 5mC and 5hmC calls, which are made by independent models) are no longer discarded. The offending positions are dropped (or resolved with `--ignore`/`--convert` when given, e.g. `--ignore h` keeps the 5mC call), the rest of the record is used, and the number of dropped positions is reported. `modbam check-tags` still reports these records as `conflict-explicit-prob-greater-than-one`. + - [pileup] A record without MM/ML tags no longer aborts the processing of the whole interval, which silently produced no output for that interval. + - [pileup] Base modification calls on the opposite strand of the primary sequence base (e.g. PacBio 6mA `T-a.` calls) no longer make every record fail in the optimized workers and the threshold estimation; the general workers are used automatically when such calls are found. + ## [v0.6.4] ### Adds - [bedmethyl] Adds `--min-samples` (an integer or `all`) and `--min-sample-coverage` to `bedmethyl merge` to require a position to be present in (and optionally covered to a minimum valid depth in) multiple inputs, enabling an inner join across replicates. Omitting them preserves the original outer-join behaviour. diff --git a/book/src/troubleshooting.md b/book/src/troubleshooting.md index 77fa2643..4abbb749 100644 --- a/book/src/troubleshooting.md +++ b/book/src/troubleshooting.md @@ -71,3 +71,30 @@ contains CG positions. However, it will not include positions for which the pass is zero (see [the column descriptions](./intro_pileup.md#description-of-bedmethyl-output)). This is to be expected. + +## PacBio (Jasmine) 5mC and 5hmC calls: "conflict-explicit-prob-greater-than-one" + +PacBio Jasmine (>= 26.1.3) calls 5mC and 5hmC with two independent models and writes them as +separate `C+m?` and `C+h?` (and `G-h?`) sub-tags. Because the models are independent, the two +probabilities at a single cytosine can sum to more than 1.0, which the SAM specification does not +allow. Older versions of `modkit` discarded the entire record when this happened (visible as +`conflict-explicit-prob-greater-than-one` in `modkit modbam check-tags` and in the debug log), which +could remove more than half of the reads of a sample. + +`modkit` now keeps the record and only drops the positions where the probabilities cannot be +reconciled; a summary of how many positions were dropped is logged at the end of the run. To keep +the 5mC calls at those positions, pass `--ignore h` (the 5hmC probability is removed at the +conflicting positions and the 5mC probability is left untouched) or `--convert h m` (the +probabilities are summed, saturating at 1.0), for example: + +```bash +modkit adjust-mods --ignore h ${pacbio_bam} ${out_bam} +modkit pileup ${out_bam} ${out_bed} --cpg --ref ${ref} +``` + +`modkit modbam check-tags` still reports these records as invalid so that non-conformant tags can +be detected. + +PacBio HiFi reads also carry 6mA calls on both strands (`A+a.` and `T-a.`). `modkit pileup` detects +the opposite-strand calls and uses the general pileup workers for these files, which is slower than +the optimized workers used for ONT data but produces the same output. diff --git a/modkit-core/src/adjust.rs b/modkit-core/src/adjust.rs index 61e3f2ba..f38b3e33 100644 --- a/modkit-core/src/adjust.rs +++ b/modkit-core/src/adjust.rs @@ -120,7 +120,7 @@ fn adjust_mod_probs<'a>( sequence_motifs: &Option>, discard_motifs: bool, ) -> MkResult { - let mod_base_info = ModBaseInfo::new_from_record(&record)?; + let mod_base_info = ModBaseInfo::new_from_record_with(&record, methods)?; let mm_style = mod_base_info.mm_style; let ml_style = mod_base_info.ml_style; @@ -302,6 +302,7 @@ pub fn adjust_modbam( spinner.finish_and_clear(); info!("done, {} records processed", total,); + crate::mod_bam::report_conflict_summary(); if !error_counts.is_empty() { info!("error/skip counts:"); diff --git a/modkit-core/src/extract/subcommand.rs b/modkit-core/src/extract/subcommand.rs index 42b4d6ed..23dc7dff 100644 --- a/modkit-core/src/extract/subcommand.rs +++ b/modkit-core/src/extract/subcommand.rs @@ -354,6 +354,7 @@ impl EntryExtractFull { n_skipped.finish_and_clear(); n_used.finish_and_clear(); n_rows.finish_and_clear(); + crate::mod_bam::report_conflict_summary(); info!( "processed {} reads, {} rows, skipped ~{} reads, failed ~{} reads", writer.num_reads(), @@ -864,6 +865,7 @@ impl EntryExtractCalls { n_skipped.finish_and_clear(); n_used.finish_and_clear(); n_rows.finish_and_clear(); + crate::mod_bam::report_conflict_summary(); info!( "processed {} reads, {} rows, skipped ~{} reads, failed ~{} reads", writer.num_reads(), diff --git a/modkit-core/src/extract/util.rs b/modkit-core/src/extract/util.rs index e628a9b5..6787ed4c 100644 --- a/modkit-core/src/extract/util.rs +++ b/modkit-core/src/extract/util.rs @@ -571,8 +571,13 @@ fn process_records_to_chan<'a, T: Read>( message: &'static str, kmer_size: usize, ) -> (usize, usize) { - let mut mod_iter = - TrackingModRecordIter::new(records, false, allow_non_primary); + let resolvers = collapse_method.cloned().into_iter().collect(); + let mut mod_iter = TrackingModRecordIter::new( + records, + false, + allow_non_primary, + resolvers, + ); let pb = multi_pb.add(get_ticker()); pb.set_message(format!("{message}records processed")); for (record, read_id, mod_base_info) in &mut mod_iter { diff --git a/modkit-core/src/mod_bam.rs b/modkit-core/src/mod_bam.rs index ef6a80e7..85f97c17 100644 --- a/modkit-core/src/mod_bam.rs +++ b/modkit-core/src/mod_bam.rs @@ -2,11 +2,13 @@ use std::cmp::Ordering; use std::collections::{HashMap, HashSet}; use std::fmt::{Debug, Display, Formatter}; use std::hash::Hash; +use std::sync::atomic::{AtomicUsize, Ordering as AtomicOrdering}; use anyhow::bail; use derive_new::new; use itertools::{Itertools, PeekingNext}; -use log::debug; +use log::{debug, warn}; +use log_once::warn_once; use nom::bytes::complete::tag; use nom::character::complete::{digit1, multispace0}; use nom::multi::separated_list1; @@ -20,14 +22,100 @@ use crate::errs::{ConflictError, MkError, MkResult}; use crate::mod_base_code::{DnaBase, ModCodeRepr, ParseChar}; use crate::motifs::iupac::nt_bytes; use crate::util::{ - get_forward_sequence, get_tag, record_is_not_primary, Strand, + get_forward_sequence, get_query_name_string, get_tag, + record_is_not_primary, Strand, }; -const MAX_PROB: f32 = 1.01f32; +/// Maximum allowed sum of base modification probabilities at a single +/// position, slightly above 1.0 to allow for rounding when converting from the +/// 8-bit ML encoding. +pub(crate) const MAX_PROB: f32 = 1.01f32; + +/// Number of positions dropped because the explicit probabilities coming from +/// separate MM sub-tags summed to more than one (e.g. PacBio Jasmine 5mC and +/// 5hmC calls, which are made by independent models). +static CONFLICT_POSITIONS: AtomicUsize = AtomicUsize::new(0); +/// Number of records that had at least one such position dropped. +static CONFLICT_RECORDS: AtomicUsize = AtomicUsize::new(0); + +/// Record that `n_positions` positions of a record were dropped because the +/// modification probabilities at those positions summed to more than one. +#[inline] +pub(crate) fn note_conflict_positions(n_positions: usize) { + if n_positions > 0 { + CONFLICT_POSITIONS.fetch_add(n_positions, AtomicOrdering::Relaxed); + CONFLICT_RECORDS.fetch_add(1, AtomicOrdering::Relaxed); + warn_once!( + "found position(s) where the base modification probabilities from \ + separate MM sub-tags sum to more than 1.0 (non-conformant tags, \ + e.g. PacBio Jasmine 5mC and 5hmC calls made by independent \ + models). These positions are dropped and the records are kept, \ + use --ignore or --convert to resolve the conflicting \ + modification code instead. This warning is only shown once, \ + per-record details are logged at DEBUG level and totals are \ + reported at the end." + ); + } +} + +/// Record that a single position was skipped in a streaming parser (where +/// records are not tracked). +#[inline] +pub(crate) fn note_conflict_position() { + CONFLICT_POSITIONS.fetch_add(1, AtomicOrdering::Relaxed); + warn_once!( + "found position(s) where the base modification probabilities from \ + separate MM sub-tags sum to more than 1.0 (non-conformant tags, e.g. \ + PacBio Jasmine 5mC and 5hmC calls made by independent models). These \ + positions are skipped and the records are kept. This warning is only \ + shown once, totals are reported at the end." + ); +} + +/// Returns (number of records, number of positions) dropped due to base +/// modification probabilities summing to more than one. +pub fn conflict_counts() -> (usize, usize) { + ( + CONFLICT_RECORDS.load(AtomicOrdering::Relaxed), + CONFLICT_POSITIONS.load(AtomicOrdering::Relaxed), + ) +} + +/// Log a summary of how many positions/records were affected by base +/// modification probabilities summing to more than one, if any. +pub fn report_conflict_summary() { + let (n_records, n_positions) = conflict_counts(); + if n_positions > 0 { + let records_message = if n_records > 0 { + format!(" in {n_records} record(s)") + } else { + String::new() + }; + warn!( + "dropped {n_positions} position(s){records_message} where the \ + base modification probabilities summed to more than 1.0, the \ + records were otherwise used. Use --ignore or --convert to \ + resolve the conflicting modification code, or `modkit modbam \ + check-tags` to inspect the modBAM." + ); + } +} + +/// Check whether the sum of a collection of base modification probabilities +/// exceeds the maximum allowed. +#[inline] +pub(crate) fn probs_exceed_max<'a>( + probs: impl Iterator, +) -> bool { + probs.sum::() > MAX_PROB +} pub(crate) struct TrackingModRecordIter<'a, T: bam::Read> { records: bam::Records<'a, T>, skip_unmapped: bool, allow_non_primary: bool, + /// collapse methods used to resolve positions where probabilities sum to + /// more than one, see `ModBaseInfo::new_from_record_with`. + resolvers: Vec, pub(crate) num_used: usize, pub(crate) num_skipped: usize, pub(crate) num_failed: usize, @@ -38,11 +126,13 @@ impl<'a, T: bam::Read> TrackingModRecordIter<'a, T> { records: bam::Records<'a, T>, skip_unmapped: bool, allow_non_primary: bool, + resolvers: Vec, ) -> Self { Self { records, skip_unmapped, allow_non_primary, + resolvers, num_used: 0, num_skipped: 0, num_failed: 0, @@ -80,7 +170,10 @@ impl<'a, T: bam::Read> Iterator for &mut TrackingModRecordIter<'a, T> { self.num_failed += 1; continue; } else { - match ModBaseInfo::new_from_record(&record) { + match ModBaseInfo::new_from_record_with( + &record, + &self.resolvers, + ) { Ok(modbase_info) => { if modbase_info.is_empty() { self.num_skipped += 1; @@ -127,13 +220,19 @@ pub(crate) struct ModBaseInfoRecordTracker< // total: usize, num_errors: usize, records: I, + /// collapse methods used to resolve positions where probabilities sum to + /// more than one, see `ModBaseInfo::new_from_record_with`. + resolvers: Vec, } pub(crate) trait WithModBaseInfos< I: Iterator>, > { - fn with_mod_base_info(self) -> ModBaseInfoRecordTracker + fn with_mod_base_info( + self, + resolvers: Vec, + ) -> ModBaseInfoRecordTracker where Self: Iterator> + Sized, @@ -142,6 +241,7 @@ pub(crate) trait WithModBaseInfos< // total: 0, num_errors: 0, records: self, + resolvers, } } } @@ -163,7 +263,10 @@ impl>> if record_is_not_primary(&record) || record.seq_len() == 0 { continue; } - match ModBaseInfo::new_from_record(&record) { + match ModBaseInfo::new_from_record_with( + &record, + &self.resolvers, + ) { Ok(modbase_info) => { if modbase_info.is_empty() { continue; @@ -640,6 +743,47 @@ impl BaseModProbs { self.check() } + /// Attempt to resolve a position whose probabilities sum to more than one + /// using the user-requested collapse methods. Ignoring a modification code + /// (either method) removes that code's probability outright (the usual + /// re-normalization is undefined when the canonical probability is + /// negative); converting sums the probabilities into the target code, + /// saturating at 1.0. Returns true when the position is valid afterwards. + pub(crate) fn resolve_conflict( + &mut self, + resolvers: &[CollapseMethod], + ) -> bool { + if self.check().is_ok() { + return true; + } + for method in resolvers { + match method { + CollapseMethod::ReNormalize(code) + | CollapseMethod::ReDistribute(code) => { + self.probs.remove(code); + } + CollapseMethod::Convert { from, to } => { + let mut converted = 0f32; + let mut any = false; + for code in from { + if let Some(p) = self.probs.remove(code) { + converted += p; + any = true; + } + } + if any { + let q = self.probs.entry(*to).or_insert(0f32); + *q = (*q + converted).min(1f32); + } + } + } + if self.check().is_ok() { + return true; + } + } + false + } + fn check(&self) -> MkResult<()> { let x = self.probs.values().sum::(); if x > MAX_PROB { @@ -915,6 +1059,13 @@ fn parse_int_list<'a>(input: &'a str) -> IResult<&'a str, Vec> { } impl MmTagInfo { + /// Strand of the modification relative to the primary sequence base, + /// `-` (Negative) means the modification is on the opposite strand, e.g. + /// `T-a` for 6mA on the complement of a thymine. + pub(crate) fn strand(&self) -> Strand { + self.strand + } + pub(crate) fn from_record(record: &bam::Record) -> MkResult> { let raw_mod_tags = parse_raw_mod_tags(record)?; Self::parse_mm_tag(&raw_mod_tags.raw_mm) @@ -1065,23 +1216,42 @@ impl MmTagInfo { } } +/// Combine the base modification probabilities from another MM sub-tag into +/// `agg`. When the probabilities at a position sum to more than one (the +/// sub-tags come from independent models, e.g. PacBio Jasmine 5mC and 5hmC) +/// the position is resolved with `resolvers` if possible, otherwise it is +/// dropped. Returns the number of dropped positions. fn combine_positions_to_probs( agg: &mut SeqPosBaseModProbs, to_add: SeqPosBaseModProbs, -) -> MkResult<()> { + resolvers: &[CollapseMethod], +) -> MkResult { if agg.skip_mode != to_add.skip_mode { agg.skip_mode = SkipMode::ImplicitUnmodified; } + let mut n_dropped = 0usize; for (position, base_mod_probs) in to_add.pos_to_base_mod_probs.into_iter() { - if let Some(probs) = agg.pos_to_base_mod_probs.get_mut(&position) { - probs.combine_checked(base_mod_probs)?; - } else { - agg.pos_to_base_mod_probs.insert(position, base_mod_probs); + let drop_position = + if let Some(probs) = agg.pos_to_base_mod_probs.get_mut(&position) { + match probs.combine_checked(base_mod_probs) { + Ok(()) => false, + Err(MkError::Conflict( + ConflictError::ProbaGreaterThanOne, + )) => !probs.resolve_conflict(resolvers), + Err(e) => return Err(e), + } + } else { + agg.pos_to_base_mod_probs.insert(position, base_mod_probs); + false + }; + if drop_position { + agg.pos_to_base_mod_probs.remove(&position); + n_dropped += 1; } } - Ok(()) + Ok(n_dropped) } // pub type SeqPosBaseModProbs = HashMap; @@ -1233,6 +1403,7 @@ pub fn extract_mod_probs( combine_positions_to_probs( &mut positions_to_probs, seq_pos_base_mod_probs, + &[], )?; } pointer += mm_tag_info.delta_list.len() * mm_tag_info.stride(); @@ -1509,22 +1680,55 @@ pub struct ModBaseInfo { converters: HashMap, pub mm_style: &'static str, pub ml_style: &'static str, + /// Number of positions that were dropped because the base modification + /// probabilities from separate MM sub-tags summed to more than one. + pub n_conflict_positions: usize, } impl ModBaseInfo { pub fn new_from_record(record: &bam::Record) -> MkResult { + Self::new_from_record_with(record, &[]) + } + + /// Parse the base modification information from a record, using + /// `resolvers` (the collapse methods requested by the user, e.g. + /// `--ignore h`) to resolve positions where the probabilities from separate + /// MM sub-tags sum to more than one. Positions that cannot be resolved are + /// dropped and counted, the record is still returned. + pub fn new_from_record_with( + record: &bam::Record, + resolvers: &[CollapseMethod], + ) -> MkResult { let raw_mod_tags = parse_raw_mod_tags(record)?; let forward_sequence = get_forward_sequence(record); let mm_tag_infos = MmTagInfo::parse_mm_tag(&raw_mod_tags.raw_mm)?; - Self::new(&mm_tag_infos, &raw_mod_tags, &forward_sequence) + let mod_base_info = Self::new( + &mm_tag_infos, + &raw_mod_tags, + &forward_sequence, + resolvers, + )?; + if mod_base_info.n_conflict_positions > 0 { + note_conflict_positions(mod_base_info.n_conflict_positions); + let read_id = get_query_name_string(record) + .unwrap_or_else(|_| "'UTF-8 decode failure'".to_string()); + debug!( + "{read_id}: dropped {} position(s) where base modification \ + probabilities summed to more than 1.0", + mod_base_info.n_conflict_positions + ); + } + Ok(mod_base_info) } pub fn new( tag_infos: &[MmTagInfo], raw_mod_tags: &RawModTags, forward_seq: &[u8], + resolvers: &[CollapseMethod], ) -> MkResult { let raw_ml = &raw_mod_tags.raw_ml; + let mut n_conflict_positions = 0usize; // todo make these DnaBase keys.. let mut pos_seq_base_mod_probs = @@ -1562,7 +1766,8 @@ impl ModBaseInfo { let agg = seq_base_mod_probs.entry(base).or_insert_with(|| { SeqPosBaseModProbs::new_empty(mm_tag_info.mode) }); - combine_positions_to_probs(agg, to_add)? + n_conflict_positions += + combine_positions_to_probs(agg, to_add, resolvers)? } pointer += mm_tag_info.delta_list.len() * mm_tag_info.stride(); @@ -1607,6 +1812,7 @@ impl ModBaseInfo { converters, mm_style: raw_mod_tags.mm_style, ml_style: raw_mod_tags.ml_style, + n_conflict_positions, }) } @@ -2161,6 +2367,132 @@ mod mod_bam_tests { ); } + #[test] + fn test_mod_base_info_conflicting_probs() { + // PacBio Jasmine-style tags: 5mC and 5hmC are called by independent + // models, so their probabilities can sum to more than one at a + // position. The C+h? track is sparse and only lists confident calls. + // 0123456789 + let dna = "ACGTCGACGT"; + // C positions 1, 4, 7; the h call is on the second C (delta 1) + let tag = "C+h?,1;C+m?,0,0,0;"; + let quals = vec![178u16, 255, 255, 255]; + let raw_mod_tags = RawModTags::new(tag, &quals, true); + let mm_tag_infos = MmTagInfo::parse_mm_tag(tag).unwrap(); + + // without a resolver the conflicting position is dropped, the record + // and the other positions are kept + let mbi = + ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes(), &[]) + .unwrap(); + assert_eq!(mbi.n_conflict_positions, 1); + let c_probs = &mbi + .pos_seq_base_mod_probs + .get(&DnaBase::C) + .unwrap() + .pos_to_base_mod_probs; + assert_eq!(c_probs.len(), 2); + assert!(c_probs.contains_key(&1)); + assert!(!c_probs.contains_key(&4)); + assert!(c_probs.contains_key(&7)); + + // a resolver for a code that is not involved does not help + let resolvers = [CollapseMethod::ReDistribute('a'.into())]; + let mbi = ModBaseInfo::new( + &mm_tag_infos, + &raw_mod_tags, + dna.as_bytes(), + &resolvers, + ) + .unwrap(); + assert_eq!(mbi.n_conflict_positions, 1); + + // the expected probabilities when only the m calls are present + let m_tag = "C+m?,0,0,0;"; + let m_quals = vec![255u16, 255, 255]; + let m_raw_mod_tags = RawModTags::new(m_tag, &m_quals, true); + let m_mm_tag_infos = MmTagInfo::parse_mm_tag(m_tag).unwrap(); + let expected = ModBaseInfo::new( + &m_mm_tag_infos, + &m_raw_mod_tags, + dna.as_bytes(), + &[], + ) + .unwrap(); + let expected_c_probs = &expected + .pos_seq_base_mod_probs + .get(&DnaBase::C) + .unwrap() + .pos_to_base_mod_probs; + + // ignoring h (with either method) removes the h probability at the + // conflicting position and leaves the m call untouched + for resolver in [ + CollapseMethod::ReDistribute('h'.into()), + CollapseMethod::ReNormalize('h'.into()), + ] { + let mbi = ModBaseInfo::new( + &mm_tag_infos, + &raw_mod_tags, + dna.as_bytes(), + &[resolver], + ) + .unwrap(); + assert_eq!(mbi.n_conflict_positions, 0); + let c_probs = &mbi + .pos_seq_base_mod_probs + .get(&DnaBase::C) + .unwrap() + .pos_to_base_mod_probs; + assert_eq!(c_probs.len(), 3); + assert_eq!(c_probs.get(&4), expected_c_probs.get(&4)); + assert_eq!(c_probs.get(&1), expected_c_probs.get(&1)); + assert_eq!(c_probs.get(&7), expected_c_probs.get(&7)); + } + + // converting h to m sums the probabilities, saturating at 1.0 + let resolvers = [CollapseMethod::Convert { + from: HashSet::from(['h'.into()]), + to: 'm'.into(), + }]; + let mbi = ModBaseInfo::new( + &mm_tag_infos, + &raw_mod_tags, + dna.as_bytes(), + &resolvers, + ) + .unwrap(); + assert_eq!(mbi.n_conflict_positions, 0); + let c_probs = &mbi + .pos_seq_base_mod_probs + .get(&DnaBase::C) + .unwrap() + .pos_to_base_mod_probs; + assert_eq!(c_probs.len(), 3); + assert_eq!( + c_probs.get(&4).unwrap(), + &BaseModProbs::new_init('m', 1.0f32) + ); + + // rounding: two probabilities that sum to slightly more than one + // (as happens with 8-bit encoding) are still accepted + let tag = "C+h.,0;C+m.,0;"; + let quals = vec![1u16, 255]; + let raw_mod_tags = RawModTags::new(tag, &quals, true); + let mm_tag_infos = MmTagInfo::parse_mm_tag(tag).unwrap(); + let mbi = + ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes(), &[]) + .unwrap(); + assert_eq!(mbi.n_conflict_positions, 0); + let c_probs = &mbi + .pos_seq_base_mod_probs + .get(&DnaBase::C) + .unwrap() + .pos_to_base_mod_probs; + assert_eq!(c_probs.len(), 3); + assert_eq!(c_probs.get(&1).unwrap().iter_probs().count(), 2); + } + #[test] fn test_format_mm_ml_tags() { let canonical_base = FundamentalBase::C; @@ -2441,7 +2773,7 @@ mod mod_bam_tests { let raw_mod_tags = RawModTags::new(tag, &quals, true); let mm_tag_infos = MmTagInfo::parse_mm_tag(tag).unwrap(); let obs_mod_base_info = - ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes()) + ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes(), &[]) .unwrap(); assert_eq!( obs_mod_base_info.pos_seq_base_mod_probs.get(&DnaBase::C).unwrap(), @@ -2476,7 +2808,7 @@ mod mod_bam_tests { let quals = vec![1, 1, 1, 100, 100, 100, 200, 200, 200]; let raw_mod_tags = RawModTags::new(tag, &quals, true); let obs_mod_base_info = - ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes()) + ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes(), &[]) .unwrap(); assert_eq!( obs_mod_base_info.pos_seq_base_mod_probs.get(&DnaBase::C).unwrap(), @@ -2512,7 +2844,7 @@ mod mod_bam_tests { let quals = vec![1, 100, 1, 100, 1, 100, 200, 200, 200]; let raw_mod_tags = RawModTags::new(tag, &quals, true); let obs_mod_base_info = - ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes()) + ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes(), &[]) .unwrap(); assert_eq!( obs_mod_base_info.pos_seq_base_mod_probs.get(&DnaBase::C).unwrap(), @@ -2550,8 +2882,8 @@ mod mod_bam_tests { let mm_tag_infos = MmTagInfo::parse_mm_tag(tag).unwrap(); let quals = vec![100, 100, 100, 1, 1, 1, 150, 150, 150, 2, 2, 2]; let tags = RawModTags::new(tag, &quals, true); - let info = - ModBaseInfo::new(&mm_tag_infos, &tags, dna.as_bytes()).unwrap(); + let info = ModBaseInfo::new(&mm_tag_infos, &tags, dna.as_bytes(), &[]) + .unwrap(); let inferred = false; let (_converters, iterator) = info.into_iter_base_mod_probs(); for (c, strand, probs) in iterator { @@ -2612,7 +2944,7 @@ mod mod_bam_tests { let quals = vec![100, 100, 100, 1, 1, 1, 150, 150, 150, 2, 2, 2]; let raw_mod_tags = RawModTags::new(tag, &quals, true); let obs_mod_base_info = - ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes()) + ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes(), &[]) .unwrap(); let top_strand_mods = @@ -2660,7 +2992,7 @@ mod mod_bam_tests { let quals = vec![1, 1, 1, 200, 200, 200, 100, 100, 100]; let raw_mod_tags = RawModTags::new(tag, &quals, true); let mut obs_mod_base_info = - ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes()) + ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes(), &[]) .unwrap(); let c_seq_base_mod_probs = obs_mod_base_info .pos_seq_base_mod_probs @@ -2691,7 +3023,7 @@ mod mod_bam_tests { // trim larger than read let edge_filter = EdgeFilter::new(50, 50, false); let mut obs_mod_base_info = - ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes()) + ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes(), &[]) .unwrap(); let c_seq_base_mod_probs = obs_mod_base_info .pos_seq_base_mod_probs @@ -2703,7 +3035,7 @@ mod mod_bam_tests { // trim with mod call _at_ the position to be trimmed let mut obs_mod_base_info = - ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes()) + ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes(), &[]) .unwrap(); let c_seq_base_mod_probs = obs_mod_base_info .pos_seq_base_mod_probs @@ -2744,14 +3076,14 @@ mod mod_bam_tests { let mm_tag_infos = MmTagInfo::parse_mm_tag(tag).unwrap(); let raw_mod_tags = RawModTags::new(tag, &quals, true); let obs_mod_base_info = - ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes()) + ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes(), &[]) .unwrap(); assert!(obs_mod_base_info.is_empty()); let tag = "C+h.;C+m.;"; let mm_tag_infos = MmTagInfo::parse_mm_tag(tag).unwrap(); let raw_mod_tags = RawModTags::new(tag, &quals, true); let obs_mod_base_info = - ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes()) + ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes(), &[]) .unwrap(); assert!(!obs_mod_base_info.is_empty()); // g c CG c gg CG CG @@ -2762,14 +3094,14 @@ mod mod_bam_tests { let mm_tag_infos = MmTagInfo::parse_mm_tag(tag).unwrap(); let raw_mod_tags = RawModTags::new(tag, &quals, true); let obs_mod_base_info = - ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes()) + ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes(), &[]) .unwrap(); assert!(obs_mod_base_info.is_empty()); let tag = "C+h.;C+m.;G-h.;G-m.;"; let mm_tag_infos = MmTagInfo::parse_mm_tag(tag).unwrap(); let raw_mod_tags = RawModTags::new(tag, &quals, true); let obs_mod_base_info = - ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes()) + ModBaseInfo::new(&mm_tag_infos, &raw_mod_tags, dna.as_bytes(), &[]) .unwrap(); assert!(!obs_mod_base_info.is_empty()); assert_eq!(obs_mod_base_info.pos_seq_base_mod_probs.len(), 1); @@ -2828,7 +3160,7 @@ mod mod_bam_tests { record.set(b"test", None, dna.as_bytes(), &vec![255; dna.len()]); let raw_mm_tags = RawModTags::new(mm, &vec![255u16; 7], true); let modbase_info = - ModBaseInfo::new(&mm_tag_infos, &raw_mm_tags, dna.as_bytes()) + ModBaseInfo::new(&mm_tag_infos, &raw_mm_tags, dna.as_bytes(), &[]) .unwrap(); let t_converter = modbase_info.converters.get(&DnaBase::T).unwrap(); let expected = vec![0, 0, 0, 0, 0, 1, 2, 3, 3, 4, 4, 4, 4, 5, 6, 7, 7]; @@ -2847,7 +3179,7 @@ mod mod_bam_tests { record.set(b"test", None, dna.as_bytes(), &vec![255; dna.len()]); let raw_mm_tags = RawModTags::new(mm, &vec![255u16; 8], true); let modbase_info = - ModBaseInfo::new(&mm_tag_infos, &raw_mm_tags, dna.as_bytes()) + ModBaseInfo::new(&mm_tag_infos, &raw_mm_tags, dna.as_bytes(), &[]) .unwrap(); let t_converter = modbase_info.converters.get(&DnaBase::T).unwrap(); let expected = vec![0, 0, 0, 0, 0, 1, 2, 3, 3, 4, 4, 4, 4, 5, 6, 7, 7]; @@ -2867,14 +3199,31 @@ mod mod_bam_tests { let mut record = bam::Record::new(); record.set(b"test", None, dna.as_bytes(), &vec![255; dna.len()]); let raw_mm_tags = RawModTags::new(mm, &vec![255u16; 9], true); + // the explicit probabilities at the first C sum to more than one, the + // position is dropped but the record is still valid let modbase_info = - ModBaseInfo::new(&mm_tag_infos, &raw_mm_tags, dna.as_bytes()); - assert!(modbase_info.is_err()); + ModBaseInfo::new(&mm_tag_infos, &raw_mm_tags, dna.as_bytes(), &[]) + .unwrap(); + assert_eq!(modbase_info.n_conflict_positions, 1); + let c_probs = &modbase_info + .pos_seq_base_mod_probs + .get(&DnaBase::C) + .unwrap() + .pos_to_base_mod_probs; + // the only C with calls was the conflicting one + assert!(c_probs.is_empty()); + // the calls on the other bases are kept + let t_probs = &modbase_info + .pos_seq_base_mod_probs + .get(&DnaBase::T) + .unwrap() + .pos_to_base_mod_probs; + assert!(t_probs.contains_key(&5)); let mm = "C+m.;N+b?,1,3,0,0,1,3,0,0;"; let mm_tag_infos = MmTagInfo::parse_mm_tag(mm).unwrap(); let raw_mm_tags = RawModTags::new(mm, &vec![255u16; 8], true); let modbase_info = - ModBaseInfo::new(&mm_tag_infos, &raw_mm_tags, dna.as_bytes()); + ModBaseInfo::new(&mm_tag_infos, &raw_mm_tags, dna.as_bytes(), &[]); assert!(modbase_info.is_err()); // todo add test for specific error // if let Err(e) = modbase_info { @@ -2891,7 +3240,7 @@ mod mod_bam_tests { let mut record = bam::Record::new(); record.set(b"test", None, dna.as_bytes(), &vec![255; dna.len()]); let _modbase_info = - ModBaseInfo::new(&mm_tag_infos, &raw_mm_tags, dna.as_bytes()) + ModBaseInfo::new(&mm_tag_infos, &raw_mm_tags, dna.as_bytes(), &[]) .unwrap(); // let c_probs = modbase_info.pos_seq_base_mod_probs.get(&'C').unwrap(); // dbg!(&c_probs.pos_to_base_mod_probs); @@ -2906,7 +3255,7 @@ mod mod_bam_tests { let mut record = bam::Record::new(); record.set(b"test", None, dna.as_bytes(), &vec![255; dna.len()]); let parse_result = - ModBaseInfo::new(&mm_tag_infos, &raw_mm_tags, dna.as_bytes()); + ModBaseInfo::new(&mm_tag_infos, &raw_mm_tags, dna.as_bytes(), &[]); match parse_result { Err(MkError::Conflict(ConflictError::ExplicitConflictInferred)) => { diff --git a/modkit-core/src/modbam_util/check_tags.rs b/modkit-core/src/modbam_util/check_tags.rs index 358dc4af..e363bb05 100644 --- a/modkit-core/src/modbam_util/check_tags.rs +++ b/modkit-core/src/modbam_util/check_tags.rs @@ -333,7 +333,23 @@ fn extract_mm_tag_info(record: &bam::Record) -> TagState { match MmTagInfo::parse_mm_tag(&raw_mod_tags.raw_mm) { Ok(mm_tags) => { let forward_sequence = get_forward_sequence(record); - match ModBaseInfo::new(&mm_tags, &raw_mod_tags, &forward_sequence) { + match ModBaseInfo::new( + &mm_tags, + &raw_mod_tags, + &forward_sequence, + &[], + ) { + // check-tags is strict, a record with positions where the + // probabilities sum to more than one is reported as invalid + // even though the other commands will now use the record. + Ok(mod_base_info) if mod_base_info.n_conflict_positions > 0 => { + TagState::ValidTagsInvalidInfo { + mmtag_infos: mm_tags, + modbase_info_err: MkError::Conflict( + crate::errs::ConflictError::ProbaGreaterThanOne, + ), + } + } Ok(mod_base_info) => { TagState::Valid { mmtag_infos: mm_tags, mod_base_info } } diff --git a/modkit-core/src/modbam_util/subcommands.rs b/modkit-core/src/modbam_util/subcommands.rs index 509b270b..f912a5ea 100644 --- a/modkit-core/src/modbam_util/subcommands.rs +++ b/modkit-core/src/modbam_util/subcommands.rs @@ -2181,6 +2181,7 @@ impl Update { spinner.finish_and_clear(); info!("done, {} records processed", total); + crate::mod_bam::report_conflict_summary(); if !error_counts.is_empty() { info!("error/skip counts:"); diff --git a/modkit-core/src/pileup/base_mods_adapter.rs b/modkit-core/src/pileup/base_mods_adapter.rs index 3dbf3d5a..fd3f5f38 100644 --- a/modkit-core/src/pileup/base_mods_adapter.rs +++ b/modkit-core/src/pileup/base_mods_adapter.rs @@ -9,17 +9,25 @@ use rust_htslib::bam::{ }; use crate::{ - errs::MkResult, + errs::{MkError, MkResult}, + mod_bam::note_conflict_position, mod_base_code::{DnaBase, ModCodeRepr}, util::qual_to_prob, }; +/// Maximum allowed sum of the ML values at a single position, slightly more +/// than 256 to allow for rounding (see `crate::mod_bam::MAX_PROB`). +const MAX_TOTAL_MOD_QUAL: u16 = 258; + #[derive(Debug, Copy, Clone, new)] pub(crate) struct ModState { pub mod_position: usize, pub modified: bool, pub filtered: bool, pub mod_code: ModCodeRepr, + /// The base the modification is called on. When the modification is on + /// the opposite strand of the primary sequence base (`neg_strand`), this + /// is the complement of the base in the read. pub primary_base: DnaBase, pub inferred: bool, pub mod_qual: u8, @@ -34,7 +42,9 @@ pub(crate) struct BaseModsAdapter<'a, const SIZE: usize = 16> { canonical_bases: [u8; SIZE], mm_pos: [usize; SIZE], ml_pos: [usize; SIZE], - // strands: [u8; SIZE], // could be bitvec + /// true when the modification is on the opposite strand of the primary + /// sequence base (`-` in the MM tag) + neg_strands: [bool; SIZE], implicits: [bool; SIZE], ml_strides: [usize; SIZE], n_codes: usize, @@ -90,7 +100,7 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { let mut ml_start = 0usize; let mut mod_codes = [ModCodeRepr::Code('N'); SIZE]; let mut canonical_bases = [0u8; SIZE]; - // let mut strands = [0u8; SIZE]; + let mut neg_strands = [false; SIZE]; let mut ml_strides = [0usize; SIZE]; let mut implicits = [false; SIZE]; let mut mm_pos = [0usize; SIZE]; @@ -102,10 +112,11 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { assert!(n_codes < SIZE); let base = mm[i]; i += 1; - let strand = mm[i]; - if strand != b'+' { - bail!("duplex data not currently supported") - } + let neg_strand = match mm[i] { + b'+' => false, + b'-' => true, + _ => bail!("invalid strand in MM tag"), + }; i += 1; let (mods_in_rec, offset) = parse_mod_code(&mm[i..], &mut mod_codes, n_codes); @@ -175,7 +186,7 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { for j in 0..mods_in_rec { mm_pos[j + n_codes] = mm_idx; canonical_bases[j + n_codes] = base; - // strands[j + n_codes] = strand; + neg_strands[j + n_codes] = neg_strand; ml_strides[j + n_codes] = mods_in_rec; implicits[j + n_codes] = implicit_mode; mm_next[j + n_codes] = delta; @@ -194,10 +205,10 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { mod_codes, canonical_bases, n_codes, + neg_strands, implicits, mm_pos, ml_pos, - // strands, ml_strides, mm_next, num_explicit_positions, @@ -264,34 +275,86 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { } let mod_state = if let Some(mod_pos) = mod_pos { let mut mod_qual = 0u8; - let mut total_mod_qual = 0u8; + let mut total_mod_qual = 0u16; let base = if self.reverse { base_complement(self.seq[mod_pos]) } else { self.seq[mod_pos] }; - let mut mod_code = ModCodeRepr::Code(base as char); + let mut mod_code = None; + // strand of the calls at this position, None until a call (or + // implicit code) for this base is found + let mut neg_strand: Option = None; for i in 0..self.n_codes { if self.canonical_bases[i] == base && self.mm_next[i].map(|x| x == 0).unwrap_or(false) { let q = self.ml.get(self.ml_pos[i]).unwrap(); - if q > mod_qual { - mod_code = self.mod_codes[i]; + if q > mod_qual || mod_code.is_none() { + mod_code = Some(self.mod_codes[i]); mod_qual = q; } - total_mod_qual = total_mod_qual.saturating_add(q); + total_mod_qual += q as u16; + match neg_strand { + Some(s) if s != self.neg_strands[i] => { + return Err(MkError::InvalidMm( + "modification calls on both strands at the \ + same position are not supported" + .to_string(), + )) + } + _ => neg_strand = Some(self.neg_strands[i]), + } + } + } + if mod_code.is_none() { + // no explicit calls here, this is an inferred canonical + // position; take the strand from the implicit codes for this + // base. + for i in 0..self.n_codes { + if self.canonical_bases[i] == base && self.implicits[i] { + match neg_strand { + Some(s) if s != self.neg_strands[i] => { + return Err(MkError::InvalidMm( + "implicit modification calls on both \ + strands for the same base are not \ + supported" + .to_string(), + )) + } + _ => neg_strand = Some(self.neg_strands[i]), + } + } } } - let canonical_qual = 255u8.checked_sub(total_mod_qual).unwrap(); - let primary_base = DnaBase::parse(base as char).unwrap(); + let neg_strand = neg_strand.unwrap_or(false); + if total_mod_qual > MAX_TOTAL_MOD_QUAL { + // the probabilities at this position (from separate sub-tags) + // sum to more than one, e.g. independent 5mC and 5hmC models. + // Skip the position and continue with the next one. + note_conflict_position(); + self.move_forward(mod_pos, base); + return self + .next_modified_position(filter_thresholds, mod_thresholds); + } + let canonical_qual = 255u16.saturating_sub(total_mod_qual) as u8; + let primary_base = { + let b = DnaBase::parse(base as char).unwrap(); + if neg_strand { + b.complement() + } else { + b + } + }; + let mod_code = + mod_code.unwrap_or(ModCodeRepr::Code(primary_base.char())); let threshold = filter_thresholds[primary_base as usize]; let mod_state = if canonical_qual > mod_qual { Some(ModState::new( mod_pos, false, qual_to_prob(canonical_qual as i32) < threshold, - ModCodeRepr::Code(base as char), + ModCodeRepr::Code(primary_base.char()), primary_base, inferred, canonical_qual, @@ -329,9 +392,23 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> { Ok(mod_state) } + /// true when any of the modification codes are called on the opposite + /// strand of the primary sequence base (e.g. PacBio `T-a.` for 6mA). + pub fn has_negative_strand_mods(&self) -> bool { + self.neg_strands.iter().take(self.n_codes).any(|x| *x) + } + + /// Bit set of the bases modifications are called on (A, C, G, T in the + /// low 4 bits). For modifications on the opposite strand, the base is + /// the complement of the primary sequence base. pub fn primary_bases_in_record(&self) -> u8 { let mut bs = 0u8; - for &raw_can_base in self.canonical_bases.iter().take(self.n_codes) { + for i in 0..self.n_codes { + let raw_can_base = if self.neg_strands[i] { + base_complement(self.canonical_bases[i]) + } else { + self.canonical_bases[i] + }; match raw_can_base { b'A' => bs.view_bits_mut::().set(0, true), b'C' => bs.view_bits_mut::().set(1, true), diff --git a/modkit-core/src/pileup/pileup_processor.rs b/modkit-core/src/pileup/pileup_processor.rs index 262371a1..c39d5268 100644 --- a/modkit-core/src/pileup/pileup_processor.rs +++ b/modkit-core/src/pileup/pileup_processor.rs @@ -14,8 +14,8 @@ use crate::{ errs::{MkError, MkResult}, interval_chunks::{ChromCoordinates, FocusPositions2}, mod_bam::{ - validate_mn_tag_on_record, BaseModCall, BaseModProbs, EdgeFilter, - MmTagInfo, + note_conflict_position, probs_exceed_max, validate_mn_tag_on_record, + BaseModCall, BaseModProbs, EdgeFilter, MmTagInfo, }, mod_base_code::{ DnaBase, ModCodeRepr, ANY_CYTOSINE, HYDROXY_METHYL_CYTOSINE, @@ -206,6 +206,14 @@ impl< erred_records = erred_records.saturating_add(1); continue 'records; }; + if modbase_iter.has_negative_strand_mods() { + // modification calls on the opposite strand of the primary + // sequence base (e.g. PacBio `T-a.`) are not supported by the + // optimized workers, the pileup command detects these up + // front and uses the general workers instead. + erred_records = erred_records.saturating_add(1); + continue 'records; + } let Ok(mut mod_state) = modbase_iter.next_modified_position( self.filter_thresholds, @@ -641,7 +649,10 @@ impl PileupWorker for GenericPileupWorker { continue 'records; }; - let mm_tag_infos = MmTagInfo::from_record(&record)?; + let Ok(mm_tag_infos) = MmTagInfo::from_record(&record) else { + erred_records = erred_records.saturating_add(1); + continue 'records; + }; let (bases_to_codes, implicit_bases) = get_base_codes_and_implicits(&mm_tag_infos); @@ -803,7 +814,7 @@ impl PileupWorker for GenericPileupWorker { if q > pos { continue 'overran; } else { - let probs = codes + let probs: FxHashMap = codes .iter() .map(|hts_code| { let mod_code = ModCodeRepr::from( @@ -813,6 +824,12 @@ impl PileupWorker for GenericPileupWorker { (mod_code, prob) }) .collect(); + if probs_exceed_max(probs.values()) { + // probabilities from separate sub-tags + // sum to more than one, skip the position + note_conflict_position(); + continue 'overran; + } let Ok(can_base) = DnaBase::try_from(codes[0].canonical_base) else { @@ -2458,7 +2475,7 @@ fn update_mods_iter2<'a>( reverse: bool, caller: &MultipleThresholdModCaller, ) -> MkResult<()> { - if let Some(res) = modbase_iter.next() { + while let Some(res) = modbase_iter.next() { match res { Ok((pos, codes)) => { let pos = pos as usize; @@ -2473,7 +2490,7 @@ fn update_mods_iter2<'a>( Strand::Positive => can_base, Strand::Negative => can_base.complement(), }; - let base_mod_probs = codes + let base_mod_probs: FxHashMap = codes .iter() .map(|hts_code| { ( @@ -2482,6 +2499,12 @@ fn update_mods_iter2<'a>( ) }) .collect(); + if probs_exceed_max(base_mod_probs.values()) { + // probabilities from separate sub-tags sum to more than + // one, skip this position and move to the next one + note_conflict_position(); + continue; + } let base_mod_probs = BaseModProbs::new(base_mod_probs, false); *pos_base_mod_call = @@ -2489,7 +2512,8 @@ fn update_mods_iter2<'a>( *mod_pos = Some(pos); *canonical_base = Some(can_base); let strand = duplex_aware_strand(tag_mod_strand, reverse)?; - *mod_strand = Some(strand) + *mod_strand = Some(strand); + break; } Err(e) => { return Err(MkError::HtsLibError(e)); diff --git a/modkit-core/src/pileup/subcommand.rs b/modkit-core/src/pileup/subcommand.rs index 4503f99e..3d1c0f57 100644 --- a/modkit-core/src/pileup/subcommand.rs +++ b/modkit-core/src/pileup/subcommand.rs @@ -26,7 +26,7 @@ use crate::fasta::MotifLocationsLookup; use crate::interval_chunks::{ ChromCoordinatesFeeder, ReferenceIntervalBatchesFeeder, TotalLength, }; -use crate::mod_bam::{CollapseMethod, EdgeFilter}; +use crate::mod_bam::{CollapseMethod, EdgeFilter, MmTagInfo}; use crate::mod_base_code::{ DnaBase, ModCodeRepr, ModifiedBasesOptions, ANY_ADENINE, ANY_CYTOSINE, ANY_GUANINE, ANY_THYMINE, METHYL_CYTOSINE, SIX_METHYL_ADENINE, @@ -50,7 +50,8 @@ use crate::sample_probs::{ use crate::util::{ create_out_directory, filter_reference_records, get_master_progress_bar, get_master_progress_bar_fancy, get_subroutine_progress_bar, get_targets, - get_ticker, reader_is_bam, reader_is_cram, Region, + get_ticker, reader_is_bam, reader_is_cram, record_is_not_primary, + unindexed_reader_is_cram, Region, Strand, }; use crate::writers::{ BedMethylWriter, BedMethylWriter2, MultipleMotifBedmethylWriter, @@ -1023,6 +1024,26 @@ impl ModBamPileup { }; let (preset, regex_motifs) = self.determine_preset()?; + // The optimized workers do not support modification calls on the + // opposite strand of the primary sequence base (e.g. PacBio 6mA + // `T-a.` sub-tags), check the first records and fall back to the + // general workers when such calls are found. + let preset = if preset.is_some() + && bam_has_negative_strand_mods( + &self.in_bam, + self.reference_fasta.as_ref(), + 1_000, + )? { + info!( + "found base modification calls on the opposite strand of the \ + primary sequence base (e.g. PacBio 6mA 'T-a.'), these are \ + not supported by the optimized pileup workers, using general \ + workers" + ); + None + } else { + preset + }; let (pileup_options, combine_strands) = match &preset { Some(preset) => match preset { @@ -1623,6 +1644,7 @@ impl ModBamPileup { drop(records_tx); }); + let mut n_failed_intervals = 0usize; for result in records_rx.into_iter() { match result { Ok(mod_base_pileup) => { @@ -1633,7 +1655,10 @@ impl ModBamPileup { write_progress.inc(rows_written); } Err(message) => { - debug!("unexpected error {message}"); + n_failed_intervals += 1; + master_progress.suspend(|| { + error!("failed to process interval, {message}"); + }); } } } @@ -1643,7 +1668,15 @@ impl ModBamPileup { if n_failed_reads > 0 { master_progress.suspend(|| { - error!("~{n_failed_reads} failed processing"); + error!("~{n_failed_reads} records failed processing"); + }); + } + if n_failed_intervals > 0 { + master_progress.suspend(|| { + error!( + "{n_failed_intervals} interval(s) failed processing, \ + output is incomplete" + ); }); } @@ -1651,6 +1684,7 @@ impl ModBamPileup { erred_reads.finish_and_clear(); master_progress.suspend(|| { info!("Done, processed {rows_processed} rows."); + crate::mod_bam::report_conflict_summary(); }); tid_progress.finish_and_clear(); @@ -1664,6 +1698,42 @@ impl ModBamPileup { } } +/// Check the first `n_records` primary records of a modBAM for base +/// modification calls on the opposite strand of the primary sequence base +/// (`-` in the MM tag, e.g. PacBio `T-a.`). +fn bam_has_negative_strand_mods( + in_bam: &PathBuf, + reference_fasta: Option<&PathBuf>, + n_records: usize, +) -> anyhow::Result { + let mut reader = bam::Reader::from_path(in_bam)?; + if unindexed_reader_is_cram(&reader) { + if let Some(reference_fp) = reference_fasta { + reader.set_reference(reference_fp)?; + } else { + bail!("CRAM input requires reference") + } + } + let mut record = bam::Record::new(); + let mut n_checked = 0usize; + while let Some(result) = reader.read(&mut record) { + result?; + if record_is_not_primary(&record) { + continue; + } + if let Ok(mm_tag_infos) = MmTagInfo::from_record(&record) { + if mm_tag_infos.iter().any(|x| x.strand() == Strand::Negative) { + return Ok(true); + } + } + n_checked += 1; + if n_checked >= n_records { + break; + } + } + Ok(false) +} + #[derive(Clone, Debug)] enum Presets { /// CpG-special, combine strands, maybe combine mods @@ -2372,6 +2442,7 @@ impl DuplexModBamPileup { "Done, processed {rows_processed} rows. Processed \ ~{n_processed_reads} reads and skipped {n_skipped_message}." ); + crate::mod_bam::report_conflict_summary(); Ok(()) } } diff --git a/modkit-core/src/read_cache.rs b/modkit-core/src/read_cache.rs index 88f48d65..077a07db 100644 --- a/modkit-core/src/read_cache.rs +++ b/modkit-core/src/read_cache.rs @@ -110,7 +110,9 @@ impl<'a> ReadCache<'a> { fn add_record(&mut self, record: &bam::Record) -> MkResult<()> { let record_name = util::get_query_name_string(record)?; - let mod_base_info = ModBaseInfo::new_from_record(record)?; + let resolvers = self.method.map(std::slice::from_ref).unwrap_or(&[]); + let mod_base_info = + ModBaseInfo::new_from_record_with(record, resolvers)?; if mod_base_info.is_empty() { return Err(MkError::NoModifiedBaseInformation); } diff --git a/modkit-core/src/read_ids_to_base_mod_probs.rs b/modkit-core/src/read_ids_to_base_mod_probs.rs index 39513aa2..221a8939 100644 --- a/modkit-core/src/read_ids_to_base_mod_probs.rs +++ b/modkit-core/src/read_ids_to_base_mod_probs.rs @@ -284,8 +284,9 @@ impl RecordProcessor for ReadIdsToBaseModProbs { } else { None }; + let resolvers = collapse_method.cloned().into_iter().collect(); let mod_base_info_iter = records - .with_mod_base_info() + .with_mod_base_info(resolvers) .filter(|(record, _)| { if only_mapped || edge_filter.is_some() { !record.is_unmapped() @@ -991,8 +992,13 @@ impl RecordProcessor for ReadsBaseModProfile { cut: Option, kmer_size: Option, ) -> anyhow::Result { - let mut mod_iter = - TrackingModRecordIter::new(records, false, allow_non_primary); + let resolvers = collapse_method.cloned().into_iter().collect(); + let mut mod_iter = TrackingModRecordIter::new( + records, + false, + allow_non_primary, + resolvers, + ); let mut agg = Vec::new(); let mut seen = HashSet::new(); let pb = if with_progress { Some(get_ticker()) } else { None }; diff --git a/modkit-core/src/util.rs b/modkit-core/src/util.rs index 7323a1a7..1853a35f 100644 --- a/modkit-core/src/util.rs +++ b/modkit-core/src/util.rs @@ -866,6 +866,13 @@ pub(crate) fn reader_is_cram(reader: &bam::IndexedReader) -> bool { } } +pub(crate) fn unindexed_reader_is_cram(reader: &bam::Reader) -> bool { + unsafe { + (*reader.htsfile()).format.format + == rust_htslib::htslib::htsExactFormat_cram + } +} + pub(crate) const KMER_SIZE: usize = 50; #[derive(Copy, Clone)] diff --git a/modkit/tests/test_pacbio_style_tags.rs b/modkit/tests/test_pacbio_style_tags.rs new file mode 100644 index 00000000..5dae7259 --- /dev/null +++ b/modkit/tests/test_pacbio_style_tags.rs @@ -0,0 +1,222 @@ +//! Tests with PacBio Jasmine-style tags (`A+a.;C+h?;C+m?;T-a.`), see +//! tests/make_pacbio_style_tags.py. The 5mC and 5hmC probabilities at the +//! first C of each record sum to more than one, one record has no MM/ML tags +//! and 6mA is called on both strands. +use std::path::PathBuf; + +use common::run_modkit; +use rust_htslib::{bam, bam::Read}; + +mod common; + +const INPUT: &str = "../tests/resources/pacbio_style_tags.bam"; +const INPUT_NO_UNTAGGED: &str = + "../tests/resources/pacbio_style_tags_no_untagged.bam"; +const REFERENCE: &str = "../tests/resources/CGI_ladder_3.6kb_ref.fa"; + +fn count_records(fp: &PathBuf) -> usize { + let mut reader = bam::Reader::from_path(fp).unwrap(); + reader.records().count() +} + +fn read_lines(fp: &PathBuf) -> Vec { + std::fs::read_to_string(fp) + .unwrap() + .lines() + .map(|l| l.to_string()) + .collect() +} + +#[test] +fn test_adjust_mods_keeps_records_with_conflicting_probs() { + let out_bam = std::env::temp_dir().join("test_pacbio_style_adjust.bam"); + run_modkit(&[ + "adjust-mods", + "--ignore", + "h", + INPUT, + out_bam.to_str().unwrap(), + ]) + .unwrap(); + // 10 records in, the record without MM/ML tags is dropped, the 9 records + // with a position where P(5mC) + P(5hmC) > 1 are all kept. + assert_eq!(count_records(&out_bam), 9); + // the output has no conflicts left, check-tags exits 0 + run_modkit(&[ + "modbam", + "check-tags", + out_bam.to_str().unwrap(), + "--interval-size", + "20", + ]) + .unwrap(); + + // converting h to m also keeps every record + let out_bam = std::env::temp_dir().join("test_pacbio_style_convert.bam"); + run_modkit(&[ + "adjust-mods", + "--convert", + "h", + "m", + INPUT, + out_bam.to_str().unwrap(), + ]) + .unwrap(); + assert_eq!(count_records(&out_bam), 9); +} + +#[test] +fn test_check_tags_reports_conflicting_probs() { + // check-tags stays strict, the conflicting records make it exit 1 + assert!(run_modkit(&[ + "modbam", + "check-tags", + INPUT, + "--interval-size", + "20" + ]) + .is_err()); + + let tmp_dir = std::env::temp_dir().join("test_pacbio_style_check_tags"); + run_modkit(&[ + "modbam", + "check-tags", + INPUT, + "--interval-size", + "20", + "--permissive", + "--force", + "--out-dir", + tmp_dir.to_str().unwrap(), + ]) + .unwrap(); + let error_counts = read_lines(&tmp_dir.join("error_counts.tsv")); + let get_count = |error: &str| -> usize { + error_counts + .iter() + .find_map(|line| { + let parts = line.split('\t').collect::>(); + if parts[0] == error { + Some(parts[1].parse::().unwrap()) + } else { + None + } + }) + .unwrap_or_else(|| panic!("{error} not in {error_counts:?}")) + }; + assert_eq!(get_count("conflict-explicit-prob-greater-than-one"), 9); + assert_eq!(get_count("MM-tag-missing"), 1); +} + +#[test] +fn test_pileup_pacbio_style_tags() { + // The general workers, a record without MM/ML tags must not abort the + // interval. + let out_general = + std::env::temp_dir().join("test_pacbio_style_pileup_general.bed"); + run_modkit(&[ + "pileup", + INPUT, + out_general.to_str().unwrap(), + "--motif", + "CG", + "0", + "--ref", + REFERENCE, + "--no-filtering", + ]) + .unwrap(); + let general_lines = read_lines(&out_general); + assert!(!general_lines.is_empty()); + + // Same input without the record lacking MM/ML tags gives the same output. + let out_no_untagged = + std::env::temp_dir().join("test_pacbio_style_pileup_no_untagged.bed"); + run_modkit(&[ + "pileup", + INPUT_NO_UNTAGGED, + out_no_untagged.to_str().unwrap(), + "--motif", + "CG", + "0", + "--ref", + REFERENCE, + "--no-filtering", + ]) + .unwrap(); + assert_eq!(general_lines, read_lines(&out_no_untagged)); + + // Requesting a preset (--cpg --modified-bases) would use the optimized + // workers, which do not support the opposite-strand `T-a.` calls; pileup + // must detect this and fall back to the general workers. With + // --modified-bases 5mC the other modifications are counted as N_other, so + // only the 5mC and canonical counts are compared with the run above. + let run_preset = |input: &str, out_fn: &str| -> Vec { + let out_preset = std::env::temp_dir().join(out_fn); + run_modkit(&[ + "pileup", + input, + out_preset.to_str().unwrap(), + "--cpg", + "--ref", + REFERENCE, + "--modified-bases", + "5mC", + "--no-filtering", + ]) + .unwrap(); + read_lines(&out_preset) + }; + let preset_lines = run_preset(INPUT, "test_pacbio_style_pileup_preset.bed"); + assert!(!preset_lines.is_empty()); + let m_counts = |lines: &[String]| -> Vec> { + lines + .iter() + .map(|l| l.split('\t').map(|x| x.to_string()).collect::>()) + .filter(|parts| parts[3] == "m") + // chrom, start, strand, N_mod, N_canonical + .map(|parts| { + vec![ + parts[0].clone(), + parts[1].clone(), + parts[5].clone(), + parts[11].clone(), + parts[12].clone(), + ] + }) + .collect() + }; + assert!(!m_counts(&preset_lines).is_empty()); + assert_eq!(m_counts(&general_lines), m_counts(&preset_lines)); + // and the record without MM/ML tags does not change the output + assert_eq!( + preset_lines, + run_preset( + INPUT_NO_UNTAGGED, + "test_pacbio_style_pileup_preset_no_untagged.bed" + ) + ); +} + +#[test] +fn test_extract_pacbio_style_tags() { + let out_tsv = std::env::temp_dir().join("test_pacbio_style_extract.tsv"); + run_modkit(&[ + "extract", + "full", + INPUT, + out_tsv.to_str().unwrap(), + "--force", + ]) + .unwrap(); + let lines = read_lines(&out_tsv); + // header plus rows from all 9 tagged records + let mut read_ids = lines + .iter() + .skip(1) + .map(|l| l.split('\t').next().unwrap().to_string()) + .collect::>(); + read_ids.sort(); + read_ids.dedup(); + assert_eq!(read_ids.len(), 9); +} diff --git a/tests/make_pacbio_style_tags.py b/tests/make_pacbio_style_tags.py new file mode 100644 index 00000000..04ce3e8d --- /dev/null +++ b/tests/make_pacbio_style_tags.py @@ -0,0 +1,79 @@ +"""Create tests/resources/pacbio_style_tags{,_all_tagged}.bam from +tests/resources/bc_anchored_10_reads.sorted.bam. + +The records mimic PacBio HiFi reads processed with Jasmine >= 26.1.3: +`A+a.;C+h?;C+m?;T-a.` where + + * the sparse `C+h?` track comes from a model that is independent of the + `C+m?` model, so at a position the two probabilities can sum to more than + one (here the first C has m=255 and h=178), + * 6mA is called on both strands (`A+a.` and `T-a.`, implicit mode). + +In `pacbio_style_tags.bam` the last record has its MM/ML tags removed, +`pacbio_style_tags_no_untagged.bam` is the same file without that record. + +Requires pysam. Run from the repository root: + python tests/make_pacbio_style_tags.py +""" +import array + +import pysam + +SRC = "tests/resources/bc_anchored_10_reads.sorted.bam" +OUT = "tests/resources/pacbio_style_tags.bam" +OUT_NO_UNTAGGED = "tests/resources/pacbio_style_tags_no_untagged.bam" + + +def parse_mm(mm, ml): + subs = [s for s in mm.rstrip(";").split(";") if s] + parsed = {} + i = 0 + for s in subs: + parts = s.split(",") + head, deltas = parts[0], parts[1:] + parsed[head] = (deltas, list(ml[i:i + len(deltas)])) + i += len(deltas) + assert i == len(ml) + return parsed + + +def transform(rec): + parsed = parse_mm(rec.get_tag("MM"), rec.get_tag("ML")) + h_deltas, _h_ml = parsed["C+h?"] + m_deltas, m_ml = parsed["C+m?"] + assert h_deltas == m_deltas + # confident 5mC call on the first C ... + m_ml[0] = 255 + # ... and an independent, confident 5hmC call at the same position + h_deltas, h_ml = [h_deltas[0]], [178] + seq = rec.get_forward_sequence() + tracks = [] + if "A" in seq: + tracks.append(("A+a.", ["0"], [10])) + tracks.append(("C+h?", h_deltas, h_ml)) + tracks.append(("C+m?", m_deltas, m_ml)) + if "T" in seq: + tracks.append(("T-a.", ["0"], [10])) + new_mm = ";".join(f"{head}," + ",".join(d) for head, d, _ in tracks) + ";" + new_ml = [x for _, _, ml in tracks for x in ml] + rec.set_tag("MM", new_mm, "Z") + rec.set_tag("ML", array.array("B", new_ml)) + return rec + + +if __name__ == "__main__": + src = pysam.AlignmentFile(SRC) + records = [transform(rec) for rec in src] + n = len(records) + with pysam.AlignmentFile(OUT_NO_UNTAGGED, "wb", template=src) as fh: + for rec in records[:-1]: + fh.write(rec) + with pysam.AlignmentFile(OUT, "wb", template=src) as fh: + for i, rec in enumerate(records): + if i == n - 1: + rec.set_tag("MM", None) + rec.set_tag("ML", None) + fh.write(rec) + pysam.index(OUT) + pysam.index(OUT_NO_UNTAGGED) + print(f"wrote {OUT} ({n} records) and {OUT_NO_UNTAGGED} ({n - 1} records)") diff --git a/tests/resources/pacbio_style_tags.bam b/tests/resources/pacbio_style_tags.bam new file mode 100644 index 0000000000000000000000000000000000000000..76a85010979dd70988b891965d125d75bc7526c9 GIT binary patch literal 5345 zcmV<76dvmziwFb&00000{{{d;LjnM91C5l+j@vdAhM6u}XwYUB?W!vy+rk`PMCt+r z@W_Lv0X!bJo)kql1Z8PDMl8vYG(n55`vmQ>i#|lx-S-jl5?SO4LQ@xuj4~6z0?7Y6 zzDU?)IxW+iv`q6&I)U&^{nWlXn?HLlb}~uQQc%cW#zLv^Myl{6eYI7w zo4Z7YNbU1tT9?~QUnG<1HmTmrBo0I<0_$Jvu1UGu{*pHxX4`z5$vrvrNFZ^Yi*=~> z``^0QN@jIY?5acXalPvA&{M_jaJqm=Rh0E4-Uyj*Zi_Na1V*}&3dS-!nf;)zW-RcL z8^k!Kgu9IM7;_TeAj+ugre3P7hN`3r9);035xmTH-@?k)yhOJa$BTUCSs&JrB>E^CS5L$ z`|ofkvpA80PSyXn&cy%T^p30etotfOm6iuriVf z001A02m}BC000301^_}s0svJJwOV_uBiCJ?x!2FJ$DSFFeSQ6&@%Y|Jvc8?C9~;!J z?WybqHl$hqfv(Wq_(_o1CCxqrict0%OI$!e)ussX5S3CuAR$#CMG>t?f{0LPDWz(O z+WwOks;CsvD2i07mWJQ4$1`Wf_m=*%;~k&zJkFfC>+k!0-*^0byM(jUqE|T@+Gk;X z&tA~B(r8U)hSsR7T6VYH(4w^6Q*e$VcE3lu% zK3{>=f^{%BU5fI6JCn|Y0XN7l<>jC&oSs<=xziOa@w6*QgO1d-28+egVmlJo$%+eW zFkmgqy0RMxg~P#MFc1t;G#v^@BJd$XMNv!B$!sp4uawJzP_9%2p;9asinRmjK&*-f z4^_*RT&~C^X(q}f5`|ngyPvD%_T%AjAQ*v;V(~~M90~?Pkw`QgiNQkg8J}vbYA`17dyA_9e?=l=+>QI z_{AA<^Imsi!v9>5BuJ9*E%qLA)spCi-t66^-PotDUf89pvqL9_F~W%_1EvQpZV-po_^swR4lN*IcjP@ zag{@de;-38+E(;tV?H&T_33;*uN%{Lqux?Az1|*a`bcdwnsAtqU__tTT_f)m*zO`k z5^F@F&w`CiX@y9Dx&y>gRFa3Y3xMllpip>i1tppndmO-`WE0M^~WrfNmSW9yB=LWWozHl61cM4p7!x z+k&CZ}JB9l3l1nhz=gT z%0;357)A$ZTTv$MX>;%F=sIHL+P16*>H=IZ$DOmOk2cjre4@I#eNRoV@A0F4` zyM_Iq*`%N3WPuk*IvD>HZ{S1a}&Zej{SQheN zfff04;mxBG&x@?YuXQ)iTRg|hA{-I2-N7z^I~LYFCd_w-TBfM(4Su` zb5H&AV%xKa^kd0ip>XNDk-u2oB)u)(hl@=dZypRg)PhThixv_q-&WfTO3{G5O}|Xu?5;zMw&7*LHIYplh=(=%c!?zlTm#(n<-{PSB#d{UJ3uleQbeP zLbEIWN&yh={*5(f?epIat(amV$jdW%<*=!P4_xi-*HfM~1fB>9Uf%uq>g=o(jvs{hghiU9_ME976<%C5~!+ z5}5%Sn1e}YvcN5+Oc|C^Nj|7nYsZHNhldBUBnq&v25eC~C|87hng^B%A(9G)kY@;? z$N;QIAe1ZD4j&QV@?yCP5(KnHu|&uAQfZJ69Fs`K;Hn@+(OIEdtxAUv9ahCksj#07 zqfTH=h8ZpqgYGCQ!SJx`<;xWj{7X8M$?{B`qQem|Ca|!`_rz&B!6aie%U8w2YPBTP zvV0^Ar=UO~0oiny2)~hs2Ax>J=Q(it!70|YIley&iC==B=13jyQ+lY~QG4g)D(=(Z zWbNi2m{vNs;`m4t^HL;s4))h3v1h@Gg!A`t4dMTrsOYOq1%1>mtv{m|qFi;gsy@=) zUz2Y5tv{vc`+xkiaq-vwV)yR7&nKVzwz5mTZVR>uqF{2JihdQNV6?3&^GOrnSDz~l zFk9N(thcqM34(1hHT2d*osP|w+1gO=r*_uVyAf!g+6Xs{nm-v_r=p8f+2&h$?Y(4- zaYc(?0B>@OWrveG(n9u1GqJmryGv`y4wk^Nr-IOtVOMtv)Y@GvL;*UI8ES&%Q$U0n7+0w#R;3J^PdYa=>*>dr)S}K#uErUcCD908Lc^7#j``LIaE6h#Hd502$jCk15F6T+IRrG8UtQ^F@jmYeLt@@ zF6+i+%{*$VO_x`ncqxH+^=riodG)dJzan1!`fq)8i&tMJZC>@}I;5Xf9XNVz9qE@r z1dc&Jdt!>vh>tUIHuIV=mln222dBEo|=j} zn$K4O@hhYo9E;Gty_w?i>B#$b%MiY)AXo*oAFfAkmh**68-}yfB?y<0s+2|0i>%e< z<*sOT1RmmLR_<^e))FkCBM4Ve!~-)ykT3xUfz}JtLCBr6d%0YGKf|Y)Bts_>JeSI6 zxJ(LoEgFkbF<`h@JQ|9^oKZ9m6bO70Ph?Xp1N0Jyud?R<{oEGmsN@A=&H=fn^ly;p zvx|Q2@xsmHxk!S%C2sH~0^^5kSBw7LdmP7DJ)5{j&5f5F`)lMKu8cEqF_L-SJdq9xxhcO}XZ|57d#*wM1E~a0<_umN9Kb5}_)BR^&MVP+3_{S|w-?cO3 zA9KNAf>(O?OHaDa2ZIM*R7o6W1YUmea_yAliwiMHVp3xx21uMi@W93`&Gv&3q`y0o z@LcDLtF6bue3xVOHly|UK}V)IN@>GFK3+i0O=#?TuddBZYs<>eYieME|1nI3-+t%& z&Rh|KZY_iC3qyPq;|8=1UqbjQPJtHM`mAou+GDUUx>ldInytyq)JGG&O}uMYbcFR! z*gxt}Z0uk0<;?^fzq955uou=Mn_b{3k;4IK%GiS z{<{GFPz2Ilm@JCGdl(9K9N7{o2F_w{Kbz(`PAXq|_>n^t198zPOfPZiIL!5g&`(1U zmLcj1LowgW?=x_5<|xOd*|-RFk)#3jEC&q(Yho{7C}vB=Qn?7pQX;ZScvi9kB59Bc zkTRzL{t=YYw81&3&6A;~pA4avYpn+sD*q>5&3DhU5W0F~XjrKH>OZ}TkpFo1Lgat^ zo39|`Q*Ufn^Y?B2+=FzW;ef+U5wLiIbg+hl-bIUBX7Px{qs{S3i)#Ul2}rD*dvXK} zdDtCY$ho7$3N~(JHZt5T3-~NMLGXCAa zSscdVIz`#oeILRux5n-Q!!FuZ8tSY$ZjJ%Adb4iKm3DpH95w4R2!0#Nyw#rSQ{s~_ z&MQ0a53g8zu4uVfg4gZV3Y@pHUO9X zey88*b^2iVd%a#C!qrae0pK(OIGqRF$)D>x=TGy3>{c5tP+zNj1b}*V0(iR+s=xBk z=MYri>AblG)o&2%`5zLmBG{MOYp1tnIIxp^1{Y6FNiVv+Su2*?Vbv773xU5mwbOCL znAdiVGbRlPk|h5R(~<(+pC|wTABzYC000000RIL6LPG)o8vp|U0000000000v63pv literal 0 HcmV?d00001 diff --git a/tests/resources/pacbio_style_tags.bam.bai b/tests/resources/pacbio_style_tags.bam.bai new file mode 100644 index 0000000000000000000000000000000000000000..65c09ced096dfb888119283f8ec70b54eddca758 GIT binary patch literal 360 zcmZ>A^kh_GU|?VZVoxCk21X#wz>v=jrWlTjfM~Bt5P5WQE{GyHiK-VaJjwuu006Fc B1wQ}) literal 0 HcmV?d00001 diff --git a/tests/resources/pacbio_style_tags_no_untagged.bam b/tests/resources/pacbio_style_tags_no_untagged.bam new file mode 100644 index 0000000000000000000000000000000000000000..f41f7985fc060ef63d56231aa82bab08b9c055f6 GIT binary patch literal 5030 zcmV;X6ItvZiwFb&00000{{{d;LjnM91C5l+j@vdAhM6u}XwYUB?W!vy+rk`PMCt+r z@W_Lv0X!bJo)kql1Z8PDMl8vYG(n55`vmQ>i#|lx-S-jl5?SO4LQ@xuj4~6z0?7Y6 zzDU?)IxW+iv`q6&I)U&^{nWlXn?HLlb}~uQQc%cW#zLv^Myl{6eYI7w zo4Z7YNbU1tT9?~QUnG<1HmTmrBo0I<0_$Jvu1UGu{*pHxX4`z5$vrvrNFZ^Yi*=~> z``^0QN@jIY?5acXalPvA&{M_jaJqm=Rh0E4-Uyj*Zi_Na1V*}&3dS-!nf;)zW-RcL z8^k!Kgu9IM7;_TeAj+ugre3P7hN`3r9);035xmTH-@?k)yhOJa$BTUCSs&JrB>E^CS5L$ z`|ofkvpA80PSyXn&cy%T^p30etotfOm6iuriVf z001A02m}BC000301^_}s0stBiwOVV8PB5{PNuE)Y&sqrI`Kg?14$4(RG+$J*TlLm zLpz6kz6OT{$8dPMl;k0QRy`Ak{4leWm&2ZTdS)%;ZcntN)1FuzcB?&WxL7PLu3P21 z8EIh+hn!_u*LDISiV6mUfnbPcm=F~X!-p^(K{dl9GTB_NTq=oTsazJta-onf)b^|U zQbpQ-q*5wpvjslEun{&M&u254y=*zV7o(^^Fboq#W8p9r3I;;qaD)m+q0mBF5R#Ed zBo>WDsmE{h=^NIe+B_N==22hOkGe{~Gf?~5KpUvq(ZC>{92}Zvr>z*BWA#|qnughI z9;;35c&_NCZWy!HY|>PfX;W$6n%)`TdRA2$Mnk70qLqmlfrzq7c)2 z^WIY*f7i2jZoZtEn}@e&&;H76T3=U445)5Hbw$0Rs#g?vY3O)tRA2QWZu~;8ABDN-e|(4g*e5>s%nRS5qk;9!QC0id z>pZ&rCm1r(xuUn4^QqaYPv`S_-I#Wo^|q?%_0Cw+$7-|Lg3Ck|R`lVWHS%78^BzJZ zu|_2NEZWGdt`G^3JQEi3DL`={oc%aZc4b+%x|THL<*Kzj?ONTg*gKU(p(}_PerO>G zg~Na`mZ1YXP=cYIKp+$hhNuVyB|?Ql;Rq87phD-10M%kXUno>0u?**BQ7V-H!(2Qb zXIP427?$I>B*$`Po{w{pa4--6!0({x0I@Wk5TJJ+a9zlk^0}-a@X<(wrekoc5FHBf zNg;wdMghE00I%6MwWGe)J{k=A1FbnQ2deTAP#d{(a`Us4jaqW@bic6$T3x?t>J3#_ zT+n{wy~Mi#v@awt1a0MeFC);t{hM!Wf%g3?(0V|(4PXx%oL@5G1sX{@KYa%%>oanF zjU+txAw9%EI<2GwMjvW%xNfeykasTFQ`-x|D1gMCf^dOOk}g4Po{7Y$0RLtP>@i4Z zq`T&GzB%_UBuVt>^dKT z`U4mppmRl;bfz6;TvrVR=%5YJx^p~H>&M4tOI7C+W2TLX&p^iqHXZzt9hjf6=>UBm z!7&`3WrnhJR+Yu!Y4!B1cV;c6OlCQhmVzYoWI^JsC2v(N9tCEBlPrEHxhx<{B4-Ib z?#y}s7!5@tp;#;$3SPPtpn?!x!(oPDXn-#jk0)p-d-;7yk}3z1R2HjMd0#&`tjmXo zhqYQ&gs`1U(R2}FICQDi_YX_OeBtmg!6Xte8wCJoQt2eiL}Ajh{OBXqat^|IW^Xsi z#UnHVKjN1#P(du7OhoAz;GaVHADIf^-yAdt%1BcmQo#OgJYetm19s=C*}7^f4O8v7 z*dP7>XHmfZo$p_W{nLN@D8hdBy`R{^{x5E@@5Wag8opmMxw(iVVFObf5Su{fg4A=x zGaZH~cZ3sN{oo@!fb*=x?M6Cy+fWSYJ{{?y80RrD!Xik5@T`O+Jvg5H2fRAaRiO2pP#%2ohu>*^_wVUP0~mG_*LeiJKfthy&Xut~Z|mb}{kUy(Ae5`q zdS|A0>Pp9)1E{C)E1VGlSkzN+Jc&ZKEw;aCWA`arTq}0{jNQ&+*ye4DI)k#1hYKvo zp9^mRg?K^Ys=``#3xXx^f-J!mG1D9F0JteC3}G0edw>ST90EK>nOH2&MxZUoC)2rH zE?X+qsu0Y}CHYYhWb0UIxFG-gpVx{OK5h) zUnu~>ZQrEkoZbKJ(26M*qP#qlS2_Ah3oN#i0gFSaK(E>xf@W}XfhF*SWXV#mTI~U^ zEUUmI!lK%fhZ3apLeCQ5a@pzP*KWeM_^R63mz*ci=FVMt7%a47GIF-)9IL>W$~NC%ZlQLJTz zFa=$Zp%4e|(o~nZnn(S?0L%qt1iaM*3ez8wMEO4LW+ylQ;e1b_)97fJr|{q3))8|( zSH6(BK9GG0F_-b+^)2SoFRdiX`4q6#j7teH1psuekmsX*5B>iyZ_0=^IunX=+|t)7Dg0Iu3gd3V-$?eRb@VD0sQK7r3uNF zHaF`Xt!08>n@kP8JyEB}<|^6RQ13@Bt*Lh-&_1#eZWuNHRdDT!E>2~eZ{>CFC1Z>$ zS;7K%lV2>me6?F$$X;nC^_FsPX)U?o5;*o$6uUAs^_D=by~RQjVWO(sTS#4$pn-0q zXvok&t^w78qZ*EL99SpuSd5_|Q`=j4qbWg1#h{c5xeVlMvE5WjEELPTsZ=VRVxe1z zf*dP^#30JIkS+udHO>i1Ai^}Pt5}mN)ijsjd2nW62%BKw23(3`)2U>N<{vPyNwV=6n@lD%X&y`=HjZYD zAzDW32?QR$(N7?c)*P+u*^$;7sof@G)DcpI%E0J@CIn(_J_KHk0c&;gY_!yt%d1bm6i2-Jg~Elr`uOqxAYT2_?|puYS6?P=UiJDqq@PtCIC|4M(l3Ju z9D{!L#1x^+FK&y&D3~3D&9_;Ik27&L^Cn>~Eo_kvPMw=cU7~igI_wOmNBs64lgerF zsNWx3aOF%o1z!f?wawg+wI1!`IBm^wh|pAH#bqwb>F6C!J*391NlnQ9+4awT<-cD1 z;&6u=O~59tbCnMK}z`h7@47S<#;+%rs}M&9P3TclZsaFv|2{p9INW_)Kt{* ze7-V>Um;!NSc3Y^%@&VOM?SD?mhf!_!6Bk{7>}}9-e)dt7|wQ=z+6JAQkFn3a#l}} zdy>@^1+bSnxyyGsOSHtUC|*Mr59|bi!vq`zT90Ca;5%n_v)SBUT1c@8mWjs&KAB7N z=_K%4BpRWkz;Mx6Bou)?qeu)W5cnh(&m=h(=p_cuaXOW&6bm%cj`4ghzne`9Nsi~k z&_HofI!dQfNzi$FdwU=Xb5ILFk&z%`9tIgqpx{QTc{BiXni%Q)@*n;y@|)6sc_HcC{@d49Uh@zBbBlCT@&YmEfZP-MH%RoQ zi+=9$!p)O&kpzh~d99<3Xd|$x31SO-1imD%K*T&RoK!rD#FcCqFSvJeNBae)BO84% zZ?0k9T%o$XQY$T_V?5A~t)X3d{p0<2z4JqY&e#$qjZiIJ=aT@n4`HZA=Zc|d6J2kE zL{Zhcp|r;JPJ7&}tH%Im&6u_(>U^~kd(B3(Lb!ggYJZlntZ*Y;p0Zisft^iW1Bx*a znnioNsC#C~U11>%A)%A`D&`fVlu!2Nh^y_Q6`qlX0XXv+u^hA|MGzcg*gXoMBr67) z-0ofh_IUF7EWi%HjW&D4qFC9lfWSBaiwLgn<+G_ImqgwXY~fN!EJQQe1o%qe{vPt+ zMMr6nCb>*H%|SW_B@T%h_{oWIXf+0TPAm&N5(=+&e<8a|Bv_6YY>P=XM&{9|4?9SJ z@D`+IU^2l+?IA&WKw(1Sf6q_QnhjG^T};1r|CD5VVx`FH_Ux1qA*y}C|j+FPj%v!(_%_~lWOdh@OC zIeSG2y0tV)Us&S(7&oAE_z}WSaSF81(PwpI);WgsMc3-ncB?&^nfiF5cZj#`h>o!S zVf&9d6dU{3e10hk{_fo7-byaJrG@Q(Q~}x=1j9dX9&Rg(JS3%NH`mLa|f;XDJ?DIXo*_0hTmK z1#p>@0RJ#@X*!TOsI8-srXP(UmTT>YW-5OTkLEk)r4YK>Ff`0me)T_IMaX}mcOmjW z^p#f-^6A&NqxsKl{kaF}K*Irt>mp$G1nFQ62fd3HH_YM|i$|N|krvki7!#0KIQRGn z81k??x{z~Mi3Mytk=e*_HA^kh_GU|?VZVoxCk21X#wz>v=jrWm?}LA2K-h&;MDCqxmPMAZuy9%TSS000<0 B1jYaW literal 0 HcmV?d00001 From 6e0afa2551cce0c398d3462c4a265a047b0ccff1 Mon Sep 17 00:00:00 2001 From: ljwharbers Date: Fri, 4 Sep 2026 11:55:51 +0200 Subject: [PATCH 2/2] pileup: support --phased and --modified-bases in the general workers The general pileup workers (used for multiple motifs, --duplex and, since the previous commit, automatically for modBAMs with opposite-strand calls such as PacBio 6mA `T-a.`) ignored two options that the optimized workers honour: - `--phased` read no `HP` tag, so the hp1/hp2 outputs were empty and only the combined file was written. The worker now keeps three tallies (combined, HP=1, HP=2) like the optimized workers and fills `phased_feature_counts`; without `--phased` the output is unchanged. - `--modified-bases` did not restrict the output rows, every code present in the modBAM produced a row (an `h` row next to every `m` row on PacBio data). The requested codes are now passed to the worker and rows with other codes are dropped after tallying, so those calls still count in `N_other` of the remaining rows, as in the optimized workers. With `--combine-mods` all codes are summed as before. Tests: the PacBio-style fixtures now carry HP tags and a regular 5hmC call (tests/make_pacbio_style_tags.py); new integration tests cover the phased partition and the row filter, and a snapshot of the previous general-worker output guards the no-flag behaviour. Co-Authored-By: Claude Fable 5.1 --- CHANGELOG.md | 2 + modkit-core/src/pileup/pileup_processor.rs | 108 ++++++++--- modkit-core/src/pileup/subcommand.rs | 20 ++ modkit/tests/test_pacbio_style_tags.rs | 179 +++++++++++++++++- tests/make_pacbio_style_tags.py | 25 ++- tests/resources/pacbio_style_tags.bam | Bin 5345 -> 5371 bytes tests/resources/pacbio_style_tags.bam.bai | Bin 360 -> 360 bytes .../pacbio_style_tags_no_untagged.bam | Bin 5030 -> 5051 bytes .../pacbio_style_tags_no_untagged.bam.bai | Bin 360 -> 360 bytes .../pacbio_style_tags_pileup_general.bed | 36 ++++ 10 files changed, 338 insertions(+), 32 deletions(-) create mode 100644 tests/resources/pacbio_style_tags_pileup_general.bed diff --git a/CHANGELOG.md b/CHANGELOG.md index afda1224..8f47af1f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - [all] Records where the base modification probabilities from separate MM sub-tags sum to more than one at a position (e.g. PacBio Jasmine 5mC and 5hmC calls, which are made by independent models) are no longer discarded. The offending positions are dropped (or resolved with `--ignore`/`--convert` when given, e.g. `--ignore h` keeps the 5mC call), the rest of the record is used, and the number of dropped positions is reported. `modbam check-tags` still reports these records as `conflict-explicit-prob-greater-than-one`. - [pileup] A record without MM/ML tags no longer aborts the processing of the whole interval, which silently produced no output for that interval. - [pileup] Base modification calls on the opposite strand of the primary sequence base (e.g. PacBio 6mA `T-a.` calls) no longer make every record fail in the optimized workers and the threshold estimation; the general workers are used automatically when such calls are found. + - [pileup] `--phased` now partitions the counts on the `HP` tag in the general workers as well (used with multiple motifs, `--duplex` and for modBAMs with opposite-strand calls). Previously the `hp1` and `hp2` outputs were empty and only `combined` was written. + - [pileup] `--modified-bases` restricts the output rows to the requested modification codes in the general workers as well, matching the optimized workers; calls with other codes on the same primary base are counted in `N_other`. Previously every modification code present in the modBAM produced a row. ## [v0.6.4] ### Adds diff --git a/modkit-core/src/pileup/pileup_processor.rs b/modkit-core/src/pileup/pileup_processor.rs index c39d5268..6872a0c7 100644 --- a/modkit-core/src/pileup/pileup_processor.rs +++ b/modkit-core/src/pileup/pileup_processor.rs @@ -498,6 +498,13 @@ pub struct GenericPileupWorker { pileup_numeric_options: PileupNumericOptions, combine_strands: bool, max_depth: u16, + /// Partition the counts on the `HP` tag into hp1/hp2 in addition to the + /// combined counts. + phased: bool, + /// Modification codes requested with `--modified-bases`. When `Some` + /// only rows with these codes are emitted, calls with other codes on the + /// same primary base are still counted in `N_other`. + mod_codes: Option>, } impl GenericPileupWorker { @@ -511,6 +518,8 @@ impl GenericPileupWorker { pileup_numeric_options: PileupNumericOptions, combine_strands: bool, max_depth: u16, + phased: bool, + mod_codes: Option>, ) -> anyhow::Result { let mut reader = bam::IndexedReader::from_path(in_bam_fp)?; if reader_is_cram(&reader) { @@ -546,6 +555,9 @@ impl GenericPileupWorker { per_mod_thresholds, default_threshold, ); + let mod_codes = mod_codes.map(|codes| { + codes.into_iter().map(|(_, code)| code).collect::>() + }); Ok(Self { reader, motif_bases, @@ -553,6 +565,8 @@ impl GenericPileupWorker { pileup_numeric_options, combine_strands, max_depth, + phased, + mod_codes, }) } } @@ -581,20 +595,32 @@ impl PileupWorker for GenericPileupWorker { start_pos as i64, end_pos as i64, ))?; + let phased = self.phased; let records = self .reader .records() .filter_ok(|record| record.tid() >= 0i32) .filter_ok(|record| record_is_primary(record)) - // TODO: Capture errors here. Also check for partition-tag - .filter_map(|res| res.ok()); + // TODO: Capture errors here. + .filter_map(|res| res.ok()) + .map(|record| { + // 0 means the record is not haplotagged + let hp = if phased { + get_haplotype_tag(&record, &HPTAG).unwrap_or(0u8) + } else { + 0u8 + }; + (record, hp) + }); - let mut chrom_features: FxHashMap = - FxHashMap::default(); + // [combined, hp1, hp2], hp1 and hp2 stay empty unless `phased` + let mut chrom_features: [FxHashMap; 3] = + [FxHashMap::default(), FxHashMap::default(), FxHashMap::default()]; let mut add_to_tally = |position_strand: PositionStrand, call: Call, motif_info: Option<&MotifInfo>, - motif_idxs: u8| { + motif_idxs: u8, + hp: u8| { let call = match motif_info.map(|x| x.primary_base) { Some(x) if call.matches_dna_base(&x) => call, Some(_x) => match call { @@ -627,7 +653,13 @@ impl PileupWorker for GenericPileupWorker { _ => position_strand, }; - chrom_features + if phased && (hp == 1 || hp == 2) { + chrom_features[hp as usize] + .entry(position_strand) + .or_insert_with(Tally2::default) + .add_call(call.clone(), motif_idxs, self.max_depth); + } + chrom_features[0] .entry(position_strand) .or_insert_with(Tally2::default) .add_call(call, motif_idxs, self.max_depth); @@ -637,7 +669,7 @@ impl PileupWorker for GenericPileupWorker { let mut canonical_base = Option::::None; let mut pos_base_mod_call = Option::::None; let mut mod_strand = Option::::None; - 'records: for record in records { + 'records: for (record, hp) in records { let reverse = record.is_reverse(); let record_strand = if record.is_reverse() { Strand::Negative @@ -749,6 +781,7 @@ impl PileupWorker for GenericPileupWorker { call, ref_base, motif_idxs, + hp, ); let Ok(_) = update_mods_iter2( @@ -790,6 +823,7 @@ impl PileupWorker for GenericPileupWorker { }, ref_base, motif_idxs, + hp, ); } else { add_to_tally( @@ -797,6 +831,7 @@ impl PileupWorker for GenericPileupWorker { Call::NoCall(base), ref_base, motif_idxs, + hp, ); } } @@ -869,6 +904,7 @@ impl PileupWorker for GenericPileupWorker { ), ref_base, motif_idxs, + hp, ); mod_pos = Some(pos); canonical_base = Some(can_base); @@ -899,6 +935,7 @@ impl PileupWorker for GenericPileupWorker { }, ref_base, motif_idxs, + hp, ); } else { add_to_tally( @@ -906,6 +943,7 @@ impl PileupWorker for GenericPileupWorker { Call::NoCall(base), ref_base, motif_idxs, + hp, ); } mod_pos = Some(pos); @@ -924,6 +962,7 @@ impl PileupWorker for GenericPileupWorker { Call::Delete, ref_base, motif_idxs, + hp, ); } (Some(q), None) => { @@ -948,6 +987,7 @@ impl PileupWorker for GenericPileupWorker { }, ref_base, motif_idxs, + hp, ); } else { add_to_tally( @@ -955,6 +995,7 @@ impl PileupWorker for GenericPileupWorker { Call::NoCall(base), ref_base, motif_idxs, + hp, ); } } @@ -967,22 +1008,43 @@ impl PileupWorker for GenericPileupWorker { PileupNumericOptions::Combine => true, _ => false, }; - let position_feature_counts = chrom_features - .into_iter() - .sorted_by(|((x, a), _), ((y, b), _)| match x.cmp(y) { - Ordering::Equal => a.cmp(b), - o @ _ => o, - }) - .flat_map(|((ref_pos, strand), tally)| { - tally.into_counts( - start_pos, - ref_pos, - combine_mods, - if self.combine_strands { '.' } else { strand.to_char() }, - ) - }) - .filter(|x| x.is_valid()) - .collect::>(); + // with --combine-mods all codes are summed into one row per primary + // base, so there is nothing to filter + let mod_codes = + if combine_mods { None } else { self.mod_codes.as_ref() }; + let into_feature_counts = + |features: FxHashMap| { + features + .into_iter() + .sorted_by(|((x, a), _), ((y, b), _)| match x.cmp(y) { + Ordering::Equal => a.cmp(b), + o @ _ => o, + }) + .flat_map(|((ref_pos, strand), tally)| { + tally.into_counts( + start_pos, + ref_pos, + combine_mods, + if self.combine_strands { + '.' + } else { + strand.to_char() + }, + ) + }) + .filter(|x| x.is_valid()) + .filter(|x| { + mod_codes + .map_or(true, |codes| codes.contains(&x.mod_code)) + }) + .collect::>() + }; + let [combined, hp1, hp2] = chrom_features; + let position_feature_counts = into_feature_counts(combined); + if phased { + pileup_space.phased_feature_counts = + [into_feature_counts(hp1), into_feature_counts(hp2)]; + } pileup_space.chrom_name = chrom_name; pileup_space.interval_width = width; diff --git a/modkit-core/src/pileup/subcommand.rs b/modkit-core/src/pileup/subcommand.rs index 3d1c0f57..ac4070fb 100644 --- a/modkit-core/src/pileup/subcommand.rs +++ b/modkit-core/src/pileup/subcommand.rs @@ -539,6 +539,24 @@ impl ModBamPileup { Ok(()) } + /// The modification codes requested with `--modified-bases`, used by the + /// general workers to restrict the output rows to these codes (the + /// optimized workers do this through their count matrices). `None` when + /// no codes were requested or when `--combine-mods` sums all codes + /// together anyway. + fn requested_mod_codes(&self) -> Option> { + if self.combine_mods { + return None; + } + self.modified_bases.as_ref().map(|modified_bases| { + modified_bases + .iter() + .map(|x| (x.primary_base, x.mod_code)) + .sorted() + .collect() + }) + } + fn determine_preset( &self, ) -> anyhow::Result<(Option, Option>)> { @@ -1534,6 +1552,8 @@ impl ModBamPileup { pileup_options.clone(), self.combine_strands, self.max_depth, + self.phased, + self.requested_mod_codes(), ) }) .collect::>>()?; diff --git a/modkit/tests/test_pacbio_style_tags.rs b/modkit/tests/test_pacbio_style_tags.rs index 5dae7259..186f3833 100644 --- a/modkit/tests/test_pacbio_style_tags.rs +++ b/modkit/tests/test_pacbio_style_tags.rs @@ -1,8 +1,9 @@ //! Tests with PacBio Jasmine-style tags (`A+a.;C+h?;C+m?;T-a.`), see //! tests/make_pacbio_style_tags.py. The 5mC and 5hmC probabilities at the -//! first C of each record sum to more than one, one record has no MM/ML tags -//! and 6mA is called on both strands. -use std::path::PathBuf; +//! first C of each record sum to more than one, the second C carries a +//! regular 5hmC call, one record has no MM/ML tags, 6mA is called on both +//! strands and the records are haplotagged (HP=1, HP=2 or untagged). +use std::{collections::HashMap, path::PathBuf}; use common::run_modkit; use rust_htslib::{bam, bam::Read}; @@ -27,6 +28,21 @@ fn read_lines(fp: &PathBuf) -> Vec { .collect() } +/// bedMethyl rows split into columns: [3] mod code, [5] strand, +/// [9] N_valid_cov, [11] N_mod, [12] N_canonical, [13] N_other_mod. +fn parse_rows(lines: &[String]) -> Vec> { + lines + .iter() + .map(|l| l.split('\t').map(|x| x.to_string()).collect()) + .collect() +} + +type RowKey = (String, String, String); + +fn row_key(row: &[String]) -> RowKey { + (row[0].clone(), row[1].clone(), row[5].clone()) +} + #[test] fn test_adjust_mods_keeps_records_with_conflicting_probs() { let out_bam = std::env::temp_dir().join("test_pacbio_style_adjust.bam"); @@ -128,6 +144,14 @@ fn test_pileup_pacbio_style_tags() { .unwrap(); let general_lines = read_lines(&out_general); assert!(!general_lines.is_empty()); + // regression: identical to the output of the general workers before the + // --phased and --modified-bases changes (snapshot made with 697de7b) + assert_eq!( + general_lines, + read_lines(&PathBuf::from( + "../tests/resources/pacbio_style_tags_pileup_general.bed" + )) + ); // Same input without the record lacking MM/ML tags gives the same output. let out_no_untagged = @@ -220,3 +244,152 @@ fn test_extract_pacbio_style_tags() { read_ids.dedup(); assert_eq!(read_ids.len(), 9); } + +#[test] +fn test_pileup_pacbio_style_modified_bases_filters_rows() { + // The general workers (opposite-strand 6mA calls). Without + // --modified-bases every code in the BAM gets a row, with + // --modified-bases 5mC only the m rows are written and the h calls are + // counted in N_other of the m rows. + let run = |extra: &[&str], out_fn: &str| -> Vec> { + let out = std::env::temp_dir().join(out_fn); + let mut args = vec![ + "pileup", + INPUT, + out.to_str().unwrap(), + "--motif", + "CG", + "0", + "--ref", + REFERENCE, + "--no-filtering", + ]; + args.extend_from_slice(extra); + run_modkit(&args).unwrap(); + parse_rows(&read_lines(&out)) + }; + let all_rows = run(&[], "test_pacbio_style_pileup_all_codes.bed"); + let m_rows = run( + &["--modified-bases", "5mC"], + "test_pacbio_style_pileup_5mc_only.bed", + ); + assert!(all_rows.iter().any(|r| r[3] == "h")); + assert!(all_rows.iter().any(|r| r[3] == "m")); + assert!(!m_rows.is_empty()); + assert!(m_rows.iter().all(|r| r[3] == "m")); + // the m rows themselves are not changed by the filter + let expected_m_rows = all_rows + .iter() + .filter(|r| r[3] == "m") + .cloned() + .collect::>>(); + assert_eq!(m_rows, expected_m_rows); + // N_other of every m row is N_mod of the h row at the same position + let h_n_mod = all_rows + .iter() + .filter(|r| r[3] == "h") + .map(|r| (row_key(r), r[11].parse::().unwrap())) + .collect::>(); + let mut n_other_total = 0usize; + for row in &m_rows { + let n_other = row[13].parse::().unwrap(); + assert_eq!(Some(&n_other), h_n_mod.get(&row_key(row)), "{row:?}"); + n_other_total += n_other; + } + // the fixture has a regular 5hmC call on the second C of each record + assert!(n_other_total > 0); + + // --combine-mods sums all codes into one row per position instead. (With + // --combine-mods a single-base `C` motif is added next to `CG`, so the + // rows carry multiple-motif labels; only the CG rows are compared.) + let combined_rows = run( + &["--modified-bases", "5mC", "--combine-mods"], + "test_pacbio_style_pileup_combine_mods.bed", + ) + .into_iter() + .filter(|r| r[3] == "C,CG,0") + .collect::>>(); + assert_eq!(combined_rows.len(), m_rows.len()); + for (combined, m) in combined_rows.iter().zip(m_rows.iter()) { + assert_eq!(row_key(combined), row_key(m)); + let n_mod_combined = combined[11].parse::().unwrap(); + let n_mod = m[11].parse::().unwrap(); + let n_other = m[13].parse::().unwrap(); + assert_eq!(n_mod_combined, n_mod + n_other); + } +} + +#[test] +fn test_pileup_pacbio_style_phased() { + // The general workers (opposite-strand 6mA calls) with --phased: the + // counts are partitioned on the HP tag, untagged records only count + // towards the combined output. + let out_dir = std::env::temp_dir().join("test_pacbio_style_pileup_phased"); + run_modkit(&[ + "pileup", + INPUT, + out_dir.to_str().unwrap(), + "--cpg", + "--ref", + REFERENCE, + "--modified-bases", + "5mC", + "--no-filtering", + "--phased", + "--prefix", + "pb", + ]) + .unwrap(); + let combined = + parse_rows(&read_lines(&out_dir.join("pb_combined.bedmethyl"))); + let hp1 = parse_rows(&read_lines(&out_dir.join("pb_hp1.bedmethyl"))); + let hp2 = parse_rows(&read_lines(&out_dir.join("pb_hp2.bedmethyl"))); + assert!(!hp1.is_empty()); + assert!(!hp2.is_empty()); + assert_ne!(hp1, hp2); + for rows in [&combined, &hp1, &hp2] { + assert!(rows.iter().all(|r| r[3] == "m")); + } + + // the combined output is the same as without --phased + let out_unphased = + std::env::temp_dir().join("test_pacbio_style_pileup_unphased.bed"); + run_modkit(&[ + "pileup", + INPUT, + out_unphased.to_str().unwrap(), + "--cpg", + "--ref", + REFERENCE, + "--modified-bases", + "5mC", + "--no-filtering", + ]) + .unwrap(); + assert_eq!(combined, parse_rows(&read_lines(&out_unphased))); + + // per position the haplotype coverages sum to at most the combined + // coverage, and in total to less because of the untagged records + let coverage = |rows: &[Vec]| -> HashMap { + rows.iter() + .map(|r| (row_key(r), r[9].parse::().unwrap())) + .collect() + }; + let combined_cov = coverage(&combined); + let hp1_cov = coverage(&hp1); + let hp2_cov = coverage(&hp2); + for key in hp1_cov.keys().chain(hp2_cov.keys()) { + assert!(combined_cov.contains_key(key), "{key:?}"); + } + let mut total_haplotagged = 0usize; + let mut total_combined = 0usize; + for (key, cov) in &combined_cov { + let haplotagged = hp1_cov.get(key).copied().unwrap_or(0) + + hp2_cov.get(key).copied().unwrap_or(0); + assert!(haplotagged <= *cov, "{key:?}"); + total_haplotagged += haplotagged; + total_combined += cov; + } + assert!(total_haplotagged > 0); + assert!(total_haplotagged < total_combined); +} diff --git a/tests/make_pacbio_style_tags.py b/tests/make_pacbio_style_tags.py index 04ce3e8d..c701dfe6 100644 --- a/tests/make_pacbio_style_tags.py +++ b/tests/make_pacbio_style_tags.py @@ -6,8 +6,12 @@ * the sparse `C+h?` track comes from a model that is independent of the `C+m?` model, so at a position the two probabilities can sum to more than - one (here the first C has m=255 and h=178), - * 6mA is called on both strands (`A+a.` and `T-a.`, implicit mode). + one (here the first C has m=255 and h=178; the second C has a regular + 5hmC call, m=10 and h=230), + * 6mA is called on both strands (`A+a.` and `T-a.`, implicit mode), + * the records are haplotagged like Longphase/WhatsHap output: record i gets + `HP:i:1` when i % 3 == 0, `HP:i:2` when i % 3 == 1 and stays untagged + otherwise (all tagged records share `PS:i:1`). In `pacbio_style_tags.bam` the last record has its MM/ML tags removed, `pacbio_style_tags_no_untagged.bam` is the same file without that record. @@ -37,15 +41,24 @@ def parse_mm(mm, ml): return parsed -def transform(rec): +def transform(i, rec): + if i % 3 == 0: + rec.set_tag("HP", 1, "i") + rec.set_tag("PS", 1, "i") + elif i % 3 == 1: + rec.set_tag("HP", 2, "i") + rec.set_tag("PS", 1, "i") parsed = parse_mm(rec.get_tag("MM"), rec.get_tag("ML")) h_deltas, _h_ml = parsed["C+h?"] m_deltas, m_ml = parsed["C+m?"] assert h_deltas == m_deltas + assert len(m_deltas) >= 2 # confident 5mC call on the first C ... m_ml[0] = 255 - # ... and an independent, confident 5hmC call at the same position - h_deltas, h_ml = [h_deltas[0]], [178] + # ... and an independent, confident 5hmC call at the same position, plus + # a non-conflicting confident 5hmC call on the second C (m low, h high) + m_ml[1] = 10 + h_deltas, h_ml = h_deltas[:2], [178, 230] seq = rec.get_forward_sequence() tracks = [] if "A" in seq: @@ -63,7 +76,7 @@ def transform(rec): if __name__ == "__main__": src = pysam.AlignmentFile(SRC) - records = [transform(rec) for rec in src] + records = [transform(i, rec) for i, rec in enumerate(src)] n = len(records) with pysam.AlignmentFile(OUT_NO_UNTAGGED, "wb", template=src) as fh: for rec in records[:-1]: diff --git a/tests/resources/pacbio_style_tags.bam b/tests/resources/pacbio_style_tags.bam index 76a85010979dd70988b891965d125d75bc7526c9..07be36eeb1c1df2cb84714c24116499fc1d0c62b 100644 GIT binary patch delta 4439 zcmV-d5vcCrDf=m~e**z;lYj#if3feqlkEC-o_=hqc5P2(FR)>o^&jX8-Ho5561$|? zhd>cxud&1h1Vn8rAs(Vq3Kd956-ZG^*wt*+j65dndxeyf3B$QSzT{6 z<@$K4>2;+snKoOq>1d>D#Pg^Jk{~!}zO-ex#Hy{peir+D8CDb4!Qf;mN(1h+b}9_G zL3SxE2VLRh)LclNu3(BMU7~@{OmB^DJfp~m`k``T);d&g=(0S1`u2^dMw83; zv_s=iB|dtSK?M|lA#dL~Qe{=G%dNW78Y=oFrEv-MqZ=1w`I0OX)4TKTldpZ?nOoOi z&d!aao3m$rVK%L=n!~4mX=)I+Z_cOBey)3aIQ&{Njwbd9OP$ z;D0Vi5+q5iE%qMLYDx4$Z}x7|ZtPR77j~(27U`bmoZ(!5wCc&(Y4=KEH}>vf-_&YJ zZA>RcEdUX&l} zU3lt)7xjw@Y}~XqE-K1J8D1({9~;)Kvxuud)9uCKLh7&n^F8WwU-;_NFMNk823Dg% z&EEU3a4~3qe;?lDDM)uO z>7sKD#yEbR;~odN2kfrx0b!g#VoyPM$R|mUF`Q>2F>b*BG6eP*q6|hKM-~@y_4=7Ozt;*|j``2;#9TTol@mV+4V= zWqHznp0?#tUD0JAh!zy|_IRSy$77?ZDD#OvQ%A&?;l%qkK|E*I?rSzdpwB~C2ZPh> zKoU=Dk}x=_ot$=0&83*lE(hY07x}Kli=4US%$mue+Rd|~$qhu07bH<&O}@*Xn)d;< zp=dM|kH(#p;15R)ktj+Sd+_>eA8C(cWH7 zfXbdvQ&b7cJe;yu-#;ps3dN(NB%Mq`Z()EuozA2fItHCqq(>jARq{~rv%5PfHW8%| z9U#Vu06D}Hsbq|b10#eHBMc207{TZ_`tnd!9#T#J4ZNmruh;bUC8K%CkPi)|?Xkgs z@S9&jHT^e#a4s91{QXOa4RRm;!Uh}stiuLg3C6*)_WdTEgE%%E;NlQ*YM}>P-w|JT z7{5$&jkv`IK9}}Z6f+PsvNJ!F$T>FXeaXMgz*BW>= zsN@!QVZ%eOd{H~QzQ6dx|7(8lFR1r_2Qc_1uW$%`zl*^aZObEV-qJ?X`naXHp}s5A zdV8j|>vG$e1Fom=MVt`<7}t}qK7s1EEy>@s;rpa5!DYLB+TPD#@aAmZI)!B+4Hj7C zKNI0Rs{6di*7%h==XsOocu9mKLbf~D27HIZ5va{j@&gp;>JUIOM#tj`CJOt1f?O(- z&*yXH@?H%pe5EWs3bMUYu8Or{QDEpeolJ*9bRtGkaA_F1 zK)j*aGQg4rgty353O_XE0gh7`nCY0EC2V8J55tqA7{p zTCEFAv#bHP@QYek8i){o1@K*yhr_jw1Z}s|Woy#HG=(mEA}j~@x3;#n(SjN<77?JD zI1&I!hEXT5Cc_Muh(UK0m0);ScJk$l2sS63$z*vZPSN2A_!n4MWRK!BonVqNn&qqF zLA6>E_Og5=45y$VA^~isDjt_Lj(hz+gc9-)n5+S$radIx@&n2@1>v*jBMD z;-9~(AufBiaxRyDeLVMB#AW*ZH#WFT+p=`a*`%^#p$xRi#g8W+<( zXM*vJXJSZ#I7{?k6hL^4gwqp#BOfukam1Q6!3hZ`UC#}z{aHx-68tn*|M;KMNA|7+ zJR^_so(89HFE7Ca)R{GCjYcucMPlb*e{B+b7OZTze;?0(HMnt6yNiV= zKu0yHyAV6b#(}g8Q{dx3(E;~DG#g288|Hz1}7FmdJ(xYBpo4)O|X0lm@xx?<0|cm)mny4avUT&a0!#7;RbA) zWishhnqf1AQYBv~h#x9JZd5K7Dxx4(_eB9nKFx8dY`J()EtN|Ldpx8_EabXWG@4=( zaVC{YWiuRvNlXIu7)Lyf#uErUezliGVXrZ?a^9iZ94eg#qSzr4jdEY_foue(Z9D{` zjRDbrz8Jv_;_f=pHZE!UCDk}=Dou}QpZIJ75$)HD=MwGX<9|U!`|rQ~wGEU8Qc zx99z^Eu1g$+xx;TcaKjXKXAdtHR&$+FbIEt+~(D=)zxb>r#du;I8D(lI`fE5L+`NZ z<2BCVwcHBTp{_YN7PzG+8HX~buhD_k0kp!L z&W@oVg+fYG5h0F3QCG%*Z&jZ*C(7KK<-KgfTPD0DS+PHh zn00jI%#$`P+_x>~W#Akg!CSC@r=B{erqtmV`~X~^#MQ8r8TmrVZzQd$Ckx<6*8sk` ztWXKsG9w6CGt5K`gCsf`30c>msET2KfLKBiYib!R zU^2=A|+ywH+3?g+dj<0afJ)8`jwlEmbP zBM?ZOQt+I|E7$ge5afD)cO>Dv&J$c4kAs;r*AF_Q^_fOjy144ep`{1Y^R1ECL&Tcl%b6l-_ z@zDzhC_Lh#QJAUX(s7tJ3ZWmOAaFxK6ozuYliy|F;>=-=OS5qi_#+8(PYITThJiJ) zlP?srrDCaEggh#L5wUV1OL#!E4PpY)>J;!n1Z6vIuop`6c&KW}L#Xy@>!F3y|Ap84 z?X!G_rX1=j7EZtVkFO#sxYj+F3O@0zR}d9YZ*124_ig;RgmkguqQ%P|uuy|^!H2Ve z2OqBx^ zK(5xT>vOqXA2&zM`V11|hCFYzr`nYG6b$w9miJ39)}D6Svg2fKO#?hn*fkwWdd;jY zPn=mL?AIiJmSl%0coj@Q&;UUOw}hFd6OQiyZnIFlLlnRc>iIT650!GjP%5>vlSv?% zp3P>!&?hs<=f^lEnTW+0hK*8DDhGytw^qsWFfmM5j_uI(-NRdc9sBQrS-HAwV?(P@M<7 zd8F%GXHOb~Hdh)RR$r?;2Uxu_0mz+;)?az#b%fS;I&W^E^&7-$Dv87^3--PC+wP6Y z66~ZkBaFwbq#xk^kQNL4unLRag}~q3dg?l7%#(Wt8CuI;Jp1vaEv|ADU1-S|n6*d@(A z1d34h8cSS2K-H!Q@eq|#Kp-JiAVm?aNP>t^Xep&?iQ4{?7OJQe(I|>ks+NY|vBxuK z#`l)~v*R6~@jT9)x$E!yecyNdd%J|Q)S_288ro-Jea~Ldw$f-#W`@?Nf2&%1RySHr zr9PhOMqO=8rp?xDIvN=|@jU8*BnS?gFYMSYv2H7{pT$03fz^U_FgRU`@_;*&&V&Is z$S&pOpevl7Sqr(-6)f?zD@cQm)U^hS#nNIs64%L!3u`c7Ez7#H8wiEN!C){D3{f;4 z3P&RFAwoq_OVi10E}yTIC(D9Ru2ck}QY;pVwFBuutcnK@Rm+uJuE-^6CdwodgvnyPkO3;(ouTtzS`6JT3^*p`Udet|JXF!EyZY$)v>NM z471S~s|{^DS9DW1j9GIwX{gGyp|oyHZ;fs|qbf(nk$Pj+I?`?!iZXus_Kl}TlPmZ1 zBlAcjK76}J-MW1}6uNV}N8J47USRUn^}i)=-#OM4e@&|^t-9J8s>WrtaT)bvn2#vR zWkn&Tcjw(FUwi*Ex30gOotwuuXV3inY+7G8M^FFK)Ff`-oKK(qO!xM5Uj5n^JGbv0 zfB5d`)}3GY#Tjw)UUy=`|6GtHNRsd^_8xN8lIVrr?A@f@*r%>u*rltpNawW88JmmN zJ#BW{f4!2}joo|LH(j-)J52WuZXejFN&6oDn{F;x4-I?j&BJak#h=_Ay?xobxpwWi z^T_DvlJaQpk*7Xz$$;i2G?&y%s(ML*mxk6yhxOGt#I>L9_F~W%_1EvQpZV-po_^sw zR4lN*IcjP@ag{@de;-38+E(;tV?H&T_33;*f3F+UcB9@>HND;*Y5GWQG@5XjkYGfg z*j*#<71-_~L=tO6qR)bjOlgHkfaDp!kWT@M3;yg!fwCjZvemJ~0Vhk=^0Z@hIzsnU z6!;D=WVwNbAQXuJ#u%Ck?7|WZ?FIs&U@#PphGB_@L!n5N4h7Ic6-ofrQlU^RRz;x# zf7=y7ESCYpY$B1MnQ)k<8J1;JEW=hfF2P14!9V~2zl*v9#8Om}hqDWS>tdl?$me*T zi$$Xp6^C1es8EPY@lkYQ48R)$@R~hSJLzeylYYO~*BX7ZuPP4!wb5I*uYWmgqn4aJ z)jQe(t*&1-^&?eRT+sg4hlmdVXkSQOe+b&@_g_Y!efL-2+5+uISD^KPZX3WJG&tX6 z!V5H#biVoyP}W=I`WQ)g?n8Qrfpqqg4j6r?#o@X+?n2%<-3kRYo zRu4t7B1n>apdTLB<>TYyT1^t5Z0FMyRe~B0r_|~P$K_I?czm3slS$|;3;?IonG{3E zpwo){@Iz8159K_&zn5YYQ3`<{@k<1#Af8AiV^kdQA4d2enhN0G>^J(#e^65&P{ICn zykPJ83wHam*}QBjN2c0#u|NFxFQbC}8$Y-Z`=|f#X@vdUhd#T7{h!%j->t7WG<@G? za&r(z!Um=|AU1{01*zwVr#lQ$?iEgT^@ES_0M0WK_cYSM+lFFD_vuIv#W;_N5f(ub zgl8lq>A|V}MEE#2(1n)|f4KS@xH{m~fvy63&ne}>7k~Ia&F}ps^>6^gZt^OJp!Yi% zcG0#n(&sIGG_8+YMjJ}GI<2>7db_T)%{hR23g5yR5r9EG3F{N6WZPoPm7V=<$75Q`F&7%^}i>$=2bvMsjJjcr-91*hJf59$*I~BJJf

I2yE98XBA((F1jGe%5MS6+LU+0PQiaE}A9v5-mAW7dCe{lBNjE3;Dg}G)6bHf!T z&<51}^Vs6KKVp9Jt}w`P-3 zeKwmZW8ETPMu z3d_O$ot>Rsw4eqYLj;H=j%s}pnE@J@gGpzyz%8Xr8J1E>KB!h}$A<@phX=AG3b3yR zY*9NXSA=|;2bKvTk_v{9X9%Il0IWwKlq=T`9}(d4e`2``5(KnHu|&uAQfZJ69Fs`K z;Hn@+(OIEdtxAUv9ahCksj#07qfTH=h8ZpqgYGCQ!SJx`<;xWj{7X8M$?{B`qQem| zCa|!`_rz&B!6aie%U8w2YPBTPvV0^Ar=UO~0o zj>A$e=Ck0};(O_`P%KsU(&=<2&A@4)e=vBg5E6sPZy{X>8ES&%Q$U0n7+0w#R;3J^ zPdYa=>*>dr)S}K#uErUcCD908Lc^7#j``LIaE6hf5fOm zqzIM1(F08g#M*cOycz>ueKCSrzkNTiHZJSNWz9Tls!f+ypLi*Oc=c<=3wiai@xLNo z{rYcxb&FSDCT(8z<~pRGRUJ5bZ5`>CK?IILKYL<|(B&6*kHaWf9EA0^8HkTFaW?as zFqam#NC&6RO{XqVyRHt~!s!vef4wK9a`t%C?~hz?e zLQ{OofBSaJ5WcA( zSOv5nu19W`^My+rhO^Tp2$ztmlts{stkvb^u4r`x9^z$I?r5@fLnjhEm&#|jObU1{8jDgfV7OR38j8Z4Q8W$|2z(Mx zWK%2y^b&{UIFrs-OGOH4f5${2U)alK_!P@=5!exCV^oYvr&FNw_V@Qe6y~Ajfg+

ud?R<{oEGmsN@A=&H=fn^ly;pvx|Q2@xsmHxk!SHOtQc87 z+3OKk(?u&hBMk#^=CVR1XiJJ9B*rj%6hI~`4o+@wzX)?Yg+dNs2jE7NJwi#S9#lbK z9703{NB0Z4bc#)(>T|vbiMWN|62@bC5;HD3B)ke{3ehf*pe;0mcmS6Ra78Oz=^AK#=Z-F(L79=Nq)f zk*TRJreDAJ-w4w`mA??v{byc9n7+IC$1P0XwKL=&bHQMOS9GFK3+i0O=#?TuddBZYs<>eYieME|1nI3-+t%&&Rh|KZY_iC3qyPq;|8=1UqbjQ zPJtHM`mAou+GDUUx>ldInytyq)JGG&O}uMYbcFR!f7n0jP;Bg9@#W0~9KW;X0I(O< zBAZ>{DUrhgXv)}<2I691EqUI<0i0}gVg6VUgj0^=0sgxH{!j$cU6?G2zPL={U^wgwRhz5SAh82}3d8%kMLAapow;rP;U$ zbdjV1fAuT}4FhXpFJCBTOT|*T2+2|+vPyVXvH~J$kP478rvUyDl+v`pIjGH(p{Ab< zp_Xf{2No*-Ctl5W&$AG^dSqx=sQl_by^4_kc=tl&fBc)TAmmeTY*+L5ZT;MXbfDpY z!%Y#ec!G4WhJ)Tki(6*#h{dDL@k)zp0gMSqf2^E)as&){*d1NSxue7iHg04#GTbfY zKEdzc+m~frl=teF+N;#TSj%^SZet6$%o@1g{)1gI{@uS>9LC~0McLSWAHpuT#_j^c zF4|Ta>a00#jsdoMvu@0lc75C&HS047ejCcX)t>58;*&7WD?9EFuULDoXva>4IbXeh zf1a@GHVe}NlPG+>kmj5Yw5{eGw4>2>;G_jB_20yv!q+{vHoJLgaHg6visE>K^qd<1}cbpm+15URiO(B}|T-|4)$1=Vj5>-irN zuOir&+H0q`W;n2ud diff --git a/tests/resources/pacbio_style_tags_no_untagged.bam b/tests/resources/pacbio_style_tags_no_untagged.bam index f41f7985fc060ef63d56231aa82bab08b9c055f6..079f8519859b6827858378fed8f33aef90c042b1 100644 GIT binary patch delta 4170 zcmV-Q5Vh~7C%Y%Ge*=Fn5Ve|Xj3e1q$E&(+mtB5UmCJ6o-&MA|Yj@hFo_@?m$=I&Q z9$=Qqjz6Fl>RCV8NX)R=o*5yn5YxsI4G4(xkPr_MX@LkNqzI%{#768cB0?b{1eCB! z@+A<$DuBppwNi)@vYczzs%Dt7h!$Z zTF|!KXia9i+NddNdsfq1O}RFnYI;p+Os37&Y&sh08u2mI14$4ZG+)}WT4LRnU_Xm} zz6Pra>tJxU6r}-oUO5*A+#tP_mV>Tvc5W`DPFFC+v#wAXbShnQuvjcjwo~CcX>nl= z2CQkC*LHmYit>N^{XTzyjz$AiFbE%lbO^Q4Xgr zBVp?48$J4ld8{-}hPrXmQ?!$g+-vuhp4wOYih9!5iD!TM$A-~v$$ERNj5W2P8;! zPJ8s-(_4YSz0)3X=NEdu$#b{np>@>N+TruXOF zXFvPi7w&)FdNn;aj_=G~_}SUCwr=Xr|NPV-PVda8FMg(bI-Qrl{!g9L2ge_}Ke~JG zgFiPT?mXyD4EUe(lLSc;8;iAvT(u;6p*L$cX*JfVs~2|Z>MYVZ%|64vXx)>&)9RJP zZtUK}y6LJV-C??SaQnbUOs%Px-^Qehwq>o^m`{ynZ91RNYx=a^sI?STtF=d}Hc}dmCLAUz zFrc1aM6MG!R`i^ZZ0MMa|w%d!cUVM~7;7h^*~zt0B**hSreglRg?!`V3ydOlao zWivd_g+n2lj=-$~bRfVb_z*fV3~CO8nvI^Jp7hk#Nx$FgtBtyNDpK(~zxd{sia)Yau}6j5iuQlV z!})0wUgeRbqXkDf>n(D9j3hkwAw8rZoxP-^&J7vk`n8X{3~&xuU0Vafm_TArLAc5% zNtYP*Gm#h@@Nb5|9)onQbjN(bH|O4kB*{(MIsMr~uyz0TYml111=)j1EhY(?# z5O4q0ukkwliWSQP2oO(Q=R(ka6a#+*+Lq->d)k&qHAR;J5G@Gk?eRpZjmJh)QRWkU zrjCd&!-)?qKz!JW-8U^jpwB~C2ZQtUKoZX@k}x={oSk>i&83)5F9+h17x}Kli=4US z%!HWO~8w=402Z%U9fI}pfh==J25P?EO7#cDV!RR;o@=#SC zQ%wIAJf`n%#`N}8qj}Yk>xO^Qc4;vD=2uWm|J5H{N`te%c?HoR^MTKA(cq_S8n_`C z7t6-)o3sz&+OUO-OT>LzxRKq^m}}meB;+S zq`u$8qiJp2(%TT~HS ztYyJvt9{DJlrD z8A85~245WjNrt15NQ{38!9G8iNM*CxOtDz4K!h(9r6<9*mx^VvlFtiFG!l&`1A%BP zOw({_I23`{?WZHq5k-Z;E$sSMT?RvOW;gGHp3QfkZqHU1yj=7;(fUC#3kpE zu|pT!S6sG`Sac^`j9_NM9y@{GiuAH;Gx!r1MR*)}o?PRiL6RF|zyW5D2BPPq(gVYZ z{f12s13-eZh6{iB@q5gVzxl;D{bB!l1Qg;&{W^!_{2h!RXj>ks^2mVvpb5ErO&_)N zT5C2L)n>DqY&0}kZ>jB7esFVxABY>OWdl~WfN=NQDLreQ|4v}Vj)fpC&!tsvzcL7m z?VQ16QO3}%bO&H0oRnb+{cpWAkDG@Y~dG`t~7rTp%lP(O&$(cIuf+qPM57n z3)2+3?3u6}JlxsY*+mO#Ah8GnY9h!0#8J)=1(q<$R2pbgNEKlz6r_W4xq5tfaCmqi zRYU>yRe>MXgJMa@CVAja0Ab4?K&eFll^P&J@dYx)>fsduF3%Us;7h=Omc4ATBtn{#Or_F1 z6QQG25b_sTSd<<`qR|)=4@X(PEFPB21)-YegA|;C3K21&nX0%fYn=4@eJD!ELm*iL zpr$=0-|_?6O=52S*~PVpcD=1)wTOTCu7+6lV(EWUmVGqyX~Z)9;Tu~l({@(AS^29@sk;@Pc}b`!qPd-FgXTpe~F-8yLk>E)qKj`?g8!S+KC- zd_R9~-bntxY1T$1MqAZl&`0%X{TaOwgO$~)Y*093Q~1q4qN5Lg|EDA3FaOEz{Rdx+ z|I&ZTyYv?_8+g&U~0XnKk-G$geX&hKPiiSK6938M1O0&Tj%R3jCxCy9 zDHvCwDwZoLHqLQS(Sb{tcoc5HCRrwxOe7gLl`E98xt#d^0@RI)`CLg9#PWeC0OXS# zmq-`$hvh=Ca9HJ`io`-)mkxyzOf14A5{Yz*gJKdBLp??ir_p$PzNc^W;;7hb3|DpD zq1qfOod!be5Sd1~ulK+<0%;qM0knT%0PV{`EJ57g1a0H0re9T!dQ)k-pnc}kF$CJL z=Pw2A)8l_cp#7I$|JoL4UnMQjdh0>tCT1OddXqnL(#H{~r`ni9OjGm~&Ro!G=p8maT;pu6Wq$YiSHJc5FMs2EZjnXYy2;)#S0(0FNoFwRW?q*Nl;(jHiaIjOO!IZTU6QtvQR(zGY7rZFtLv zR?83@(@L-kXg^$!a$atu3T1yo&MHGe6$&{`Ndyncnq6M%ie^XPp)$=%9j?Qgf+=(a z;Toz1!K@OLq5wnS|EQ=RO4jMUOeVXZ;*(6AiN<0)m&m5LR07Zz3Ww-0KrS2!1wt_Y z6p8=^0V$DKI>9nPmy%zrPPg zF$*mZlnnV1e4?ll1+H$W7$<$GP;2luG~lhHz&{|6(tb=1?=8Ozp5FT}OjfnJ-oUu? zGr#rMs4$iKgG;${=aZjXm8t*t3tQY#$V&{KqjQfPu*JBMRBqzOt2j?)MiL~p%Uffz zjU0h>O%Pl1Bk&=4O(TC+9K*TdaVf5i!_eT~&AoclbUIe@k(p$});(;1%v&o5(?s6J zqxGn*^*27-d+)nHBu7q4>((B8)}nWJr4m(_`;wZOJ0N==uWqguN)YSff5 z=v&pN&51H!&GKHegd)NA)pj zyiZ}~;-4?h;u>{BRa}vO?!kW}k^e;YQjzyxcpZuS{^D=9M1J3@KL3adrW-u^yDA{v z^a#@n9%xBicLW}i@euBW>C1yNNn-ZHHV7n+6ui#khHHQ8K?riQJCg8R=Nhi9$H9^_ zCm*y&>q#0N>*Ba48%vMOgP=aW28C=n0%e`?w6^k{n9U0RzYP=ATW^2gnRY^2uBK3q z!w?_AV1c&bR|&t(DR82#&1(9rJ%$`dQ)|;!vo)C++GwJ+iFd4+kHr58>rY5z3kKIV z>ed7tzqfye0^}g(BAs5~wIv4%Xg=Dk48+C2T=KjJ3Y=tiVRBj!gfougfd{)lgFq0f z#4zs^gm(%;dIaSzbQsc)z5R5O=eSDo=*cUGsCdMMLNHUsB_l9x6hME7g2W94A_~I& zUUr{>i&OOsmt-R%@FNa$PcfE*hJiJQfdxP)}H;;6-~Jzzx*(h(m{0are5AjS(b9s`fZV%#lYXdtnG@2NguQpE1) zYS2C%B^J{0Y{&NMZm{3T3LL?2R6w>W#hP6x){fe%>V&@b_@Lg#mUdZNyT7wbM!x$O Ui^Fh4qyG=4j<4YYCXpEWjCHG0$qx?9JyyC`53-cbc~+M*BSQ51ck0o(^#pg_<9hS3DIT{I}%v}qb7 zHIn>E+5|=(NDRY366E20SE9&EY9`5_y;|y0m-k*?YIV+czH@))k9P>C)FRh;2I}YF zc*id2Txqr^Gec|ERjo6t8|{`-Kc4DFU2RULt@dm>9veFGK{NwN5Ij_$x@6bHx-LUI zhkd>VhXu!Qc)FD2A%9jq6NmgTvy_*^o_Km@E#z)bw4~FXSRHn&J!`mFEG@2E<+~Ya zVGW0zWm(sD0wI5j3I>CLV2Eaz5ETx?hcF#MHNzw_*<7w%Dv4sLTo%Q0p^z`s_N)6+ zMcRL)QYvS&1wO&B5jGyrXET|-Y&p9Zqo_bI3=>6T;V=~n214O*gbGKY&_Y@el95Ox z7L7%z$8YrM8`h!PJQ^A1QD4=Mx=Oz@Q2W|I8>rgRz#xB~92}Zvr>z*BWA#|qnughI z9;;35c&_NCZWy!HY|>PfX;W$6n%)`TdRA2$Mnk70qLqmlfrzq7c)2 z^WIY*f7gGrcW%C%nVW~VXV3o1Y+7Gejc0y+Y7!^6=hNpt);pQbD_{6@_vHTJ`|geJ z-2K&Gn-RC~_a-L%&jm?>Bne-!TgX*Sq8ECzn@PK}yRKf?&{Z$e>9*;y8MK~h(`e61 zVmEg0Vc&FBlkPIzJGgUTt0wJx_;0!yuwEMW)SG{Y-3-N_++DprY+c>BaoD{)Zd_3w z>tBBQ-B%2#ZbNlNy`rjD6nJUqcx+T(^&xKjLa!f%x#)j|_03UL z``PO}y8I^?GSRuBx0>^**{V zXdwuN!+J3#_T+n{wy~Mi#v@awt1a0MeFC);t{hM!Wf%g3?(0V|(4PXx%oL@5G1sX{@KYf1( zDC;wFeT^hM_aQyRKsv3Y14bWeaky@-yO4J-*i+jJ!YF{mo`P_JPLeJ`Y@Uh4r~vvjn^0LPl8gfdP!r^zZ?_=7TT2^)J6G81zeapbjBA zc-;qWlQBoZ(i1psGK=_JcUVbZev=p)r~ z4#Ih6Z#T)sBQyd(;+HT`K`fq3MClmdpF;Q_nF`?F95e^YNK+qD!2WGKVDI5H2fRAaRiO2pP#%2ohu>*^_wVUP0~mG_*LeiJKfr&mi_Vp?K5y&e zY5llubRd+g(|Tv7cj`*VoCBz*@GG1V0a(;ia6E}Zwk@{5Xk+&&TU;x4{fyntV%X+w ziaLX`kcSH_$e#;u0fl%$;;O=0cMF0g@PaJC6*1Er?f|$cDhy#5qI-Y_#T)`WMwwVF z&PJdu$S2deTrOKG)vA9G%*!SDQ4r+iQbnp33L?wIm_#ZRV&YMnhN00&48nAfj=@9} z6$N3i6Ie|djwIN<0+2yXt^p_1YGN^+ib7ww5{7>O?J#r(T(f{{y|1>9MuWj z%F#fFe*kryrm`d9zvQwZp1 zcR#ZQ^bf9;xhMX4vF&Lg{aErxC|vq(d4X zpL^qz3Htfq`hI^DVvELg9)b5;7+auoWvnS<6NG;Ya`L({ZX5OXY%;FTW;4ZX>Wb0U zIxFG-gpVx{OK5h)Unu~>ZQrEkoZbKJ(26M*qP#qlS2_Ah3oN#i0gFSaK(E>xf@W}X zfhF*SWXV#mTI~U^EUUmI!lK%fhZ3apLeCQ5adaYIoMP>1Kd(fm!K4j@_waKJ3QDwIM|n~k_c@zV2j#*sVwGF z0ji}}3_g(iU`L9u*13e(YaoE4z#=1OG=@|RRPoe|g=%}`-TOrUTm?};%? zoJ~X-PN+x+l}b^pWrQ#VU67#=2kz2Tm${lp{lNgt1!V-h)dULDACg4*KJ8{FH~-;$ zPodN3Xqcz)-`~~|b3IqSkhwmPeF-s_@!<6>=F)#Jtt86%6tLBdOB7GQc4Io$&rJ{< z0>cH2gJ4g}fw794osa}^p1#2VPq_SplM#L+4-GoDf-mR5=?AA+H_h>VEhIh!@8(Dy z?^AlH-BEkz_$uzx;AHLk9+*}-H{$q66U(JYtPl2Wli0IhLBjcd+=TEaIx70gQ$ZiK zE9-yH=!Ga(T^*{AbobVzOTF>eH1psuekmsX*5B>iyZ_0=^IunX=+|t)7Dg0Iu3gd3 zV-$?eRb@VD0sQK7r3uNFHaF`Xt!08>n@kP8JyEB}<|^6RQ13@Bt*Lh-&_1#eZWuNH zRdDT!E>2~eZ{>CFC1Z>$S;7K%lV2>me6@dDUC3T(CiRwbZ)q*L;SxCZR1~{1H1(E1 zt-ZxU5@DjM+*?Rpl%RobqiD#`K&}DRf}isSNg%>BtgBd)D%CWX;CXOnUVj&r9zqg8x%0lQAP($U zIuc2;@fe#-CNpUsOd&RoW{e?PM(YU#9>38~Adl7@t?b#6)*7kZCSuePQiRID=z}H% zVr@PIUX22;z8J<*zdb*%Hm~Z&Rn32FwA7Z%t53cZN4)xl!iBv0`0@WBUj5SVeSV8q zUnXr{_4+!bpH&?=deb`6FM|jigMRkJ6rsy6Zi~Yxm>q=8w^@jfGjTTaCSfiuY>^I5 zotsHrqIR=7>F6C!J*391NlnQ9+4awT<-cD1;&6u=O~59tbCnMK}z`h7@47S<#;+%rs}M& z9P3TclZsaFv|2{p9INW_)Kt{*e7-V>Um;!NSc3Y^%@&VOM?SD?mhf!_!6Bk{7>}}9 z-e)dt7|wQ=z+6JAQkFn3a#nv&kb9EV6$P-DIl0SsIZL#}t|(qZ77y$Mfx`qG1X_<` zg5W!6cC*>sURp@836_b+1wNTe^XVk;S|l2wqrh;{SR@pIJ)=krC=mE07SAL(7U(4g z&T%@Gs}u_~(vIi3rUXW!_Yu+Q94SeQc2KxdwY8z3Ug2kK#_luAYvW{8BCzy zMyh!<0GnBdw~+~NT^0TT{ivOXq~_iXx~$WG&$(%+-Y}XN>HP8^{wwmE(tmj&>D>O? z*H&Kh5B_tDbX4*JG3S8X6Z$tu^ree_?(xFSlXH;-i8Xnxqm5`Iu&D`R3ws2p}M_ND=nmBJkXAvuilJx|U2lU#QPsMkw8r&Ld)%z6#{g%|n6@VBe6*H%|SW_B@T%h_{oWIXf+0TPAm&N z5(=+&e<8a|Bv^lr7i^13HAd#qs1G|xfbbTiW?(YGN9`d&dO%@9;(yOi(3%ZXQ(a8I zcK@3Q(?6fP5YvNaUqzU{xA^xhOy9FTg z&XWlnc^s=4GZ~r z0adr5vg5tFPG;I$sSLBG1~&NRQIdM|t?xN|MF_gJG)iAs;{6ympmX>U!cTDuw9wIK zbz{~!hV(_(>eF_sJ(-#Mc%pZRx9y0Iu>N8Dk2(|^``3JaGXdA{t~mhG3u}?dEO3{| z;Q+K{Y*l}U(qd>W1;N7sylnMg|5y~oQ=S(9{yPBvP#D}@*enXeI|TtdhSCx`3Yo?3 zUM3~*e6@7t(aQ(O2I38mTucJGNHBnU zj)#hcBfgu<7c#{{u~YzODIQ)qJS$lNmNZBOaG8IT0RJ#@X*!TOsI8-srXP(UmTT>Y zW-5OTkLEk)r4YK>Ff`0me)T_IMaX}mcOmjW^p#f-^6A&NqxsKl{kaF}K*Irt>mp$G z1nFQ62fd3HH_YM|i$|N|krvki7!#0KIQRGn81k??x{z~Mi3Mytk=e*_HYBgtXv0UE#CpU%`M=vYv6wUPj|@JxBqEz5RK{dj{q6