Methylation (long-read)
2 videos in this tutorial
The tutorials target the JBrowse v5 beta, and the v4.3.0 release
on the download page lacks some of what they show. To install the
beta of JBrowse Web, run npm install -g @jbrowse/cli@next, then
jbrowse create jbrowse2 --branch v5.0.0-beta.13. Desktop beta builds are coming soon.
Nanopore and PacBio basecallers record each read's DNA methylation in
base-modification tags (MM and ML), and JBrowse draws those calls straight
from the BAM. We follow one dataset, HG002 nanopore reads over the SNRPN
imprinting center, where one parent's copy is methylated and the other's is not,
from per-read calls to an aggregate profile to the two parental alleles pulled
apart.
Prerequisites
- for your own data, long reads whose BAM or CRAM already has
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) - modkit 0.6 or later, Oxford Nanopore's single-binary tool for tallying per-read calls into per-CpG methylation fractions, for the aggregate section only
- WhatsHap, to haplotag reads of your own
that have no
HPtag
Where the data comes from
The files are region slices of public
ONT open data on the ont-open-data S3
bucket.
Nothing to download: the track configs below load region slices we host, cut from these files.
The files
- the HP1 bedMethyl (per-CpG methylation fractions) 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
Two parental alleles at the SNRPN imprinting center
Every view below is at chr15:24,948,000-24,962,000, the SNRPN imprinting
center; type it into the location box once the tracks are loaded. The reads and
the aggregate profile should agree on which allele is methylated.
Loading GRCh38
The track's assemblyNames has to name the assembly the reads were aligned to,
GRCh38.
Goes in the assemblies array of config.json. See Assemblies.
{
"name": "hg38",
"uri": "https://jbrowse.org/genomes/GRCh38/fasta/hg38.prefix.fa.gz",
"refNameAliases": {
"uri": "https://s3.amazonaws.com/jbrowse.org/genomes/GRCh38/hg38_aliases.txt"
},
"cytobands": "https://jbrowse.org/genomes/GRCh38/cytoBand.txt"
}jbrowse add-assembly https://jbrowse.org/genomes/GRCh38/fasta/hg38.prefix.fa.gz \
--name hg38 \
--refNameAliases https://s3.amazonaws.com/jbrowse.org/genomes/GRCh38/hg38_aliases.txt \
--config '{"cytobands":"https://jbrowse.org/genomes/GRCh38/cytoBand.txt"}'hg38 is one of the genomes JBrowse Desktop hosts, with gene tracks already set up: on the start screen click Show all available genomes and pick it. To load the files of this config instead:
In JBrowse Desktop, Open new genome on the start screen (or File → Open genome... in a session), then Open from a URL and paste, one per line:
https://jbrowse.org/genomes/GRCh38/fasta/hg38.prefix.fa.gz
https://jbrowse.org/genomes/GRCh38/fasta/hg38.prefix.fa.gz.fai
https://jbrowse.org/genomes/GRCh38/fasta/hg38.prefix.fa.gz.gziJBrowse reads the format off the file name. Then fill in:
- Genome name:
hg38 - refName aliases (under More options):
https://s3.amazonaws.com/jbrowse.org/genomes/GRCh38/hg38_aliases.txt - cytobands (under More options):
https://jbrowse.org/genomes/GRCh38/cytoBand.txt
The figures also show UCSC's CpG islands and the NCBI RefSeq genes. The hosted hg38 hub config, https://jbrowse.org/ucsc/hg38/config.json, has both:
- CpG Islands reads https://jbrowse.org/ucsc/hg38/cpgIslandExt.bed.gz with a CSI index
- NCBI RefSeq - RefSeq All reads https://jbrowse.org/ucsc/hg38/ncbiRefSeq.gff.gz
Per-read methylation from the alignments
Load the modBAM as an alignments track. For your own reads, swap uri for a
sorted modBAM or CRAM with the .bai or .crai beside it, aligned to the same
assembly:
Goes in the tracks array of config.json. See Tracks.
{
"trackId": "HG002_snrpn_5mC_reads",
"name": "HG002 ONT reads (5mC, haplotagged)",
"uri": "https://jbrowse.org/demos/methylation/HG002_SNRPN_5mC_haplotagged.bam",
"assemblyNames": ["hg38"]
}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 hg38In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track... and, in Add a track from file or URL, enter:
- Main file:
https://jbrowse.org/demos/methylation/HG002_SNRPN_5mC_haplotagged.bam
Click Next. JBrowse reads the adapter and track type off the file name. Then fill in:
- Track name:
HG002 ONT reads (5mC, haplotagged) - Assembly:
hg38
Click Add.
Open this track in JBrowse Web ↗
Open this track in JBrowse Desktop ↗ (JBrowse Desktop 5.0+)
File → Open track... also opens the file by URL and infers the .bai.
Color by... → Modifications in the track menu paints each read with its 5mC calls, in one of two modes:
- One color per modification type paints the positions the MM tag reports as modified
- One color per type, plus low-probability & unmodified in blue (IGV's "2-color" scheme) also fills in every CpG the tag left implicit, so an unmethylated region is solid blue
The alignments track guide covers both modes, the probability threshold, and the cytosine-context submenu.
Aggregate methylation with modkit bedMethyl
modkit pileup collapses the per-read calls into a bedMethyl file with one row per CpG per modification type, each holding the fraction of reads that were modified. The file is much smaller than the reads.1
# --cpg: count at CpG sites only
# --combine-strands: one row per CpG, the two strands summed
# --combine-mods with --modified-bases C: 5mC and 5hmC summed into one fraction
modkit pileup sample.bam output.bedmethyl --ref reference.fa \
--cpg --combine-strands --combine-mods --modified-bases C
bgzip output.bedmethyl
tabix -p bed output.bedmethyl.gzTwo modkit pileup options shape the output:
--combine-modscollapses 5mC and 5hmC into a single fraction (bisulfite-equivalent). Pass--modified-bases 5mC 5hmCwithout it to keep separate rows per modification type (mfor 5mC,hfor 5hmC)--phasedwrites one file per haplotype from the reads'HPtags, and this dataset uses it
With --phased the output argument is a directory, which gets hp1.bedmethyl,
hp2.bedmethyl and combined.bedmethyl:
modkit pileup sample.bam hp_bedmethyl/ --ref reference.fa \
--cpg --combine-strands --combine-mods --modified-bases C --phasedbedMethyl is a BED file with a numeric score column, so it loads as a
MultiQuantitativeTrack (see the
multi-quantitative track config guide).
JBrowse reads the modification type from the name column and draws one
subtrack per type, with a vertical bar per CpG on a percent-methylation axis:
Goes in the tracks array of config.json. See Tracks.
{
"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 MultiQuantitativeTrackIn JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add track from pasted JSON, and paste:
{
"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"
}
}Open this track in JBrowse Web ↗
Open this track in JBrowse Desktop ↗ (JBrowse Desktop 5.0+)
Splitting the alleles apart
Each long read is one DNA molecule, so reads with an HP haplotype tag (from
WhatsHap, HiPhase or ONT's wf-human-variation) separate by allele. In
grammar-of-graphics terms, Group by... sets a facet, which splits the
display into stacked sections, one per value of a field.
Pick Group by... → Tag... in the track menu and enter HP. The pileup then
stacks into:
- one band per haplotype, one methylated over the island and the other unmethylated
- an
HP: noneband for the reads the haplotagging could not assign
The dialog also offers to color reads by the tag; with methylation coloring on,
that box starts unchecked, and leaving it unchecked keeps the coloring.
Advanced → Edit plot... shows the resulting facet as text.
The two per-haplotype bedMethyl files from --phased load as one track above
the reads, with a row per file. Pinning the axis at 0 to 100 puts both rows on
one scale, so an unmethylated row stays flat:
Goes in the tracks array of config.json. See Tracks.
{
"type": "MultiQuantitativeTrack",
"trackId": "HG002_snrpn_modkit_multi",
"name": "HG002 5mC by haplotype (modkit)",
"assemblyNames": ["hg38"],
"adapter": {
"type": "MultiWiggleAdapter",
"subadapters": [
{
"type": "BedTabixAdapter",
"name": "HP1",
"color": "#d62728",
"uri": "https://jbrowse.org/demos/methylation/HG002_SNRPN_hp1.modkit.bed.gz"
},
{
"type": "BedTabixAdapter",
"name": "HP2",
"color": "#1f77b4",
"uri": "https://jbrowse.org/demos/methylation/HG002_SNRPN_hp2.modkit.bed.gz"
}
]
},
"displayDefaults": {
"mark": "bar",
"scales": { "y": { "domainMin": 0, "domainMax": 100, "title": "% 5mC" } }
}
}jbrowse add-track-json '{
"type": "MultiQuantitativeTrack",
"trackId": "HG002_snrpn_modkit_multi",
"name": "HG002 5mC by haplotype (modkit)",
"assemblyNames": ["hg38"],
"adapter": {
"type": "MultiWiggleAdapter",
"subadapters": [
{
"type": "BedTabixAdapter",
"name": "HP1",
"color": "#d62728",
"uri": "https://jbrowse.org/demos/methylation/HG002_SNRPN_hp1.modkit.bed.gz"
},
{
"type": "BedTabixAdapter",
"name": "HP2",
"color": "#1f77b4",
"uri": "https://jbrowse.org/demos/methylation/HG002_SNRPN_hp2.modkit.bed.gz"
}
]
},
"displayDefaults": {
"mark": "bar",
"scales": { "y": { "domainMin": 0, "domainMax": 100, "title": "% 5mC" } }
}
}'In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add track from pasted JSON, and paste:
{
"type": "MultiQuantitativeTrack",
"trackId": "HG002_snrpn_modkit_multi",
"name": "HG002 5mC by haplotype (modkit)",
"assemblyNames": ["hg38"],
"adapter": {
"type": "MultiWiggleAdapter",
"subadapters": [
{
"type": "BedTabixAdapter",
"name": "HP1",
"color": "#d62728",
"uri": "https://jbrowse.org/demos/methylation/HG002_SNRPN_hp1.modkit.bed.gz"
},
{
"type": "BedTabixAdapter",
"name": "HP2",
"color": "#1f77b4",
"uri": "https://jbrowse.org/demos/methylation/HG002_SNRPN_hp2.modkit.bed.gz"
}
]
},
"displayDefaults": {
"mark": "bar",
"scales": { "y": { "domainMin": 0, "domainMax": 100, "title": "% 5mC" } }
}
}Open this track in JBrowse Web ↗
Open this track in JBrowse Desktop ↗ (JBrowse Desktop 5.0+)
Reads of your own without an HP tag get one from whatshap haplotag, which
reads a phased VCF. The
alignments track guide
covers the Group-by dialog:
whatshap haplotag --reference reference.fa -o haplotagged.bam phased.vcf.gz reads.bam
samtools index haplotagged.bamSee also
- Tutorial: Methylation (bisulfite)
- Tutorial: Phased trio analysis (1000 Genomes)
- Tutorial: Comparing one genome's two haplotypes (T2T-HG002)
- Tutorial: RNA-seq visualization
- User guide: Alignments track
- User guide: Alignments track
- User guide: Quantitative track
External links
Notes
-
To compare two samples, run
modkit dmron their per-sample pileups and load its BED output as aFeatureTrackbeside the bedMethyl tracks, so the differentially-methylated regions line up with the positions driving them. ↩
Feedback on this tutorial is welcome: contact us.