Methylation on PacBio HiFi — what a consumer actually holds, and whether a module can bin it¶
Question: a consumer holding PacBio HiFi samples has methylation data. In which files and fields,
at what scale, and can a just-dna annotation module point at it and bin it the way repeat_alleles.csv
bins FORMAT/REPCN? The open question from the VCF 4.5 audit (§6) was whether any real tool writes
the 4.5 M5mC/DPM5mC/ADM5mC FORMAT keys, or whether PacBio users hold bedMethyl-like BED files
instead, which a VCF pointer never reaches.
Date of analysis: 2026-09-30. This is evidence, never contract, the standing rule for
everything under docs/probes/. It proposes no design and files nothing; §6 prices what the
measurements imply.
Basis — tools, read from source and run where marked. Every repository was cloned on 2026-09-30
into data/interim/methylation/ (git-ignored); versions are the latest GitHub release on that day.
| Tool | Version read | Released | Source read | Run here |
|---|---|---|---|---|
| TRGT | v5.1.0 (HEAD b21c621) |
2026-06-10 | docs/vcf_files.md, src/trgt/writers/write_vcf.rs, src/trgt/workflows/tr.rs, src/trgt/reads/read.rs, CHANGELOG.md, tags v0.7.0–v0.9.0 |
yes, binary v5.1.0 |
| pb-CpG-tools | v3.0.0 (HEAD 119ff2b, 2025-07-09) |
2025-01-28 | README.md, CHANGELOG.md |
yes, binary v3.0.0, model and count modes |
| MethBat | v1.1.0 (HEAD 859e01c, 2026-07-15) |
2026-06-02 | docs/pileup_guide.md, docs/profile_guide.md, docs/report_guide.md, docs/migration_v1.md, CHANGELOG.md, data/report_regions/ |
yes, binary v1.1.0: pileup, profile, report |
| HiFi-human-WGS-WDL | v4.0.0 (15e82cb) |
2026-08-18 | docs/tools.md, docs/tools_containers.md, docs/singleton.md, workflows/wdl-common/wdl/tasks/methbat.wdl; docs/tools_containers.md at tags v2.1.0, v3.0.0, v3.1.0, v3.2.0 |
no |
| modkit (ONT) | v0.6.4 (HEAD 5cecc3f) |
2026-06-11 | book/src/intro_pileup.md, book/src/migrating_060.md, tests/resources/ |
yes, binary v0.6.4 (see §1.6) |
| pbsv | v2.11.0 | 2025-02-26 | whole repo grepped | no |
| sawfish | v2.2.1 | 2025-11-04 | whole repo grepped | no |
| HiFiCNV | v1.0.1 | 2024-10-22 | whole repo grepped | no |
| Paraphase | v4.1.0 | 2026-09-25 | whole repo grepped | no |
| DeepVariant | v1.10.0 (latest release 2026-03-05) | — | GitHub code search over google/deepvariant; deepvariant/channels/base_methylation_channel.cc, deepvariant/methylation_aware_phasing.h |
no |
| rastair | v2.2.0 | 2026-08-24 | src/vcf/schema.rs, tests/snapshots/call_cli__vcf_with_ml.snap |
no |
| VCF 4.5 | hts-specs 510c107 (2026-08-12) |
— | data/interim/vcf45/hts-specs/VCFv4.5.tex, the §1.6.2 block and Table 2 |
— |
Basis — data. Every file below is public and was opened here; §2 gives the URLs.
- PacBio's 2022Q4 Revio WGS pipeline outputs for the HG002 trio (TRGT 0.5.0, pb-CpG-tools 2.3.2).
- PacBio's 2026Q2 HG002 SPRQ-Nx run (TRGT 5.0.0, pb-CpG-tools 3.0.0, MethBat 0.17.0). The TRGT VCF and CpG BEDs were read by tabix range request, and the haplotagged BAM was sliced over three loci.
- PacBio's 2024Q4 PureTarget repeat-expansion panel on Coriell lines (TRGT 1.1.2), including three FMR1 expansion carriers.
The single most load-bearing measurement. Two female FMR1 carriers, each with one expanded and one normal allele, as TRGT reports them:
| sample | karyotype | MC (CGG per allele) |
AM (methylation per allele) |
|---|---|---|---|
| HM06968 | XX | 33,112 |
0.90,0.07 |
| NA06905 | XX | 23,79 |
0.03,0.88 |
| NA09237 | XY | 898 |
0.85 |
In HM06968 the expanded allele is the unmethylated one. Every element rule the format has picks a
value by the field's own magnitude, so largest over AM selects 0.90, which belongs to the
33-repeat allele. The bin a fragile X module means is "the methylation of the expanded allele", and
that means selecting an index of AM by the values of a different field, MC. No member of
VALID_ELEMENT_RULES can say that. The male record is the case the fmr1_cgg_repeat example's prose
describes, and it holds there: 898 CGG, AM=0.85.
1. What each tool emits¶
1.1 TRGT — a per-allele methylation mean in a VCF FORMAT field¶
The field is AM. The header line TRGT writes today (write_vcf.rs:39, and every file opened here from
TRGT 1.1.2 on):
docs/vcf_files.md describes it as "Mean 5mCpG methylation level per allele (. when no CpG sites or
methylation data unavailable)", one value per called allele, in TRGT allele order, which the same
page says is not reordered to match a phased GT.
What the number is, read from the source rather than the page:
- Per read,
get_tr_meth(tr.rs:260) averages the ML probability (value / 255.0) over the CpG dinucleotides that fall inside the repeat span of that read. Flanks are excluded, so this is the methylation of the tract, not of the promoter CpG island around it. - Per allele,
get_meth_by_hap(tr.rs:199) assigns each spanning read to an allele by its tract length (assign_read) and takes the mean of those per-read means. The allele assignment is by length, not by theHPhaplotag. encode_am(write_vcf.rs:401) rounds to two decimals and writes a missing value for an allele with no reads that carried a methylation profile.- A haploid call (
GT=1, a male chrX, or--karyotype XY) carries one value.
So it is a mean of means, a probability-weighted fraction, not a count of methylated reads. It is 5mC
only: get_meth keeps modifications whose canonical base is C and never asks the modification code.
That is safe for 5mC-only BAMs, and TRGT 1.5.1 fixed a bug in exactly this path for BAMs carrying
several modifications. What it does with a BAM carrying both 5mC and 5hmC calls on the same C (the
SPRQ-Nx data in §2 does) was not established.
The type and scale moved once, and the changelog does not say so. The 2022Q4 file, written by TRGT 0.5.0, declares:
##FORMAT=<ID=AM,Number=.,Type=Integer,Description="Mean methylation level per allele">
chrX 147912051 . CGGCGG…CGG CGGCGG…CGG 0 . TRID=FMR1;END=147912110;MOTIFS=CGG;STRUC=(CGG)n GT:AL:ALLR:SD:MC:MS:AP:AM 1:93:90-97:18:31:0(0-93):0.978495:50
chr1 44836 . AAATAAATAAATAAATAAATAAATAAATAAAT AAATAAATAAATAAATAAATAAATAAATAAATAAATAAATAAAT 0 . TRID=chr1_44835_44867;END=44867;MOTIFS=AAAT;STRUC=(AAAT)n GT:AL:ALLR:SD:MC:MS:AP:AM 0/1:32,44:31-32,43-44:8,3:8,11:0(0-32),0(0-44):1,1:159,185
The earliest source tag in the repository, v0.7.0 (2023-12-11), already declares Type=Float. TRGT
0.1.0–0.5.0 were released without source, and CHANGELOG.md has no entry for the change, so the
integer scale (values up to 185 are seen, which reads like the 0–255 ML byte) is inferred, not
verified. An integer AM also appears on the AAAT repeat above, which has no CpG in its span, so the
0.5.0 computation was not the one quoted from 5.1.0 either. A band authored on [0, 1] misreads every
0.5.0 file.
A current record, from the 2026Q2 pipeline run on HG002 (male, normal FMR1):
##trgtVersion=5.0.0-e9acee0
chrX 147912049 . CGCGG…GGC CGCGG…GGC . . TRID=FXS_FMR1;END=147912111;MOTIFS=CGG;STRUC=<TR> GT:AL:ALLR:SD:MC:MS:AP:AM 1:95:88-98:27:29:0(1-31)_0(34-61)_0(64-94):0.957895:0.1
The same run carries a second, wider record over the same tract
(TRID=chrX_147911948_147912141_GGC, AM=0.11), from a different catalog entry. The TRID of the
FMR1 record is FMR1 in the 2022Q4 and PureTarget files and FXS_FMR1 in the 2026Q2 one, so which
record a consumer reads for "FMR1" depends on the catalog the pipeline shipped.
Rerunning TRGT 5.1.0 here on the PureTarget BAMs (--preset targeted) and on the HG002 slice
reproduced the published values: NA09237 GT=1, MC=895, AM=0.85; HM06968 MC=32,112,
AM=0.9,0.07; HG002 MC=29, AM=0.1.
AM=. has two readings, and the file cannot tell them apart. get_tr_meth returns nothing both
when the span holds no CpG (the AAAT repeat in a modern file) and when the input BAM carries no MM/ML
tags at all. No FORMAT field says whether methylation was called. SD is a depth, but a depth of
length support, not of methylation calls.
1.2 pb-CpG-tools — per-CpG BED, one file per haplotype¶
aligned_bam_to_cpg_scores writes <prefix>.combined.bed.gz, and hap1/hap2 when the BAM is
haplotagged, each with a tabix index and a bigWig. Since v3.0.0 the BEDs are bgzipped and carry a
## header; v2 wrote plain, headerless BED (the 2022Q4 files, 1.3 GB uncompressed, are v2.3.2).
Columns 1–6 are shared by both pileup modes: chrom, 0-based begin (the C of the CpG: the reference
at chr15:24954889–24954890, 1-based, reads CG), end (one base), mod_score
(percent), type (Total, hap1, hap2), cov. model mode (the default) adds
est_mod_count, est_unmod_count, discretized_mod_score; count mode adds mod_count,
unmod_count, avg_mod_score, avg_unmod_score. One row per CpG, strands collapsed. The bigWig keeps
columns 1–4.
A real header and rows, from the 2026Q2 HG002 run, over the SNURF DMR:
##pb-cpg-tools-version=3.0.0
##cmdline=aligned_bam_to_cpg_scores --threads 16 --bam HG002.GRCh38.haplotagged.bam --ref human_GRCh38_no_alt_analysis_set.fasta --output-prefix HG002.GRCh38 --min-mapq 1 --min-coverage 4
##pileup-mode=model
##modsites-mode=denovo
##min-coverage=4
##min-mapq=1
##basemod-source=jasmine 26.1.3 (commit v26.1.3)
#chrom begin end mod_score type cov est_mod_count est_unmod_count discretized_mod_score
chr15 24954888 24954889 46.2 Total 63 30 33 47.6
chr15 24954888 24954889 86.7 hap1 32 28 4 87.5
chr15 24954888 24954889 4.1 hap2 31 1 30 3.2
The 2022Q4 v2.3.2 row shape, for comparison (no header at all):
Running v3.0.0 here on the BAM slice reproduced those three rows exactly in model mode. In count
mode the same site reads 47.6 Total, 78.1 hap1, 16.1 hap2. So the two modes of one tool
disagree by 12 points on the unmethylated haplotype at one CpG, and the mode is a header line, not a
column.
1.3 MethBat — the pileup that replaced pb-CpG-tools, plus region labels¶
MethBat 1.0.0 (2026-05-19) added methbat pileup and dropped pb-CpG-tools as an input. It writes one
BED per modification, with Total/hap1/hap2 rows interleaved in one file, 15 columns. From the run
here (methbat pileup, v1.1.0, same BAM slice):
##methbat_version=1.1.0-fbf4686
##base_modification=5mC
#chrom start end name score strand mod_score type cov mod_count unmod_count inferred_unmod_count diff_base_count avg_mod_score avg_unmod_score
chr15 24954888 24954889 m 476 . 47.6 Total 63 30 33 0 0 0.897 0.114
chr15 24954888 24954889 m 781 . 78.1 hap1 32 25 7 0 0 0.929 0.158
chr15 24954888 24954889 m 161 . 16.1 hap2 31 5 26 0 0 0.736 0.102
and in the 5hmC file, stranded:
mod_score is percent, score is the same value in tenths (0–1000), and name is the modification
code. The default is pb-CpG-tools count mode (the numbers above match it), not model mode.
Two downstream commands aggregate over regions:
methbat profilelabels each regionMethylated(mean ≥ 80 %),Unmethylated(≤ 20 %),AlleleSpecificMethylation(≥ 50 % of sites phased, Fisher p ≤ 0.01, haplotype delta ≥ 50 points),UncategorizedorNoData, withmean_hap1_methyl,mean_hap2_methyl,mean_combined_methyl, site counts and median coverages.methbat reportcompares those labels against an expected category per region and writesPASS,Inconclusive,AnomalousQcWarningorAnomalous, with QC warnings (LowPhasedSites,LowHaplotypeCoverage,WeakASM). MethBat shipsdata/report_regions/hg38_imprinting_targets.tsv, 15 imprinted DMRs from Mackay et al. 2022, Table 2, eachexpected_category=AlleleSpecificMethylationandanomalous_categories=Methylated;Unmethylated.
methbat report run here against that file:
chrom start end region_label report_summary qc_warnings expected_category summary_label mean_combined_methyl mean_meth_delta mean_hap1_methyl mean_hap2_methyl
chr11 2698717 2701029 KCNQ1OT1:TSS-DMR PASS PASS AlleleSpecificMethylation AlleleSpecificMethylation 49.4 71.1 15.3 86.4
chr15 24954856 24956829 SNURF:TSS-DMR PASS PASS AlleleSpecificMethylation AlleleSpecificMethylation 48.9 -65.5 80.1 14.6
(The other 13 regions report NoData or Inconclusive here only because the slice did not include
them.)
The scale moved at 1.0.0. docs/migration_v1.md calls it a breaking change: profile, report and
compare columns that were unit fractions are now percentages. The 2026Q2 pipeline ran MethBat 0.17.0
and its profile.tsv holds fractions:
#methbat_version=0.17.0-7c5e249
#command=methbat profile --input-prefix HG002.GRCh38.cpg_pileup --input-regions …/cpgIslandExt.sorted.hg38.tsv --output-region-profile HG002.GRCh38.methbat.profile.tsv
chr15 24954888 24955907 CpG:_77 AlleleSpecificMethylation … 0.8806080283353007 0.041076080548510176 -0.8395319477867907 … 0.48489584542693237 …
That is the same DMR in the same sample, 0.88/0.04 here and 80.1/14.6 from 1.1.0, under the
same column names.
1.4 The HiFi-human-WGS-WDL pipeline — which of these a user ends up holding¶
From docs/tools_containers.md at each release tag:
| pipeline | date | CpG pileup | region profile | TRGT |
|---|---|---|---|---|
| v2.1.0 | 2025-01-10 | pb-CpG-tools 2.3.2 | — | 1.4.1 |
| v3.0.0 | 2025-07-07 | pb-CpG-tools 3.0.0 | — | 3.0.0 |
| v3.1.0 | 2025-09-09 | pb-CpG-tools 3.0.0 | MethBat 0.15.0 | 4.0.0 |
| v3.2.0 | 2026-01-15 | pb-CpG-tools 3.0.0 | MethBat 0.15.0 | 5.0.0 |
| v4.0.0 | 2026-08-18 | MethBat 1.1.0 pileup |
MethBat 1.1.0 | 5.1.0 |
In v4.0.0 (methbat.wdl) the user holds, per sample:
<prefix>.5mC.bed.gz(+.tbi) — per-CpG,Total/hap1/hap2rows, the §1.3 layout;<prefix>.5hmC.bed.gz(+.tbi) — per-strand;<prefix>.6mA.bed.gz— declared as an output but skipped by default (skip_6mA = true);<prefix>.methbat.profile.tsv—methbat profileover themethbat_region_tsvshipped in the reference container. The 2026Q2 run usedcpgIslandExt.sorted.hg38.tsv, UCSC CpG islands;- the phased TRGT VCF, with
AM.
The pipeline runs methbat profile, not methbat report, so no pipeline output labels the
imprinted DMRs. A user holds CpG-island numbers, at island boundaries, which for SNURF (CpG:_77,
chr15:24954888–24955907) and KCNQ1OT1 (CpG:_159) sit inside the Mackay DMRs but are not them.
The 2026Q2 run's VCFs declare VCFv4.2 (TRGT, DeepVariant small variants), VCFv4.3 (MitorSaw)
and VCFv4.4 (sawfish structural variants). None declares 4.5.
1.5 The other PacBio callers: none writes a methylation field¶
- pbsv 2.11.0, HiFiCNV 1.0.1, Paraphase 4.1.0 — no occurrence of
methyl,M5mCorVCFv4.5anywhere in the repository. - sawfish 2.2.1 — two hits, both in a vendored BAM utility (
lib/rust-vc-utils/src/bam_utils/basemod.rs) that reads 5mC from MM/ML. Nothing in its VCF writer. - DeepVariant — reads 5mC as an input pileup channel (
BaseMethylationChannel) and uses it for methylation-aware phasing of unphased reads (methylation_aware_phasing.h). A code search of the repository forM5mCreturns zero hits. It consumes methylation and emits none. - GitHub,
org:PacificBiosciences M5mC— zero hits.
1.6 Beyond PacBio, for the record¶
ONT modkit 0.6.4 writes bedMethyl, 18 columns: chrom, start, end, modification code, score
(= N_valid_cov), strand, thickStart, thickEnd, colour, N_valid_cov, percent modified, N_mod,
N_canonical, N_other_mod, N_delete, N_fail, N_diff, N_nocall (book/src/intro_pileup.md). With
--phased it writes combined, hp1 and hp2 files. A fixture row from tests/resources:
The book mentions no VCF output anywhere. Run here on the PacBio SPRQ-Nx slice, modkit produced zero
rows: modkit modbam check-tags flagged 103 of 161 records conflict-explicit-prob-greater-than-one
(jasmine 26.1.3 writes 5mC and 5hmC calls on the same C), and pileup still processed nothing after
adjust-mods --ignore h and --no-filtering. That was not investigated further, and it says nothing
about ONT data.
Who writes VCF 4.5 base-modification fields at all. GitHub code search on 2026-09-30:
"ID=M5mC" 22 hits, "ID=DPM5mC" 13. Apart from spec copies and RDF fixtures, two producers:
- rastair (TAPS / 5-base, Ludwig Institute; used by JAX's
5baseTAPSpipeline) writes##fileformat=VCFv4.5withM5mC/DPM5mC/ADM5mC, but declares themNumber=., not the spec'sNumber=M. Its source says why: "Number=M… is not in the VCF grammar. htslib tolerates it, but a strict parser (noodles, and so every tool built on it) rejects the whole file rather than the line." It writes a value on both the C and the G of a CpG, which is the strand-specific form the spec allows (its INFOM5mC_Strandscounts per strand); theNumber=.is the departure. A real record from its test snapshot:
##FORMAT=<ID=M5mC,Number=.,Type=Float,Description="Methylation level at CpG sites, one value per CpG context">
chr19 6105712 . C . 99 PASS AD=10;BQ=35.5778;DP=18;MQ=60;M5mC_Strands=4,8,6,0;CPG GT:GL:GC:DP:M5mC:DPM5mC:ADM5mC 0/0:99:18:18:0.666667:12:8
- Illumina DRAGEN 5-base gVCFs, as documented by Illumina Connected Annotations: FORMAT
M5mC/DPM5mCon reference blocks, plus an INFOM5mCdeclaredNumber=R,Type=String, a per-base context string, which is not the spec's field at all. The FORMAT header line is not shown in that page, and no DRAGEN file was opened here.
So the fields exist in the wild, from short-read methylation chemistries, and neither producer writes them in the spec's shape. Nothing in the PacBio stack writes them.
2. The open samples¶
| what | URL (under https://downloads.pacbcloud.com/public/) |
opened here |
|---|---|---|
| 2022Q4 HG002 trio, pipeline outputs (README lists tool versions) | revio/2022Q4/WGS-variant-pipeline-analysis/ |
README; trgt/HG002.GRCh38.haplotagged.trgt.sorted.vcf.gz (7 MB, whole); first bytes of cpg/HG002.GRCh38.combined.bed and hap1.bed by HTTP range |
| 2024Q4 PureTarget panel, 24 Coriell lines | 2024Q4/Vega/PureTargetCoriell24/ |
README; puretarget_report.genotype.csv; TRGT_VCF_files/ for HM06968, NA06905, NA09237 (both reps), NA03697, NA13509, ND11494; PBMM2-BAM-Input-For-IGV-And-TRGT/ BAMs for NA09237_rep1 and HM06968 (≈ 0.4–0.6 MB each) |
| 2024Q4 HG002, pipeline v2.0.0-rc4 | 2024Q4/Vega/HG002/Human-WGS-variant-pipeline/ |
listed only (headerless v2-era CpG BEDs, TRGT VCF) |
| 2026Q2 HG002 SPRQ-Nx, three uses (5mC, 5hmC, 6mA called on instrument) | 2026Q2/HG002-SPRQ-Nx/Use1/analysis/ |
README, inputs.json, outputs.json, HG002.stats.txt, HG002.GRCh38.methbat.profile.tsv (5 MB, whole); tabix slices of HG002.GRCh38.trgt.sorted.vcf.gz and the three cpg_pileup.*.bed.gz; a 161-read slice of the 64.7 GB HG002.GRCh38.haplotagged.bam over SNURF, KCNQ1OT1 and FMR1 |
The PureTarget README states the data cover only each sample's target locus, for privacy. The
consumer-shaped files a module would read are all here: the TRGT VCF, the CpG BEDs in two generations,
and a MethBat profile. What was not found is any public PacBio file with a VCF 4.5
base-modification field. The searches: the full directory crawl of revio/ and of every quarterly
folder from 2024Q4 to 2026Q2 for vcf, and the GitHub code searches in §1.5–§1.6.
The runs here (§1.1–§1.3) used the reference as UCSC hg38 chr11, chr15 and chrX; the scratch
outputs are in data/interim/methylation/run/.
3. Six measured facts a design would have to meet¶
- The allele a band is about is selected by another field. HM06968 and NA06905 above. In both
female carriers one allele is about 0.9 and the other below 0.1, and which one is which does not
follow length. That is consistent with X inactivation (n = 2, not established here). The format's
consequence does not depend on the biology: on a two-allele record a band on
AMreads the value of whichever allele the rule picks, and no current rule can pick "the allele with the most CGG". - An expansion's CpGs have no reference coordinate. NA09237's reads carry about 2,700 bp of tract.
pb-CpG-tools run here on that BAM reports 88 CpG sites for the whole locus, all at reference
positions; the inserted tract has no position to report a CpG at. Per-base files cannot hold the
methylation of an expanded repeat. Only a per-allele summary like
AMcan. - The same quantity arrives on four scales across seven fields. TRGT 0.5.0
AMinteger (scale unverified), TRGT ≥ 0.7AMfraction, pb-CpG-toolsmod_scorepercent, MethBat 0.17 profile fraction, MethBat 1.x profile percent, MethBatscoretenths of a percent, and VCF 4.5M5mCfraction. Twice the scale changed under an unchanged field or column name. - The number depends on the tool and its mode. At chr15:24954888, hap2 is 4.1 % in pb-CpG-tools
modelmode and 16.1 % incountmode (MethBat's default). The mode is a header line. - Haplotype labels are not parental. In one sample the methylated allele is
hap1at SNURF andhap2at KCNQ1OT1. HiPhase haplotypes are phase-block-local, so no file here can say "the maternal allele is methylated". A band can say "one allele methylated and the other not", or bin the delta or the pooled value. - Absence has more than one cause. A missing
AMis either no CpG in the span or no MM/ML in the BAM (§1.1). A missing BED row is coverage below--min-coverage(4) or no CpG in the reads (denovosites mode), never "unmethylated". The H19/IGF2 IG-DMR in the 2026Q2 run has combined rows (mean coverage 11.9) but no haplotype rows at all, which MethBat's own README attributes to ALT contigs pulling reads away from the locus.
4. Mapping onto the format¶
4.1 FMR1 repeat methylation — reachable today, and one rule short¶
What a row would point at: the TRGT VCF, FORMAT/AM, Number=., Type=Float, one value per called
allele, on [0, 1] since TRGT 0.7. source_field: FORMAT/AM already passes the pointer grammar. The
bands the corpus example asserts in prose are a fraction of the tract methylated: NA09237 (full
mutation, male) 0.85, HG002 (normal, male) 0.10.
What it needs that the format lacks:
- a measure kind. None of the five
VALID_MEASURE_KINDSis a methylation fraction; binning it asallele_fractionwould put two quantities under one name (P5), as the VCF 4.5 audit §6 said. - a host table.
repeat_alleles.csvpinsmeasure_kind=repeat_count(RepeatAlleleRow._EXPECTED_KIND), and it is the table keyed(gene, repeat_unit)that already names the FMR1 tract. - an element rule that selects by another field. "The
AMvalue at the allele index whereMCis largest." On a hemizygous record (GT=1) any rule gives the one value, which is why the male case works. On a two-allele record none of the eight rules gives the right one (§3.1). - the unresolved sentinel, with its two causes named. The existing
unresolvedrow covers "no value"; it cannot say which of §1.1's two absences happened.
An adjacent observation, for the same module: no PacBio file opened here carries REPCN. The
example's count pointer FORMAT/REPCN is ExpansionHunter's key. TRGT carries the count as MC, a
Type=String field with _-joined per-motif counts inside each allele (18_8,25_8 at HTT), so the
count half of the module does not reach a PacBio user either, unless the author writes
FORMAT/REPCN|FORMAT/MC and the consumer can parse a motif-segmented string.
4.2 Imprinting (SNURF/SNRPN at 15q11–13, KCNQ1OT1 at 11p15 IC2, H19 at IC1) — real data, no pointer¶
What the data are: per-CpG rows in <prefix>.5mC.bed.gz (MethBat 1.x) or
<prefix>.cpg_pileup.{combined,hap1,hap2}.bed.gz (pb-CpG-tools 3.x), and a region aggregate in
profile.tsv. For HG002, measured three ways:
| DMR (Mackay 2022, GRCh38, 0-based) | pipeline pb-CpG-tools model, mean of sites | MethBat 1.1.0 report |
pipeline MethBat 0.17.0 profile (CpG island) |
|---|---|---|---|
| SNURF:TSS-DMR chr15:24954856–24956829 | Total 48.7, hap1 84.1, hap2 6.6 (113 sites) | ASM, PASS; 48.9 / 80.1 / 14.6 | CpG:_77: ASM; 0.88 / 0.04 |
| KCNQ1OT1:TSS-DMR chr11:2698717–2701029 | Total 48.5, hap1 8.9, hap2 87.3 (193 sites) | ASM, PASS; 49.4 / 15.3 / 86.4 | CpG:_159: ASM; 0.07 / 0.91 |
| H19/IGF2:IG-DMR chr11:1997581–2003510 | Total 38.6 (250 sites, mean cov 11.9), no haplotype rows | not in the slice | CpG:_19 Methylated, CpG:_27 Uncategorized, unphased |
What a module row would need to state:
- the region, as chrom, start and end on a named build. No VCF record exists for it, so
variant_keyandrsiddo not apply. - the file and column: a BED column (
mod_score) in rows of onetype, or aprofile.tsvcolumn (mean_combined_methyl).source_fieldis a VCF pointer by grammar and description, so it cannot reach either. - the aggregation from sites to region (MethBat's is the mean of site percentages over the sites
inside the region).
profile.tsvhas done it already, but only for the regions the pipeline was given, which in the 2026Q2 run were CpG islands. - the haplotype selector:
Total, or a statement about both haplotypes together (the delta, or "one high and one low"). §3.5 rules out "the maternal allele". - the unit (§3.3) and the modification (5mC; the 5hmC file sits beside it).
- the bands. MethBat's
report_regions/README.mdsays loss of imprinting "typically manifest[s] as fully methylated or unmethylated", and that where phasing is absent a combined value near 50 % may be treated as normal. That is a three-band pooled shape, with MethBat's own thresholds at 20 and 80 %, plus a haplotype-delta band for the ASM confirmation. - callability: a coverage floor (
covper site;median_hap*_coverageornum_phased_sitesper region), since a missing row is not an unmethylated one (§3.6).
MethBat's hg38_imprinting_targets.tsv is itself close to such a module: region, label, expected
category, anomalous categories. What it lacks for this format is a citation per row and bands a
consumer applies without running MethBat.
5. The design questions, against the data¶
- Per-base vs region. Settled by the data: both clinical questions are regional. TRGT has already
aggregated over the tract, per allele; MethBat
profile/reportaggregate over a region. The per-base BED is what a consumer holds, not what a band reads. A per-base element rule, the VCF 4.5 audit's "the value at this base on this strand", has no customer in the PacBio stack. - Pooled vs per-haplotype. Dissolved for FMR1:
AMis per allele, selected by length, not by haplotag. For imprinting both exist (Totalandhap1/hap2rows) and both are meaningful: the pooled value carries the three-band answer, and the haplotype values confirm allele-specificity, but never with a parent attached (§3.5). - Modification type as an axis. Real. MethBat writes one file per modification and a
namecolumn, the SPRQ-Nx instrument calls 5mC, 5hmC and 6mA, and VCF 4.5 names each (M5mC,M5hmC,M6mA). TRGT and pb-CpG-tools are 5mC only. If it is ever a kind, P5 wants the modification as its own axis rather than folded into several kind names. - Tissue. In no file opened here. HG002 and the Coriell samples are cell lines. MethBat's
data/README.mdwarns that its HPRC cell-line background gives more false calls on blood. It would be a column, asHeteroplasmyRow.tissueis, and nothing in the data supplies it. - The element rule. Sharpened. What is missing is not a per-base rule but a cross-field one: pick an index of one field by the values of another (§3.1).
- Callability via a depth field. Partly answered. BED rows carry
cov,mod_count/unmod_count, and profiles carry median coverages and phased-site counts. TRGT carries none for methylation, and its missing value is ambiguous (§1.1). - Whether
source_fieldreaches the data. Yes forFORMAT/AM, and nothing else. Every per-base or per-region methylation number on PacBio is in a BED or TSV. VCF 4.5M5mCwould be reachable by grammar, but no PacBio tool writes it, and the two producers found write it withNumber=..
6. What it would cost¶
Read against docs/CONSTITUTION.md in full, legality first (P3, P6, P8), then price (P9). Everything
below is additive, so everything is minor-legal and nothing is patch scope.
| addition | legality | P9 cost |
|---|---|---|
a measure_kind member for a methylation fraction |
minor: an additive vocabulary member (P3, P6) | small but authored, and a one-way door under P5, so the name wants the audit P5 asks for. The data say it should be a fraction on [0, 1] with the unit in the description, and consumers of MethBat 1.x or pb-CpG-tools divide by 100 |
hosting AM bins in repeat_alleles.csv |
minor, if the table accepts a second kind and the kind joins the bin-group key; a relaxed validation accepts more and invalidates nothing (P8) | no new column, but a table whose rows then mean two quantities; the alternative, a new table kind keyed like repeat_alleles.csv, is full cost |
| a cross-field element rule | minor: a new VALID_ELEMENT_RULES member is additive |
a member alone would have to hard-code MC, a TRGT convention, into the vocabulary. Naming the selecting field needs a new optional pointer column, full cost and a new shape of rule. The one-way door here is larger than §6 of the VCF 4.5 audit priced |
| a region-level methylation table (imprinting) | minor: a new optional table kind (P3) | full cost: region, build, file and column, row selector, aggregation, unit, modification, coverage floor, bands. P9's "one concern per table" gate is met only if it stays that one concern |
widening source_field to name a BED/TSV column |
grammar widening is minor-legal, since every value that passed still passes | cheap in columns, dear in meaning: one field would then name either a VCF key or a column in some other file, which is the overloading P5 forbids. The priced alternative is a separate pointer column |
The FMR1 half is the cheaper one and the only one with a live pointer: a kind member plus the
cross-field selection. The imprinting half is a new table at full cost, against no consumer request.
The RELEASE_CYCLE rule that a commit on main is patch scope only puts either half on the open minor
branch.
7. Verdict¶
Three answers, not one.
- FMR1 repeat methylation: real. TRGT writes
FORMAT/AMin every version whose output or source was opened (0.5.0, 0.7.0 to 5.1.0); every PacBio pipeline output opened here carries it; open samples exist, including a male full mutation (NA09237, 898 CGG,AM=0.85), reproduced here with TRGT 5.1.0. The pointer grammar reaches it today. A module cannot yet bin it correctly on a two-allele record. - Imprinting: the data are real, the pointer is not. HG002 shows allele-specific methylation at
SNURF and KCNQ1OT1 in pipeline outputs and in MethBat run here. All of it is in BED and TSV files,
and
source_fieldcannot name a column in either. - VCF 4.5
M5mCfrom PacBio: not real yet. No PacBio tool writes it, no PacBio VCF opened declares 4.5, and the two producers found anywhere (rastair, DRAGEN 5-base) are short-read and write it in a shape that departs from the spec'sNumber=M.
8. Could not verify¶
- The TRGT 0.5.0 integer
AMscale. Source for 0.1.0–0.5.0 is not in the repository; values up to 185 suggest the 0–255 ML byte. Also unverified: which release between 0.5.0 and 0.7.0 moved it to Float, sinceCHANGELOG.mdis silent. - TRGT on a BAM carrying 5mC and 5hmC calls on the same C.
get_methdoes not check the modification code. The HG002 SPRQ-Nx value (0.1) is low, as a normal male's should be, but whether 5hmC probabilities leak intoAMwas not tested. - The default
methbat_region_tsvin pipeline v4.0.0. It sits in a 3.2 GB reference-data tar (Zenodo 21517827) that was not downloaded. The 2026Q2 run usedcpgIslandExt.sorted.hg38.tsv. - Which pipeline release produced the 2026Q2 outputs.
inputs.jsonnamesGRCh38.ref_map.v3p1p0, and the tool versions in the file headers (pb-CpG-tools 3.0.0, MethBat 0.17.0, TRGT 5.0.0) match no taggedtools_containers.mdexactly. - X inactivation as the reason for the female
AMpattern. n = 2, cell lines, no XCI assay. - modkit on PacBio data. It produced no rows on the SPRQ-Nx slice; the cause was not diagnosed. modkit on ONT data was not run; its columns are quoted from its book and a test fixture.
- DRAGEN's actual VCF bytes and FORMAT header. Only Illumina's documentation page was read.
- DeepVariant's VCF writer. Known to consume 5mC; the claim that it emits no methylation field rests
on a code search for
M5mCand a read of two headers, not on the writer itself. - Clinical band values. The 20 and 80 % imprinting thresholds are MethBat's defaults, not a guideline. No paper was opened for a methylation threshold at FMR1 or at any DMR, so any bin values quoted here are observations from samples, not proposed boundaries.