diff --git a/CHANGELOG.md b/CHANGELOG.md index d74fa17f..8f47af1f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,14 @@ 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. + - [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 - [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..6872a0c7 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, @@ -490,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 { @@ -503,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) { @@ -538,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, @@ -545,6 +565,8 @@ impl GenericPileupWorker { pileup_numeric_options, combine_strands, max_depth, + phased, + mod_codes, }) } } @@ -573,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 { @@ -619,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); @@ -629,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 @@ -641,7 +681,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); @@ -738,6 +781,7 @@ impl PileupWorker for GenericPileupWorker { call, ref_base, motif_idxs, + hp, ); let Ok(_) = update_mods_iter2( @@ -779,6 +823,7 @@ impl PileupWorker for GenericPileupWorker { }, ref_base, motif_idxs, + hp, ); } else { add_to_tally( @@ -786,6 +831,7 @@ impl PileupWorker for GenericPileupWorker { Call::NoCall(base), ref_base, motif_idxs, + hp, ); } } @@ -803,7 +849,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 +859,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 { @@ -852,6 +904,7 @@ impl PileupWorker for GenericPileupWorker { ), ref_base, motif_idxs, + hp, ); mod_pos = Some(pos); canonical_base = Some(can_base); @@ -882,6 +935,7 @@ impl PileupWorker for GenericPileupWorker { }, ref_base, motif_idxs, + hp, ); } else { add_to_tally( @@ -889,6 +943,7 @@ impl PileupWorker for GenericPileupWorker { Call::NoCall(base), ref_base, motif_idxs, + hp, ); } mod_pos = Some(pos); @@ -907,6 +962,7 @@ impl PileupWorker for GenericPileupWorker { Call::Delete, ref_base, motif_idxs, + hp, ); } (Some(q), None) => { @@ -931,6 +987,7 @@ impl PileupWorker for GenericPileupWorker { }, ref_base, motif_idxs, + hp, ); } else { add_to_tally( @@ -938,6 +995,7 @@ impl PileupWorker for GenericPileupWorker { Call::NoCall(base), ref_base, motif_idxs, + hp, ); } } @@ -950,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; @@ -2458,7 +2537,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 +2552,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 +2561,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 +2574,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..ac4070fb 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, @@ -538,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>)> { @@ -1023,6 +1042,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 { @@ -1513,6 +1552,8 @@ impl ModBamPileup { pileup_options.clone(), self.combine_strands, self.max_depth, + self.phased, + self.requested_mod_codes(), ) }) .collect::>>()?; @@ -1623,6 +1664,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 +1675,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 +1688,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 +1704,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 +1718,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 +2462,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..186f3833 --- /dev/null +++ b/modkit/tests/test_pacbio_style_tags.rs @@ -0,0 +1,395 @@ +//! 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, 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}; + +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() +} + +/// 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"); + 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()); + // 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 = + 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); +} + +#[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 new file mode 100644 index 00000000..c701dfe6 --- /dev/null +++ b/tests/make_pacbio_style_tags.py @@ -0,0 +1,92 @@ +"""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; 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. + +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(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, 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: + 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(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]: + 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 00000000..07be36ee Binary files /dev/null and b/tests/resources/pacbio_style_tags.bam differ diff --git a/tests/resources/pacbio_style_tags.bam.bai b/tests/resources/pacbio_style_tags.bam.bai new file mode 100644 index 00000000..6bf9f178 Binary files /dev/null and b/tests/resources/pacbio_style_tags.bam.bai differ 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 00000000..079f8519 Binary files /dev/null and b/tests/resources/pacbio_style_tags_no_untagged.bam differ diff --git a/tests/resources/pacbio_style_tags_no_untagged.bam.bai b/tests/resources/pacbio_style_tags_no_untagged.bam.bai new file mode 100644 index 00000000..a1d01ce4 Binary files /dev/null and b/tests/resources/pacbio_style_tags_no_untagged.bam.bai differ diff --git a/tests/resources/pacbio_style_tags_pileup_general.bed b/tests/resources/pacbio_style_tags_pileup_general.bed new file mode 100644 index 00000000..02c7c6fc --- /dev/null +++ b/tests/resources/pacbio_style_tags_pileup_general.bed @@ -0,0 +1,36 @@ +oligo_1512_adapters 19 20 h 4 + 19 20 255,0,0 4 100.00 4 0 0 0 0 0 2 +oligo_1512_adapters 19 20 m 4 + 19 20 255,0,0 4 0.00 0 0 4 0 0 0 2 +oligo_1512_adapters 63 64 h 6 + 63 64 255,0,0 6 33.33 2 1 3 0 0 0 0 +oligo_1512_adapters 63 64 m 6 + 63 64 255,0,0 6 50.00 3 1 2 0 0 0 0 +oligo_1512_adapters 64 65 h 1 - 64 65 255,0,0 1 0.00 0 1 0 1 0 1 0 +oligo_1512_adapters 64 65 m 1 - 64 65 255,0,0 1 0.00 0 1 0 1 0 1 0 +oligo_1512_adapters 69 70 h 5 + 69 70 255,0,0 5 0.00 0 0 5 1 0 0 0 +oligo_1512_adapters 69 70 m 5 + 69 70 255,0,0 5 100.00 5 0 0 1 0 0 0 +oligo_1512_adapters 70 71 h 3 - 70 71 255,0,0 3 0.00 0 0 3 0 0 0 0 +oligo_1512_adapters 70 71 m 3 - 70 71 255,0,0 3 100.00 3 0 0 0 0 0 0 +oligo_1512_adapters 72 73 h 6 + 72 73 255,0,0 6 0.00 0 2 4 0 0 0 0 +oligo_1512_adapters 72 73 m 6 + 72 73 255,0,0 6 66.67 4 2 0 0 0 0 0 +oligo_1512_adapters 73 74 h 3 - 73 74 255,0,0 3 0.00 0 0 3 0 0 0 0 +oligo_1512_adapters 73 74 m 3 - 73 74 255,0,0 3 100.00 3 0 0 0 0 0 0 +oligo_1512_adapters 90 91 h 5 + 90 91 255,0,0 5 0.00 0 3 2 0 0 0 1 +oligo_1512_adapters 90 91 m 5 + 90 91 255,0,0 5 40.00 2 3 0 0 0 0 1 +oligo_1512_adapters 91 92 h 3 - 91 92 255,0,0 3 0.00 0 3 0 0 0 0 0 +oligo_1512_adapters 91 92 m 3 - 91 92 255,0,0 3 0.00 0 3 0 0 0 0 0 +oligo_1512_adapters 93 94 h 4 + 93 94 255,0,0 4 0.00 0 2 2 1 0 1 0 +oligo_1512_adapters 93 94 m 4 + 93 94 255,0,0 4 50.00 2 2 0 1 0 1 0 +oligo_1512_adapters 94 95 h 2 - 94 95 255,0,0 2 0.00 0 2 0 1 0 0 0 +oligo_1512_adapters 94 95 m 2 - 94 95 255,0,0 2 0.00 0 2 0 1 0 0 0 +oligo_1512_adapters 100 101 h 6 + 100 101 255,0,0 6 0.00 0 1 5 0 0 0 0 +oligo_1512_adapters 100 101 m 6 + 100 101 255,0,0 6 83.33 5 1 0 0 0 0 0 +oligo_1512_adapters 101 102 h 2 - 101 102 255,0,0 2 0.00 0 1 1 0 0 1 0 +oligo_1512_adapters 101 102 m 2 - 101 102 255,0,0 2 50.00 1 1 0 0 0 1 0 +oligo_1512_adapters 124 125 h 3 + 124 125 255,0,0 3 0.00 0 0 3 0 0 0 1 +oligo_1512_adapters 124 125 m 3 + 124 125 255,0,0 3 100.00 3 0 0 0 0 0 1 +oligo_1512_adapters 125 126 h 3 - 125 126 255,0,0 3 33.33 1 0 2 0 0 0 0 +oligo_1512_adapters 125 126 m 3 - 125 126 255,0,0 3 66.67 2 0 1 0 0 0 0 +oligo_1512_adapters 135 136 h 2 + 135 136 255,0,0 2 0.00 0 2 0 0 0 0 0 +oligo_1512_adapters 135 136 m 2 + 135 136 255,0,0 2 0.00 0 2 0 0 0 0 0 +oligo_1512_adapters 136 137 h 2 - 136 137 255,0,0 2 50.00 1 1 0 0 0 0 1 +oligo_1512_adapters 136 137 m 2 - 136 137 255,0,0 2 0.00 0 1 1 0 0 0 1 +oligo_1512_adapters 146 147 h 1 - 146 147 255,0,0 1 100.00 1 0 0 0 0 0 1 +oligo_1512_adapters 146 147 m 1 - 146 147 255,0,0 1 0.00 0 0 1 0 0 0 1