v2.9.0: FIX reverse-strand MM/ML parsing bug
⚠️ CORRECTNESS FIX — reverse-strand MM/ML parsing
What was wrong
The SAM spec requires MM tag skip-counts to walk through base positions in the original sequencing direction. For reverse-aligned reads, this is the reverse complement of the stored SEQ field. Our `parse_mm_tag_query_positions` walked in SEQ direction regardless of strand, producing modification positions that landed at A bases in SEQ (for an "A+a" MM entry) instead of at the correct T bases in SEQ (since A in original = T in SEQ for reverse reads).
On real hia5 PacBio data (~47% reverse-aligned), every modification on reverse reads landed ~24 bp (median) from where it should have been — far enough to shift TF footprint calls, close enough that the HMM still produced "plausible" output. Footprint counts on reverse reads were ~7% too high from spurious calls generated by the shifted modification pattern.
Impact per code path
| Code path | Bug live since | Status |
|---|---|---|
| `fiberhmm-apply` (nucleosome/MSP) | v2.7.0 (Apr 16, 2026) — ~1 day | Now FIXED |
| `fiberhmm-recall-tfs` (TF calls) | v2.6.0 (Apr 13, 2026) — 4 days | Now FIXED |
| `fiberhmm-train` / `fiberhmm-probs` / `fiberhmm-utils` | v2.0.0 (Feb 22, 2026) — ~2 months | Now FIXED |
| `fiberhmm-extract` m6a/m5c BED | never affected (uses `pysam.modified_bases`) | unchanged |
| DAF with IUPAC `R`/`Y` in SEQ | never affected (different code path) | unchanged |
Bundled models were trained on output from the buggy parser. The HMM is robust enough that calls are still usable, but models are technically miscalibrated on reverse-read data — retraining is recommended long-term.
The fix
Three lines in `parse_mm_tag_query_positions`: for reverse-aligned reads, reverse-complement SEQ before searching for target bases, walk there, then flip positions back to SEQ frame via `q_len - 1 - pos`. Same transformation pysam does internally. Zero performance change.
Correctness validation
New regression test (`tests/test_mm_parser_vs_pysam.py`) validates byte-for-byte against `pysam.modified_bases` on every commit:
| Test | Result |
|---|---|
| Synthetic forward-strand read | PASS (byte-identical to pysam) |
| Synthetic reverse-strand read | PASS (byte-identical to pysam) |
| Synthetic DAF MM/ML reverse-strand | PASS |
| Empty MM, empty ML, overshoot skip | PASS |
| 500 real hia5 PacBio reads (mixed forward/reverse) | PASS — byte-identical on every (base, mod_code) group |
Full suite: 235 existing + 10 new = 245 tests pass, no regressions.
End-to-end behavior on real hia5 BAM
| Metric | v2.8.2 (buggy) | v2.9.0 (fixed) |
|---|---|---|
| Forward-read mean footprints/read | 65.8 | 65.8 (unchanged) |
| Reverse-read mean footprints/read | 71.4 | 66.5 (spurious calls removed) |
| Forward-read mean MSPs/read | 64.5 | 64.5 (unchanged) |
| Reverse-read mean MSPs/read | 62.8 | 65.1 |
Reverse-read statistics now track forward-read statistics symmetrically, as expected biologically.
How the bug fell through the cracks
Self-referential validation — every internal correctness check compared fiberhmm output to fiberhmm output (`fiberhmm-call` fused vs `fiberhmm-apply | fiberhmm-recall-tfs`, etc.). All internal paths had the same bug so they all agreed. Comparing to `pysam.modified_bases` catches it immediately; that comparison is now part of the standard test suite.
For users
If you have output BAMs from v2.0.0 through v2.8.2:
- Forward-read calls are correct (always were)
- Reverse-read modifications in `ns`/`nl`/`as`/`al` (apply output, since v2.7.0) and `MA`/`AQ` (recall-tfs output, since v2.6.0) are shifted by median ~24 bp
- Re-running with v2.9.0 gives correct output for all reads; no code changes needed
Install
```bash
pip install --upgrade fiberhmm
```
🤖 Generated with Claude Code