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

[ljwharbers]

Fixes #612 (PacBio BAMs give "~N failed processing" / 0 rows in 0.6.x pileup; the MM-tag-missing reads noted there are the second point below). Related: #545.

Problem

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. At the sparse 5hmC positions the two probabilities sum to well over 1.0 (PacBio documents that they "may sum to more than 100%"). modkit rejects the whole record with conflict-explicit-prob-greater-than-one, which removed 32–65% of the reads in our samples in adjust-mods, extract, summary, call-mods and the pre-0.6.1 pileup paths, with the loss only visible at DEBUG level.

Two further issues surfaced on the same BAMs in the 0.6.1+ pileup:

• a single record without MM/ML tags aborts the whole interval in the general workers (the aggregator only logs at DEBUG), so --region runs returned 0 rows;
• BaseModsAdapter bails on any opposite-strand sub-tag (T-a., present on every PacBio HiFi read), so the optimized workers fail every record and threshold sampling sees no probabilities. This is the failure reported in V0.6.1 does not work on pacbio samples, while v0.5.0 worked fine (no problems with ONT) #612.

Changes

• mod_bam: a position whose probabilities from separate sub-tags sum to more than one is dropped instead of the 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). Dropped positions are counted and reported at the end of the run; modbam check-tags still reports the records as invalid.
• pileup: the same positions are skipped in the general worker and the optimized adapter; a record that fails to parse no longer aborts the interval, failed intervals are reported as errors.
• pileup: BaseModsAdapter accepts - strand sub-tags (fixes threshold estimation) and pileup falls back to the general workers when such calls are present.
• Tests: synthetic PacBio-style BAMs (tests/make_pacbio_style_tags.py) with unit and integration tests; CHANGELOG and troubleshooting docs.

Verification

On a chr20 slice of an affected sample: adjust-mods --ignore h keeps 2,687 of 2,687 tagged primary records (was 1,035); pileup --cpg --modified-bases 5mC gives 21,958 CpG rows (was 0), matching modkit 0.4.4 on a BAM with the 5hmC track stripped; pileup --region on the full BAM gives 1.14 M rows (was 0); ONT output is byte-identical to 0.6.2. cargo test passes.

🤖 Generated with Claude Code

[ljwharbers]

Pushed a second commit (6e0afa2) that fixes two gaps in the general pileup workers, which this PR routes PacBio HiFi BAMs (opposite-strand T-a. 6mA calls) to:

• --phased was ignored by GenericPileupWorker: the HP tag was never read and phased_feature_counts stayed empty, so the hp1/hp2 outputs were empty files and only combined was written. The worker now keeps three tallies (combined, HP=1, HP=2) like DnaPileupWorker does. Without --phased the output is byte-identical to before (guarded by a snapshot test).
• --modified-bases did not restrict the output rows in the general path; every code present in the BAM produced a row (an h row next to every m row on PacBio data, unlike the optimized workers). 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. --combine-mods is unaffected. Note that this also applies to the other users of the general workers (multiple motifs, --duplex).

Tests: the PacBio-style fixtures now carry HP tags and a regular 5hmC call; test_pacbio_style_tags.rs gained tests for the phased partition (hp1 + hp2 coverage <= combined per position, combined == unphased run) and for the row filter (N_other of each m row equals N_mod of the dropped h row). Full cargo test is green.

Verified on a haplotagged PacBio chr20 subset (8,604 reads, 86 % with HP): pileup --cpg --modified-bases 5mC --phased now gives 74,937 m-only rows in combined (was 149,874 rows with h) and 72,204 / 73,178 rows in hp1 / hp2 (were 0), with no position where hp1 + hp2 coverage exceeds combined.