In PhanstielLab/BentoBox: Coordinate-Based Genomic Visualization Package for R

knitr::opts_chunk$set(
    fig.align = "center",
    fig.width = 3,
    fig.height = 3,
    collapse = TRUE,
    comment = "#>",
    warning = FALSE,
    message = FALSE
)
library(grid)
library(plotgardener)
library(plotgardenerData)

plotgardener makes it easy to create reproducible, publication-quality figures from multi-omic data. Since each plot can be placed in exactly the desired location, users can stack multiple types of genomic data so that their axes and data are correctly aligned. In this section we will show some examples of plotting multi-omic data and how the pgParams object and "below" y-coordinate can facilitate this process.

In the following example, we plot the same genomic region (i.e. chr21:28000000-30300000) represented in Hi-C data, loop annotations, signal track data, GWAS data, all along a common gene track and genome label axis:

## Load example data
library(plotgardenerData)
data("IMR90_HiC_10kb")
data("IMR90_DNAloops_pairs")
data("IMR90_ChIP_H3K27ac_signal")
data("hg19_insulin_GWAS")

## Create a plotgardener page
pageCreate(
    width = 3, height = 5, default.units = "inches",
    showGuides = FALSE, xgrid = 0, ygrid = 0
)

## Plot Hi-C data in region
plotHicSquare(
    data = IMR90_HiC_10kb,
    chrom = "chr21", chromstart = 28000000, chromend = 30300000,
    assembly = "hg19",
    x = 0.5, y = 0.5, width = 2, height = 2,
    just = c("left", "top"), default.units = "inches"
)

## Plot loop annotations
plotPairsArches(
    data = IMR90_DNAloops_pairs,
    chrom = "chr21", chromstart = 28000000, chromend = 30300000,
    assembly = "hg19",
    x = 0.5, y = 2.5, width = 2, height = 0.25,
    just = c("left", "top"), default.units = "inches",
    fill = "black", linecolor = "black", flip = TRUE
)

## Plot signal track data
plotSignal(
    data = IMR90_ChIP_H3K27ac_signal,
    chrom = "chr21", chromstart = 28000000, chromend = 30300000,
    assembly = "hg19",
    x = 0.5, y = 2.75, width = 2, height = 0.5,
    just = c("left", "top"), default.units = "inches"
)

## Plot GWAS data
plotManhattan(
    data = hg19_insulin_GWAS,
    chrom = "chr21", chromstart = 28000000, chromend = 30300000,
    assembly = "hg19",
    ymax = 1.1, cex = 0.20,
    x = 0.5, y = 3.5, width = 2, height = 0.5,
    just = c("left", "top"), default.units = "inches"
)

## Plot gene track
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(org.Hs.eg.db)
plotGenes(
    chrom = "chr21", chromstart = 28000000, chromend = 30300000,
    assembly = "hg19",
    x = 0.5, y = 4, width = 2, height = 0.5,
    just = c("left", "top"), default.units = "inches"
)

## Plot genome label
plotGenomeLabel(
    chrom = "chr21", chromstart = 28000000, chromend = 30300000,
    assembly = "hg19",
    x = 0.5, y = 4.5, length = 2, scale = "Mb",
    just = c("left", "top"), default.units = "inches"
)

Using the `pgParams` object

The pgParams() function creates a pgParams object that can contain any argument from plotgardener functions.

We can recreate and simplify the multi-omic plot above by saving the genomic region, left-based x-coordinate, and width into a pgParams object:

params <- pgParams(
    chrom = "chr21", chromstart = 28000000, chromend = 30300000,
    assembly = "hg19",
    x = 0.5, just = c("left", "top"),
    width = 2, length = 2, default.units = "inches"
)

Since these values are the same for each of the functions we are using to build our multi-omic figure, we can now pass the pgParams object into our functions so we don't need to write the same parameters over and over again:

## Load example data
data("IMR90_HiC_10kb")
data("IMR90_DNAloops_pairs")
data("IMR90_ChIP_H3K27ac_signal")
data("hg19_insulin_GWAS")

## Create a plotgardener page
pageCreate(
    width = 3, height = 5, default.units = "inches",
    showGuides = FALSE, xgrid = 0, ygrid = 0
)

## Plot Hi-C data in region
plotHicSquare(
    data = IMR90_HiC_10kb,
    params = params,
    y = 0.5, height = 2
)

## Plot loop annotations
plotPairsArches(
    data = IMR90_DNAloops_pairs,
    params = params,
    y = 2.5, height = 0.25,
    fill = "black", linecolor = "black", flip = TRUE
)

## Plot signal track data
plotSignal(
    data = IMR90_ChIP_H3K27ac_signal,
    params = params,
    y = 2.75, height = 0.5
)

## Plot GWAS data
plotManhattan(
    data = hg19_insulin_GWAS,
    params = params,
    ymax = 1.1, cex = 0.20,
    y = 3.5, height = 0.5
)

## Plot gene track
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(org.Hs.eg.db)
plotGenes(
    params = params,
    y = 4, height = 0.5
)

## Plot genome label
plotGenomeLabel(
    params = params,
    y = 4.5, scale = "Mb"
)

The pgParams object also simplifies the code for making complex multi-omic figures when we want to change the genomic region of our plots. If we want to change the region for the figure above, we can simply put it into the pgParams object and re-run the code to generate the figure:

params <- pgParams(
    chrom = "chr21", chromstart = 29000000, chromend = 30000000,
    assembly = "hg19",
    x = 0.5, just = c("left", "top"),
    width = 2, length = 2, default.units = "inches"
)

## Load example data
data("IMR90_HiC_10kb")
data("IMR90_DNAloops_pairs")
data("IMR90_ChIP_H3K27ac_signal")
data("hg19_insulin_GWAS")

## Create a plotgardener page
pageCreate(
    width = 3, height = 5, default.units = "inches",
    showGuides = FALSE, xgrid = 0, ygrid = 0
)

## Plot Hi-C data in region
plotHicSquare(
    data = IMR90_HiC_10kb,
    params = params,
    y = 0.5, height = 2
)

## Plot loop annotations
plotPairsArches(
    data = IMR90_DNAloops_pairs,
    params = params,
    y = 2.5, height = 0.25,
    fill = "black", linecolor = "black", flip = TRUE
)

## Plot signal track data
plotSignal(
    data = IMR90_ChIP_H3K27ac_signal,
    params = params,
    y = 2.75, height = 0.5
)

## Plot GWAS data
plotManhattan(
    data = hg19_insulin_GWAS,
    params = params,
    ymax = 1.1, cex = 0.20,
    y = 3.5, height = 0.5
)

## Plot gene track
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(org.Hs.eg.db)
plotGenes(
    params = params,
    y = 4, height = 0.5
)

## Plot genome label
plotGenomeLabel(
    params = params,
    y = 4.5, scale = "Mb"
)

Alternatively, if we want to plot around a particular gene rather than a genomic region we can use pgParams() to specify gene and geneBuffer. If geneBuffer is not included, the default buffer adds (gene length) / 2 base pairs to the ends of the gene coordinates.

params <- pgParams(
    gene = "LINC00113", geneBuffer = 100000, assembly = "hg19",
    x = 0.5, just = c("left", "top"),
    width = 2, length = 2, default.units = "inches"
)

## Load example data
data("IMR90_HiC_10kb")
data("IMR90_DNAloops_pairs")
data("IMR90_ChIP_H3K27ac_signal")
data("hg19_insulin_GWAS")

## Create a plotgardener page
pageCreate(
    width = 3, height = 5, default.units = "inches",
    showGuides = FALSE, xgrid = 0, ygrid = 0
)

## Plot Hi-C data in region
plotHicSquare(
    data = IMR90_HiC_10kb,
    params = params,
    y = 0.5, height = 2
)

## Plot loop annotations
plotPairsArches(
    data = IMR90_DNAloops_pairs,
    params = params,
    y = 2.5, height = 0.25,
    fill = "black", linecolor = "black", flip = TRUE
)

## Plot signal track data
plotSignal(
    data = IMR90_ChIP_H3K27ac_signal,
    params = params,
    y = 2.75, height = 0.5
)

## Plot GWAS data
plotManhattan(
    data = hg19_insulin_GWAS,
    params = params,
    ymax = 1.1, cex = 0.20,
    y = 3.5, height = 0.5
)

## Plot gene track
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(org.Hs.eg.db)
plotGenes(
    params = params,
    y = 4, height = 0.5
)

## Plot genome label
plotGenomeLabel(
    params = params,
    y = 4.5, scale = "bp", fontsize = 7
)

The "below" y-coordinate

Since multi-omic plots often involve vertical stacking, the placement of multi-omic plots can be facilitated with the "below" y-coordinate. Rather than providing a numeric value or unit object to the y parameter in plotting functions, we can place a plot below the previously drawn plotgardener plot with a character value consisting of the distance below the last plot, in page units, and "b". For example, on a page made in inches, y = "0.1b" will place a plot 0.1 inches below the last drawn plot.

We can further simplify the placement code of our multi-omic figure above by using the "below" y-coordinate to easily stack our plots:

## Load example data
data("IMR90_HiC_10kb")
data("IMR90_DNAloops_pairs")
data("IMR90_ChIP_H3K27ac_signal")
data("hg19_insulin_GWAS")

## pgParams
params <- pgParams(
    chrom = "chr21", chromstart = 28000000, chromend = 30300000,
    assembly = "hg19",
    x = 0.5, just = c("left", "top"),
    width = 2, length = 2, default.units = "inches"
)

## Create a plotgardener page
pageCreate(
    width = 3, height = 5, default.units = "inches",
    showGuides = FALSE, xgrid = 0, ygrid = 0
)

## Plot Hi-C data in region
plotHicSquare(
    data = IMR90_HiC_10kb,
    params = params,
    y = 0.5, height = 2
)

## Plot loop annotations
plotPairsArches(
    data = IMR90_DNAloops_pairs,
    params = params,
    y = "0b",
    height = 0.25,
    fill = "black", linecolor = "black", flip = TRUE
)

## Plot signal track data
plotSignal(
    data = IMR90_ChIP_H3K27ac_signal,
    params = params,
    y = "0b",
    height = 0.5
)

## Plot GWAS data
plotManhattan(
    data = hg19_insulin_GWAS,
    params = params,
    ymax = 1.1, cex = 0.20,
    y = "0.25b",
    height = 0.5
)

## Plot gene track
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(org.Hs.eg.db)
plotGenes(
    params = params,
    y = "0b",
    height = 0.5
)

## Plot genome label
plotGenomeLabel(
    params = params,
    y = "0b",
    scale = "Mb"
)

Plotting and comparing multiple signal tracks

In many multi-omic visualizations, multiple signal tracks are often aligned and stacked to compare different kinds of signal data and/or signals from different samples. plotgardener does not normalize signal data based on variables like read depth, but it is possible to scale plotgardener signal plots to the same y-axis.

To determine the appropriate y-axis range, we first must get the maximum signal score from all of our datasets to be compared:

library(plotgardenerData)
data("IMR90_ChIP_H3K27ac_signal")
data("GM12878_ChIP_H3K27ac_signal")

maxScore <- max(c(IMR90_ChIP_H3K27ac_signal$score, 
                    GM12878_ChIP_H3K27ac_signal$score))
print(maxScore)

In each of our signal plotting calls, we will then use the range parameter to set the range of both our y-axes to c(0, maxScore). Here we can do this with our pgParams object:

params <- pgParams(
    chrom = "chr21",
    chromstart = 28000000, chromend = 30300000,
    assembly = "hg19",
    range = c(0, maxScore)
)

We are now ready to plot, align, and compare our signal plots along the genomic x-axis and the score y-axis:

## Create a page
pageCreate(width = 7.5, height = 2.1, default.units = "inches",
            showGuides = FALSE, xgrid = 0, ygrid = 0)

## Plot and place signal plots
signal1 <- plotSignal(
    data = IMR90_ChIP_H3K27ac_signal, params = params,
    x = 0.5, y = 0.25, width = 6.5, height = 0.65,
    just = c("left", "top"), default.units = "inches"
)

signal2 <- plotSignal(
    data = GM12878_ChIP_H3K27ac_signal, params = params,
    linecolor = "#7ecdbb",
    x = 0.5, y = 1, width = 6.5, height = 0.65,
    just = c("left", "top"), default.units = "inches"
)

## Plot genome label
plotGenomeLabel(
    chrom = "chr21",
    chromstart = 28000000, chromend = 30300000,
    assembly = "hg19",
    x = 0.5, y = 1.68, length = 6.5,
    default.units = "inches"
)

## Add text labels
plotText(
    label = "IMR90", fonsize = 10, fontcolor = "#37a7db",
    x = 0.5, y = 0.25, just = c("left", "top"),
    default.units = "inches"
)
plotText(
    label = "GM12878", fonsize = 10, fontcolor = "#7ecdbb",
    x = 0.5, y = 1, just = c("left", "top"),
    default.units = "inches"
)