Skip to main content

Methylation (bisulfite)

TL;DR: align WGBS/EM-seq short reads with bwameth, load the plain BAM, and JBrowse colors per-read methylation straight from the C→T conversion, with CpG, CHG, and CHH each selectable. No MM/ML tags and no methylation caller.

Prerequisites

Where the data comes from

TAIR10 (RefSeq GCF_000001735.4) and one wild-type Col-0 WGBS run from the European Nucleotide Archive, DRR029742 (paired-end 150 bp).

What bisulfite data looks like

Bisulfite sequencing (WGBS) and its enzymatic cousin EM-seq read DNA methylation without any long-read basecaller. A chemical (sodium bisulfite) or enzymatic (APOBEC) step converts every unmethylated cytosine to uracil, which reads as T, while a methylated cytosine is protected and still reads as C. Methylation is therefore recoverable from ordinary short Illumina reads by comparing each read to the reference: a C→T change at a cytosine means it was unmethylated, a retained C means it was methylated.

JBrowse 2 makes that comparison itself, per read, at render time.

Plants methylate in three sequence contexts: CpG, CHG, and CHH (H is A, C, or T). JBrowse restricts the coloring to any one of them, so all three read off the same pileup. Everything below runs on real Arabidopsis thaliana data, from SRA reads to a colored browser view.

Producing the BAM

The reproduce script runs the whole pipeline, from the TAIR10 reference and one wild-type Col-0 WGBS run (DRR029742, paired-end 150 bp) through Trim Galore to a sorted BAM. Two of its steps decide what JBrowse can read.

The aligner is bwameth, which handles bisulfite reads by in-silico C→T converting both reads and reference and then running bwa mem. It emits an ordinary BAM carrying the original read sequences, so the C→T signal survives for JBrowse to compare against the reference at render time. JBrowse reads Bismark's BAMs the same way.

Trimming and alignment are four commands, on any pair of WGBS or EM-seq FASTQs:

trim_galore --paired R1.fastq.gz R2.fastq.gz
# index once per reference: bwameth aligns against a C->T copy of it
bwameth.py index tair10.fa
# no methylation flags anywhere: the BAM keeps the original read sequences,
# and JBrowse makes the comparison at render time
bwameth.py --reference tair10.fa -t 8 R1_val_1.fq.gz R2_val_2.fq.gz \
  | samtools sort -o arabidopsis_wgbs.bam -
samtools index arabidopsis_wgbs.bam

bwameth.py index writes a C→T converted copy of the reference next to the original, so the reference directory has to be writable.

Check the conversion rate

An unconverted cytosine is indistinguishable from a methylated one, so the library's conversion rate is worth having before reading anything off the track. The chloroplast is unmethylated, which makes it the control:

MethylDackel extract --CHH -r NC_000932.1 -o conversion \
  tair10.fa arabidopsis_wgbs.bam

Column 4 of the resulting conversion_CHH.bedGraph is the methylated percentage at each cytosine, and on an unmethylated sequence that is the fraction the conversion missed. Modern libraries convert above 99%. An organism with no plastid uses whatever the library spiked in for this, usually unmethylated lambda or pUC19, added to the reference as an extra contig. The reproduce script prints the rate for this run.

Aggregate methylation, optionally

A per-position methylation fraction across the whole genome is complementary to the per-read coloring, and MethylDackel calls one in all three plant contexts. It writes a bedGraph per context, which becomes a bigWig once the header line is dropped and the percentage column kept:

# CpG is emitted always; --CHG --CHH add the two plant contexts.
# -o fixes the output prefix, which the loop below reads back.
MethylDackel extract --CHG --CHH -o arabidopsis_wgbs \
  tair10.fa arabidopsis_wgbs.bam
samtools faidx tair10.fa
cut -f1,2 tair10.fa.fai > tair10.chrom.sizes  # two columns, not the .fai itself
for ctx in CpG CHG CHH; do
  # tail drops MethylDackel's track line; cut keeps the percentage column
  tail -n +2 arabidopsis_wgbs_${ctx}.bedGraph | cut -f1-4 |
    sort -k1,1 -k2,2n > ${ctx}.bg
  bedGraphToBigWig ${ctx}.bg tair10.chrom.sizes arabidopsis_wgbs_${ctx}.bw
done

Group the resulting bigWigs into a single MultiQuantitativeTrack, a subadapter per context, and they render as three labeled rows, the Aggregate methylation track in the figures below. This is the same mechanism as the DNA methylation tutorial's aggregate section.

{
  "type": "MultiQuantitativeTrack",
  "trackId": "arabidopsis_wgbs_methyldackel",
  "name": "Aggregate methylation (MethylDackel)",
  "assemblyNames": ["tair10"],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "subadapters": [
      {
        "type": "BigWigAdapter",
        "name": "CpG",
        "color": "red",
        "uri": "arabidopsis_wgbs_CpG.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "CHG",
        "color": "green",
        "uri": "arabidopsis_wgbs_CHG.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "CHH",
        "color": "blue",
        "uri": "arabidopsis_wgbs_CHH.bw"
      }
    ]
  }
}

Load the tracks

The figures below use three tracks: the TAIR10 assembly, its gene models, and the arabidopsis_wgbs.bam produced above.

Set up the assembly from the same tair10.fa reference (the CLI indexes and bgzips it for you):

jbrowse add-assembly tair10.fa --name tair10 --load copy

The gene models come with the reference (datasets download genome accession GCF_000001735.4 --include gff3). Sort, compress, and index the GFF3, then add it as a FeatureTrack:

jbrowse sort-gff genomic.gff | bgzip > tair10.gff.gz
tabix -p gff tair10.gff.gz
jbrowse add-track tair10.gff.gz --name "TAIR10 genes" --load copy

The alignments track's displayDefaults decides which context it opens on, and the track menu switches it afterwards:

{
  "type": "AlignmentsTrack",
  "trackId": "arabidopsis_wgbs",
  "name": "Arabidopsis WGBS (bwameth)",
  "assemblyNames": ["tair10"],
  "adapter": {
    "type": "BamAdapter",
    "uri": "arabidopsis_wgbs.bam"
  },
  "displayDefaults": {
    "colorBy": {
      "type": "bisulfite",
      "modifications": { "cytosineContext": "CG" }
    }
  }
}
jbrowse add-track arabidopsis_wgbs.bam \
  --trackId arabidopsis_wgbs \
  --name "Arabidopsis WGBS (bwameth)" \
  --assemblyNames tair10 \
  --displayDefaults '{"colorBy":{"type":"bisulfite","modifications":{"cytosineContext":"CG"}}}' \
  --load copy

cytosineContext takes CG, CHG, CHH or all.

See the assemblies configuration guide for the equivalent assembly JSON.

The Aggregate methylation row in the figures is the optional MethylDackel track from the section above.

Every step here works identically on JBrowse Desktop, which opens tair10.fa, the BAM, and the bigWigs straight from local disk with no web server.

Color the reads

In the alignments track menu, pick Color by... → Bisulfite / EM-seq, then a cytosine context: CpG, CHG, CHH, or All cytosines. Methylated cytosines paint red. Once a context is set, that same submenu carries a Show unmethylated (blue) checkbox below the contexts, which paints the converted sites blue, separating an unmethylated cytosine from a position with no cytosine in it. The mode is reference-based, and applies to bisulfite and EM-seq libraries.

The figure and clip below leave Show unmethylated off.

Two methylation regimes

Plants run two unrelated methylation programs, and the three contexts are what tell them apart:

ProgramCpGCHGCHHEffect on the locus
Gene body methylation (gbM)yesnonoNone obvious, the gene stays transcribed. Maintained by MET1, depleted at both ends of the gene
Transposon and repeat silencingyesyesyesHeterochromatin, transcriptionally off. CMT3 maintains CHG, RdDM and CMT2 maintain CHH

So red in the CpG row alone is gene body methylation, and red in all three rows is silencing.

Type NC_003070.9:4,398,000-4,412,000 into the location box to reach a window on chromosome 1 that carries one of each: the expressed gene AT1G12930 on the left, and a transposon on the right. The gene body is methylated in CpG only; the transposon is methylated in all three contexts. The reproduce script prints the fraction per context for both regions.

The RepeatMasker lane names the element: META1_LTR#LTR/Copia, an LTR retrotransposon, which TAIR10's own transposable-element annotation calls AT1TE14315, family META1. The gene track carries a pseudogene, AT1G12935, over the same interval. The lane comes from UCSC's GenArk hub for TAIR10, whose sequence names are the RefSeq accessions this assembly already uses, so it loads with no aliasing.

TAIR10 genes, the RepeatMasker lane, the aggregate MethylDackel track, and three copies of the same WGBS pileup colored by CpG, CHG and CHH. AT1G12930 is red in CpG only; the LTR/Copia element on the right is red in all three. Open this view in JBrowse ↗

The steps behind the figure above, run against your own data.

  1. Open your genome: on the JBrowse Desktop start screen click Open new genome (or Show all available genomes to pick a hosted one).This figure uses arabidopsis, loaded from an indexed FASTA (.fa + .fai).
  2. Add your own track: File → Open track..., then paste a URL or choose a local file. This one needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “Genes (TAIR10 / RefSeq)”.
  3. Add your own track: File → Open track..., then paste a URL or choose a local file. This one needs a BigBed file (.bb).This figure uses “RepeatMasker (TAIR10)”.
  4. Track menu → Filter by... → add the expressionThis figure filters on `feature.end-feature.start>200`.“RepeatMasker (TAIR10)”.
  5. Drag the bar at the bottom edge of the track to resize it.“RepeatMasker (TAIR10)”.
  6. Add your own track: File → Open track..., then paste a URL or choose a local file.This figure uses “Aggregate methylation (MethylDackel)”.
  7. Track menu → Plot type → Multi-row → XY plot“Aggregate methylation (MethylDackel)”.
  8. Track menu → Score → Set min/max score...Sets the score-axis minimum (0 here).“Aggregate methylation (MethylDackel)”.
  9. Track menu → Score → Set min/max score...Sets the score-axis maximum (100 here).“Aggregate methylation (MethylDackel)”.
  10. Drag the bar at the bottom edge of the track to resize it.“Aggregate methylation (MethylDackel)”.
  11. Add your own track: File → Open track..., then paste a URL or choose a local file. This one needs a CRAM file (.cram + .crai).This figure uses “Per-read WGBS — CpG context”.
  12. Track menu → Color by... → Bisulfite / EM-seq → CpG“Per-read WGBS — CpG context”.
  13. Track menu → Show... → Show coverage (unchecked)“Per-read WGBS — CpG context”.
  14. Track menu → Read height → Track sizing → Fixed read height + fixed track height“Per-read WGBS — CpG context”.
  15. Track menu → Read height → Compact“Per-read WGBS — CpG context”.
  16. Drag the bar at the bottom edge of the track to resize it.“Per-read WGBS — CpG context”.
  17. Add your own track: File → Open track..., then paste a URL or choose a local file. This one needs a CRAM file (.cram + .crai).This figure uses “Per-read WGBS — CHG context”.
  18. Track menu → Color by... → Bisulfite / EM-seq → CHG“Per-read WGBS — CHG context”.
  19. Track menu → Show... → Show coverage (unchecked)“Per-read WGBS — CHG context”.
  20. Track menu → Read height → Track sizing → Fixed read height + fixed track height“Per-read WGBS — CHG context”.
  21. Track menu → Read height → Compact“Per-read WGBS — CHG context”.
  22. Drag the bar at the bottom edge of the track to resize it.“Per-read WGBS — CHG context”.
  23. Add your own track: File → Open track..., then paste a URL or choose a local file. This one needs a CRAM file (.cram + .crai).This figure uses “Per-read WGBS — CHH context”.
  24. Track menu → Color by... → Bisulfite / EM-seq → CHH“Per-read WGBS — CHH context”.
  25. Track menu → Show... → Show coverage (unchecked)“Per-read WGBS — CHH context”.
  26. Track menu → Read height → Track sizing → Fixed read height + fixed track height“Per-read WGBS — CHH context”.
  27. Track menu → Read height → Compact“Per-read WGBS — CHH context”.
  28. Drag the bar at the bottom edge of the track to resize it.“Per-read WGBS — CHH context”.
  29. Type your region of interest into the location box and press Enter.This figure is at NC_003070.9:4,398,000-4,412,000.

Open this view in JBrowse Desktop ↗

Opens JBrowse Desktop (JBrowse Desktop 5.0+) at this view and saves it as a reopenable session — swap in your own files afterwards.

Nothing happens? Paste this link into Desktop's Open JBrowse Web link... (start screen, or File → Session):

https://jbrowse.org/code/jb2/main/?config=test_data%2Farabidopsis_methylation%2Fconfig_emseq_bisulfite.json&session=spec-%7B%22sessionTracks%22%3A%5B%7B%22type%22%3A%22AlignmentsTrack%22%2C%22trackId%22%3A%22arabidopsis_wgbs_cg%22%2C%22name%22%3A%22Per-read+WGBS+%E2%80%94+CpG+context%22%2C%22assemblyNames%22%3A%5B%22arabidopsis%22%5D%2C%22adapter%22%3A%7B%22type%22%3A%22CramAdapter%22%2C%22uri%22%3A%22https%3A%2F%2Fjbrowse.org%2Fdemos%2Fbisulfite%2Farabidopsis_wgbs_bisulfite.cram%22%7D%7D%2C%7B%22type%22%3A%22AlignmentsTrack%22%2C%22trackId%22%3A%22arabidopsis_wgbs_chg%22%2C%22name%22%3A%22Per-read+WGBS+%E2%80%94+CHG+context%22%2C%22assemblyNames%22%3A%5B%22arabidopsis%22%5D%2C%22adapter%22%3A%7B%22type%22%3A%22CramAdapter%22%2C%22uri%22%3A%22https%3A%2F%2Fjbrowse.org%2Fdemos%2Fbisulfite%2Farabidopsis_wgbs_bisulfite.cram%22%7D%7D%2C%7B%22type%22%3A%22AlignmentsTrack%22%2C%22trackId%22%3A%22arabidopsis_wgbs_chh%22%2C%22name%22%3A%22Per-read+WGBS+%E2%80%94+CHH+context%22%2C%22assemblyNames%22%3A%5B%22arabidopsis%22%5D%2C%22adapter%22%3A%7B%22type%22%3A%22CramAdapter%22%2C%22uri%22%3A%22https%3A%2F%2Fjbrowse.org%2Fdemos%2Fbisulfite%2Farabidopsis_wgbs_bisulfite.cram%22%7D%7D%5D%2C%22views%22%3A%5B%7B%22type%22%3A%22LinearGenomeView%22%2C%22assembly%22%3A%22arabidopsis%22%2C%22loc%22%3A%22NC_003070.9%3A4%2C398%2C000-4%2C412%2C000%22%2C%22tracks%22%3A%5B%7B%22trackId%22%3A%22arabidopsis_genes%22%7D%2C%7B%22trackId%22%3A%22arabidopsis_rmsk%22%2C%22type%22%3A%22LinearBasicDisplay%22%2C%22jexlFiltersSetting%22%3A%5B%22jexl%3Afeature.end-feature.start%3E200%22%5D%2C%22height%22%3A50%7D%2C%7B%22trackId%22%3A%22arabidopsis_methyldackel%22%2C%22type%22%3A%22MultiLinearWiggleDisplay%22%2C%22defaultRendering%22%3A%22multirowxy%22%2C%22minScore%22%3A0%2C%22maxScore%22%3A100%2C%22height%22%3A170%7D%2C%7B%22trackId%22%3A%22arabidopsis_wgbs_cg%22%2C%22type%22%3A%22LinearAlignmentsDisplay%22%2C%22colorBy%22%3A%7B%22type%22%3A%22bisulfite%22%2C%22modifications%22%3A%7B%22cytosineContext%22%3A%22CG%22%7D%7D%2C%22showCoverage%22%3Afalse%2C%22heightMode%22%3A%22fixed%22%2C%22featureHeight%22%3A3%2C%22height%22%3A90%7D%2C%7B%22trackId%22%3A%22arabidopsis_wgbs_chg%22%2C%22type%22%3A%22LinearAlignmentsDisplay%22%2C%22colorBy%22%3A%7B%22type%22%3A%22bisulfite%22%2C%22modifications%22%3A%7B%22cytosineContext%22%3A%22CHG%22%7D%7D%2C%22showCoverage%22%3Afalse%2C%22heightMode%22%3A%22fixed%22%2C%22featureHeight%22%3A3%2C%22height%22%3A90%7D%2C%7B%22trackId%22%3A%22arabidopsis_wgbs_chh%22%2C%22type%22%3A%22LinearAlignmentsDisplay%22%2C%22colorBy%22%3A%7B%22type%22%3A%22bisulfite%22%2C%22modifications%22%3A%7B%22cytosineContext%22%3A%22CHH%22%7D%7D%2C%22showCoverage%22%3Afalse%2C%22heightMode%22%3A%22fixed%22%2C%22featureHeight%22%3A3%2C%22height%22%3A90%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+methylation%2Farabidopsis_wgbs_contexts

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 https://jbrowse.org/demos/bisulfite/arabidopsis_wgbs_bisulfite.cram \
  --trackId arabidopsis_wgbs_cg \
  --name "Per-read WGBS — CpG context" \
  --assemblyNames arabidopsis

jbrowse add-track https://jbrowse.org/demos/bisulfite/arabidopsis_wgbs_bisulfite.cram \
  --trackId arabidopsis_wgbs_chg \
  --name "Per-read WGBS — CHG context" \
  --assemblyNames arabidopsis

jbrowse add-track https://jbrowse.org/demos/bisulfite/arabidopsis_wgbs_bisulfite.cram \
  --trackId arabidopsis_wgbs_chh \
  --name "Per-read WGBS — CHH context" \
  --assemblyNames arabidopsis

Run these where the config.json is, or add --out <dir>, and point each uri at your own file. The config needs arabidopsis 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 Session spec tab, or jbrowse set-default-session.

This session spec pastes after &session=spec- on any JBrowse Web instance.

{
  "sessionTracks": [
    {
      "type": "AlignmentsTrack",
      "trackId": "arabidopsis_wgbs_cg",
      "name": "Per-read WGBS — CpG context",
      "assemblyNames": [
        "arabidopsis"
      ],
      "adapter": {
        "type": "CramAdapter",
        "uri": "https://jbrowse.org/demos/bisulfite/arabidopsis_wgbs_bisulfite.cram"
      }
    },
    {
      "type": "AlignmentsTrack",
      "trackId": "arabidopsis_wgbs_chg",
      "name": "Per-read WGBS — CHG context",
      "assemblyNames": [
        "arabidopsis"
      ],
      "adapter": {
        "type": "CramAdapter",
        "uri": "https://jbrowse.org/demos/bisulfite/arabidopsis_wgbs_bisulfite.cram"
      }
    },
    {
      "type": "AlignmentsTrack",
      "trackId": "arabidopsis_wgbs_chh",
      "name": "Per-read WGBS — CHH context",
      "assemblyNames": [
        "arabidopsis"
      ],
      "adapter": {
        "type": "CramAdapter",
        "uri": "https://jbrowse.org/demos/bisulfite/arabidopsis_wgbs_bisulfite.cram"
      }
    }
  ],
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "arabidopsis",
      "loc": "NC_003070.9:4,398,000-4,412,000",
      "tracks": [
        {
          "trackId": "arabidopsis_genes"
        },
        {
          "trackId": "arabidopsis_rmsk",
          "type": "LinearBasicDisplay",
          "jexlFiltersSetting": [
            "jexl:feature.end-feature.start>200"
          ],
          "height": 50
        },
        {
          "trackId": "arabidopsis_methyldackel",
          "type": "MultiLinearWiggleDisplay",
          "defaultRendering": "multirowxy",
          "minScore": 0,
          "maxScore": 100,
          "height": 170
        },
        {
          "trackId": "arabidopsis_wgbs_cg",
          "type": "LinearAlignmentsDisplay",
          "colorBy": {
            "type": "bisulfite",
            "modifications": {
              "cytosineContext": "CG"
            }
          },
          "showCoverage": false,
          "heightMode": "fixed",
          "featureHeight": 3,
          "height": 90
        },
        {
          "trackId": "arabidopsis_wgbs_chg",
          "type": "LinearAlignmentsDisplay",
          "colorBy": {
            "type": "bisulfite",
            "modifications": {
              "cytosineContext": "CHG"
            }
          },
          "showCoverage": false,
          "heightMode": "fixed",
          "featureHeight": 3,
          "height": 90
        },
        {
          "trackId": "arabidopsis_wgbs_chh",
          "type": "LinearAlignmentsDisplay",
          "colorBy": {
            "type": "bisulfite",
            "modifications": {
              "cytosineContext": "CHH"
            }
          },
          "showCoverage": false,
          "heightMode": "fixed",
          "featureHeight": 3,
          "height": 90
        }
      ]
    }
  ]
}

Loaded against config: test_data/arabidopsis_methylation/config_emseq_bisulfite.json

The same view with jbrowse-anywidget. Point the adapter at your own file.

from jbrowse_anywidget import LinearGenomeView, make_assembly

view = LinearGenomeView(
    assembly=make_assembly("arabidopsis", "https://your-server/your-genome.fa.gz"),
    location="NC_003070.9:4,398,000-4,412,000",
)

view.add_track({
    "type": "FeatureTrack",
    "trackId": "my_track_1",
    "name": "My track 1",
    "assemblyNames": ["arabidopsis"],
    "adapter": {"type": "Gff3TabixAdapter", "uri": "https://your-server/your-file"},
})
view.add_track({
    "type": "FeatureTrack",
    "trackId": "my_track_2",
    "name": "My track 2",
    "assemblyNames": ["arabidopsis"],
    "adapter": {"type": "BigBedAdapter", "uri": "https://your-server/your-file"},
})
view.add_track({
    "type": "MultiQuantitativeTrack",
    "trackId": "my_track_3",
    "name": "My track 3",
    "assemblyNames": ["arabidopsis"],
    "adapter": {"type": "MultiWiggleAdapter", "uri": "https://your-server/your-file"},
})
view.add_track({
    "type": "AlignmentsTrack",
    "trackId": "my_track_4",
    "name": "My track 4",
    "assemblyNames": ["arabidopsis"],
    "adapter": {"type": "CramAdapter", "uri": "https://your-server/your-file"},
})
view.add_track({
    "type": "AlignmentsTrack",
    "trackId": "my_track_5",
    "name": "My track 5",
    "assemblyNames": ["arabidopsis"],
    "adapter": {"type": "CramAdapter", "uri": "https://your-server/your-file"},
})
view.add_track({
    "type": "AlignmentsTrack",
    "trackId": "my_track_6",
    "name": "My track 6",
    "assemblyNames": ["arabidopsis"],
    "adapter": {"type": "CramAdapter", "uri": "https://your-server/your-file"},
})

view  # display the widget

Hand this to a coding agent. It rebuilds the figure above headlessly, and the session file is then the thing to edit — swap an adapter uri for your own file and rerun.

cat > session.json <<'JSON'
{
  "sessionTracks": [
    {
      "type": "AlignmentsTrack",
      "trackId": "arabidopsis_wgbs_cg",
      "name": "Per-read WGBS — CpG context",
      "assemblyNames": [
        "arabidopsis"
      ],
      "adapter": {
        "type": "CramAdapter",
        "uri": "https://jbrowse.org/demos/bisulfite/arabidopsis_wgbs_bisulfite.cram"
      }
    },
    {
      "type": "AlignmentsTrack",
      "trackId": "arabidopsis_wgbs_chg",
      "name": "Per-read WGBS — CHG context",
      "assemblyNames": [
        "arabidopsis"
      ],
      "adapter": {
        "type": "CramAdapter",
        "uri": "https://jbrowse.org/demos/bisulfite/arabidopsis_wgbs_bisulfite.cram"
      }
    },
    {
      "type": "AlignmentsTrack",
      "trackId": "arabidopsis_wgbs_chh",
      "name": "Per-read WGBS — CHH context",
      "assemblyNames": [
        "arabidopsis"
      ],
      "adapter": {
        "type": "CramAdapter",
        "uri": "https://jbrowse.org/demos/bisulfite/arabidopsis_wgbs_bisulfite.cram"
      }
    }
  ],
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "arabidopsis",
      "loc": "NC_003070.9:4,398,000-4,412,000",
      "tracks": [
        {
          "trackId": "arabidopsis_genes"
        },
        {
          "trackId": "arabidopsis_rmsk",
          "type": "LinearBasicDisplay",
          "jexlFiltersSetting": [
            "jexl:feature.end-feature.start>200"
          ],
          "height": 50
        },
        {
          "trackId": "arabidopsis_methyldackel",
          "type": "MultiLinearWiggleDisplay",
          "defaultRendering": "multirowxy",
          "minScore": 0,
          "maxScore": 100,
          "height": 170
        },
        {
          "trackId": "arabidopsis_wgbs_cg",
          "type": "LinearAlignmentsDisplay",
          "colorBy": {
            "type": "bisulfite",
            "modifications": {
              "cytosineContext": "CG"
            }
          },
          "showCoverage": false,
          "heightMode": "fixed",
          "featureHeight": 3,
          "height": 90
        },
        {
          "trackId": "arabidopsis_wgbs_chg",
          "type": "LinearAlignmentsDisplay",
          "colorBy": {
            "type": "bisulfite",
            "modifications": {
              "cytosineContext": "CHG"
            }
          },
          "showCoverage": false,
          "heightMode": "fixed",
          "featureHeight": 3,
          "height": 90
        },
        {
          "trackId": "arabidopsis_wgbs_chh",
          "type": "LinearAlignmentsDisplay",
          "colorBy": {
            "type": "bisulfite",
            "modifications": {
              "cytosineContext": "CHH"
            }
          },
          "showCoverage": false,
          "heightMode": "fixed",
          "featureHeight": 3,
          "height": 90
        }
      ]
    }
  ]
}
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/code/jb2/main/test_data/arabidopsis_methylation/config_emseq_bisulfite.json \
  --session 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.

One WGBS pileup recolored CpG, then CHG, then CHH from the track menu, under the TAIR10 genes, the RepeatMasker lane and the aggregate MethylDackel rows: red holds over the LTR element through all three contexts and drains from the gene body. Open the session this video starts in ↗

Reproduce it end to end

The whole pipeline is wrapped in one script, build_arabidopsis_wgbs.sh:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_arabidopsis_wgbs.sh
bash build_arabidopsis_wgbs.sh          # builds ./arabidopsis_wgbs_build/jbrowse2
npx --yes serve arabidopsis_wgbs_build/jbrowse2 # then open the printed URL

It downloads the TAIR10 reference and the DRR029742 WGBS run, trims and bisulfite-aligns them with bwameth, downloads JBrowse, and writes a config.json with the assembly, the gene models, and the per-read pileup pre-colored Bisulfite / CpG, opening on the window above.

With MethylDackel on PATH it also prints the conversion rate and the per-context fraction over both regions, and skips both with a warning if it is not. The aggregate bigWig track is left out either way.

On Debian/Ubuntu, apt install wget samtools tabix covers several of those. bwameth, Trim Galore, and the NCBI datasets CLI install from their own instructions, and node comes from nodejs.org. The alignment step downloads a full WGBS run, so allow time and disk for it.

See also