Keep records with conflicting mod probabilities; fix pileup on PacBio (Jasmine) tags

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.

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.

modkit-core/src/adjust.rs

@@ -120,7 +120,7 @@ fn adjust_mod_probs<'a>(
     sequence_motifs: &Option<SequenceMotifs<'a>>,
     discard_motifs: bool,
 ) -> MkResult<bam::Record> {
-    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:");

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(),

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 {