Skip to main content

Methylation (long-read)

2 videos in this tutorial

JBrowse v5 beta

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.17. 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. At the SNRPN imprinting center, one parent's copy is methylated and the other's is not. In HG002 nanopore reads there we:

  • color each read by its modification calls
  • tally the calls into an aggregate track with modkit
  • group the reads by haplotype to pull the two parental alleles apart
  • smooth each allele into a moving average across the imprinted domain

Prerequisites

  • for your own data, long reads whose BAM or CRAM already has MM/ML modification 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 HP tag
pip install whatshap

The data

The files are region slices of public ONT open data on the ont-open-data S3 bucket.

The track configs below load region slices we host, cut from these files, so no download is needed. Reads of your own with MM and ML tags load by swapping uri, as Per-read methylation from the alignments shows.

The files
What it isFile
The HP1 bedMethyl (per-CpG methylation fractions) from the wf-human-variation sup run on HG002, restricted to m (5mC) rows and cut to the SNRPN locus, and to chr15:22,000,000-27,000,000 for the last sectionSAMPLE.wf_mods.1.bedmethyl.gz
The HP2 bedMethyl from that same runSAMPLE.wf_mods.2.bedmethyl.gz
The reads, from the HG002 sup basecalls, sliced to the same locus and haplotagged with whatshap haplotag against the phased SNP calls from that same wf-human-variation runPAO83395.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.gzi

JBrowse 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:

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 hg38

In 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 ↗

File → Open track... also opens the file by URL and infers the .bai.

Color by... → Per-base coloring → 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.

The modBAM opened by URL and colored from its track menu: Color by..., Modifications, and the two-color mode.

Open the session the video starts in ↗, then take the steps it shows:

  1. From File, Open track, then paste the BAM URL
  2. The form infers the .bai index and the adapter; Add
  3. Then color the reads by their modification calls
  4. Each read painted by its 5mC calls: red methylated, blue not
Open this view in JBrowse Web ↗ ·

The short route: open the view in JBrowse Web, where File → Open track... adds your own files beside its hosted genome and tracks.

From scratch with your own data, the steps behind the figure:

  1. Open your genome: on the JBrowse Desktop start screen click Open new genome (or Show all available genomes to pick a hosted one). For JBrowse Web, jbrowse add-assembly puts it in your config.json.This figure uses hg38, loaded from a bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  2. Add a track: check one your config has in the track selector, or File → Open track... to paste a URL or choose a local file.This figure uses “UCSC CpG islands”.
  3. Add the next track the same way. This one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “NCBI RefSeq (GFF3Tabix)”.
    1. Track menu → Gene glyph → Representative transcript
    2. Track menu → Feature height → Compact
  4. Type your region of interest into the location box and press Enter.This figure is at chr15:24,953,500-24,957,500.

This session spec draws the figure. It goes after &session=spec- on a JBrowse Web link that loads https://jbrowse.org/code/jb2/main/test_data/config_demo.json, which is what the Web button above opens.

{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,953,500-24,957,500",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        }
      ]
    }
  ]
}

Hand this to a coding agent. It rebuilds the figure above headlessly, and the session file is then the thing to edit: change the location or a track's settings, or add your own track under sessionTracks, and rerun.

cat > session.json <<'JSON'
{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,953,500-24,957,500",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        }
      ]
    }
  ]
}
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/code/jb2/main/test_data/config_demo.json \
  --width 1920 --height 740 \
  --spec session.json -o figure.png

Using JBrowse with AI agents covers the rest: hosted genomes to point it at, and why waiting for the render is the part that goes wrong.

HG002 ONT reads over the SNRPN CpG island in both modification color modes. The bottom mode fills every unmarked CpG in blue.
Open in JBrowse Web ↗

Or rebuild it from your own data:

The short route: open the view in JBrowse Web, where File → Open track... adds your own files beside its hosted genome and tracks.

From scratch with your own data, the steps behind the figure:

  1. Open your genome: on the JBrowse Desktop start screen click Open new genome (or Show all available genomes to pick a hosted one). For JBrowse Web, jbrowse add-assembly puts it in your config.json.This figure uses hg38, loaded from a bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  2. Add a track: check one your config has in the track selector, or File → Open track... to paste a URL or choose a local file.This figure uses “UCSC CpG islands”.
  3. Add the next track the same way. This one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “NCBI RefSeq (GFF3Tabix)”.
    1. Track menu → Gene glyph → Representative transcript
    2. Track menu → Feature height → Compact
  4. Add the next track the same way. This one needs a BAM file (.bam + .bai).This figure uses “HG002 ONT reads (5mC, haplotagged)”.
    1. Click Force load in the track's "Zoom in to see features or force load" message.Loads the region even past the byte-size limit, which can be slow.
    2. Track menu → Color by... → ModificationsNeeds MM/ML modification tags in your BAM/CRAM.
    3. Track menu → Color by... → Per-base coloring → Modifications → One color per modification type
  5. Type your region of interest into the location box and press Enter.This figure is at chr15:24,948,000-24,962,000.

This session spec draws the figure. It goes after &session=spec- on a JBrowse Web link that loads https://jbrowse.org/code/jb2/main/test_data/config_demo.json, which is what the Web button above opens.

{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,948,000-24,962,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_snrpn_5mC_reads",
          "height": 320,
          "forceLoad": true,
          "pileup": {
            "bases": {
              "show": true,
              "encoding": {
                "color": "modifications"
              },
              "unmodified": {
                "show": false
              }
            }
          }
        }
      ]
    }
  ]
}

Hand this to a coding agent. It rebuilds the figure above headlessly, and the session file is then the thing to edit: change the location or a track's settings, or add your own track under sessionTracks, and rerun.

The figure stacks 2 frames: one command draws each, and ImageMagick stacks them.

Frame 1: Modified only

cat > session-1.json <<'JSON'
{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,948,000-24,962,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_snrpn_5mC_reads",
          "height": 320,
          "forceLoad": true,
          "pileup": {
            "bases": {
              "show": true,
              "encoding": {
                "color": "modifications"
              },
              "unmodified": {
                "show": false
              }
            }
          }
        }
      ]
    }
  ]
}
JSON

cat > callouts-1.json <<'JSON'
[
  {
    "type": "text",
    "text": "Modified positions only",
    "anchor": {
      "text": "HG002 ONT reads",
      "alignX": "right",
      "dx": 90
    }
  }
]
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/code/jb2/main/test_data/config_demo.json \
  --width 1500 --height 742 \
  --timeout 90000 \
  --spec session-1.json --annotations callouts-1.json -o frame-1.png

Frame 2: Every CpG

cat > session-2.json <<'JSON'
{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,948,000-24,962,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_snrpn_5mC_reads",
          "height": 320,
          "forceLoad": true,
          "pileup": {
            "bases": {
              "show": true,
              "encoding": {
                "color": "modifications"
              },
              "unmodified": {
                "show": true
              }
            }
          }
        }
      ]
    }
  ]
}
JSON

cat > callouts-2.json <<'JSON'
[
  {
    "type": "text",
    "text": "Every CpG filled in",
    "anchor": {
      "text": "HG002 ONT reads",
      "alignX": "right",
      "dx": 90
    }
  }
]
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/code/jb2/main/test_data/config_demo.json \
  --width 1500 --height 742 \
  --timeout 90000 \
  --spec session-2.json --annotations callouts-2.json -o frame-2.png
magick frame-1.png frame-2.png -append figure.png

Using JBrowse with AI agents covers the rest: hosted genomes to point it at, and why waiting for the render is the part that goes wrong.

Open in JBrowse Web ↗

Or rebuild it from your own data:

The short route: open the view in JBrowse Web, where File → Open track... adds your own files beside its hosted genome and tracks.

From scratch with your own data, the steps behind the figure:

  1. Open your genome: on the JBrowse Desktop start screen click Open new genome (or Show all available genomes to pick a hosted one). For JBrowse Web, jbrowse add-assembly puts it in your config.json.This figure uses hg38, loaded from a bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  2. Add a track: check one your config has in the track selector, or File → Open track... to paste a URL or choose a local file.This figure uses “UCSC CpG islands”.
  3. Add the next track the same way. This one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “NCBI RefSeq (GFF3Tabix)”.
    1. Track menu → Gene glyph → Representative transcript
    2. Track menu → Feature height → Compact
  4. Add the next track the same way. This one needs a BAM file (.bam + .bai).This figure uses “HG002 ONT reads (5mC, haplotagged)”.
    1. Click Force load in the track's "Zoom in to see features or force load" message.Loads the region even past the byte-size limit, which can be slow.
    2. Track menu → Color by... → ModificationsNeeds MM/ML modification tags in your BAM/CRAM.
    3. Track menu → Color by... → Per-base coloring → Modifications → One color per type, plus low-probability & unmodified in blue
  5. Type your region of interest into the location box and press Enter.This figure is at chr15:24,948,000-24,962,000.

This session spec draws the figure. It goes after &session=spec- on a JBrowse Web link that loads https://jbrowse.org/code/jb2/main/test_data/config_demo.json, which is what the Web button above opens.

{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,948,000-24,962,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_snrpn_5mC_reads",
          "height": 320,
          "forceLoad": true,
          "pileup": {
            "bases": {
              "show": true,
              "encoding": {
                "color": "modifications"
              },
              "unmodified": {
                "show": true
              }
            }
          }
        }
      ]
    }
  ]
}

Hand this to a coding agent. It rebuilds the figure above headlessly, and the session file is then the thing to edit: change the location or a track's settings, or add your own track under sessionTracks, and rerun.

The figure stacks 2 frames: one command draws each, and ImageMagick stacks them.

Frame 1: Modified only

cat > session-1.json <<'JSON'
{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,948,000-24,962,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_snrpn_5mC_reads",
          "height": 320,
          "forceLoad": true,
          "pileup": {
            "bases": {
              "show": true,
              "encoding": {
                "color": "modifications"
              },
              "unmodified": {
                "show": false
              }
            }
          }
        }
      ]
    }
  ]
}
JSON

cat > callouts-1.json <<'JSON'
[
  {
    "type": "text",
    "text": "Modified positions only",
    "anchor": {
      "text": "HG002 ONT reads",
      "alignX": "right",
      "dx": 90
    }
  }
]
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/code/jb2/main/test_data/config_demo.json \
  --width 1500 --height 742 \
  --timeout 90000 \
  --spec session-1.json --annotations callouts-1.json -o frame-1.png

Frame 2: Every CpG

cat > session-2.json <<'JSON'
{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,948,000-24,962,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_snrpn_5mC_reads",
          "height": 320,
          "forceLoad": true,
          "pileup": {
            "bases": {
              "show": true,
              "encoding": {
                "color": "modifications"
              },
              "unmodified": {
                "show": true
              }
            }
          }
        }
      ]
    }
  ]
}
JSON

cat > callouts-2.json <<'JSON'
[
  {
    "type": "text",
    "text": "Every CpG filled in",
    "anchor": {
      "text": "HG002 ONT reads",
      "alignX": "right",
      "dx": 90
    }
  }
]
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/code/jb2/main/test_data/config_demo.json \
  --width 1500 --height 742 \
  --timeout 90000 \
  --spec session-2.json --annotations callouts-2.json -o frame-2.png
magick frame-1.png frame-2.png -append figure.png

Using JBrowse with AI agents covers the rest: hosted genomes to point it at, and why waiting for the render is the part that goes wrong.

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.gz

Two modkit pileup options shape the output:

  • --combine-mods collapses 5mC and 5hmC into a single fraction (bisulfite-equivalent). Pass --modified-bases 5mC 5hmC without it to keep separate rows per modification type (m for 5mC, h for 5hmC)
  • --phased writes one file per haplotype from the reads' HP tags, 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 --phased

bedMethyl 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 MultiQuantitativeTrack

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_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 ↗

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: none band 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.

Group by... → Tag... with HP, splitting the pileup by haplotype.

Open the session the video starts in ↗, then take the steps it shows:

  1. Both haplotypes' reads interleave; red is a methylated CpG, blue is not
  2. Group the reads by their HP haplotype tag
  3. HP 1 is methylated at the SNRPN promoter; HP 2 is not
Open this view in JBrowse Web ↗ ·

The short route: open the view in JBrowse Web, where File → Open track... adds your own files beside its hosted genome and tracks.

From scratch with your own data, the steps behind the figure:

  1. Open your genome: on the JBrowse Desktop start screen click Open new genome (or Show all available genomes to pick a hosted one). For JBrowse Web, jbrowse add-assembly puts it in your config.json.This figure uses hg38, loaded from a bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  2. Add a track: check one your config has in the track selector, or File → Open track... to paste a URL or choose a local file.This figure uses “UCSC CpG islands”.
  3. Add the next track the same way. This one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “NCBI RefSeq (GFF3Tabix)”.
    1. Track menu → Gene glyph → Representative transcript
    2. Track menu → Feature height → Compact
  4. Add the next track the same way. This one needs a BAM file (.bam + .bai).This figure uses “HG002 ONT reads (5mC, haplotagged)”.
    1. Click Force load in the track's "Zoom in to see features or force load" message.Loads the region even past the byte-size limit, which can be slow.
    2. Track menu → Color by... → ModificationsNeeds MM/ML modification tags in your BAM/CRAM.
    3. Track menu → Color by... → Per-base coloring → Modifications → One color per type, plus low-probability & unmodified in blue
  5. Type your region of interest into the location box and press Enter.This figure is at chr15:24,948,000-24,962,000.

This session spec draws the figure. It goes after &session=spec- on a JBrowse Web link that loads https://jbrowse.org/code/jb2/main/test_data/config_demo.json, which is what the Web button above opens.

{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,948,000-24,962,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_snrpn_5mC_reads",
          "height": 320,
          "forceLoad": true,
          "pileup": {
            "bases": {
              "show": true,
              "encoding": {
                "color": "modifications"
              },
              "unmodified": {
                "show": true
              }
            }
          }
        }
      ]
    }
  ]
}

Hand this to a coding agent. It rebuilds the figure above headlessly, and the session file is then the thing to edit: change the location or a track's settings, or add your own track under sessionTracks, and rerun.

cat > session.json <<'JSON'
{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,948,000-24,962,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_snrpn_5mC_reads",
          "height": 320,
          "forceLoad": true,
          "pileup": {
            "bases": {
              "show": true,
              "encoding": {
                "color": "modifications"
              },
              "unmodified": {
                "show": true
              }
            }
          }
        }
      ]
    }
  ]
}
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/code/jb2/main/test_data/config_demo.json \
  --width 1920 --height 740 \
  --spec session.json -o figure.png

Using JBrowse with AI agents covers the rest: hosted genomes to point it at, and why waiting for the render is the part that goes wrong.

HG002 ONT reads over the SNRPN CpG island colored by 5mC. Top: file order. Bottom: grouped by the HP tag.
Open in JBrowse Web ↗

Or rebuild it from your own data:

The short route: open the view in JBrowse Web, where File → Open track... adds your own files beside its hosted genome and tracks.

From scratch with your own data, the steps behind the figure:

  1. Open your genome: on the JBrowse Desktop start screen click Open new genome (or Show all available genomes to pick a hosted one). For JBrowse Web, jbrowse add-assembly puts it in your config.json.This figure uses hg38, loaded from a bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  2. Add a track: check one your config has in the track selector, or File → Open track... to paste a URL or choose a local file.This figure uses “UCSC CpG islands”.
  3. Add the next track the same way. This one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “NCBI RefSeq (GFF3Tabix)”.
    1. Track menu → Gene glyph → Representative transcript
    2. Track menu → Feature height → Compact
  4. Add the next track the same way. This one needs a BAM file (.bam + .bai).This figure uses “HG002 ONT reads (5mC, haplotagged)”.
    1. Click Force load in the track's "Zoom in to see features or force load" message.Loads the region even past the byte-size limit, which can be slow.
    2. Track menu → Color by... → ModificationsNeeds MM/ML modification tags in your BAM/CRAM.
    3. Track menu → Color by... → Per-base coloring → Modifications → One color per type, plus low-probability & unmodified in blue
  5. Type your region of interest into the location box and press Enter.This figure is at chr15:24,948,000-24,962,000.

This session spec draws the figure. It goes after &session=spec- on a JBrowse Web link that loads https://jbrowse.org/code/jb2/main/test_data/config_demo.json, which is what the Web button above opens.

{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,948,000-24,962,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_snrpn_5mC_reads",
          "height": 320,
          "forceLoad": true,
          "pileup": {
            "bases": {
              "show": true,
              "encoding": {
                "color": "modifications"
              },
              "unmodified": {
                "show": true
              }
            }
          }
        }
      ]
    }
  ]
}

Hand this to a coding agent. It rebuilds the figure above headlessly, and the session file is then the thing to edit: change the location or a track's settings, or add your own track under sessionTracks, and rerun.

The figure stacks 2 frames: one command draws each, and ImageMagick stacks them.

Frame 1: Ungrouped

cat > session-1.json <<'JSON'
{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,948,000-24,962,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_snrpn_5mC_reads",
          "height": 320,
          "forceLoad": true,
          "pileup": {
            "bases": {
              "show": true,
              "encoding": {
                "color": "modifications"
              },
              "unmodified": {
                "show": true
              }
            }
          }
        }
      ]
    }
  ]
}
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/code/jb2/main/test_data/config_demo.json \
  --width 1500 --height 730 \
  --timeout 90000 \
  --spec session-1.json -o frame-1.png

The figure was taken after clicks this command does not make, so its picture can differ.

It leaves out 2 callouts anchored to what the clicks bring up.

Frame 2: Grouped by HP

cat > session-2.json <<'JSON'
{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,948,000-24,962,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_snrpn_5mC_reads",
          "height": 320,
          "forceLoad": true,
          "facet": "tags.HP",
          "pileup": {
            "bases": {
              "show": true,
              "encoding": {
                "color": "modifications"
              },
              "unmodified": {
                "show": true
              }
            }
          }
        }
      ]
    }
  ]
}
JSON

cat > callouts-2.json <<'JSON'
[
  {
    "type": "text",
    "text": "HP1",
    "anchor": {
      "text": "HP: 1",
      "alignX": "right",
      "dx": 90,
      "dy": 18
    }
  },
  {
    "type": "text",
    "text": "HP2",
    "anchor": {
      "text": "HP: 2",
      "alignX": "right",
      "dx": 90
    }
  },
  {
    "type": "text",
    "text": "HP unknown",
    "anchor": {
      "text": "HP: none",
      "alignX": "right",
      "dx": 90
    }
  }
]
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/code/jb2/main/test_data/config_demo.json \
  --width 1500 --height 730 \
  --timeout 90000 \
  --spec session-2.json --annotations callouts-2.json -o frame-2.png
magick frame-1.png frame-2.png -append figure.png

Using JBrowse with AI agents covers the rest: hosted genomes to point it at, and why waiting for the render is the part that goes wrong.

Open in JBrowse Web ↗

Or rebuild it from your own data:

The short route: open the view in JBrowse Web, where File → Open track... adds your own files beside its hosted genome and tracks.

From scratch with your own data, the steps behind the figure:

  1. Open your genome: on the JBrowse Desktop start screen click Open new genome (or Show all available genomes to pick a hosted one). For JBrowse Web, jbrowse add-assembly puts it in your config.json.This figure uses hg38, loaded from a bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  2. Add a track: check one your config has in the track selector, or File → Open track... to paste a URL or choose a local file.This figure uses “UCSC CpG islands”.
  3. Add the next track the same way. This one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “NCBI RefSeq (GFF3Tabix)”.
    1. Track menu → Gene glyph → Representative transcript
    2. Track menu → Feature height → Compact
  4. Add the next track the same way. This one needs a BAM file (.bam + .bai).This figure uses “HG002 ONT reads (5mC, haplotagged)”.
    1. Click Force load in the track's "Zoom in to see features or force load" message.Loads the region even past the byte-size limit, which can be slow.
    2. Track menu → Group by... → Tag... → enter "HP"HP is the haplotype tag written by phasing tools like WhatsHap or Longphase.
    3. Track menu → Color by... → ModificationsNeeds MM/ML modification tags in your BAM/CRAM.
    4. Track menu → Color by... → Per-base coloring → Modifications → One color per type, plus low-probability & unmodified in blue
  5. Type your region of interest into the location box and press Enter.This figure is at chr15:24,948,000-24,962,000.

This session spec draws the figure. It goes after &session=spec- on a JBrowse Web link that loads https://jbrowse.org/code/jb2/main/test_data/config_demo.json, which is what the Web button above opens.

{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,948,000-24,962,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_snrpn_5mC_reads",
          "height": 320,
          "forceLoad": true,
          "facet": "tags.HP",
          "pileup": {
            "bases": {
              "show": true,
              "encoding": {
                "color": "modifications"
              },
              "unmodified": {
                "show": true
              }
            }
          }
        }
      ]
    }
  ]
}

Hand this to a coding agent. It rebuilds the figure above headlessly, and the session file is then the thing to edit: change the location or a track's settings, or add your own track under sessionTracks, and rerun.

The figure stacks 2 frames: one command draws each, and ImageMagick stacks them.

Frame 1: Ungrouped

cat > session-1.json <<'JSON'
{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,948,000-24,962,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_snrpn_5mC_reads",
          "height": 320,
          "forceLoad": true,
          "pileup": {
            "bases": {
              "show": true,
              "encoding": {
                "color": "modifications"
              },
              "unmodified": {
                "show": true
              }
            }
          }
        }
      ]
    }
  ]
}
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/code/jb2/main/test_data/config_demo.json \
  --width 1500 --height 730 \
  --timeout 90000 \
  --spec session-1.json -o frame-1.png

The figure was taken after clicks this command does not make, so its picture can differ.

It leaves out 2 callouts anchored to what the clicks bring up.

Frame 2: Grouped by HP

cat > session-2.json <<'JSON'
{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,948,000-24,962,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_snrpn_5mC_reads",
          "height": 320,
          "forceLoad": true,
          "facet": "tags.HP",
          "pileup": {
            "bases": {
              "show": true,
              "encoding": {
                "color": "modifications"
              },
              "unmodified": {
                "show": true
              }
            }
          }
        }
      ]
    }
  ]
}
JSON

cat > callouts-2.json <<'JSON'
[
  {
    "type": "text",
    "text": "HP1",
    "anchor": {
      "text": "HP: 1",
      "alignX": "right",
      "dx": 90,
      "dy": 18
    }
  },
  {
    "type": "text",
    "text": "HP2",
    "anchor": {
      "text": "HP: 2",
      "alignX": "right",
      "dx": 90
    }
  },
  {
    "type": "text",
    "text": "HP unknown",
    "anchor": {
      "text": "HP: none",
      "alignX": "right",
      "dx": 90
    }
  }
]
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/code/jb2/main/test_data/config_demo.json \
  --width 1500 --height 730 \
  --timeout 90000 \
  --spec session-2.json --annotations callouts-2.json -o frame-2.png
magick frame-1.png frame-2.png -append figure.png

Using JBrowse with AI agents covers the rest: hosted genomes to point it at, and why waiting for the render is the part that goes wrong.

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",
    "encoding": {
      "y": {
        "scale": { "domainMin": 0, "domainMax": 100 },
        "axis": { "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",
    "encoding": {
      "y": {
        "scale": { "domainMin": 0, "domainMax": 100 },
        "axis": { "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",
    "encoding": {
      "y": {
        "scale": { "domainMin": 0, "domainMax": 100 },
        "axis": { "title": "% 5mC" }
      }
    }
  }
}

Open this track in JBrowse Web ↗

The SNRPN CpG island: the per-haplotype profile on top, the reads grouped by HP below.
Open in JBrowse Web ↗

Or rebuild it from your own data:

The short route: open the view in JBrowse Web, where File → Open track... adds your own files beside its hosted genome and tracks.

From scratch with your own data, the steps behind the figure:

  1. Open your genome: on the JBrowse Desktop start screen click Open new genome (or Show all available genomes to pick a hosted one). For JBrowse Web, jbrowse add-assembly puts it in your config.json.This figure uses hg38, loaded from a bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  2. Add a track: check one your config has in the track selector, or File → Open track... to paste a URL or choose a local file.This figure uses “UCSC CpG islands”.
  3. Add the next track the same way. This one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “NCBI RefSeq (GFF3Tabix)”.
    1. Track menu → Gene glyph → Representative transcript
    2. Track menu → Feature height → Compact
  4. Add the next track the same way.This figure uses “HG002 5mC by haplotype (modkit)”.
  5. Add the next track the same way. This one needs a BAM file (.bam + .bai).This figure uses “HG002 ONT reads (5mC, haplotagged)”.
    1. Click Force load in the track's "Zoom in to see features or force load" message.Loads the region even past the byte-size limit, which can be slow.
    2. Track menu → Group by... → Tag... → enter "HP"HP is the haplotype tag written by phasing tools like WhatsHap or Longphase.
    3. Track menu → Color by... → ModificationsNeeds MM/ML modification tags in your BAM/CRAM.
    4. Track menu → Color by... → Per-base coloring → Modifications → One color per type, plus low-probability & unmodified in blue
  6. Type your region of interest into the location box and press Enter.This figure is at chr15:24,948,000-24,962,000.

A figure adds its tracks to one session. The jbrowse CLI writes the same tracks into a config.json instead, where every session that opens it has them: add-track where flags cover the whole track, add-track-json where they do not.

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",
    "encoding": {
      "y": {
        "scale": {
          "domainMin": 0,
          "domainMax": 100
        },
        "axis": {
          "title": "% 5mC"
        }
      }
    }
  }
}'

Run these where the config.json is, or add --out <dir>, and point each uri at your own file. The config needs hg38 already — jbrowse add-assembly adds yours. The location and the settings the steps carry are session state rather than track config — that half is the Spec tab, or jbrowse set-default-session.

This session spec draws the figure. It goes after &session=spec- on a JBrowse Web link that loads https://jbrowse.org/code/jb2/main/test_data/config_demo.json, which is what the Web button above opens.

{
  "sessionTracks": [
    {
      "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",
        "encoding": {
          "y": {
            "scale": {
              "domainMin": 0,
              "domainMax": 100
            },
            "axis": {
              "title": "% 5mC"
            }
          }
        }
      }
    }
  ],
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,948,000-24,962,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_snrpn_modkit_multi",
          "height": 170
        },
        {
          "trackId": "HG002_snrpn_5mC_reads",
          "height": 460,
          "forceLoad": true,
          "facet": "tags.HP",
          "pileup": {
            "bases": {
              "show": true,
              "encoding": {
                "color": "modifications"
              },
              "unmodified": {
                "show": true
              }
            }
          }
        }
      ]
    }
  ]
}

Hand this to a coding agent. It rebuilds the figure above headlessly, and the session file is then the thing to edit: change the location or a track's settings, or add your own track under sessionTracks, and rerun.

cat > session.json <<'JSON'
{
  "sessionTracks": [
    {
      "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",
        "encoding": {
          "y": {
            "scale": {
              "domainMin": 0,
              "domainMax": 100
            },
            "axis": {
              "title": "% 5mC"
            }
          }
        }
      }
    }
  ],
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,948,000-24,962,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_snrpn_modkit_multi",
          "height": 170
        },
        {
          "trackId": "HG002_snrpn_5mC_reads",
          "height": 460,
          "forceLoad": true,
          "facet": "tags.HP",
          "pileup": {
            "bases": {
              "show": true,
              "encoding": {
                "color": "modifications"
              },
              "unmodified": {
                "show": true
              }
            }
          }
        }
      ]
    }
  ]
}
JSON

cat > callouts.json <<'JSON'
[
  {
    "type": "text",
    "text": "HP1",
    "anchor": {
      "text": "HP: 1",
      "alignX": "right",
      "dx": 90,
      "dy": 18
    }
  },
  {
    "type": "text",
    "text": "HP2",
    "anchor": {
      "text": "HP: 2",
      "alignX": "right",
      "dx": 90
    }
  },
  {
    "type": "text",
    "text": "HP unknown",
    "anchor": {
      "text": "HP: none",
      "alignX": "right",
      "dx": 90
    }
  }
]
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/code/jb2/main/test_data/config_demo.json \
  --width 1500 --height 1085 \
  --timeout 90000 \
  --spec session.json --annotations callouts.json -o figure.png

Using JBrowse with AI agents covers the rest: hosted genomes to point it at, and why waiting for the render is the part that goes wrong.

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.bam

Smoothing across the imprinted domain

A single CpG's fraction rests on the few reads covering it, so per-site values scatter from 0 to 100 and a difference spread over many sites is hard to see. A window step averages each CpG with a fixed number of neighbouring CpGs, so a CpG island weighs as many sites as it holds. We'll load both haplotypes around SNRPN, each in a row with its raw values under a moving average:

Goes in the tracks array of config.json. See Tracks.

{
  "type": "MultiQuantitativeTrack",
  "trackId": "HG002_imprinted_5mC_smoothed",
  "name": "HG002 5mC by haplotype, moving average",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "subadapters": [
      {
        "type": "BedTabixAdapter",
        "name": "HP1",
        "color": "#d62728",
        "uri": "https://jbrowse.org/demos/methylation/HG002_chr15_imprinted_hp1.modkit.bed.gz"
      },
      {
        "type": "BedTabixAdapter",
        "name": "HP2",
        "color": "#1f77b4",
        "uri": "https://jbrowse.org/demos/methylation/HG002_chr15_imprinted_hp2.modkit.bed.gz"
      }
    ]
  },
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "displayId": "HG002_imprinted_5mC_smoothed-LinearMarkDisplay",
      "rows": "source",
      "legend": { "show": false },
      "fetchSizeLimit": 5000000,
      "encoding": {
        "y": {
          "scale": { "domainMin": 0, "domainMax": 100 },
          "axis": { "title": "% 5mC" }
        }
      },
      "marks": [
        {
          "mark": "point",
          "encoding": {
            "y": "score",
            "color": { "value": "#bbbbbb" },
            "size": 2
          }
        },
        {
          "mark": { "type": "line", "interpolate": "linear" },
          "transform": [
            {
              "type": "window",
              "groupby": ["source"],
              "frame": [-50, 50],
              "ops": [
                {
                  "op": "mean",
                  "field": "score",
                  "weight": "n_valid_cov",
                  "as": "meth"
                }
              ]
            }
          ],
          "encoding": {
            "y": "meth",
            "color": {
              "field": "source",
              "scale": {
                "domain": ["HP1", "HP2"],
                "range": ["#d62728", "#1f77b4"]
              }
            },
            "strokeWidth": 2
          }
        }
      ]
    }
  ]
}
jbrowse add-track-json '{
  "type": "MultiQuantitativeTrack",
  "trackId": "HG002_imprinted_5mC_smoothed",
  "name": "HG002 5mC by haplotype, moving average",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "subadapters": [
      {
        "type": "BedTabixAdapter",
        "name": "HP1",
        "color": "#d62728",
        "uri": "https://jbrowse.org/demos/methylation/HG002_chr15_imprinted_hp1.modkit.bed.gz"
      },
      {
        "type": "BedTabixAdapter",
        "name": "HP2",
        "color": "#1f77b4",
        "uri": "https://jbrowse.org/demos/methylation/HG002_chr15_imprinted_hp2.modkit.bed.gz"
      }
    ]
  },
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "displayId": "HG002_imprinted_5mC_smoothed-LinearMarkDisplay",
      "rows": "source",
      "legend": { "show": false },
      "fetchSizeLimit": 5000000,
      "encoding": {
        "y": {
          "scale": { "domainMin": 0, "domainMax": 100 },
          "axis": { "title": "% 5mC" }
        }
      },
      "marks": [
        {
          "mark": "point",
          "encoding": {
            "y": "score",
            "color": { "value": "#bbbbbb" },
            "size": 2
          }
        },
        {
          "mark": { "type": "line", "interpolate": "linear" },
          "transform": [
            {
              "type": "window",
              "groupby": ["source"],
              "frame": [-50, 50],
              "ops": [
                {
                  "op": "mean",
                  "field": "score",
                  "weight": "n_valid_cov",
                  "as": "meth"
                }
              ]
            }
          ],
          "encoding": {
            "y": "meth",
            "color": {
              "field": "source",
              "scale": {
                "domain": ["HP1", "HP2"],
                "range": ["#d62728", "#1f77b4"]
              }
            },
            "strokeWidth": 2
          }
        }
      ]
    }
  ]
}'

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_imprinted_5mC_smoothed",
  "name": "HG002 5mC by haplotype, moving average",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "subadapters": [
      {
        "type": "BedTabixAdapter",
        "name": "HP1",
        "color": "#d62728",
        "uri": "https://jbrowse.org/demos/methylation/HG002_chr15_imprinted_hp1.modkit.bed.gz"
      },
      {
        "type": "BedTabixAdapter",
        "name": "HP2",
        "color": "#1f77b4",
        "uri": "https://jbrowse.org/demos/methylation/HG002_chr15_imprinted_hp2.modkit.bed.gz"
      }
    ]
  },
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "displayId": "HG002_imprinted_5mC_smoothed-LinearMarkDisplay",
      "rows": "source",
      "legend": { "show": false },
      "fetchSizeLimit": 5000000,
      "encoding": {
        "y": {
          "scale": { "domainMin": 0, "domainMax": 100 },
          "axis": { "title": "% 5mC" }
        }
      },
      "marks": [
        {
          "mark": "point",
          "encoding": {
            "y": "score",
            "color": { "value": "#bbbbbb" },
            "size": 2
          }
        },
        {
          "mark": { "type": "line", "interpolate": "linear" },
          "transform": [
            {
              "type": "window",
              "groupby": ["source"],
              "frame": [-50, 50],
              "ops": [
                {
                  "op": "mean",
                  "field": "score",
                  "weight": "n_valid_cov",
                  "as": "meth"
                }
              ]
            }
          ],
          "encoding": {
            "y": "meth",
            "color": {
              "field": "source",
              "scale": {
                "domain": ["HP1", "HP2"],
                "range": ["#d62728", "#1f77b4"]
              }
            },
            "strokeWidth": 2
          }
        }
      ]
    }
  ]
}

Open this track in JBrowse Web ↗

The display draws two marks over the same rows:

  • the point mark plots each CpG's score, its percent methylated, in grey
  • the line mark plots meth, which the window step writes: each CpG's mean over the 50 CpGs on either side, weighted by n_valid_cov, the read depth modkit counted there
  • with groupby: ["source"] the step averages each haplotype over that haplotype's CpGs, and with rows: "source" the display gives each haplotype a row whose label names it, so legend.show is off
  • with fetchSizeLimit raised, the track reads the whole slice when zoomed out

Type chr15:24,800,000-25,200,000 into the location box to see the SNRPN promoter and the stretch of SNHG14 downstream of it:

HG002 5mC by haplotype over SNRPN and SNHG14, each CpG in grey under its moving average.
Open in JBrowse Web ↗

Or rebuild it from your own data:

The short route: open the view in JBrowse Web, where File → Open track... adds your own files beside its hosted genome and tracks.

From scratch with your own data, the steps behind the figure:

  1. Open your genome: on the JBrowse Desktop start screen click Open new genome (or Show all available genomes to pick a hosted one). For JBrowse Web, jbrowse add-assembly puts it in your config.json.This figure uses hg38, loaded from a bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  2. Add a track: check one your config has in the track selector, or File → Open track... to paste a URL or choose a local file.This figure uses “UCSC CpG islands”.
  3. Add the next track the same way. This one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “NCBI RefSeq (GFF3Tabix)”.
    1. Track menu → Gene glyph → Representative transcript
    2. Track menu → Feature height → Compact
  4. Add the next track the same way.This figure uses “HG002 5mC by haplotype, moving average”.
    1. Click Force load in the track's "Zoom in to see features or force load" message.Loads the region even past the byte-size limit, which can be slow.
  5. Type your region of interest into the location box and press Enter.This figure is at chr15:24,800,000-25,200,000.

A figure adds its tracks to one session. The jbrowse CLI writes the same tracks into a config.json instead, where every session that opens it has them: add-track where flags cover the whole track, add-track-json where they do not.

jbrowse add-track-json '{
  "type": "MultiQuantitativeTrack",
  "trackId": "HG002_imprinted_5mC_smoothed",
  "name": "HG002 5mC by haplotype, moving average",
  "assemblyNames": [
    "hg38"
  ],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "subadapters": [
      {
        "type": "BedTabixAdapter",
        "name": "HP1",
        "color": "#d62728",
        "uri": "https://jbrowse.org/demos/methylation/HG002_chr15_imprinted_hp1.modkit.bed.gz"
      },
      {
        "type": "BedTabixAdapter",
        "name": "HP2",
        "color": "#1f77b4",
        "uri": "https://jbrowse.org/demos/methylation/HG002_chr15_imprinted_hp2.modkit.bed.gz"
      }
    ]
  },
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "displayId": "HG002_imprinted_5mC_smoothed-LinearMarkDisplay",
      "rows": "source",
      "legend": {
        "show": false
      },
      "fetchSizeLimit": 5000000,
      "encoding": {
        "y": {
          "scale": {
            "domainMin": 0,
            "domainMax": 100
          },
          "axis": {
            "title": "% 5mC"
          }
        }
      },
      "marks": [
        {
          "mark": "point",
          "encoding": {
            "y": "score",
            "color": {
              "value": "#bbbbbb"
            },
            "size": 2
          }
        },
        {
          "mark": {
            "type": "line",
            "interpolate": "linear"
          },
          "transform": [
            {
              "type": "window",
              "groupby": [
                "source"
              ],
              "frame": [
                -50,
                50
              ],
              "ops": [
                {
                  "op": "mean",
                  "field": "score",
                  "weight": "n_valid_cov",
                  "as": "meth"
                }
              ]
            }
          ],
          "encoding": {
            "y": "meth",
            "color": {
              "field": "source",
              "scale": {
                "domain": [
                  "HP1",
                  "HP2"
                ],
                "range": [
                  "#d62728",
                  "#1f77b4"
                ]
              }
            },
            "strokeWidth": 2
          }
        }
      ]
    }
  ]
}'

Run these where the config.json is, or add --out <dir>, and point each uri at your own file. The config needs hg38 already — jbrowse add-assembly adds yours. The location and the settings the steps carry are session state rather than track config — that half is the Spec tab, or jbrowse set-default-session.

This session spec draws the figure. It goes after &session=spec- on a JBrowse Web link that loads https://jbrowse.org/code/jb2/main/test_data/config_demo.json, which is what the Web button above opens.

{
  "sessionTracks": [
    {
      "type": "MultiQuantitativeTrack",
      "trackId": "HG002_imprinted_5mC_smoothed",
      "name": "HG002 5mC by haplotype, moving average",
      "assemblyNames": [
        "hg38"
      ],
      "adapter": {
        "type": "MultiWiggleAdapter",
        "subadapters": [
          {
            "type": "BedTabixAdapter",
            "name": "HP1",
            "color": "#d62728",
            "uri": "https://jbrowse.org/demos/methylation/HG002_chr15_imprinted_hp1.modkit.bed.gz"
          },
          {
            "type": "BedTabixAdapter",
            "name": "HP2",
            "color": "#1f77b4",
            "uri": "https://jbrowse.org/demos/methylation/HG002_chr15_imprinted_hp2.modkit.bed.gz"
          }
        ]
      },
      "displays": [
        {
          "type": "LinearMarkDisplay",
          "displayId": "HG002_imprinted_5mC_smoothed-LinearMarkDisplay",
          "rows": "source",
          "legend": {
            "show": false
          },
          "fetchSizeLimit": 5000000,
          "encoding": {
            "y": {
              "scale": {
                "domainMin": 0,
                "domainMax": 100
              },
              "axis": {
                "title": "% 5mC"
              }
            }
          },
          "marks": [
            {
              "mark": "point",
              "encoding": {
                "y": "score",
                "color": {
                  "value": "#bbbbbb"
                },
                "size": 2
              }
            },
            {
              "mark": {
                "type": "line",
                "interpolate": "linear"
              },
              "transform": [
                {
                  "type": "window",
                  "groupby": [
                    "source"
                  ],
                  "frame": [
                    -50,
                    50
                  ],
                  "ops": [
                    {
                      "op": "mean",
                      "field": "score",
                      "weight": "n_valid_cov",
                      "as": "meth"
                    }
                  ]
                }
              ],
              "encoding": {
                "y": "meth",
                "color": {
                  "field": "source",
                  "scale": {
                    "domain": [
                      "HP1",
                      "HP2"
                    ],
                    "range": [
                      "#d62728",
                      "#1f77b4"
                    ]
                  }
                },
                "strokeWidth": 2
              }
            }
          ]
        }
      ]
    }
  ],
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,800,000-25,200,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_imprinted_5mC_smoothed",
          "height": 300,
          "forceLoad": true
        }
      ]
    }
  ]
}

Hand this to a coding agent. It rebuilds the figure above headlessly, and the session file is then the thing to edit: change the location or a track's settings, or add your own track under sessionTracks, and rerun.

cat > session.json <<'JSON'
{
  "sessionTracks": [
    {
      "type": "MultiQuantitativeTrack",
      "trackId": "HG002_imprinted_5mC_smoothed",
      "name": "HG002 5mC by haplotype, moving average",
      "assemblyNames": [
        "hg38"
      ],
      "adapter": {
        "type": "MultiWiggleAdapter",
        "subadapters": [
          {
            "type": "BedTabixAdapter",
            "name": "HP1",
            "color": "#d62728",
            "uri": "https://jbrowse.org/demos/methylation/HG002_chr15_imprinted_hp1.modkit.bed.gz"
          },
          {
            "type": "BedTabixAdapter",
            "name": "HP2",
            "color": "#1f77b4",
            "uri": "https://jbrowse.org/demos/methylation/HG002_chr15_imprinted_hp2.modkit.bed.gz"
          }
        ]
      },
      "displays": [
        {
          "type": "LinearMarkDisplay",
          "displayId": "HG002_imprinted_5mC_smoothed-LinearMarkDisplay",
          "rows": "source",
          "legend": {
            "show": false
          },
          "fetchSizeLimit": 5000000,
          "encoding": {
            "y": {
              "scale": {
                "domainMin": 0,
                "domainMax": 100
              },
              "axis": {
                "title": "% 5mC"
              }
            }
          },
          "marks": [
            {
              "mark": "point",
              "encoding": {
                "y": "score",
                "color": {
                  "value": "#bbbbbb"
                },
                "size": 2
              }
            },
            {
              "mark": {
                "type": "line",
                "interpolate": "linear"
              },
              "transform": [
                {
                  "type": "window",
                  "groupby": [
                    "source"
                  ],
                  "frame": [
                    -50,
                    50
                  ],
                  "ops": [
                    {
                      "op": "mean",
                      "field": "score",
                      "weight": "n_valid_cov",
                      "as": "meth"
                    }
                  ]
                }
              ],
              "encoding": {
                "y": "meth",
                "color": {
                  "field": "source",
                  "scale": {
                    "domain": [
                      "HP1",
                      "HP2"
                    ],
                    "range": [
                      "#d62728",
                      "#1f77b4"
                    ]
                  }
                },
                "strokeWidth": 2
              }
            }
          ]
        }
      ]
    }
  ],
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:24,800,000-25,200,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_imprinted_5mC_smoothed",
          "height": 300,
          "forceLoad": true
        }
      ]
    }
  ]
}
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/code/jb2/main/test_data/config_demo.json \
  --width 1500 --height 725 \
  --timeout 90000 \
  --spec session.json -o figure.png

Using JBrowse with AI agents covers the rest: hosted genomes to point it at, and why waiting for the render is the part that goes wrong.

Edit marks... in the track menu opens the marks. In the line's window step, set rows before and rows after to 500 for a 1,001-CpG average, then type chr15:22,200,000-26,900,000 to take in the whole imprinted domain:

The same track across the Prader-Willi and Angelman domain with a wider moving average.
Open in JBrowse Web ↗

Or rebuild it from your own data:

The short route: open the view in JBrowse Web, where File → Open track... adds your own files beside its hosted genome and tracks.

From scratch with your own data, the steps behind the figure:

  1. Open your genome: on the JBrowse Desktop start screen click Open new genome (or Show all available genomes to pick a hosted one). For JBrowse Web, jbrowse add-assembly puts it in your config.json.This figure uses hg38, loaded from a bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
  2. Add a track: check one your config has in the track selector, or File → Open track... to paste a URL or choose a local file.This figure uses “UCSC CpG islands”.
  3. Add the next track the same way. This one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “NCBI RefSeq (GFF3Tabix)”.
    1. Track menu → Gene glyph → Representative transcript
    2. Track menu → Feature height → Compact
  4. Add the next track the same way.This figure uses “HG002 5mC by haplotype, moving average”.
    1. Click Force load in the track's "Zoom in to see features or force load" message.Loads the region even past the byte-size limit, which can be slow.
    2. Track menu → Edit marks... → Add markThis figure draws point marks. The form picks a mark's type and maps its channels to fields.
  5. Type your region of interest into the location box and press Enter.This figure is at chr15:22,200,000-26,900,000.

A figure adds its tracks to one session. The jbrowse CLI writes the same tracks into a config.json instead, where every session that opens it has them: add-track where flags cover the whole track, add-track-json where they do not.

jbrowse add-track-json '{
  "type": "MultiQuantitativeTrack",
  "trackId": "HG002_imprinted_5mC_smoothed",
  "name": "HG002 5mC by haplotype, moving average",
  "assemblyNames": [
    "hg38"
  ],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "subadapters": [
      {
        "type": "BedTabixAdapter",
        "name": "HP1",
        "color": "#d62728",
        "uri": "https://jbrowse.org/demos/methylation/HG002_chr15_imprinted_hp1.modkit.bed.gz"
      },
      {
        "type": "BedTabixAdapter",
        "name": "HP2",
        "color": "#1f77b4",
        "uri": "https://jbrowse.org/demos/methylation/HG002_chr15_imprinted_hp2.modkit.bed.gz"
      }
    ]
  },
  "displays": [
    {
      "type": "LinearMarkDisplay",
      "displayId": "HG002_imprinted_5mC_smoothed-LinearMarkDisplay",
      "rows": "source",
      "legend": {
        "show": false
      },
      "fetchSizeLimit": 5000000,
      "encoding": {
        "y": {
          "scale": {
            "domainMin": 0,
            "domainMax": 100
          },
          "axis": {
            "title": "% 5mC"
          }
        }
      },
      "marks": [
        {
          "mark": "point",
          "encoding": {
            "y": "score",
            "color": {
              "value": "#bbbbbb"
            },
            "size": 2
          }
        },
        {
          "mark": {
            "type": "line",
            "interpolate": "linear"
          },
          "transform": [
            {
              "type": "window",
              "groupby": [
                "source"
              ],
              "frame": [
                -50,
                50
              ],
              "ops": [
                {
                  "op": "mean",
                  "field": "score",
                  "weight": "n_valid_cov",
                  "as": "meth"
                }
              ]
            }
          ],
          "encoding": {
            "y": "meth",
            "color": {
              "field": "source",
              "scale": {
                "domain": [
                  "HP1",
                  "HP2"
                ],
                "range": [
                  "#d62728",
                  "#1f77b4"
                ]
              }
            },
            "strokeWidth": 2
          }
        }
      ]
    }
  ]
}'

Run these where the config.json is, or add --out <dir>, and point each uri at your own file. The config needs hg38 already — jbrowse add-assembly adds yours. The location and the settings the steps carry are session state rather than track config — that half is the Spec tab, or jbrowse set-default-session.

This session spec draws the figure. It goes after &session=spec- on a JBrowse Web link that loads https://jbrowse.org/code/jb2/main/test_data/config_demo.json, which is what the Web button above opens.

{
  "sessionTracks": [
    {
      "type": "MultiQuantitativeTrack",
      "trackId": "HG002_imprinted_5mC_smoothed",
      "name": "HG002 5mC by haplotype, moving average",
      "assemblyNames": [
        "hg38"
      ],
      "adapter": {
        "type": "MultiWiggleAdapter",
        "subadapters": [
          {
            "type": "BedTabixAdapter",
            "name": "HP1",
            "color": "#d62728",
            "uri": "https://jbrowse.org/demos/methylation/HG002_chr15_imprinted_hp1.modkit.bed.gz"
          },
          {
            "type": "BedTabixAdapter",
            "name": "HP2",
            "color": "#1f77b4",
            "uri": "https://jbrowse.org/demos/methylation/HG002_chr15_imprinted_hp2.modkit.bed.gz"
          }
        ]
      },
      "displays": [
        {
          "type": "LinearMarkDisplay",
          "displayId": "HG002_imprinted_5mC_smoothed-LinearMarkDisplay",
          "rows": "source",
          "legend": {
            "show": false
          },
          "fetchSizeLimit": 5000000,
          "encoding": {
            "y": {
              "scale": {
                "domainMin": 0,
                "domainMax": 100
              },
              "axis": {
                "title": "% 5mC"
              }
            }
          },
          "marks": [
            {
              "mark": "point",
              "encoding": {
                "y": "score",
                "color": {
                  "value": "#bbbbbb"
                },
                "size": 2
              }
            },
            {
              "mark": {
                "type": "line",
                "interpolate": "linear"
              },
              "transform": [
                {
                  "type": "window",
                  "groupby": [
                    "source"
                  ],
                  "frame": [
                    -50,
                    50
                  ],
                  "ops": [
                    {
                      "op": "mean",
                      "field": "score",
                      "weight": "n_valid_cov",
                      "as": "meth"
                    }
                  ]
                }
              ],
              "encoding": {
                "y": "meth",
                "color": {
                  "field": "source",
                  "scale": {
                    "domain": [
                      "HP1",
                      "HP2"
                    ],
                    "range": [
                      "#d62728",
                      "#1f77b4"
                    ]
                  }
                },
                "strokeWidth": 2
              }
            }
          ]
        }
      ]
    }
  ],
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:22,200,000-26,900,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_imprinted_5mC_smoothed",
          "height": 300,
          "forceLoad": true,
          "marks": [
            {
              "mark": "point",
              "encoding": {
                "y": "score",
                "color": {
                  "value": "#bbbbbb"
                },
                "size": 2
              }
            },
            {
              "mark": {
                "type": "line",
                "interpolate": "linear"
              },
              "transform": [
                {
                  "type": "window",
                  "groupby": [
                    "source"
                  ],
                  "frame": [
                    -500,
                    500
                  ],
                  "ops": [
                    {
                      "op": "mean",
                      "field": "score",
                      "weight": "n_valid_cov",
                      "as": "meth"
                    }
                  ]
                }
              ],
              "encoding": {
                "y": "meth",
                "color": {
                  "field": "source",
                  "scale": {
                    "domain": [
                      "HP1",
                      "HP2"
                    ],
                    "range": [
                      "#d62728",
                      "#1f77b4"
                    ]
                  }
                },
                "strokeWidth": 2
              }
            }
          ]
        }
      ]
    }
  ]
}

Hand this to a coding agent. It rebuilds the figure above headlessly, and the session file is then the thing to edit: change the location or a track's settings, or add your own track under sessionTracks, and rerun.

cat > session.json <<'JSON'
{
  "sessionTracks": [
    {
      "type": "MultiQuantitativeTrack",
      "trackId": "HG002_imprinted_5mC_smoothed",
      "name": "HG002 5mC by haplotype, moving average",
      "assemblyNames": [
        "hg38"
      ],
      "adapter": {
        "type": "MultiWiggleAdapter",
        "subadapters": [
          {
            "type": "BedTabixAdapter",
            "name": "HP1",
            "color": "#d62728",
            "uri": "https://jbrowse.org/demos/methylation/HG002_chr15_imprinted_hp1.modkit.bed.gz"
          },
          {
            "type": "BedTabixAdapter",
            "name": "HP2",
            "color": "#1f77b4",
            "uri": "https://jbrowse.org/demos/methylation/HG002_chr15_imprinted_hp2.modkit.bed.gz"
          }
        ]
      },
      "displays": [
        {
          "type": "LinearMarkDisplay",
          "displayId": "HG002_imprinted_5mC_smoothed-LinearMarkDisplay",
          "rows": "source",
          "legend": {
            "show": false
          },
          "fetchSizeLimit": 5000000,
          "encoding": {
            "y": {
              "scale": {
                "domainMin": 0,
                "domainMax": 100
              },
              "axis": {
                "title": "% 5mC"
              }
            }
          },
          "marks": [
            {
              "mark": "point",
              "encoding": {
                "y": "score",
                "color": {
                  "value": "#bbbbbb"
                },
                "size": 2
              }
            },
            {
              "mark": {
                "type": "line",
                "interpolate": "linear"
              },
              "transform": [
                {
                  "type": "window",
                  "groupby": [
                    "source"
                  ],
                  "frame": [
                    -50,
                    50
                  ],
                  "ops": [
                    {
                      "op": "mean",
                      "field": "score",
                      "weight": "n_valid_cov",
                      "as": "meth"
                    }
                  ]
                }
              ],
              "encoding": {
                "y": "meth",
                "color": {
                  "field": "source",
                  "scale": {
                    "domain": [
                      "HP1",
                      "HP2"
                    ],
                    "range": [
                      "#d62728",
                      "#1f77b4"
                    ]
                  }
                },
                "strokeWidth": 2
              }
            }
          ]
        }
      ]
    }
  ],
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr15:22,200,000-26,900,000",
      "tracks": [
        {
          "trackId": "cpgisland_ucsc_hg38",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "ncbi_refseq_109_hg38_latest",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 90
        },
        {
          "trackId": "HG002_imprinted_5mC_smoothed",
          "height": 300,
          "forceLoad": true,
          "marks": [
            {
              "mark": "point",
              "encoding": {
                "y": "score",
                "color": {
                  "value": "#bbbbbb"
                },
                "size": 2
              }
            },
            {
              "mark": {
                "type": "line",
                "interpolate": "linear"
              },
              "transform": [
                {
                  "type": "window",
                  "groupby": [
                    "source"
                  ],
                  "frame": [
                    -500,
                    500
                  ],
                  "ops": [
                    {
                      "op": "mean",
                      "field": "score",
                      "weight": "n_valid_cov",
                      "as": "meth"
                    }
                  ]
                }
              ],
              "encoding": {
                "y": "meth",
                "color": {
                  "field": "source",
                  "scale": {
                    "domain": [
                      "HP1",
                      "HP2"
                    ],
                    "range": [
                      "#d62728",
                      "#1f77b4"
                    ]
                  }
                },
                "strokeWidth": 2
              }
            }
          ]
        }
      ]
    }
  ]
}
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/code/jb2/main/test_data/config_demo.json \
  --width 1500 --height 725 \
  --timeout 90000 \
  --spec session.json -o figure.png

Using JBrowse with AI agents covers the rest: hosted genomes to point it at, and why waiting for the render is the part that goes wrong.

See also

Notes

  1. To compare two samples, 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. ↩