Repository navigation
v2.10.4: auto-skip --deam on fiber-seq BAMs
Bug fix
fiberhmm-extract --all (default mode) now auto-skips --deam when the input BAM looks like fiber-seq.
What was happening
On vg-surjected fiber-seq CRAMs (pangenome-graph aligned, projected onto linear reference):
--allmode ran every track including--deam— even on fiber-seq data that has no deaminations (Hia5 is a methyltransferase).--deampath-3 (MD walker) fell through on every read because paths 1 (MM/MLucode) and 2 (R/Y) had nothing to match.- Surjected BAMs frequently have MDs that disagree with CIGAR on ~97% of reads (documented upstream issue in
vg surject). - pysam's
get_aligned_pairs(with_seq=True)on bad MD raises AssertionError AND corrupts malloc state before propagating. - A later allocation in the worker segfaults. Python's
finallyblock doesn't run on segfault, so buffered m6a/m5c BED writes are lost. - Result:
--deamis empty, andm6a/m5cfiles are also empty or missing entirely.
Fix
Auto-skip --deam on BAMs where:
- MM/ML carries m6A (
A+a/T-a) or 5mC (C+m) subtypes, AND - no dU code (
+u/+55797), AND - no R/Y IUPAC in the query sequence
When triggered, prints:
Auto-skipping --deam: BAM looks like fiber-seq (MM/ML has m6A/5mC subtypes,
no dU code, no R/Y). Hia5 is a methyltransferase, not a deaminase --
no deamination calls to extract. Pass --deam explicitly to override.
Explicit --deam is still honored — the auto-skip only fires in --all / default mode.
Verified
Tested on a 5 Mb chr1 slice of a real surjected fiber-seq CRAM: 29,679 reads → 22.8M m6a + 4.3M m5c + 1.9M footprint + 1.9M msp positions, zero crashes, all bigBeds produced. Pre-fix this would have segfaulted mid-region and lost m6a/m5c writes.
Tests
294 pass. 8 new cover the fiber-seq signature detection across PacBio modern MM/ML, SAM 1.7 suffix variants (A+a., C+m?), DAF IUPAC, modkit dU, ChEBI 55797, and mixed cases.
Install
pip install --upgrade fiberhmm