Skip to content

v2.9.3: fiberhmm-extract --ry (DAF IUPAC deamination calls)

Choose a tag to compare

View all tags
@mtcicero26 mtcicero26 released this 21 Apr 04:31
· 122 commits to main since this release
v2.9.3
40612ac

New: fiberhmm-extract --ry

The DAF-seq analogue of --m6a for BAMs where deamination calls live as IUPAC R/Y codes in the query sequence (written by fiberhmm-daf-encode) rather than in MM/ML tags. Scans the stored sequence for R (G→A) and Y (C→T) positions, maps them to reference coordinates, and emits one BED12 row per read with each call as a 1 bp block.

Why this exists

fiberhmm-daf-encode converts DAF deamination mismatches into IUPAC ambiguity codes in the query sequence and does not produce MM/ML tags. Previously there was no way to get per-deamination positions back out of a FiberHMM-processed DAF BAM into a bigBed for FiberBrowser / UCSC / pyft-style visualization — --m6a and --m5c rely on MM/ML. --ry closes that gap.

Output

column content
1–12 classic BED12 (one row per read, each R/Y as a 1 bp block)
score (col 5) constant 255 — deamination calls from IUPAC encoding are deterministic, not probabilistic
strand (col 6) read alignment strand (+ or -)
13 (with --block-scores) int[blockCount] blockMod — per-block byte, 0 = R (G→A), 1 = Y (C→T)

autoSQL schema (table fiberhmm_ry) is embedded in the bigBed via bedToBigBed -as= so downstream tools can introspect it with bigBedInfo -as file.bb.

CLI

# DAF BAM (IUPAC-encoded) -> BED12 + bigBed
fiberhmm-extract -i recalled.bam -o out/ --ry

# With per-block R/Y disambiguation
fiberhmm-extract -i recalled.bam -o out/ --ry --block-scores

# Default "all types" now includes --ry; no-op on non-DAF BAMs (0 rows, skipped)
fiberhmm-extract -i recalled.bam -o out/ --all

End-to-end verification

chr3L slice of a Drosophila DddB recalled BAM:

  • 4,373 reads, 1.27 M individual R/Y deamination calls
  • 1.7 MB bigBed
  • All rows pass blockCount == len(blockSizes) == len(blockStarts) == len(blockMod)
  • autoSQL correctly reports definedFieldCount: 12, extraFieldCount: 1

Tests

6 new unit tests cover the R/Y scan, blockMod disambiguation, indel handling, empty sequences, and schema shape. Full suite (276 tests) passes.

Install

pip install --upgrade fiberhmm
Assets 3
Loading