Pangenome (hosting your own graph)

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.

To let people browse a pangenome graph in JBrowse, from a whole chromosome down to single nodes, the graph has to be cut into small indexed files that JBrowse reads one window at a time. One command writes those files and the config that puts them on tracks. We run it on HPRC release 2, and:

load the plugin and run the command
open the track as a graph and as a row of segments
check the index against the graph
add two optional layers: haplotypes per segment, and haplotype routes
convert the haplotype routes to PAF for a whole-genome synteny view

Experimental

The graph view is a beta plugin. We welcome your feedback.

Prerequisites

the GraphGenomeView plugin
htslib (bgzip, tabix) and bcftools
gfa-to-tabix 1.4.0 or later, for the graph's files and the walk index (walk format 5, the only one the plugin reads)
gfatools, for an rGFA's bubbles
minigraph, for each assembly's path through the graph
vg, for a plain GFA's bubbles
for the walks as PAF: python3 and the JBrowse CLI for make-pif

# Debian/Ubuntu
sudo apt install tabix bcftools
# or conda, which has every tool above
conda install -c bioconda -c conda-forge htslib bcftools gfatools minigraph \
  vg
# any platform
cargo install gfa-to-tabix
npm install -g @jbrowse/cli

Download the data

The input is HPRC release 2's Minigraph-Cactus graph. gfa-to-tabix build reads the rGFA, and the walk index reads the base-level GFA, both as local files. Fetch the rGFA first and the base-level GFA when you reach the walks:

B=https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/release2/minigraph-cactus/v2.1/hprc-v2.1-mc-grch38/hprc-v2.1-mc-grch38
# SV-resolution rGFA, 0.8 GB
wget -c $B.sv.gfa.gz
# base-level GFA with every haplotype's walk, 64 GB
wget -c $B.gfa.gz

Other files (no download needed)

What it is	File
The finished files, for comparison	README.txt

The GraphGenomeView plugin

Install GraphGenomeView from the plugin store, at Tools → Plugin store in JBrowse Web or Global plugins... on the JBrowse Desktop start screen (Installing plugins). It needs JBrowse 5. The config the command writes names it with the plugins entry below, in which storePlugin picks the build for the JBrowse reading the config, and esmUrl is the fallback when the store cannot be read (configuring plugins).

{
  "plugins": [
    {
      "name": "GraphGenomeView",
      "storePlugin": "GraphGenomeView",
      "esmUrl": "https://jbrowse.org/plugins/jbrowse-plugin-graphgenomeviewer/latest/dist/jbrowse-plugin-graphgenomeviewer.esm.js"
    }
  ]
}

Indexing a graph with gfa-to-tabix build

Install gfa-to-tabix 0.12.0 or later, whose build command writes every file below (its releases have Linux and macOS binaries if you have no Rust toolchain):

cargo install gfa-to-tabix

An rGFA, from minigraph or the minigraph stage of Minigraph-Cactus, needs an output prefix and the assembly name your config uses. Its bubbles come from gfatools bubble, so put gfatools on PATH:

gfa-to-tabix build graph.rgfa.gz -o out --assembly hg38

A plain GFA, from pggb, odgi, vg or base-level Minigraph-Cactus, also needs the backbone sample and its bubbles. The bubbles come from a snarl VCF (vg's word for a bubble is snarl):

vg deconstruct -a writes one record per snarl against the reference path, tagged LV for its level in the snarl tree, and pggb -V writes the same file. gfa-to-tabix bubbles keeps the LV=0 records, so the tag has to be there.
The VCF's CHROM must be the assembly's refName, so rename the PanSN path (sample#haplotype#contig) as the pggb tutorial does:

# -p: the reference path to decompose against
vg deconstruct -a -p K12#1#chr graph.gfa > graph.snarls.vcf
printf 'K12#1#chr\tchr\n' > rename_chrs.tsv
bcftools annotate --rename-chrs rename_chrs.tsv graph.snarls.vcf \
  | bcftools sort -Oz -o graph.snarls.vcf.gz

gfa-to-tabix build graph.gfa -o out --reference K12 --assembly K12 --snarls graph.snarls.vcf.gz

At human-chromosome scale, index the SV-resolution rGFA; a pggb graph's index does not finish there.

The command writes these files beside the prefix:

file	what it holds
`.segs.bed.gz`	one row per segment, at its reference coordinate
`.links.bed.gz`	one row per link per endpoint, both ends stated in full
`.bubbles.bed.gz`	where haplotypes diverge and rejoin, with each bubble's shortest and longest allele
`.fold10000.segs.bed.gz`	the graph with variants under 10 kb folded into the reference, so a whole chromosome draws
`.fold10000.links.bed.gz`	the tier's links
`.contig.segs.bed.gz`, `.contig.links.bed.gz`	the segment and link rows under each segment's own coordinate, which the allele inventory reads
`.alleles.bed.gz`	one row per allele, with a CIGAR that states its size
`.graph.json`	a manifest naming the files above and the tier's fold size
`.config.json`	the tracks below, with the plugin entry

The coarse tier for whole-chromosome views

The coarse tier is a copy of the graph with every variant under a threshold folded into the reference. It keeps each allele whose length, or the reference it replaces, reaches the threshold, and the way back to the reference from both of its ends. The graph track folds what it draws the same way at ten of the linear view's pixels, so zooming out onto the tier keeps every loop the fine cut drew. The threshold is in the file name, 10,000 bp by default for an rGFA; --tier sets it, and a plain GFA defaults to 50, since most of its variants are single bases.

Configuring the graph track

The config's first track is the graph. Its adapter holds:

uri, the prefix
coarse, which names the tier and its threshold, foldBelowBp; the track draws the tier past a tenth of that in bp per pixel
assemblyNameToPanSN, which maps your assembly name to the graph's PanSN sample; the command writes it only when the two differ

Config file

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

{
  "type": "GraphTrack",
  "trackId": "hprc_graph",
  "name": "hprc graph",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "RgfaTabixAdapter",
    "uri": "hprc",
    "assemblyNameToPanSN": { "hg38": "GRCh38" },
    "coarse": { "uri": "hprc.fold10000", "foldBelowBp": 10000 }
  },
  "displayDefaults": { "showLabels": "none" },
  "displays": [
    { "type": "LinearGraphDisplay" },
    { "type": "LinearBasicDisplay" }
  ]
}

CLI (add-track-json)

jbrowse add-track-json '{
  "type": "GraphTrack",
  "trackId": "hprc_graph",
  "name": "hprc graph",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "RgfaTabixAdapter",
    "uri": "hprc",
    "assemblyNameToPanSN": { "hg38": "GRCh38" },
    "coarse": { "uri": "hprc.fold10000", "foldBelowBp": 10000 }
  },
  "displayDefaults": { "showLabels": "none" },
  "displays": [
    { "type": "LinearGraphDisplay" },
    { "type": "LinearBasicDisplay" }
  ]
}'

In the app

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": "GraphTrack",
  "trackId": "hprc_graph",
  "name": "hprc graph",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "RgfaTabixAdapter",
    "uri": "hprc",
    "assemblyNameToPanSN": { "hg38": "GRCh38" },
    "coarse": { "uri": "hprc.fold10000", "foldBelowBp": 10000 }
  },
  "displayDefaults": { "showLabels": "none" },
  "displays": [
    { "type": "LinearGraphDisplay" },
    { "type": "LinearBasicDisplay" }
  ]
}

hprc, hprc.fold10000 are relative to a config.json. Replace each with its URL or its path on this computer.

Turned on, the track draws as a graph. Display types → Feature display in its track menu draws the same segments as a row.

The HPRC graph's segment index over C4 on hg38.

The command writes out.config.json. Merge its plugins and tracks entries into your own config, then run jbrowse validate config.json, which reports a misspelled or undeclared slot JBrowse otherwise ignores. Without a config, the plugin's Add pangenome graph track workflow reads out.graph.json beside out.segs.bed.gz, so a track opened from that url gets the coarse tier too.

The other three tracks in the config draw the bubbles as a row and as a curve, and the allele inventory, one row per alternative path, as an alignments track. Turn them on and switch the graph back with Display types → Graph.

The four tracks the command writes, over the C4 region on hg38: the bubbles as a row and as a curve, the allele inventory, and the graph track at the bottom.

A node's right-click menu offers Open in the haplotype its rGFA id names (NA20809#2#CM094351.1 opens NA20809#2), when the session holds an assembly named or aliased sample#haplotype. The HPRC tutorial takes that route.

Checking the index against the graph

Query a locus out of the index, by the graph's PanSN name for the reference contig:

tabix hprc.segs.bed.gz 'GRCh38#0#chr1:103,690,000-103,700,000' | head -3

Columns one to three are the contig and span, four is the segment id, and five is its rank (build order). Ask the graph about one of those segments:

gfatools view -l s12829 -r 0 hprc-v2.1-mc-grch38.sv.gfa.gz

The S-line's SN and SO tags match the row's first two columns, and SR matches the fifth.

An empty result over a tiled reference means the query used the wrong name: chr1 finds nothing where GRCh38#0#chr1 finds every segment.
Backbone rows with no alleles mark a place where the graph collapsed, which minigraph does to near-identical segmental duplications.

Recording which haplotypes pass through each segment

The rank in an rGFA is build order, so it names the first assembly a segment came from. The full list of haplotypes whose paths pass through each segment comes from one of two sources.

With the assemblies, map each one back through the graph, reference first, and read the path it takes:

# --call: the path this sample takes through every bubble, one line per
#   `gfatools bubble` line, in the same order for every sample
# -xasm: the assembly-to-graph preset
# -c: base-level alignment, which --call reads
minigraph -cxasm --call -t 8 graph.rgfa.gz sample.fa > sample.call.bed

gfa-to-tabix paths joins the call files, the reference's first, into one tabix-indexed row per bubble per sample, drawn as one row per haplotype:

gfa-to-tabix paths -o graph.paths ref.call.bed sample.call.bed

With a plain GFA, the command records the haplotypes whose paths visit each segment as an SM:Z: tag while it reads the paths. The track below is written for a build like the pggb graph of five E. coli strains, indexed against K12.

The node panel lists them as samples; a track reads them as feature.samples and their count as feature.sampleCount.
Color by... → Attribute... with sampleCount gives each count a separate color.
For more than a handful of haplotypes, Edit plot... in the same dialog takes a ramp, here red for a segment on one haplotype to grey for a segment on most:

Config file

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

{
  "type": "FeatureTrack",
  "trackId": "graph_haplotypes_per_segment",
  "name": "graph: haplotypes per segment",
  "assemblyNames": ["K12"],
  "adapter": {
    "type": "RgfaTabixAdapter",
    "uri": "graph"
  },
  "displayDefaults": {
    "encoding": {
      "color": {
        "field": "sampleCount",
        "scale": {
          "type": "linear",
          "domainMin": 1,
          "range": ["#e31a1c", "#bdbdbd"]
        },
        "legend": { "title": "Haplotypes" }
      }
    }
  }
}

CLI (add-track-json)

jbrowse add-track-json '{
  "type": "FeatureTrack",
  "trackId": "graph_haplotypes_per_segment",
  "name": "graph: haplotypes per segment",
  "assemblyNames": ["K12"],
  "adapter": {
    "type": "RgfaTabixAdapter",
    "uri": "graph"
  },
  "displayDefaults": {
    "encoding": {
      "color": {
        "field": "sampleCount",
        "scale": {
          "type": "linear",
          "domainMin": 1,
          "range": ["#e31a1c", "#bdbdbd"]
        },
        "legend": { "title": "Haplotypes" }
      }
    }
  }
}'

In the app

In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add pangenome graph track, and fill in:

File type: rGFA segments (tabix BED pair)
Path: graph.segs.bed.gz
Track name: graph: haplotypes per segment
Assembly: K12

graph is relative to a config.json. Replace it with its URL or its path on this computer.

The encoding.color block is a color scale that maps the sampleCount field linearly from scale.domainMin onto the two colors of scale.range. The count is per haplotype (HG002.1), so a diploid sample's two copies count separately.

Indexing haplotype walks with gfa-to-tabix

A walk is one haplotype's route through the graph. gfa-to-tabix --walks (1.1.0 or later) writes every walk to one tabix-indexed file per reference, which the browser reads for the graph track and for the haplotype lanes. It cuts each walk into pieces, one for each stretch the walk spends in a 64 kb chunk of a reference, and files each piece under its chunk, so the browser reads one range of the file per window. A chunk holds the reference's path through it once, each distinct route a haplotype takes off that path once, and every other walk's piece as the list of routes it takes.

The base-level GFA holds each walk as a W line. Index it against both of HPRC's references:¹

# --refs: one walk file per reference sample, in its coordinates
gfa-to-tabix hprc-v2.1-mc-grch38.gfa.gz --walks --refs GRCh38,CHM13 \
  -o hprc-v2.1-mc-grch38

vg convert -f graph.gbz writes the same GFA for a graph published only as a GBZ.

The tool writes one file and its .tbi index for each sample in --refs, named <prefix>.<sample>.walks.bed.gz, so a track on GRCh38 reads an index that covers GRCh38 alone:

file	what it holds
`.GRCh38.walks.bed.gz`	every walk through each 64 kb chunk of GRCh38
`.CHM13.walks.bed.gz`	every walk through each 64 kb chunk of T2T-CHM13

The file opens with header lines, which tabix -H prints:

tabix -H hprc-v2.1-mc-grch38.CHM13.walks.bed.gz | head -4

#walks	format:i:5	chunk:i:65536	maxnode:i:1024	stay:i:311	abridge:i:1000000	abridgebytes:i:3000000
#reference	CHM13
#haplotype	GRCh38#0
#haplotype	HG00097#1

The first line gives the file's format, the chunk size, the longest node in bp, the --stay the tool ran with, and the --abridge-steps and --abridge-bytes budgets. Where a chunk's walks cost more than those budgets, usually in a satellite array, gfa-to-tabix moves the costliest walks to detail rows that a window read skips. The file then lists its reference sample and every other haplotype with rows in it, the other reference included, once each.

A track's walksUri names one reference's set by its prefix, <prefix>.<sample>. The track below draws the graph and, as a second display, one lane per haplotype. Each haplotype in assemblyNames is an assembly whose aliases include its sample#haplotype name, as the HPRC tutorial declares one; assemblyNameToPanSN covers the reference, which has none:

Config file

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

{
  "type": "GraphTrack",
  "trackId": "my_graph_lanes",
  "name": "My graph: haplotypes vs the reference",
  "assemblyNames": ["hg38", "HG00097.1", "HG00099.1"],
  "adapter": {
    "type": "WalkTabixSyntenyAdapter",
    "walksUri": "https://example.com/graphs/my_graph.GRCh38",
    "assemblyNames": ["hg38"],
    "assemblyNameToPanSN": { "hg38": "GRCh38#0" }
  },
  "displays": [
    { "type": "MultiWaySyntenyDisplay", "height": 600 },
    { "type": "LinearGraphDisplay" }
  ]
}

CLI (add-track-json)

jbrowse add-track-json '{
  "type": "GraphTrack",
  "trackId": "my_graph_lanes",
  "name": "My graph: haplotypes vs the reference",
  "assemblyNames": ["hg38", "HG00097.1", "HG00099.1"],
  "adapter": {
    "type": "WalkTabixSyntenyAdapter",
    "walksUri": "https://example.com/graphs/my_graph.GRCh38",
    "assemblyNames": ["hg38"],
    "assemblyNameToPanSN": { "hg38": "GRCh38#0" }
  },
  "displays": [
    { "type": "MultiWaySyntenyDisplay", "height": 600 },
    { "type": "LinearGraphDisplay" }
  ]
}'

In the app

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": "GraphTrack",
  "trackId": "my_graph_lanes",
  "name": "My graph: haplotypes vs the reference",
  "assemblyNames": ["hg38", "HG00097.1", "HG00099.1"],
  "adapter": {
    "type": "WalkTabixSyntenyAdapter",
    "walksUri": "https://example.com/graphs/my_graph.GRCh38",
    "assemblyNames": ["hg38"],
    "assemblyNameToPanSN": { "hg38": "GRCh38#0" }
  },
  "displays": [
    { "type": "MultiWaySyntenyDisplay", "height": 600 },
    { "type": "LinearGraphDisplay" }
  ]
}

The lanes align each walk to the reference's on the nodes both visit, so the track needs no sequence beyond the walk file.

A graph track with no lanes names the same set through RgfaTabixAdapter, here T2T-CHM13's:

Config file

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

{
  "type": "GraphTrack",
  "trackId": "hprc_walks_chm13",
  "name": "HPRC v2.1 haplotype walks",
  "assemblyNames": ["hs1"],
  "adapter": {
    "type": "RgfaTabixAdapter",
    "walksUri": "hprc-v2.1-mc-grch38.CHM13",
    "assemblyNameToPanSN": { "hs1": "CHM13" },
    "defaultHaplotypes": [
      "HG00097#1",
      "HG00099#1",
      "HG00128#1",
      "HG00133#1",
      "HG01109#1",
      "HG01123#1",
      "HG01960#1",
      "HG02055#1"
    ]
  }
}

CLI (add-track-json)

jbrowse add-track-json '{
  "type": "GraphTrack",
  "trackId": "hprc_walks_chm13",
  "name": "HPRC v2.1 haplotype walks",
  "assemblyNames": ["hs1"],
  "adapter": {
    "type": "RgfaTabixAdapter",
    "walksUri": "hprc-v2.1-mc-grch38.CHM13",
    "assemblyNameToPanSN": { "hs1": "CHM13" },
    "defaultHaplotypes": [
      "HG00097#1",
      "HG00099#1",
      "HG00128#1",
      "HG00133#1",
      "HG01109#1",
      "HG01123#1",
      "HG01960#1",
      "HG02055#1"
    ]
  }
}'

In the app

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": "GraphTrack",
  "trackId": "hprc_walks_chm13",
  "name": "HPRC v2.1 haplotype walks",
  "assemblyNames": ["hs1"],
  "adapter": {
    "type": "RgfaTabixAdapter",
    "walksUri": "hprc-v2.1-mc-grch38.CHM13",
    "assemblyNameToPanSN": { "hs1": "CHM13" },
    "defaultHaplotypes": [
      "HG00097#1",
      "HG00099#1",
      "HG00128#1",
      "HG00133#1",
      "HG01109#1",
      "HG01123#1",
      "HG01960#1",
      "HG02055#1"
    ]
  }
}

assemblyNameToPanSN maps the assembly name to the set's reference sample, here hs1 to CHM13
defaultHaplotypes lists the haplotypes the track draws until the reader picks others, as sample#haplotype, or a bare sample for both of its haplotypes

The HPRC tracks built this way open at chr22:20,000,000-20,100,000 on hs1, from the CHM13 set and hg38, from the GRCh38 set. The HPRC page reads the GRCh38 set for its Haplotypes launch, for the graph it draws under a clicked row, and for its BandageJS link.

The track menu's Haplotypes submenu switches the track between The track's default haplotypes, Every haplotype in the graph and Chosen in Settings..., whose dialog offers the names from the walk file's header. A graph view opened from the track draws the same set. Where the graph collapses the copies of a repeat onto one set of nodes, the track draws a walk through them where the reference places those nodes.

Two budgets in the adapter's config refuse a window with a notice to zoom in. walkByteBudget caps the compressed bytes that the three indexes estimate a window would fetch. The estimate is the same for any set of haplotypes, because the browser downloads every haplotype's rows and drops the unwanted ones by name. walkStepBudget caps the walk steps a window keeps, counting only the haplotypes the track draws, so fewer haplotypes draw a wider window.

Whole-genome synteny from the walks, as PAF

The walk files serve one window at a time. For a whole-genome synteny view, a GFA's walks also convert to PAF against the reference walk, one record per stretch of nodes the two walks share. Convert the walks you want:

curl -fO https://raw.githubusercontent.com/cmdcolin/gfa-to-pairwise-paf/v1.0.0/gfa_to_pairwise_paf.py
# --contig-lengths: each assembly's .fai, since walks omit contig lengths
gzip -dc graph.gfa.gz | python3 gfa_to_pairwise_paf.py --reference GRCh38#0 \
  --queries HG01109#1,HG00099#1 --contig-lengths contigs.fai > graph.paf

Index the PAF for JBrowse:

jbrowse make-pif graph.paf --csi --out graph.pif.gz

build_hprc_multiway_synteny.sh runs both on the release graph for eight haplotypes, with each haplotype's contig lengths and CAT genes, and its output is hosted, one lane per haplotype under GRCh38.

Reproduce it end to end

The one command builds the graph track, the bubbles, the tier and the allele inventory, with the tools under Prerequisites. It:

places every segment on a genome: an rGFA holds each position in tags on the segment, and for a plain GFA the command follows the backbone's paths first, then places each remaining segment on the first other haplotype whose path visits it
finds the bubbles, with gfatools bubble on an rGFA or from the snarl VCF on a plain GFA, and builds the tier from them
reads each allele out of the links, following it from where it leaves the backbone to where it rejoins. The reference between those two points and the sequence the allele passes through give the CIGAR its size

gfa-to-tabix build hprc-v2.1-mc-grch38.sv.gfa.gz -o hprc --assembly hg38

gfa-to-tabix build reads the graph once and writes the index, the tier, the alleles and the bubbles concurrently. The coarse tier is also gfa-to-tabix fold, which runs alone:

gfa-to-tabix fold hprc-v2.1-mc-grch38.sv.gfa.gz --below 10000 -o hprc.fold10000

The walk files come from the command under Indexing haplotype walks.

External links

Li H. The rGFA format and gfatools: the SN/SO/SR tags and the bubble calls.
gfa-to-tabix: --walks, and the rows it writes.
gfa-to-tabix and gbz-base compared: what each holds, sizes and requests on HPRC v2.1.
ggbandage and ggtubemap: R packages that draw a region of a walk file as ggplot2 layers.

Citations

HPRC release 2, the worked example here.

Notes

On a 24-core machine gfa-to-tabix took 39 minutes and 13 GB of memory, and wrote 2.5 GB for the two references. ↩

Feedback on this tutorial is welcome: contact us.

Pangenome (hosting your own graph)

Prerequisites#

Download the data#

The GraphGenomeView plugin#

Indexing a graph with gfa-to-tabix build#

The coarse tier for whole-chromosome views#

Configuring the graph track#

Checking the index against the graph#

Recording which haplotypes pass through each segment#

Indexing haplotype walks with gfa-to-tabix#

Whole-genome synteny from the walks, as PAF#

Reproduce it end to end#

See also#

External links#

Citations#

Notes

Prerequisites

Download the data

The GraphGenomeView plugin

Indexing a graph with gfa-to-tabix build

The coarse tier for whole-chromosome views

Configuring the graph track

Checking the index against the graph

Recording which haplotypes pass through each segment

Indexing haplotype walks with gfa-to-tabix

Whole-genome synteny from the walks, as PAF

Reproduce it end to end

See also

External links

Citations