In anthony-aylward/refflat: Download and plot refFlat data

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Plotting genes with `refflat`

By Anthony Aylward, 2019-03-11

This vignette will cover a simple example use case for refflat: plotting gene locations beneath GWAS association values, to provide context for the locus.

Import the libraries we'll need.

library(RColorBrewer)
library(refflat)

Get a couple of good colors using RColorBrewer.

palette <- brewer.pal(3, "Paired")[1:2]

refflat includes an example data frame called t2d_fto, which contains P-values of GWAS association (type 2 diabetes risk) for SNPs in a region of roughly 1 Mb around the FTO gene. See help(t2d_fto) for details.

Now we can create the plot. We'll use par to divide the plot into two sections, invoke a basic plot call to plot the GWAS association data, then invoke plot_refflat to plot the gene locations.

par(mar = c(2, 4, 0, 0), mfcol = c(2, 1))
plot(
  t2d_fto[["POS"]],
  -1 * log10(t2d_fto[["PVALUE"]]),
  pch = 21,
  col = palette[2],
  bg = palette[1],
  xlab = "",
  xaxt = "n",
  ylab = "-log10 p-value"
)
plot_refflat("chr16", t2d_fto[1, "POS"], t2d_fto[nrow(t2d_fto), "POS"])

Note that the transcription start sites of RPGRIP1L and FTO are very close to one another. We can clarify things by imposing a 1 kb buffer around genes, requiring that they be drawn at different y-axis levels if their ends are less than 1 kb apart.

par(mar = c(2,4,0,0), mfcol = c(2, 1))
plot(
  t2d_fto[["POS"]],
  -1 * log10(t2d_fto[["PVALUE"]]),
  pch = 21,
  col = palette[2],
  bg = palette[1],
  xlab = "",
  xaxt = "n",
  ylab = "-log10 p-value"
)
plot_refflat(
  "chr16",
  t2d_fto[1, "POS"],
  t2d_fto[nrow(t2d_fto), "POS"],
  buffer = 1e3
)

Finally, we can do some more drawing on the lower panel. For example, to draw a vertical line at the position of the lead SNP:

lead_snp_position = t2d_fto[
  t2d_fto[["PVALUE"]] == min(t2d_fto[["PVALUE"]]),
  "POS"
]
par(mar = c(2,4,0,0), mfcol = c(2, 1))
plot(
  t2d_fto[["POS"]],
  -1 * log10(t2d_fto[["PVALUE"]]),
  pch = 21,
  col = palette[2],
  bg = palette[1],
  xlab = "",
  xaxt = "n",
  ylab = "-log10 p-value"
)
plot_refflat(
  "chr16",
  t2d_fto[1, "POS"],
  t2d_fto[nrow(t2d_fto), "POS"],
  buffer = 1e3
)
abline(v = lead_snp_position, lty = 2, lwd = 2)