Methylation (long-read)
TL;DR: JBrowse reads DNA methylation straight from the MM/ML tags nanopore and PacBio basecallers write. This tutorial follows one dataset, HG002 nanopore reads over an imprinting center, from per-read calls to an aggregate profile to the two parental alleles pulled apart.
Prerequisites
- nothing to install to read along: every figure loads hosted data
- for your own data, long reads whose BAM or CRAM already carries
MM/MLmodification tags, which modern ONT and PacBio basecallers write by default, plus a JBrowse instance to load them into (the web quickstart, or the desktop quickstart, which opens a local modBAM with no hosting step) - modkit for the aggregate section only, a single-binary download from its releases page
Where the data comes from
Both files are region slices of public
ONT open data on the ont-open-data S3
bucket.
- the HP1 bedMethyl from the
wf-human-variationsup run on HG002, restricted to the SNRPN locus and tom(5mC) rows: https://ont-open-data.s3.amazonaws.com/giab_2025.01/analysis/wf-human-variation/sup/HG002/PAW70337/output/SAMPLE.wf_mods.1.bedmethyl.gz - the HP2 bedMethyl from that same run: https://ont-open-data.s3.amazonaws.com/giab_2025.01/analysis/wf-human-variation/sup/HG002/PAW70337/output/SAMPLE.wf_mods.2.bedmethyl.gz
- the reads, from the HG002 sup basecalls, sliced to the same locus and
haplotagged with
whatshap haplotagagainst the phased SNP calls from that samewf-human-variationrun: https://ont-open-data.s3.amazonaws.com/giab_2023.05/analysis/hg002/sup/PAO83395.pass.cram - the HP1 slice the figures actually load, rehosted so they load without the S3 round trip: https://jbrowse.org/demos/methylation/HG002_SNRPN_hp1.modkit.bed.gz
- the haplotagged read slice beside it: https://jbrowse.org/demos/methylation/HG002_SNRPN_5mC_haplotagged.bam
The SNRPN imprinting center
At this locus on chr15, one parental allele is methylated and the other is not. That makes it a dataset with its own control: the views have to come out as two populations, and the reads and the aggregate profile have to agree on which allele is which.
Per-read methylation from the alignments
Load the modBAM as an AlignmentsTrack. Its assemblyNames must match an
assembly already configured in JBrowse (see the
assemblies configuration guide), and the
.bai index sits beside the file:
{
"type": "AlignmentsTrack",
"trackId": "HG002_snrpn_5mC_reads",
"name": "HG002 ONT reads (5mC, haplotagged)",
"assemblyNames": ["hg38"],
"adapter": {
"type": "BamAdapter",
"uri": "https://jbrowse.org/demos/methylation/HG002_SNRPN_5mC_haplotagged.bam"
}
}jbrowse add-track https://jbrowse.org/demos/methylation/HG002_SNRPN_5mC_haplotagged.bam \
--trackId HG002_snrpn_5mC_reads \
--name "HG002 ONT reads (5mC, haplotagged)" \
--assemblyNames hg38Set Color by... → Modifications from the track menu and each read is painted with its own 5mC calls. Two modes are offered: one paints only the positions the MM tag reports as modified, the other (IGV's "2-color" scheme) also fills in every CpG the tag left implicit, so an unmethylated region reads as solid blue. The alignments track guide covers both modes, the probability threshold, and the cytosine-context submenu.
The pileup over the CpG island is an interleaved mix of methylated and unmethylated reads. Splitting it by allele is one setting, two sections below.
Aggregate methylation with modkit bedMethyl
modkit pileup collapses the per-read calls into a bedMethyl file, one row per CpG per modification type, carrying the fraction of reads that were modified. It is the compact form of the same information, and stays fast at whole-genome zoom.
modkit pileup sample.bam output.bedmethyl --ref reference.fa --preset traditional
bgzip output.bedmethyl
tabix -p bed output.bedmethyl.gz--preset traditional collapses 5mC and 5hmC into a single 5mC fraction
(bisulfite-equivalent). Omit it to keep separate rows per modification type (m
for 5mC, h for 5hmC). Passing --partition-tag HP writes one file per
haplotype, which is what this dataset uses.
Because bedMethyl is a BED file with a numeric score column, it loads through a
BedTabixAdapter in a MultiQuantitativeTrack (see the
multi-quantitative track config guide).
JBrowse reads the modification type from the name column and gives each type
its own subtrack:
{
"type": "MultiQuantitativeTrack",
"trackId": "HG002_snrpn_modkit_hp1",
"name": "HG002 methylation HP1 (modkit)",
"assemblyNames": ["hg38"],
"adapter": {
"type": "BedTabixAdapter",
"uri": "https://jbrowse.org/demos/methylation/HG002_SNRPN_hp1.modkit.bed.gz"
}
}jbrowse add-track https://jbrowse.org/demos/methylation/HG002_SNRPN_hp1.modkit.bed.gz \
--trackId HG002_snrpn_modkit_hp1 \
--name "HG002 methylation HP1 (modkit)" \
--assemblyNames hg38 \
--trackType MultiQuantitativeTrackThe Y axis is percent methylation, each CpG a vertical bar.
Splitting the alleles apart
Each long read is a single DNA molecule, so reads that carry an HP haplotype
tag (from WhatsHap, HiPhase, or ONT's wf-human-variation) can be separated by
allele. Pick Group by... → Tag... from the track menu and enter HP. The
dialog scans the reads in view, reports the values it found, and offers to color
reads by the same tag, with that box arriving checked; uncheck it to keep
the methylation coloring. The pileup then stacks into one band per haplotype,
computed in the browser, one band methylated over the island and the other not.
Loading the two per-haplotype bedMethyl files above the reads puts the summary and its source in the same view, on one x scale.
The aggregate and the reads below it split the same way, molecule by molecule, with the same haplotype on the same side.
See the
alignments track guide for
the Group-by dialog and the phased-trio tutorial
for producing HP-tagged reads.
Aggregate for navigation, reads for detail
Keep the bedMethyl track for whole-genome navigation, since it stays quick at any zoom and is what a tumor-versus-normal comparison reads off, then drop the per-read BAM or CRAM below it once you are there, for the single-molecule and allele-specific detail only the reads carry.
To compare two samples rather than two alleles, run modkit dmr on their
per-sample pileups and load its BED output as a FeatureTrack beside the
bedMethyl tracks, so the differentially-methylated regions line up with the
positions driving them.
See also
- Alignments track
- Alignments track
- Methylation (bisulfite)
- Phased trio analysis (1000 Genomes)
- Comparing one genome's two haplotypes (T2T-HG002)
- RNA-seq visualization
- Multi-quantitative track
- modkit documentation
- Gallery: methylation and base modifications
Feedback on this tutorial is welcome: contact us.