Pangenome (hosting your own graph)

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 lane of segments
• check the index against the graph
• add two optional layers: the haplotypes on each segment, and the walk of each haplotype

Prerequisites#

• the GraphGenomeView plugin
• htslib (bgzip, tabix), bcftools, python3, sort
• gfatools and GNU awk, for an rGFA
• minigraph, for each assembly's path through the graph
• vg 1.69.0+, gbz-base and gbz-haplotype-index, for haplotype walks

Where the data comes from#

The GraphGenomeView plugin#

GraphGenomeView loads by URL, from a plugins array at the top of config.json (configuring plugins). The config the command writes holds this entry:

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

On JBrowse Desktop v5.0.0-beta.1 or later, install it once at Global plugins... → Add custom plugin: open Advanced options, paste that esmUrl into ESM build URL and leave the rest empty.

Indexing a graph with build_pangenome_graph.sh#

Fetch the script:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_pangenome_graph.sh

An rGFA, from minigraph or the minigraph stage of Minigraph-Cactus, needs an output prefix and the assembly name your config uses:

bash build_pangenome_graph.sh graph.rgfa.gz out --assembly hg38

A plain GFA, from pggb, odgi, vg or base-level Minigraph-Cactus, also needs the backbone sample and its bubbles. vg deconstruct writes the snarl VCF the bubbles come from, one record per top-level snarl against the reference path (pggb -V writes the same file). The VCF's CHROM must be the assembly's refName, so rename the PanSN path as the pggb tutorial does:

# -p: the reference path to decompose against
vg deconstruct -p K12#1#chr graph.gbz > 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

bash build_pangenome_graph.sh graph.gfa out --reference K12 --assembly K12 --snarls graph.snarls.vcf.gz

The rGFA route runs GNU awk, which is brew install gawk on macOS with gnubin first on PATH.¹ At human-chromosome scale, index the SV-resolution rGFA; a pggb graph's index does not finish there. An assembly graph from SPAdes or Flye has no reference to index against; open it in Bandage.

The command writes these files beside the prefix:

The bubble tier#

The tier draws each bubble as one node on the reference backbone, and folds bubbles under its threshold into the backbone. A bubble's size for that test is the larger of its reference span and its longest allele. The threshold is in the file name, 10,000 bp by default; --tier sets it, and a pggb graph defaults to 50, since most of its bubbles are single bases.

The graph track#

The config's first track is the graph. uri is the prefix, and coarse names the tier the track draws past aboveBpPerPx bp per pixel. assemblyNameToPanSN maps your assembly name to the graph's PanSN sample, and the command writes it only when the two differ.

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

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.tier10000", "aboveBpPerPx": 1014 }
  },
  "displayDefaults": { "showLabels": "none" },
  "displays": [
    { "type": "LinearGraphDisplay" },
    { "type": "LinearBasicDisplay" }
  ]
}'

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

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

✕

Open in JBrowse Web ↗Open in JBrowse Desktop ↗

This session loads a large dataset and takes a while to open.

Or rebuild it from your own data:

Desktop

Opening in Desktop needs JBrowse Desktop 5.0+, and saves the view as a reopenable session.

Nothing happens?

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

Copy

https://jbrowse.org/code/jb2/main/?config=test_data%2Fgraphgenomeview%2Fhprc.json&session=spec-%7B%22views%22%3A%5B%7B%22type%22%3A%22LinearGenomeView%22%2C%22assembly%22%3A%22hg38%22%2C%22loc%22%3A%22chr6%3A31%2C980%2C000-32%2C050%2C000%22%2C%22tracks%22%3A%5B%7B%22trackId%22%3A%22hg38_ncbiRefSeq_ucsc%22%2C%22type%22%3A%22LinearBasicDisplay%22%2C%22geneGlyphMode%22%3A%22longestCoding%22%2C%22displayMode%22%3A%22compact%22%2C%22height%22%3A70%7D%2C%7B%22trackId%22%3A%22hprc_minigraph_segments%22%2C%22type%22%3A%22LinearBasicDisplay%22%2C%22showLabels%22%3A%22none%22%2C%22height%22%3A110%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+pangenome%2Fprepare_graph_segments

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).This figure uses hg38.
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 “NCBI RefSeq genes (hg38)”.
   1. Track menu → Gene glyph → Representative transcript
   2. Track menu → Feature height → Compact
   3. Drag the bar at the bottom edge of the track to resize it (70px here).
3. Add your own track: File → Open track..., then paste a URL or choose a local file.This figure uses “HPRC release 2 graph (rGFA segments)”.
   1. Track menu → Show... → Labels → None
   2. Drag the bar at the bottom edge of the track to resize it (110px here).
4. Type your region of interest into the location box and press Enter.This figure is at chr6:31,980,000-32,050,000.

Web

Its genome and tracks are hosted, so File → Open track... there adds your own files beside them.

With your own data, the steps behind the figure:

1. Put your genome in your JBrowse Web's config.json: jbrowse add-assembly adds it.This figure uses hg38.
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 “NCBI RefSeq genes (hg38)”.
   1. Track menu → Gene glyph → Representative transcript
   2. Track menu → Feature height → Compact
   3. Drag the bar at the bottom edge of the track to resize it (70px here).
3. Add your own track: File → Open track..., then paste a URL or choose a local file.This figure uses “HPRC release 2 graph (rGFA segments)”.
   1. Track menu → Show... → Labels → None
   2. Drag the bar at the bottom edge of the track to resize it (110px here).
4. Type your region of interest into the location box and press Enter.This figure is at chr6:31,980,000-32,050,000.

Spec

This session spec draws the figure from test_data/graphgenomeview/hprc.json: paste it after &session=spec- on a JBrowse Web link that loads that config.

Copy

{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr6:31,980,000-32,050,000",
      "tracks": [
        {
          "trackId": "hg38_ncbiRefSeq_ucsc",
          "type": "LinearBasicDisplay",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 70
        },
        {
          "trackId": "hprc_minigraph_segments",
          "type": "LinearBasicDisplay",
          "showLabels": "none",
          "height": 110
        }
      ]
    }
  ]
}

Agent

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.

Copy

cat > session.json <<'JSON'
{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr6:31,980,000-32,050,000",
      "tracks": [
        {
          "trackId": "hg38_ncbiRefSeq_ucsc",
          "type": "LinearBasicDisplay",
          "geneGlyphMode": "longestCoding",
          "displayMode": "compact",
          "height": 70
        },
        {
          "trackId": "hprc_minigraph_segments",
          "type": "LinearBasicDisplay",
          "showLabels": "none",
          "height": 110
        }
      ]
    }
  ]
}
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/code/jb2/main/test_data/graphgenomeview/hprc.json \
  --width 1500 --height 454 \
  --timeout 180000 \
  --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.

The command writes out.config.json. Merge its plugins and tracks entries into your own config.

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

✕

Open in JBrowse Web ↗Open in JBrowse Desktop ↗

This session loads a large dataset and takes a while to open.

Or rebuild it from your own data:

Desktop

Opening in Desktop needs JBrowse Desktop 5.0+, and saves the view as a reopenable session.

Nothing happens?

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

Copy

https://jbrowse.org/code/jb2/main/?config=https%3A%2F%2Fjbrowse.org%2Fpangenome%2Fhprc-grch38%2Fconfig.json&session=spec-%7B%22views%22%3A%5B%7B%22type%22%3A%22LinearGenomeView%22%2C%22assembly%22%3A%22hg38%22%2C%22loc%22%3A%22chr6%3A31%2C980%2C000-32%2C050%2C000%22%2C%22tracks%22%3A%5B%22hprc_minigraph_bubbles%22%2C%22hprc_bubble_score%22%2C%22hprc_minigraph_alleles%22%2C%7B%22trackId%22%3A%22hprc_minigraph_segments%22%2C%22type%22%3A%22LinearGraphDisplay%22%2C%22colorScheme%22%3A%22reference-position%22%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+pangenome%2Fhost_your_own

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).This figure uses hg38, loaded from a .2bit genome.
2. Add your own track: File → Open track..., then paste a URL or choose a local file.This figure uses “HPRC release 2 bubbles”.
3. Add your own track: File → Open track..., then paste a URL or choose a local file.This figure uses “HPRC release 2 graph: variability (segments per bubble)”.
4. Add your own track: File → Open track..., then paste a URL or choose a local file. This one needs a bgzip-compressed, tabix-indexed BED (.bed.gz + .tbi).This figure uses “HPRC release 2 graph: allele inventory”.
5. Add your own track: File → Open track..., then paste a URL or choose a local file.This figure uses “HPRC release 2 graph (rGFA segments)”.
   1. Track menu → Color → Reference position
6. Type your region of interest into the location box and press Enter.This figure is at chr6:31,980,000-32,050,000.

Web

Its genome and tracks are hosted, so File → Open track... there adds your own files beside them.

With your own data, the steps behind the figure:

1. Put your genome in your JBrowse Web's config.json: jbrowse add-assembly adds it.This figure uses hg38, loaded from a .2bit genome.
2. Add your own track: File → Open track..., then paste a URL or choose a local file.This figure uses “HPRC release 2 bubbles”.
3. Add your own track: File → Open track..., then paste a URL or choose a local file.This figure uses “HPRC release 2 graph: variability (segments per bubble)”.
4. Add your own track: File → Open track..., then paste a URL or choose a local file. This one needs a bgzip-compressed, tabix-indexed BED (.bed.gz + .tbi).This figure uses “HPRC release 2 graph: allele inventory”.
5. Add your own track: File → Open track..., then paste a URL or choose a local file.This figure uses “HPRC release 2 graph (rGFA segments)”.
   1. Track menu → Color → Reference position
6. Type your region of interest into the location box and press Enter.This figure is at chr6:31,980,000-32,050,000.

Spec

This session spec draws the figure from https://jbrowse.org/pangenome/hprc-grch38/config.json: paste it after &session=spec- on a JBrowse Web link that loads that config.

Copy

{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr6:31,980,000-32,050,000",
      "tracks": [
        "hprc_minigraph_bubbles",
        "hprc_bubble_score",
        "hprc_minigraph_alleles",
        {
          "trackId": "hprc_minigraph_segments",
          "type": "LinearGraphDisplay",
          "colorScheme": "reference-position"
        }
      ]
    }
  ]
}

Agent

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.

Copy

cat > session.json <<'JSON'
{
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr6:31,980,000-32,050,000",
      "tracks": [
        "hprc_minigraph_bubbles",
        "hprc_bubble_score",
        "hprc_minigraph_alleles",
        {
          "trackId": "hprc_minigraph_segments",
          "type": "LinearGraphDisplay",
          "colorScheme": "reference-position"
        }
      ]
    }
  ]
}
JSON

npx @jbrowse/capture --instance https://jbrowse.org/code/jb2/main/ \
  --config https://jbrowse.org/pangenome/hprc-grch38/config.json \
  --width 1500 --height 1080 \
  --timeout 180000 \
  --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.

A node's right-click menu offers Open in the haplotype named in its rGFA id, such as NA20809#2#CM094351.1, when the session holds an assembly named or aliased sample#haplotype, here NA20809#2. 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. 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.

Which haplotypes walk each segment#

The rank in an rGFA is build order, so it names the first assembly a segment came from. The haplotypes whose paths walk each segment come 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

build_minigraph_paths.sh runs that per assembly and joins the output into one tabix-indexed row per bubble per sample, drawn as one lane per haplotype.

With a plain GFA, the command records the haplotypes whose paths visit each segment as an SM:Z: tag while it walks the paths. The node panel lists them as samples, and a track reads them as feature.samples and their count as feature.sampleCount. Color by... → Attribute... with sampleCount gives each count a separate colour. Past a handful of haplotypes a ramp reads better; Edit plot... in the same dialog takes one, here red for a segment on one haplotype to grey for a segment on most:

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

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

The E. coli pggb tutorial draws the count over an IS5 insertion.

The count is per haplotype (HG002.1), so a diploid sample's two copies count separately. The HPRC tutorial reads HPRC's phased VCF, which lists every haplotype's allele at every bubble.

Haplotype walks: a gbz-base database#

A .gbz is vg's indexed form of a graph, with one walk per haplotype. The browser reads it as a gbz-base database, the graph in SQLite, which three commands build.

Build the distance-index chains. vg 1.69.0 or newer reads them out of a distance index, and a top-level one (vg index --no-nested-distance) is enough:

vg chains graph.gbz graph.dist > graph.chains

Build the database, one row per node and per path, plus the chains. Without --chains, a window comes back as the reference walk alone.

gbz-base construct --chains graph.chains graph.gbz

Name the haplotypes. gbz-base reports the walks in a subgraph as unknown#1, unknown#2, and gbz-haplotype-index writes their names to a companion file. It reads the database beside the GBZ to check that the two match:

# --interval: bp between recorded GBWT positions per path; denser is bigger
#   and faster
# --anchor-spacing: bp between anchor nodes on the reference path, so a window
#   walks only the chosen lanes
gbz-haplotype-index --interval 16384 --anchor-spacing 131072 \
  graph.gbz graph.gbz.db graph.haplotype-index.db

cargo install gbz-base and cargo install gbz-haplotype-index install the two tools. The browser reads only the format 3 companion that gbz-haplotype-index 0.3.0 and later writes.²

Serve the database and the companion from URLs that answer range requests, and point the track's uri and haplotypeIndexLocation at them. 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:

{
  "type": "GraphTrack",
  "trackId": "my_graph_lanes",
  "name": "My graph: haplotypes vs the reference, read from the graph",
  "assemblyNames": ["hg38", "HG00097.1", "HG00099.1"],
  "adapter": {
    "type": "GbzBaseSyntenyAdapter",
    "uri": "https://example.com/graphs/my_graph.gbz.db",
    "haplotypeIndexLocation": {
      "uri": "https://example.com/graphs/my_graph.haplotype-index.db"
    },
    "assemblyNames": ["hg38"],
    "assemblyNameToPanSN": { "hg38": "GRCh38#0" },
    "context": 1000,
    "nodeLimit": 50000
  },
  "displays": [
    { "type": "MultiWaySyntenyDisplay", "height": 600 },
    { "type": "LinearGraphDisplay" }
  ]
}

jbrowse add-track-json '{
  "type": "GraphTrack",
  "trackId": "my_graph_lanes",
  "name": "My graph: haplotypes vs the reference, read from the graph",
  "assemblyNames": ["hg38", "HG00097.1", "HG00099.1"],
  "adapter": {
    "type": "GbzBaseSyntenyAdapter",
    "uri": "https://example.com/graphs/my_graph.gbz.db",
    "haplotypeIndexLocation": {
      "uri": "https://example.com/graphs/my_graph.haplotype-index.db"
    },
    "assemblyNames": ["hg38"],
    "assemblyNameToPanSN": { "hg38": "GRCh38#0" },
    "context": 1000,
    "nodeLimit": 50000
  },
  "displays": [
    { "type": "MultiWaySyntenyDisplay", "height": 600 },
    { "type": "LinearGraphDisplay" }
  ]
}'

{
  "type": "GraphTrack",
  "trackId": "my_graph_lanes",
  "name": "My graph: haplotypes vs the reference, read from the graph",
  "assemblyNames": ["hg38", "HG00097.1", "HG00099.1"],
  "adapter": {
    "type": "GbzBaseSyntenyAdapter",
    "uri": "https://example.com/graphs/my_graph.gbz.db",
    "haplotypeIndexLocation": {
      "uri": "https://example.com/graphs/my_graph.haplotype-index.db"
    },
    "assemblyNames": ["hg38"],
    "assemblyNameToPanSN": { "hg38": "GRCh38#0" },
    "context": 1000,
    "nodeLimit": 50000
  },
  "displays": [
    { "type": "MultiWaySyntenyDisplay", "height": 600 },
    { "type": "LinearGraphDisplay" }
  ]
}

The adapter rejects a companion built from a graph with a different path count. Past nodeLimit nodes, or 5 Mb for the Graph display, both displays ask the reader to zoom in. Haplotypes against each other draws the lanes this track produces.

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:

1. places every segment on a genome. An rGFA states each segment's sequence and offset in its own tags. For a plain GFA the command walks the backbone's paths first, so every segment they visit lands on the reference, and places each remaining segment on the first other haplotype that walks it
2. finds the bubbles, with gfatools bubble on an rGFA or from the snarl VCF on a plain GFA, and builds the tier from them
3. 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 walks give the CIGAR its size

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_pangenome_graph.sh
bash build_pangenome_graph.sh hprc-v2.1-mc-grch38.sv.gfa.gz hprc --assembly hg38

build_pangenome_graph.sh fetches and runs build_rgfa_tabix.sh or build_pggb_tabix.sh, then build_bubble_tier.sh and build_rgfa_alleles.sh, each runnable alone. HPRC publishes the gbz-base database itself, so a separate script builds only the haplotype-walk companion, from the 5.5 GB .gbz and the 10 GB database:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_hprc_gbz_index.sh
bash build_hprc_gbz_index.sh out

See also#

• Pangenome (HPRC) part 1: the graph's alleles and the haplotypes that have them
• Pangenome (HPRC) part 2: haplotypes against each other
• Pangenome (pggb)
• Pangenome (Minigraph-Cactus)
• Graph genome view

External links#

• Li H. The rGFA format and gfatools: the SN/SO/SR tags and the bubble calls.
• gbz-base: a GBZ as a SQLite database, range-requested per window.

Citations#

• HPRC release 2, the worked example here.

Feedback on this tutorial is welcome: contact us.