ggtubemap

Pangenome graphs as tube maps, the view of sequenceTubeMap, drawn with ggplot2. The layout is sequenceTubeMap's own, from @jbrowse/tubemap-core, run in V8; the drawing is ordinary ggplot2 layers, so a tube map takes scales, themes, patchwork and vector output like any other plot.

• In: a GFA, sequenceTubeMap's extracted JSON, or a region of a gfa-to-tabix walk file such as the HPRC release 2.1 graph's, read by range request
• Out: every tube, node outline, base and read mismatch as data frames, which autoplot() draws as sequenceTubeMap does

The package is not on CRAN yet, so install it from GitHub:

pak::pak("GMOD/ggtubemap")

The layout bundle is tracked in git and the V8 package runs it, so there is no Node or browser to install.

Gallery

Nine assemblies, bases printed	A tandem repeat from the HPRC walk file
HLA-A in fourteen haplotypes	Sixteen haplotypes by phase
91 haplotypes through MICB	Reads under BRCA1
43 mitochondria, a panel per sample	With ggbandage's force layout

tools/readme-figures.R draws every figure.

Usage

library(ggtubemap)
library(ggplot2)

tm <- tubemap_layout(tubemap_example("example-6"))
autoplot(tm)

tubemap_layout() takes a graph, from read_gfa() or as a list of nodes and tracks, and returns every tube, node outline, base and read mismatch as data frames in layout pixels: shapes (the tubes), nodes, outlines, bases, mismatches, panels and ruler. autoplot() draws them in sequenceTubeMap's order.

tubemap_example() loads one of four bundled graphs: example-6, brca1, micb-kir3dl1 and hprc-chrM.

g <- read_gfa("graph.gfa", reference = "GRCh38#0#chr6")
tm <- tubemap_layout(g, node_width = "compressed")

read_gfa() reads P and W lines; walks get PanSN names, sample#haplotype#contig.

From a walk file

A whole-genome pangenome is too large to load, so gfa-to-tabix indexes it once and read_gfa() cuts one region at a time. A walk file is a .walks.gz with its .walks.gz.idx directory beside it, local or a URL. It holds every haplotype's walk, named in PanSN form, sample#haplotype#contig:

vg convert -f graph.gbz > graph.gfa
gfa-to-tabix graph.gfa --walks --refs GRCh38 -o graph

That writes graph.GRCh38.walks.gz, in walk format 7. ggtubemap reads format 7 only, so rebuild an older .walks.bed.gz with gfa-to-tabix 4.0 or later. JBrowse's Hosting your own graph covers the steps, and ggbandage's Region-specific pangenomes the other kind of index. The JBrowse HPRC demo hosts the release 2.1 graph's chr6, built with --sequences so its nodes print their bases:

walks <- "https://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38-chr6-sequences.GRCh38.walks.gz"
str <- read_gfa(
  walks,
  region = "chr6:160,653,940-160,654,130",
  haplotypes = c("CHM13#0", "HG00097#1", "HG00099#1", "HG00133#2", "HG00140#1")
)
autoplot(tubemap_layout(str), xlim = c(700, 2000))

• Range requests fetch only the blocks the region needs: this cut takes under a second
• The reference sample's walk comes first; haplotypes picks the others by sample#haplotype, every one by default
• chr6 names the reference sample's sequence on it, GRCh38#0#chr6
• @jbrowse/graph-walks does the cut in V8, and R serves each byte range
• At normal width the repeat is 6,200 layout pixels across, so xlim shows part of it; node_width = "compressed" fits it whole, without bases
• The whole genome is in https://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38.GRCh38.walks.gz, 1.6 GB, with the same nodes and walks but no bases

Twelve haplotypes through 400 bp of HLA-A, a bubble at nearly every node:

samples <- c("HG00097", "HG00099", "HG00126", "HG00128", "HG00133", "HG002")
hla <- read_gfa(
  walks,
  region = "chr6:29,942,850-29,943,250",
  haplotypes = c("CHM13#0", paste0(rep(samples, each = 2), "#", 1:2))
)
autoplot(tubemap_layout(hla, node_width = "compressed"))

Colour by anything

Each shape row has its track's name, strand, mapping_quality, sample_name, read_group and freq. Add a column and map it:

micb <- tubemap_layout(tubemap_example("micb-kir3dl1"), node_width = "compressed")
micb$shapes$phase <- sub("^[^#]+#([^#]+)#.*", "\\1", micb$shapes$name)
autoplot(micb, haplotypes = aes(fill = phase)) +
  scale_fill_manual(values = c("0" = "black", "1" = "#d97706", "2" = "#0e7490"))

The same works on a cut from a walk file, here eight samples' sixteen haplotypes through the tandem repeat:

samples <- c(samples, "HG00140", "HG00146")
phased <- tubemap_layout(
  read_gfa(walks, region = "chr6:160,653,940-160,654,130",
           haplotypes = paste0(rep(samples, each = 2), "#", 1:2)),
  node_width = "compressed"
)
phased$shapes$phase <- sub("^[^#]+#([^#]+)#.*", "\\1", phased$shapes$name)
autoplot(phased, haplotypes = aes(fill = phase)) +
  scale_fill_manual(values = c("0" = "black", "1" = "#d97706", "2" = "#0e7490"))

Pass mappings through autoplot()'s haplotypes and reads arguments rather than adding another tube layer, which would draw over the nodes. A fill on the reads alone, such as reads = aes(fill = mapping_quality), takes the fill scale, and the haplotypes keep their default colours.

Reads

Reads come with the graph in sequenceTubeMap's extracted form, as in the brca1 and hprc-chrM examples. autoplot() colours them by strand, fades them by mapping quality, and marks their substitutions, deletions and insertions. xlim and ylim zoom in by layout pixel:

br <- tubemap_layout(tubemap_example("brca1"))
autoplot(br, xlim = c(250, 900))

Ruler and positions

autoplot() draws a ruler along the top in the coordinates of the first track with an index_of_first_base. read_gfa() gives every track one, taken from a walk's start or from a path named name[start-end] or name:start-end, and 0 for other paths. Draw the ruler on a custom plot with scale_x_tubemap(tm).

tubemap_x() turns positions on that track into layout x, to mark a region:

micb <- tubemap_layout(tubemap_example("micb-kir3dl1"), node_width = "compressed")
autoplot(micb) +
  annotate("rect", xmin = tubemap_x(micb, 31500700), xmax = tubemap_x(micb, 31500800),
           ymin = -Inf, ymax = Inf, fill = "gold", alpha = 0.3)

Facets

facet_by stacks one panel per sample_name, read_group or PanSN haplotype_sample, every panel sharing the same x:

chrm <- tubemap_layout(tubemap_example("hprc-chrM"),
                       node_width = "compressed", facet_by = "sample_name")
autoplot(chrm)

With ggbandage

ggbandage cuts the same walk file for its Bandage-style layouts, and patchwork stacks the two. This reuses walks and str from From a walk file:

library(patchwork)
cut <- tempfile(fileext = ".gfa")
writeLines(
  ggbandage::bandage_gfa(walks, "chr6:160,653,940-160,654,130",
                         haplotypes = c("CHM13#0", "HG00097#1", "HG00099#1", "HG00133#2", "HG00140#1"),
                         context = 0),
  cut
)
g <- ggbandage::bandage_layout(cut, walks = "HG00133#2#CM090065.1", width = 600)
autoplot(tubemap_layout(str, node_width = "compressed")) / autoplot(g, labels = TRUE)

Custom plots

autoplot() builds the plot from these layers, each taking the layout:

ggplot() +
  geom_tubemap_haplotypes(tm) +
  geom_tubemap_reads(tm) +
  geom_tubemap_nodes(tm) +
  geom_tubemap_sequence(tm) +
  geom_tubemap_mismatches(tm) +
  coord_tubemap(tm) +
  theme_tubemap()

coord_tubemap() reverses the y scale and fixes the aspect ratio at 1, because the layout's y runs downward and its curves assume square pixels.

Output size

Text and node outlines are sized in layout pixels and scale with the plot, so a map looks the same at any size. The text layers drop sequence that would print smaller than 1.5 pt, as the viewer hides it when zoomed out.

tubemap_size() gives the inches at which the map prints at a chosen pixels per inch; at the default 96 it matches sequenceTubeMap at 100%:

size <- tubemap_size(tm, legend = 0.5)
ggsave("tubemap.pdf", autoplot(tm), width = size[["width"]],
       height = size[["height"]], device = cairo_pdf)

Development

The layout bundle, the files shared with ggbandage and the update workflows are described in DEVELOPMENT.md.

See also

• @jbrowse/tubemap-core - sequenceTubeMap's layout, without the DOM
• ggbandage - Bandage-style graph figures as ggplot2 layers
• jbrowse-plugin-graphgenomeviewer - JBrowse 2 plugin that browses these graphs by locus
• sequenceTubeMap, MemPanG26 edition - our fork of the tube map app

Tutorials on jbrowse.org

• HPRC part 1: graph alleles and haplotypes