Zoom browser

Here we demonstrate how to use the zoom browser to interactively browse whole genome GWAS or eQTL data within an R session. Via a shiny/plotly interface, users can click on peaks in a whole genome manhattan plot to zoom in and out of regions and pan left or right using the interactive browser.

Example GWAS data

Example whole genome summary statistics from an SLE GWAS study can be downloaded from EBI here and loaded as follows. The main locuszoomr vignette shows how this data is formatted.

# FTP download full summary statistics (308 MB) from SLE GWAS from
# https://www.ebi.ac.uk/gwas/studies/GCST003156
library(data.table)
SLE_gwas <- fread('../bentham_2015_26502338_sle_efo0002690_1_gwas.sumstats.tsv')

The zoom browser is started using zoom(). Like the locus() function, zoom() tries to autodetect the columns for SNP id, chromosome, position and p-value. If you have complex or ambiguous column names, then the arguments chrom, pos, p and labs can be used to specify the appropriate columns in your dataset.

An appropriate hg19 or hg38 (or other species) ensembl database of genome information must be specified. In the example below, the GWAS data is hg19.

library(locuszoomr)
library(EnsDb.Hsapiens.v75)

zoom(SLE_gwas, ens_db = "EnsDb.Hsapiens.v75")

Navigation

Click on points in the whole genome or individual chromosome manhattan plots to explore a particular locus. Alternatively, enter a gene name or SNP id or chromosome position and location into the input text box. The locus plot can be navigated by zooming in/out, or moving left/right using the navigation icons. Dragging the locus plot will move the locus window. Also, the plotly selection tool can be used to select a specific region.

The icon allows users to toggle options such as showing the recombination rate (only if recomb data has been provided via recomb argument), expanding the gene tracks to show all genes, or selecting individual genes in eQTL data (if the eqtl_gene argument has been set).

Add recombination rate

Recombination rate can be shown on a secondary y axis. Whole recombination rate track data file (around 30 MB) should be downloaded from UCSC genome browser (documented here).

The download site can be accessed at https://hgdownload.soe.ucsc.edu/gbdb/hg38/recombRate/. For hg38, download recomb1000GAvg.bw.

For hg19, the link is https://hgdownload.soe.ucsc.edu/gbdb/hg19/decode/ and the default track we use is hapMapRelease24CombinedRecombMap.bw.

The .bw track file can then be loaded into R as a GRanges object using import.bw(), and then used directly by zoom().

library(rtracklayer)
recomb.hg19 <- import.bw("/../hapMapRelease24CombinedRecombMap.bw")

zoom(SLE_gwas, ens_db = "EnsDb.Hsapiens.v75", recomb = recomb.hg19)

Add LD via LDlink

The LDlink API can be queried from within the zoom browser. This requires an LDlink API access token, which can be obtained from the LDlink website. The token is then supplied via the ld_token argument when calling zoom().

zoom(SLE_gwas, ens_db = "EnsDb.Hsapiens.v75", recomb = recomb.hg19, ld_token = "my_token")

This enables extra buttons for Get LD in the settings dropdown menu. LD API requests are cached and retained during an R session, so subsequent requests for the same SNP are much faster. At each locus, the SNP with the lowest p-value is presumed to be the index SNP for LD requests. But clicking on other points can re-base the LD to alternative SNPs.

Dual GWAS browsing

Using the data2 argument, users can browse 2 GWAS at the same time. Shown below, eQTL data from synovial tissue¹ is compared to data from a rheumatoid arthritis GWAS.² eQTL data is best placed first, so that the eqtl_gene argument can be used to show which gene a SNP drives, as each SNP sometimes drives multiple genes. The beta argument is used to show up/down triangles for positive/ negative beta values. It points to the relevant column of beta values and can be applied to either dataset. But here NA is used to indicate that beta is not being applied to the 2nd GWAS dataset.

The display now has 2 whole genome manhattan plots followed by 2 optional chromosome plots for the current locus, followed by the dual layered locus plot. LD information can also be overlaid and either locus plot can be clicked to rebase the LD information.

# load ensembl Homo sapiens database v110
library(AnnotationHub)
ah <- AnnotationHub()
query(ah, c("EnsDb", "Homo sapiens"))  # shows latest versions
ensDb_v110 <- ah[["AH113665"]]

# load hg38 recombination rate data downloaded from UCSC
recomb1000G <- import.bw("/../recomb1000GAvg.bw")

# double GWAS, eQTL first
zoom(eqtlres, data2 = RA_gwas, ens_db = ensDb_v110,
     eqtl_gene = "gene", beta = c("beta", NA),
     recomb = recomb1000G, ld_token = "my_token")

Saving plots

There are 3 methods for saving plots. The easiest is to click the camera icon in the plotly modebar. This exports the plotly as an svg based on the current screen dimensions.

The 2nd method is to click the save icon. This exports the current plot to a pdf using base graphics (note: the formatting will be different).

In the settings dropdown the export filetype can be changed from pdf to plotly. This exports the current plotly locus object as an .rds file. Once loaded back into R, this can be viewed as a standard plotly by calling the object (which generates the html). This also allows the plotly to be embedded in markdown documents as the plotly object will output html when called. The plotly features which don’t rely on shiny for updates are preserved, e.g. hover information. However shiny dependent features such as gene track reordering based on display width won’t function.

Citation

If you use this package please cite as:

Lewis MJ, Wang S. (2025) locuszoomr: an R package for visualising publication-ready regional gene locus plots. Bioinformatics Advances 2025; vbaf006, doi:10.1093/bioadv/vbaf006

References

Goldmann, K. et al. Ann Rheum Dis 2023; doi:10.1136/ard-2023-224540
Ishigaki, K. et al. Nat Genet 2022; 54: 1640-1651, doi:10.1038/s41588-022-01213-w