Pangenome (Minigraph-Cactus)

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.20. Desktop beta builds for Linux and macOS are on the releases page; the Windows build waits on code signing.

We build a pangenome graph of five E. coli strains with Minigraph-Cactus, which lays the other strains onto a backbone from one reference, and turn its outputs into tracks on the K12 genome. The pggb tutorial builds the same strains with no reference and takes them through the same sections, so the two builders can be compared at the same loci. One cactus-pangenome run writes the graph, a VCF, an odgi file, a HAL alignment and short-read indexes. We:

load the synteny, variants, whole-genome alignment (MAF), depth and per-strain presence as tracks on K12, and set both builders' depth in one track
load the short reads of an isolate outside the graph, mapped through it, as an alignments track
load the graph as a graph track in the force-directed layout
set odgi viz's picture of the graph beside the presence heatmap on K12

Experimental

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

Prerequisites

docker or singularity, for the cactus image, which includes odgi, halSynteny, hal2maf, vg and samtools
htslib (bgzip, tabix) and python3
node, for the JBrowse CLI
the GraphGenomeView plugin, for drawing the graph as a graph. In your own JBrowse, install it from the plugin store
gfa-to-tabix 2.4.0 or later, for indexing the graph
for the whole build: the NCBI datasets CLI, bedGraphToBigWig (UCSC kentUtils), samtools, unzip and wget

# Debian/Ubuntu
sudo apt install samtools tabix nodejs
# or conda, which has every tool above
conda install -c bioconda -c conda-forge samtools htslib nodejs \
  ucsc-bedgraphtobigwig ncbi-datasets-cli
# or Linux x86_64 binaries for the tools apt lacks
wget -P ~/.local/bin \
  https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/bedGraphToBigWig \
  https://ftp.ncbi.nlm.nih.gov/pub/datasets/command-line/v2/linux-amd64/datasets \
  https://ftp.ncbi.nlm.nih.gov/pub/datasets/command-line/v2/linux-amd64/dataformat
chmod +x ~/.local/bin/bedGraphToBigWig ~/.local/bin/datasets ~/.local/bin/dataformat
# any platform
cargo install gfa-to-tabix
npm install -g @jbrowse/cli@next

The data

The page's commands fetch the five E. coli RefSeq assemblies by accession and build local files from them; run the build script first, or follow the same commands on your own assemblies. The script fetches the short reads too, and the graph track reads our hosted copy.

What it is	File
K12	GCF_000005845.2_ASM584v2/
Sakai	GCF_000008865.2_ASM886v2/
CFT073	GCF_000007445.1_ASM744v1/
NCTC86	GCF_002007705.1_ASM200770v1/
IAI39	GCF_000026345.1_ASM2634v1/
KTa004 short reads mapped through the finished graph, forward mate	DRR063408_1.fastq.gz
KTa004 short reads, reverse mate	DRR063408_2.fastq.gz

Read by URL (no download needed)

What it is	File
The graph's segments and links, tabix-indexed, which the graph track below reads	ecoli_cactus.graph.bed.gz

Building the graph with cactus-pangenome

Each strain's FASTA holds its chromosome as one record named chr. Fetch and trim one assembly with the NCBI datasets CLI:

# $acc: the RefSeq accession; $strain: the name the graph gives the sample
datasets download genome accession "$acc" --filename "$strain.zip"
unzip -o "$strain.zip" -d "$strain"
# keep the first record, the chromosome, renamed chr; the rest are plasmids
awk '/^>/{n++; if (n == 1) print ">chr"; next} n == 1' \
  "$strain"/ncbi_dataset/data/*/*.fna > "$strain.fa"

Minigraph-Cactus (cactus-pangenome) builds a pangenome graph reference-first. minigraph lays down a backbone from the reference you pick, and Cactus aligns every other sample onto it and normalizes the result into a graph. It takes a seqFile: one name<TAB>path line per sample.

cat > seqfile.txt <<'EOF'
K12     K12.fa
Sakai   Sakai.fa
CFT073  CFT073.fa
NCTC86  NCTC86.fa
IAI39   IAI39.fa
EOF

Contigs keep their plain names here (chr). Cactus applies PanSN sample#haplotype#contig naming to the graph internally.

Every step runs in the same cactus image (Cactus also ships a binary release for a machine with no container runtime), so wrap the docker run once:

in_cactus() {
  docker run --rm -u "$(id -u):$(id -g)" -w /data -v "$PWD":/data \
    quay.io/comparative-genomics-toolkit/cactus:v3.2.1 "$@"
}

Under singularity, singularity exec --bind "$PWD":/data --pwd /data docker://<image> replaces the wrapper body. The build script picks the runtime off PATH.

--reference K12 makes K12 the minigraph backbone and the path every output is decomposed against:

in_cactus cactus-pangenome /data/js /data/seqfile.txt \
  --outDir /data/mc --outName ecoli --reference K12 \
  --vcf --gfa --gbz --odgi --viz --draw --giraffe --consCores 8

/data/js is the Toil job store, and must not already exist on a fresh run.
--outName ecoli prefixes every output file.
--vcf deconstructs the graph into the variants track's input.
--odgi writes the .og that the depth and presence tracks read.
--viz writes the odgi raster shown at the end.
--giraffe writes the indexes the read-mapping section needs (mc/ecoli.d2.gbz and its .dist/.min/.zipcodes), built over the graph filtered to sequence at least two haplotypes have.

Load the genomes

Every projection is a track on a strain's assembly. We'll add K12, the reference the projections map to, from a bgzip-compressed, indexed FASTA:

Config file

Goes in the assemblies array of config.json. See Assemblies.

{
  "name": "K12",
  "uri": "https://jbrowse.org/demos/ecoli_pangenome/K12.fa.gz"
}

CLI (add-assembly)

jbrowse add-assembly https://jbrowse.org/demos/ecoli_pangenome/K12.fa.gz \
  --name K12

JBrowse Desktop

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/demos/ecoli_pangenome/K12.fa.gz
https://jbrowse.org/demos/ecoli_pangenome/K12.fa.gz.fai
https://jbrowse.org/demos/ecoli_pangenome/K12.fa.gz.gzi

JBrowse reads the format off the file name. Then fill in:

Genome name: K12

Sakai, CFT073, NCTC86 and IAI39 load the same way, one assembly per strain, named as in the assemblyNames of the tracks below.

Synteny between all the strains (halSynteny)

halSynteny reads the HAL's base-level alignment and emits synteny blocks per genome pair as PSL with no sample tag in its sequence names. Run it once per pair of strains, and append each pair to one PAF with PanSN sample#0#chr names:

# $q, $t: the query and target strains of one pair
in_cactus halSynteny --queryGenome "$q" --targetGenome "$t" \
  /data/mc/ecoli.full.hal "/data/hs_${t}_${q}.psl"
# PSL to PAF; halSynteny keeps the query on +, so the 2nd strand character is
# the target's
awk -v OFS='\t' -v qn="${q}#0#chr" -v tn="${t}#0#chr" \
  '{ strand = (substr($9, 2, 1) == "-") ? "-" : "+";
     print qn, $11, $12, $13, strand, tn, $15, $16, $17, $1, ($13 - $12), 255 }' \
  "hs_${t}_${q}.psl" >> ecoli_cactus_ava.paf

Index the combined PAF so a range query fetches only the region in view:

jbrowse make-pif ecoli_cactus_ava.paf --csi   # -> ecoli_cactus_ava.pif.gz (+ .csi)

Load it with an MultiGenomeIndexedPAFAdapter, whose PanSN sample# prefix on every record is how it maps a record to its strain:

Config file

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

{
  "type": "SyntenyTrack",
  "trackId": "ecoli_cactus_ava",
  "name": "MC graph: all-vs-all synteny (halSynteny)",
  "assemblyNames": ["K12", "Sakai", "CFT073", "NCTC86", "IAI39"],
  "adapter": {
    "type": "MultiGenomeIndexedPAFAdapter",
    "uri": "ecoli_cactus_ava.pif.gz",
    "csi": true,
    "assemblyNames": ["K12", "Sakai", "CFT073", "NCTC86", "IAI39"]
  }
}

CLI (add-track-json)

jbrowse add-track-json '{
  "type": "SyntenyTrack",
  "trackId": "ecoli_cactus_ava",
  "name": "MC graph: all-vs-all synteny (halSynteny)",
  "assemblyNames": ["K12", "Sakai", "CFT073", "NCTC86", "IAI39"],
  "adapter": {
    "type": "MultiGenomeIndexedPAFAdapter",
    "uri": "ecoli_cactus_ava.pif.gz",
    "csi": true,
    "assemblyNames": ["K12", "Sakai", "CFT073", "NCTC86", "IAI39"]
  }
}'

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": "SyntenyTrack",
  "trackId": "ecoli_cactus_ava",
  "name": "MC graph: all-vs-all synteny (halSynteny)",
  "assemblyNames": ["K12", "Sakai", "CFT073", "NCTC86", "IAI39"],
  "adapter": {
    "type": "MultiGenomeIndexedPAFAdapter",
    "uri": "ecoli_cactus_ava.pif.gz",
    "csi": true,
    "assemblyNames": ["K12", "Sakai", "CFT073", "NCTC86", "IAI39"]
  }
}

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

That shows it in the linear view. For the synteny view, open Add → Linear synteny view, pick the track under Quick start, and click Launch.

To stack the strains, a linear synteny view takes one panel per strain and names the track once in tracks, which puts it on every band. Add → Linear synteny view, whose Quick start fills in a row per assembly the track lists, opens the same view as this session:

Config file

Goes at the top level of config.json, replacing any defaultSession there. See Default session.

{
  "defaultSession": {
    "name": "Minigraph-Cactus strains stacked",
    "views": [
      {
        "type": "LinearSyntenyView",
        "views": [
          { "assembly": "K12" },
          { "assembly": "Sakai" },
          { "assembly": "CFT073" },
          { "assembly": "NCTC86" },
          { "assembly": "IAI39" }
        ],
        "tracks": ["ecoli_cactus_ava"],
        "filter": [{ "field": "length", "gte": 10000 }],
        "levelHeights": [110, 110, 110, 110],
        "encoding": { "color": "strand" }
      }
    ]
  }
}

CLI (set-default-session)

jbrowse set-default-session --session - << 'EOF'
{
  "name": "Minigraph-Cactus strains stacked",
  "views": [
    {
      "type": "LinearSyntenyView",
      "views": [
        { "assembly": "K12" },
        { "assembly": "Sakai" },
        { "assembly": "CFT073" },
        { "assembly": "NCTC86" },
        { "assembly": "IAI39" }
      ],
      "tracks": ["ecoli_cactus_ava"],
      "filter": [{ "field": "length", "gte": 10000 }],
      "levelHeights": [110, 110, 110, 110],
      "encoding": { "color": "strand" }
    }
  ]
}
EOF

Open in JBrowse

Open this session in JBrowse Web ↗

In grammar-of-graphics terms, encoding.color is a color scale that maps each ribbon's strand to red where the two strains run the same way and blue where one is inverted against the other.

The Minigraph-Cactus graph's halSynteny projection, the strains stacked and colored by strand.

Variants against K12 (vg deconstruct)

--vcf decomposes the graph against K12 with vg deconstruct, genotyped across the other four strains. Its CHROM is already chr, so the .gz and .tbi Cactus wrote load unchanged as a VariantTrack:

Config file

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

{
  "type": "VariantTrack",
  "trackId": "ecoli_cactus_variants",
  "name": "MC graph: pangenome variants (vs K12)",
  "uri": "mc/ecoli.vcf.gz",
  "assemblyNames": ["K12"],
  "displays": [
    {
      "type": "LinearVariantDisplay",
      "variantLayout": "columns",
      "rows": "haplotype"
    }
  ]
}

CLI (add-track-json)

jbrowse add-track-json '{
  "type": "VariantTrack",
  "trackId": "ecoli_cactus_variants",
  "name": "MC graph: pangenome variants (vs K12)",
  "uri": "mc/ecoli.vcf.gz",
  "assemblyNames": ["K12"],
  "displays": [
    {
      "type": "LinearVariantDisplay",
      "variantLayout": "columns",
      "rows": "haplotype"
    }
  ]
}'

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": "VariantTrack",
  "trackId": "ecoli_cactus_variants",
  "name": "MC graph: pangenome variants (vs K12)",
  "uri": "mc/ecoli.vcf.gz",
  "assemblyNames": ["K12"],
  "displays": [
    {
      "type": "LinearVariantDisplay",
      "variantLayout": "columns",
      "rows": "haplotype"
    }
  ]
}

mc/ecoli.vcf.gz is relative to a config.json. Replace it with its URL or its path on this computer.

rows: "haplotype" draws one row per strain, and variantLayout: "columns" one column per variant. The multi-sample variant track guide covers columns versus genomic positions, the genotype colors, and clustering samples by genotype.

vg deconstruct writes a record per snarl (vg's word for a bubble) at every nesting level, so wide records paint over finer ones. cactus-pangenome prunes that tree with vcfbub by default; --vcfbub 0 turns that off, and --vcfwave realigns the remaining records into primitive variants.

Whole-genome alignment onto K12 (MAF)

The MAF comes out of mc/ecoli.full.hal. hal2maf --refGenome K12 roots every block on K12, and the rows come out as K12.chr, Sakai.chr, and so on, the sample.contig naming the MAF display splits species on:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/maf_to_bed.py
in_cactus hal2maf --refGenome K12 --noAncestors /data/mc/ecoli.full.hal /data/ecoli_cactus.maf
python3 maf_to_bed.py ecoli_cactus.maf ecoli_cactus.maf.bed
bgzip ecoli_cactus.maf.bed
tabix -C -p bed ecoli_cactus.maf.bed.gz

maf_to_bed.py writes one line per block, which a MafTabixAdapter reads. The streaming maf2bed converts this MAF too; see producing the tabix BED.

Config file

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

{
  "type": "MafTrack",
  "trackId": "ecoli_cactus_maf",
  "name": "MC graph: whole-genome alignment (MAF, vs K12)",
  "assemblyNames": ["K12"],
  "adapter": {
    "type": "MafTabixAdapter",
    "samples": ["K12", "Sakai", "CFT073", "NCTC86", "IAI39"],
    "uri": "ecoli_cactus.maf.bed.gz",
    "csi": true
  }
}

CLI (add-track-json)

jbrowse add-track-json '{
  "type": "MafTrack",
  "trackId": "ecoli_cactus_maf",
  "name": "MC graph: whole-genome alignment (MAF, vs K12)",
  "assemblyNames": ["K12"],
  "adapter": {
    "type": "MafTabixAdapter",
    "samples": ["K12", "Sakai", "CFT073", "NCTC86", "IAI39"],
    "uri": "ecoli_cactus.maf.bed.gz",
    "csi": true
  }
}'

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": "MafTrack",
  "trackId": "ecoli_cactus_maf",
  "name": "MC graph: whole-genome alignment (MAF, vs K12)",
  "assemblyNames": ["K12"],
  "adapter": {
    "type": "MafTabixAdapter",
    "samples": ["K12", "Sakai", "CFT073", "NCTC86", "IAI39"],
    "uri": "ecoli_cactus.maf.bed.gz",
    "csi": true
  }
}

ecoli_cactus.maf.bed.gz is relative to a config.json. Replace it with its URL or its path on this computer.

The HAL projected onto K12 as a MAF, one row per strain, colored where each differs from K12.

A row with no colored columns shares its sequence with K12, and a row that stops has no alignment to K12 there. The coverage band tells the two apart.

samples names the rows and fixes their order. To order them by shared graph content instead, build a tree from odgi similarity and point nhLocation at it, as the pggb tutorial's MAF track does:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/odgi_similarity_to_newick.py
# -D '#' -p 1: group the paths by sample, the first field of each PanSN name
in_cactus odgi similarity -i /data/mc/ecoli.full.og -D '#' -p 1 > ecoli_cactus_similarity.tsv
python3 odgi_similarity_to_newick.py ecoli_cactus_similarity.tsv ecoli_cactus.nh

Then add "nhLocation": { "uri": "ecoli_cactus.nh" } to the MAF track's adapter.

Pangenome depth and per-strain presence

odgi depth counts how many paths traverse the graph under each K12 base, and odgi pav splits that per strain. Both run over mc/ecoli.full.og in 500 bp windows of K12's path, which Cactus names K12#0#chr:

# the windows, and the reference length bedGraphToBigWig needs
reflen=$(awk '!/^>/ { c += length($0) } END { print c }' K12.fa)
printf 'chr\t%s\n' "$reflen" > chrom.sizes
awk -v p="K12#0#chr" -v len="$reflen" -v w=500 \
  'BEGIN { for (s = 0; s < len; s += w) { e = s + w; if (e > len) e = len
           print p "\t" s "\t" e } }' > depth_windows.bed

in_cactus odgi depth -i /data/mc/ecoli.full.og -b /data/depth_windows.bed |
  awk -v p="K12#0#chr" -v OFS='\t' '$1 == p && $4 + 0 == $4 { print "chr", $2, $3, $4 }' |
  sort -k1,1 -k2,2n > ecoli_cactus_depth.bedgraph
bedGraphToBigWig ecoli_cactus_depth.bedgraph chrom.sizes ecoli_cactus_depth.bw

Each other strain's path name ends in a subpath tag (Sakai#0#chr#0), so the per-strain filter matches the sample prefix:

in_cactus odgi pav -i /data/mc/ecoli.full.og -b /data/depth_windows.bed > pav.tsv
for strain in Sakai CFT073 NCTC86 IAI39; do
  awk -F'\t' -v OFS='\t' -v s="$strain" \
    '$5 ~ "^" s "#" && $6 + 0 == $6 { print "chr", $2, $3, $6 }' pav.tsv |
    sort -k1,1 -k2,2n > "pav_${strain}.bedgraph"
  bedGraphToBigWig "pav_${strain}.bedgraph" chrom.sizes "ecoli_cactus_pav_${strain}.bw"
done

The four bigWigs load as one heatmap track, a MultiQuantitativeTrack with a row per strain, where a white column is a stretch that strain lacks:

Config file

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

{
  "type": "MultiQuantitativeTrack",
  "trackId": "ecoli_cactus_pav",
  "name": "MC graph: per-strain presence (odgi pav, vs K12)",
  "assemblyNames": ["K12"],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "subadapters": [
      {
        "type": "BigWigAdapter",
        "name": "Sakai",
        "uri": "ecoli_cactus_pav_Sakai.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "CFT073",
        "uri": "ecoli_cactus_pav_CFT073.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "NCTC86",
        "uri": "ecoli_cactus_pav_NCTC86.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "IAI39",
        "uri": "ecoli_cactus_pav_IAI39.bw"
      }
    ]
  },
  "displayDefaults": { "mark": "rect" }
}

CLI (add-track-json)

jbrowse add-track-json '{
  "type": "MultiQuantitativeTrack",
  "trackId": "ecoli_cactus_pav",
  "name": "MC graph: per-strain presence (odgi pav, vs K12)",
  "assemblyNames": ["K12"],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "subadapters": [
      {
        "type": "BigWigAdapter",
        "name": "Sakai",
        "uri": "ecoli_cactus_pav_Sakai.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "CFT073",
        "uri": "ecoli_cactus_pav_CFT073.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "NCTC86",
        "uri": "ecoli_cactus_pav_NCTC86.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "IAI39",
        "uri": "ecoli_cactus_pav_IAI39.bw"
      }
    ]
  },
  "displayDefaults": { "mark": "rect" }
}'

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": "MultiQuantitativeTrack",
  "trackId": "ecoli_cactus_pav",
  "name": "MC graph: per-strain presence (odgi pav, vs K12)",
  "assemblyNames": ["K12"],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "subadapters": [
      {
        "type": "BigWigAdapter",
        "name": "Sakai",
        "uri": "ecoli_cactus_pav_Sakai.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "CFT073",
        "uri": "ecoli_cactus_pav_CFT073.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "NCTC86",
        "uri": "ecoli_cactus_pav_NCTC86.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "IAI39",
        "uri": "ecoli_cactus_pav_IAI39.bw"
      }
    ]
  },
  "displayDefaults": { "mark": "rect" }
}

ecoli_cactus_pav_Sakai.bw, ecoli_cactus_pav_CFT073.bw, ecoli_cactus_pav_NCTC86.bw, ecoli_cactus_pav_IAI39.bw are relative to a config.json. Replace each with its URL or its path on this computer.

Depth counts path steps, so a repeat the graph folded onto one run of nodes reads above the strain count. seqwish, pggb's graph builder, folds the rrnC and rrnB operons onto each other; the reference-first graph keeps them apart.

To compare the two builders, load both depth curves as rows of one track on a fixed axis, and type chr:3,935,000-3,955,000, the rrnC operon (a ribosomal RNA gene copy):

Config file

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

{
  "type": "MultiQuantitativeTrack",
  "trackId": "ecoli_depth_by_builder",
  "name": "Pangenome depth over K12, by builder (odgi depth)",
  "assemblyNames": ["K12"],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "subadapters": [
      {
        "type": "BigWigAdapter",
        "name": "pggb",
        "uri": "https://jbrowse.org/demos/ecoli_pangenome/ecoli_pggb_depth.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "Minigraph-Cactus",
        "uri": "https://jbrowse.org/demos/ecoli_pangenome/ecoli_cactus_depth.bw"
      }
    ]
  },
  "displayDefaults": {
    "encoding": { "y": { "scale": { "domainMin": 0, "domainMax": 10 } } }
  }
}

CLI (add-track-json)

jbrowse add-track-json '{
  "type": "MultiQuantitativeTrack",
  "trackId": "ecoli_depth_by_builder",
  "name": "Pangenome depth over K12, by builder (odgi depth)",
  "assemblyNames": ["K12"],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "subadapters": [
      {
        "type": "BigWigAdapter",
        "name": "pggb",
        "uri": "https://jbrowse.org/demos/ecoli_pangenome/ecoli_pggb_depth.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "Minigraph-Cactus",
        "uri": "https://jbrowse.org/demos/ecoli_pangenome/ecoli_cactus_depth.bw"
      }
    ]
  },
  "displayDefaults": {
    "encoding": { "y": { "scale": { "domainMin": 0, "domainMax": 10 } } }
  }
}'

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": "MultiQuantitativeTrack",
  "trackId": "ecoli_depth_by_builder",
  "name": "Pangenome depth over K12, by builder (odgi depth)",
  "assemblyNames": ["K12"],
  "adapter": {
    "type": "MultiWiggleAdapter",
    "subadapters": [
      {
        "type": "BigWigAdapter",
        "name": "pggb",
        "uri": "https://jbrowse.org/demos/ecoli_pangenome/ecoli_pggb_depth.bw"
      },
      {
        "type": "BigWigAdapter",
        "name": "Minigraph-Cactus",
        "uri": "https://jbrowse.org/demos/ecoli_pangenome/ecoli_cactus_depth.bw"
      }
    ]
  },
  "displayDefaults": {
    "encoding": { "y": { "scale": { "domainMin": 0, "domainMax": 10 } } }
  }
}

odgi depth over the rrnC operon against each builder's graph, on one axis.

Mapping a new isolate through the graph

We'll map the short reads of a sample outside the graph through the whole pangenome and project the result onto K12. A read that follows another strain's path through a variable region has no mismatches or soft clips and maps to that region's K12 coordinates. A read over sequence K12 lacks has no K12 coordinate and stays unmapped.

--giraffe wrote the indexes during the build. vg giraffe emits a GAM, and vg surject projects it onto one path as a BAM. The reads are E. coli KTa004 (ENA DRR063408, Illumina MiSeq), a strain the graph has never seen. The first 150,000 pairs of each mate are enough for the depth to read:

mkdir -p reads
for r in 1 2; do
  # awk reads to the end; head closes the pipe early and gzip dies of SIGPIPE
  gzip -dc "DRR063408_$r.fastq.gz" | awk 'NR<=600000' | gzip > "reads/sub_$r.fastq.gz"
done

# -Z/-d/-m/-z are the graph and the three indexes --giraffe built; -f takes
# the paired FASTQs, one flag per mate; -p shows progress
in_cactus vg giraffe -p \
  -Z /data/mc/ecoli.d2.gbz -d /data/mc/ecoli.d2.dist \
  -m /data/mc/ecoli.d2.shortread.withzip.min -z /data/mc/ecoli.d2.shortread.zipcodes \
  -f /data/reads/sub_1.fastq.gz -f /data/reads/sub_2.fastq.gz > mapped.gam

# -x is the graph the GAM was mapped against; -b writes BAM instead of GAM;
# -p is the path to surject onto, so the BAM's one reference sequence is that
# PanSN name; -N/-R set the sample name and read group on every read
in_cactus vg surject -x /data/mc/ecoli.d2.gbz -b -p K12#0#chr \
  -N KTa004 -R KTa004 /data/mapped.gam > mapped.raw.bam

The BAM's one reference sequence is the PanSN path name, which matches no assembly. Rename it to the assembly's refName (chr) in the header, drop the unmapped reads, then sort and index:

samtools view -H mapped.raw.bam | sed 's|SN:K12#0#chr|SN:chr|' > reads_hdr.sam
samtools reheader reads_hdr.sam mapped.raw.bam > reads_reheader.bam
samtools view -b -F 4 reads_reheader.bam > reads_mapped.bam
samtools sort -o ecoli_cactus_reads.bam reads_mapped.bam
samtools index ecoli_cactus_reads.bam

The result loads as an ordinary alignments track:

Config file

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

{
  "trackId": "ecoli_cactus_reads",
  "name": "KTa004 reads mapped through the graph (vs K12)",
  "uri": "ecoli_cactus_reads.bam",
  "assemblyNames": ["K12"]
}

CLI (add-track)

jbrowse add-track ecoli_cactus_reads.bam \
  --trackId ecoli_cactus_reads \
  --name "KTa004 reads mapped through the graph (vs K12)" \
  --assemblyNames K12 \
  --load copy

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... and, in Add a track from file or URL, enter:

Main file: ecoli_cactus_reads.bam

Click Next. JBrowse reads the adapter and track type off the file name. Then fill in:

Track name: KTa004 reads mapped through the graph (vs K12)
Assembly: K12

Click Add.

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

Turn on K12 genes and type chr:555,000-595,000. KTa004's reads cover K12 on either side of the prophage that begins at intD and nearly vanish across it, sequence K12 carries and KTa004 lacks.

KTa004's short reads, mapped through the graph and surjected onto K12, across the prophage that begins at intD.

vg giraffe rewrites ecoli.d2.dist as it runs, leaving it newer than the .min and .zipcodes built from it, so a second run refuses to start. touch the two derived files before re-mapping.

Drawing the graph as a graph

The graph genome view plugin draws the graph as a track of the linear view. Pangenome (hosting your own graph) covers installing the plugin and the one command that indexes a graph for it.

Indexing the graph

gfa-to-tabix build indexes the base-level graph, where every SNP is a bubble. That GFA has no rGFA tags, so --reference makes the command place segments by following the reference's path lines:

gfa-to-tabix build mc/ecoli.gfa.gz -o ecoli_cactus --reference K12

Without --snarls the command writes no bubble file. The index is in the contig layout. The uri below is our hosted copy.

Config file

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

{
  "type": "GraphTrack",
  "trackId": "ecoli_cactus_segments",
  "name": "MC graph: segments (whole graph, by locus)",
  "assemblyNames": ["K12"],
  "adapter": {
    "type": "RgfaTabixAdapter",
    "uri": "https://jbrowse.org/demos/ecoli_pangenome/ecoli_cactus",
    "format": 2
  },
  "displayDefaults": {
    "labels": { "name": { "show": false }, "description": { "show": false } }
  },
  "displays": [
    { "type": "LinearGraphDisplay" },
    { "type": "LinearBasicDisplay" }
  ]
}

CLI (add-track-json)

jbrowse add-track-json '{
  "type": "GraphTrack",
  "trackId": "ecoli_cactus_segments",
  "name": "MC graph: segments (whole graph, by locus)",
  "assemblyNames": ["K12"],
  "adapter": {
    "type": "RgfaTabixAdapter",
    "uri": "https://jbrowse.org/demos/ecoli_pangenome/ecoli_cactus",
    "format": 2
  },
  "displayDefaults": {
    "labels": { "name": { "show": false }, "description": { "show": false } }
  },
  "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": "ecoli_cactus_segments",
  "name": "MC graph: segments (whole graph, by locus)",
  "assemblyNames": ["K12"],
  "adapter": {
    "type": "RgfaTabixAdapter",
    "uri": "https://jbrowse.org/demos/ecoli_pangenome/ecoli_cactus",
    "format": 2
  },
  "displayDefaults": {
    "labels": { "name": { "show": false }, "description": { "show": false } }
  },
  "displays": [
    { "type": "LinearGraphDisplay" },
    { "type": "LinearBasicDisplay" }
  ]
}

Turned on, the track draws the window as a graph on K12's coordinates; Display types → Feature display in its menu draws the segments as a row.

Adding the graph track to an empty K12 session with Open track... → Add pangenome graph track.

Video

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

Load the graph by URL from File, Open track
The graph draws under the genes, cut from the window on screen
Select the IS1 element past flhD and zoom to it
At IS1 the graph re-cuts: the element is one node the other strains skip

Steps

Open this view in JBrowse Web ↗ ·

Open in JBrowse Desktop ↗

Needs JBrowse Desktop 5.0 or newer. If nothing opens, copy the link and paste it into Desktop's Open JBrowse Web link..., on the start screen or under File → Session.

https://jbrowse.org/code/jb2/main/?config=test_data%2Fgraphgenomeview%2Fconfig.json&session=spec-%7B%22sessionTracks%22%3A%5B%7B%22type%22%3A%22FeatureTrack%22%2C%22trackId%22%3A%22K12_genes%22%2C%22name%22%3A%22K12+genes%22%2C%22assemblyNames%22%3A%5B%22K12%22%5D%2C%22adapter%22%3A%7B%22type%22%3A%22Gff3TabixAdapter%22%2C%22gffGzLocation%22%3A%7B%22uri%22%3A%22https%3A%2F%2Fjbrowse.org%2Fdemos%2Fecoli_pangenome%2FK12.gff.gz%22%7D%2C%22index%22%3A%7B%22location%22%3A%7B%22uri%22%3A%22https%3A%2F%2Fjbrowse.org%2Fdemos%2Fecoli_pangenome%2FK12.gff.gz.csi%22%7D%2C%22indexType%22%3A%22CSI%22%7D%7D%7D%5D%2C%22views%22%3A%5B%7B%22type%22%3A%22LinearGenomeView%22%2C%22assembly%22%3A%22K12%22%2C%22loc%22%3A%22chr%3A1%2C972%2C900-1%2C984%2C900%22%2C%22tracks%22%3A%5B%7B%22trackId%22%3A%22K12_genes%22%2C%22type%22%3A%22LinearBasicDisplay%22%2C%22height%22%3A70%7D%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+pangenome_cactus%2Fsubgraph_launch

Quickest: open the view, then File → Open track... for your files.

From scratch:

Open your genome: Open new genome in JBrowse Desktop, or jbrowse add-assembly for JBrowse Web.This figure uses K12, loaded from a bgzip-compressed FASTA (.fa.gz + .fai + .gzi).
Add your track: File → Open track..., or check it in the track selector. Needs a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi).This figure uses “K12 genes”.
Go to your region in the location box.This figure is at chr:1,972,900-1,984,900.

Agent

Paste into a coding agent, with your goal on the first line.

My goal: [write yours, e.g. this view of my own files, or these track settings in my config.json]

Start from this JBrowse view.
- Config it loads: https://jbrowse.org/code/jb2/main/test_data/graphgenomeview/config.json
- Tracks:
  - K12 genes (K12_genes), FeatureTrack, a bgzip-compressed, tabix-indexed GFF3 (.gff.gz + .tbi), drawn as LinearBasicDisplay

Docs: https://jbrowse.org/jb2/docs/agents_authoring/ (writing a config), https://jbrowse.org/jb2/docs/agents/ (driving a running JBrowse). Use the jbrowse-authoring, jbrowse-hosted-data and jbrowse-capture skills if you have them. Check any config with `npx @jbrowse/cli validate`, then render it and look: a config that loads is not yet a track that drew.

The session spec that draws the view:
{
  "sessionTracks": [
    {
      "type": "FeatureTrack",
      "trackId": "K12_genes",
      "name": "K12 genes",
      "assemblyNames": [
        "K12"
      ],
      "adapter": {
        "type": "Gff3TabixAdapter",
        "gffGzLocation": {
          "uri": "https://jbrowse.org/demos/ecoli_pangenome/K12.gff.gz"
        },
        "index": {
          "location": {
            "uri": "https://jbrowse.org/demos/ecoli_pangenome/K12.gff.gz.csi"
          },
          "indexType": "CSI"
        }
      }
    }
  ],
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "K12",
      "loc": "chr:1,972,900-1,984,900",
      "tracks": [
        {
          "trackId": "K12_genes",
          "type": "LinearBasicDisplay",
          "height": 70
        }
      ]
    }
  ]
}

Redraw the figure headlessly

cat > session.json <<'JSON'
{
  "sessionTracks": [
    {
      "type": "FeatureTrack",
      "trackId": "K12_genes",
      "name": "K12 genes",
      "assemblyNames": [
        "K12"
      ],
      "adapter": {
        "type": "Gff3TabixAdapter",
        "gffGzLocation": {
          "uri": "https://jbrowse.org/demos/ecoli_pangenome/K12.gff.gz"
        },
        "index": {
          "location": {
            "uri": "https://jbrowse.org/demos/ecoli_pangenome/K12.gff.gz.csi"
          },
          "indexType": "CSI"
        }
      }
    }
  ],
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "K12",
      "loc": "chr:1,972,900-1,984,900",
      "tracks": [
        {
          "trackId": "K12_genes",
          "type": "LinearBasicDisplay",
          "height": 70
        }
      ]
    }
  ]
}
JSON

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

More: Using JBrowse with AI agents.

CLI

Adds these tracks to your config.json, styled as in the figure. Point each uri at your file.

jbrowse add-track-json '{
  "type": "FeatureTrack",
  "trackId": "K12_genes",
  "name": "K12 genes",
  "assemblyNames": [
    "K12"
  ],
  "adapter": {
    "type": "Gff3TabixAdapter",
    "gffGzLocation": {
      "uri": "https://jbrowse.org/demos/ecoli_pangenome/K12.gff.gz"
    },
    "index": {
      "location": {
        "uri": "https://jbrowse.org/demos/ecoli_pangenome/K12.gff.gz.csi"
      },
      "indexType": "CSI"
    }
  },
  "displayDefaults": {
    "height": 70
  }
}'

Run beside config.json, or add --out <dir>. Needs K12 in the config (jbrowse add-assembly).

Notebook

This view in a notebook, with jbrowse-anywidget. Swap a track's uri for your file.

from jbrowse_anywidget import JBrowseApp, fetch_hub

config = fetch_hub("https://jbrowse.org/code/jb2/main/test_data/graphgenomeview/config.json")
app = JBrowseApp(
    assemblies=config["assemblies"],
    tracks=config.get("tracks", []) + [
        {
            "type": "FeatureTrack",
            "trackId": "K12_genes",
            "name": "K12 genes",
            "assemblyNames": [
                "K12"
            ],
            "adapter": {
                "type": "Gff3TabixAdapter",
                "gffGzLocation": {
                    "uri": "https://jbrowse.org/demos/ecoli_pangenome/K12.gff.gz"
                },
                "index": {
                    "location": {
                        "uri": "https://jbrowse.org/demos/ecoli_pangenome/K12.gff.gz.csi"
                    },
                    "indexType": "CSI"
                }
            }
        }
    ],
    views=[
        {
            "type": "LinearGenomeView",
            "assembly": "K12",
            "loc": "chr:1,972,900-1,984,900",
            "tracks": [
                {
                    "trackId": "K12_genes",
                    "type": "LinearBasicDisplay",
                    "height": 70
                }
            ]
        }
    ],
)
app  # display the widget

Image

jbrowse-img draws this view as an SVG, no browser needed. Swap --config or a track for your own.

cat > tracks.json <<'JSON'
[
  {
    "type": "FeatureTrack",
    "trackId": "K12_genes",
    "name": "K12 genes",
    "assemblyNames": [
      "K12"
    ],
    "adapter": {
      "type": "Gff3TabixAdapter",
      "gffGzLocation": {
        "uri": "https://jbrowse.org/demos/ecoli_pangenome/K12.gff.gz"
      },
      "index": {
        "location": {
          "uri": "https://jbrowse.org/demos/ecoli_pangenome/K12.gff.gz.csi"
        },
        "indexType": "CSI"
      }
    }
  }
]
JSON

npx @jbrowse/img --config https://jbrowse.org/code/jb2/main/test_data/graphgenomeview/config.json \
  --tracks tracks.json \
  --assembly K12 --loc chr:1,972,900-1,984,900 \
  --track K12_genes display:LinearBasicDisplay '{"height":70}' \
  --width 1920 --out figure.svg

For a PNG, --out figure.png (needs rsvg-convert: apt install librsvg2-bin, brew install librsvg).

Session

The session spec behind the figure, over this config.

{
  "sessionTracks": [
    {
      "type": "FeatureTrack",
      "trackId": "K12_genes",
      "name": "K12 genes",
      "assemblyNames": [
        "K12"
      ],
      "adapter": {
        "type": "Gff3TabixAdapter",
        "gffGzLocation": {
          "uri": "https://jbrowse.org/demos/ecoli_pangenome/K12.gff.gz"
        },
        "index": {
          "location": {
            "uri": "https://jbrowse.org/demos/ecoli_pangenome/K12.gff.gz.csi"
          },
          "indexType": "CSI"
        }
      }
    }
  ],
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "K12",
      "loc": "chr:1,972,900-1,984,900",
      "tracks": [
        {
          "trackId": "K12_genes",
          "type": "LinearBasicDisplay",
          "height": 70
        }
      ]
    }
  ]
}

The index records which strains pass through each segment as an SB:H: bitmask, and a feature track over the same files colors each segment by that count:

Config file

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

{
  "type": "FeatureTrack",
  "trackId": "ecoli_cactus_carriage",
  "name": "MC graph: strains per segment",
  "assemblyNames": ["K12"],
  "adapter": {
    "type": "RgfaTabixAdapter",
    "uri": "https://jbrowse.org/demos/ecoli_pangenome/ecoli_cactus",
    "format": 2
  },
  "displayDefaults": {
    "displayMode": "collapsed",
    "labels": { "name": { "show": false }, "description": { "show": false } },
    "encoding": {
      "color": {
        "field": "sampleCount",
        "scale": {
          "domain": ["5", "4", "3", "2", "1"],
          "range": ["#bdbdbd", "#fed976", "#feb24c", "#fd8d3c", "#e31a1c"],
          "names": [
            "All 5 strains (core)",
            "4 strains",
            "3 strains",
            "2 strains",
            "1 strain (private)"
          ]
        },
        "legend": { "title": "Strains" }
      }
    }
  }
}

CLI (add-track-json)

jbrowse add-track-json '{
  "type": "FeatureTrack",
  "trackId": "ecoli_cactus_carriage",
  "name": "MC graph: strains per segment",
  "assemblyNames": ["K12"],
  "adapter": {
    "type": "RgfaTabixAdapter",
    "uri": "https://jbrowse.org/demos/ecoli_pangenome/ecoli_cactus",
    "format": 2
  },
  "displayDefaults": {
    "displayMode": "collapsed",
    "labels": { "name": { "show": false }, "description": { "show": false } },
    "encoding": {
      "color": {
        "field": "sampleCount",
        "scale": {
          "domain": ["5", "4", "3", "2", "1"],
          "range": ["#bdbdbd", "#fed976", "#feb24c", "#fd8d3c", "#e31a1c"],
          "names": [
            "All 5 strains (core)",
            "4 strains",
            "3 strains",
            "2 strains",
            "1 strain (private)"
          ]
        },
        "legend": { "title": "Strains" }
      }
    }
  }
}'

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: Indexed GFA (gfa-to-tabix .segs.bed.gz or .graph.bed.gz)
Path: https://jbrowse.org/demos/ecoli_pangenome/ecoli_cactus.graph.bed.gz
Track name: MC graph: strains per segment
Assembly: K12

Past the flagellar operon, K12 has an IS1 element (a short insertion sequence) the other four skip. To see it:

Type chr:1,978,100-1,979,700.
Turn on the strains-per-segment track.
Show the MAF track above the graph: the strains-per-segment track counts the strains, and the MAF's rows name them.

1.6 kb of K12 past flhD as tracks above and as the graph track below. The shaded span is the IS1 element.

In the track menu, Show → Show deletion edges hides the dashed link the graph draws for the other four strains' route, which is labelled with the length of the node it skips.

Comparing odgi viz's node-order axis with K12's coordinates

Open mc/ecoli.viz/chr.full.viz.png, the odgi viz raster --viz wrote, with one row per strain and graph node order on the horizontal axis.

The graph drawn by odgi viz, one row per strain, on graph node order. The gold band marks the locus the next figure opens.

The odgi pav track contains the same paths. Drawn on K12's coordinates, in the raster's row order and colors, they differ only on the horizontal axis. The gold band marks chr:1,000,000-1,100,000 in both.

The same paths and colors on K12's coordinates.

The graph axis counts pangenome bases, so a locus where other strains have sequence K12 lacks takes up more of it. The band contains Sakai's stx2 prophage and a second Sakai-only stretch.

Reproduce it end to end

build_ecoli_pangenome_cactus.sh runs everything above in one shot. It:

downloads the five RefSeq genomes, keeping only each one's chromosome, and runs cactus-pangenome
converts the HAL, VCF, odgi depth and odgi pav into the projections above, and maps the KTa004 reads
indexes the base-level graph by following the reference's paths, as the graph-hosting tutorial describes
writes a config.json with the assemblies, per-strain gene tracks, the projection and segments tracks, the plugin declaration and a default session

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_ecoli_pangenome_cactus.sh
bash build_ecoli_pangenome_cactus.sh   # builds ./ecoli_cactus_build/jbrowse2
npx --yes serve ecoli_cactus_build/jbrowse2

The script needs the tools under Prerequisites; force a container runtime with CONTAINER=singularity.

External links

Minigraph-Cactus

Feedback on this tutorial is welcome: contact us.

Pangenome (Minigraph-Cactus)

Prerequisites#

The data#

Building the graph with cactus-pangenome#

Load the genomes#

Synteny between all the strains (halSynteny)#

Variants against K12 (vg deconstruct)#

Whole-genome alignment onto K12 (MAF)#

Pangenome depth and per-strain presence#

Mapping a new isolate through the graph#

Drawing the graph as a graph#

Indexing the graph#

Comparing odgi viz's node-order axis with K12's coordinates#

Reproduce it end to end#

See also#

External links#

Prerequisites

The data

Building the graph with cactus-pangenome

Load the genomes

Synteny between all the strains (halSynteny)

Variants against K12 (vg deconstruct)

Whole-genome alignment onto K12 (MAF)

Pangenome depth and per-strain presence

Mapping a new isolate through the graph

Drawing the graph as a graph

Indexing the graph

Comparing odgi viz's node-order axis with K12's coordinates

Reproduce it end to end

See also

External links