Choosing the chemistry¶
FiberHMM needs to know two things about your data: the enzyme that marked accessible DNA and, for Hia5, the sequencing platform. Together they pick the bundled model, the observation mode and every chemistry-dependent default.
--enzyme¶
| Your assay | --enzyme |
What FiberHMM reads |
|---|---|---|
| Fiber-seq with the m6A methyltransferase Hia5 | hia5 |
m6A calls in MM/ML |
| DAF-seq with the double-stranded deaminase DddB | dddb |
C→T / G→A conversions |
| DAF-seq with the single-strand-preferring deaminase DddA | ddda |
C→T / G→A conversions |
fiberhmm-pipeline, fiberhmm-call, fiberhmm-apply,
fiberhmm-recall-tfs/-recall-nucs and fiberhmm-posteriors take
--enzyme. Only these three presets are accepted. fiberhmm-qc also takes
--enzyme (default auto, read from the header).
For anything else, supply a model file with -m/--model (see
Training a model).
--seq¶
--seq pacbio|nanopore matters for Hia5 only:
- PacBio reports m6A on both strands of the molecule (A calls on the read
strand as
A+a, and A calls on the opposite strand asT-a). - Nanopore reports m6A only on the sequenced strand (
A+a).
The two platforms therefore have different models (hia5_pacbio.json,
hia5_nanopore.json), different observation modes (pacbio-fiber,
nanopore-fiber) and different defaults (see below).
For DddA and DddB the deamination read-out does not depend on the platform:
both use the daf mode on every platform, and --seq only sets the platform
recorded in the output header. Without --seq, fiberhmm-call takes it from
the header (@RG PL, aligner or basecaller @PG) and otherwise records
pacbio for DddA and nanopore for DddB; fiberhmm-pipeline records
nanopore for DAF reads it aligns that name no platform.
Automatic detection of --seq¶
When --seq is omitted for Hia5, fiberhmm-pipeline, fiberhmm-call,
fiberhmm-apply, fiberhmm-posteriors and
fiberhmm-recall-tfs/-recall-nucs look at the input, strongest evidence
first:
- a
FIBERHMM-CHEMISTRYdeclaration in the header (a BAM FiberHMM already called); - the
MMspecifications of the first 200 records:T-apresent means PacBio,A+aalone means Nanopore; @RG PLvalues and@PGprogram names (for examplepbmm2,ccs,jasminefor PacBio;dorado,guppy,minimap2 -x map-ontfor Nanopore).
The choice is printed:
NOTE: --seq not given; using --seq nanopore (detected from MM specs of 200 read(s) (A+a only, no T-a)).
If the sources disagree, or more than 10% of the inspected reads show the
other MM pattern, the run stops and asks for --seq:
error: fiberhmm-call: cannot infer the sequencing platform for --enzyme hia5: MM specs are mixed: ... Pass --seq pacbio or --seq nanopore.
Input that cannot be inspected (stdin) falls back to PacBio with a warning,
so pass --seq when you stream:
WARNING: --seq not given for --enzyme hia5 and the platform could not be detected (stdin (not inspected)); assuming PacBio (--seq pacbio). Pass --seq nanopore for Nanopore data.
For Hia5, an explicit --seq that the reads contradict is refused, because
the other platform's model changes the calls:
error: fiberhmm-call: --seq pacbio was given, but the input looks like nanopore (MM specs of 200 read(s) (A+a only, no T-a)). Calling nanopore reads with the pacbio model changes the calls (for example, Nanopore reads called as PacBio give ~100x more TF calls). Pass --seq nanopore, omit --seq to detect it, or add --force-seq to use --seq pacbio anyway.
Hia5 runs on reads whose first records carry no m6A calls at all
(MM A+a or T-a) are refused too, with a hint that the data may be
DAF-seq. --force-seq (in fiberhmm-call, -apply, -posteriors and
-recall-tfs/-recall-nucs) runs anyway in both cases. For DddA and DddB
a --seq that disagrees with the header only warns.
What the chemistry changes¶
| Setting | Hia5 PacBio | Hia5 Nanopore | DddB | DddA |
|---|---|---|---|---|
| Model | hia5_pacbio.json |
hia5_nanopore.json |
dddb_nanopore.json |
ddda_nuc.json (HMM), ddda_TF.json (TF recall) |
| Observation mode | pacbio-fiber |
nanopore-fiber |
daf |
daf |
ML threshold in fiberhmm-call |
128 | 248 | 128 (MM/ML dU only) | 128 (MM/ML dU only) |
| Nucleosome recall | accessible-cut split | accessible-cut split | accessible-cut split | phase-aware radial |
CpG mask (first-pass HMM and recall; --use-m5c) |
off | off | off | on |
| Adjacent-target thinning | — | — | off | runs ≥ 2, keep-one |
| Duplicate marking, SNP screen | — | — | on (file input) | on (file input) |
The full table, including the thresholds of the other tools, is in Chemistries and platforms.
Re-calling a BAM FiberHMM already called¶
fiberhmm-call, fiberhmm-apply and fiberhmm-recall-tfs/-recall-nucs write a
FIBERHMM-CHEMISTRY line into
the output header. Tools that read such a BAM take the chemistry from it:
fiberhmm-recall-tfs/-recall-nucsfill a missing--seqfrom it;fiberhmm-extractandfiberhmm-qcuse it for their ML threshold and QC reference;fiberhmm-consensusandfiberhmm-transferuse it to pick the emission model.
A run whose chemistry contradicts the declaration stops with a one-line fix:
error: fiberhmm-recall-tfs: the input BAM declares chemistry [assay=daf enzyme=dddb platform=nanopore mode=daf] but this run would declare [assay=fiber-seq enzyme=hia5 platform=nanopore mode=nanopore-fiber]. Pass --enzyme/--seq matching the input (with --model for a custom table), or --replace-chemistry to re-declare the output deliberately.
A custom --model given without --enzyme inherits the declared enzyme and
platform when its observation mode matches, and then uses that enzyme's
defaults:
NOTE: fiberhmm-recall-tfs: custom model on an input declaring [assay=daf enzyme=dddb platform=nanopore mode=daf]; using the defaults of --enzyme dddb --seq nanopore with the given model file(s).
--replace-chemistry drops the input's declaration and writes this run's
(for a custom model, enzyme=custom;platform=unknown). Use it only when you
are deliberately re-calling with a different chemistry.