R/link_eqtl.R

In locuszoomr: Gene Locus Plot with Gene Annotations

#' Obtain GTEx eQTL data via LDlinkR
#' 
#' Adds eQTL (expression quantitative trait loci) information from GTEx
#' (https://gtexportal.org/) to a 'locus' class object. It queries LDlink
#' (https://ldlink.nci.nih.gov/) via the `LDlinkR` package to retrieve GTEx eQTL
#' information on a reference SNP.
#' 
#' @param loc Object of class 'locus' generated by [locus()]
#' @param pop A 1000 Genomes Project population, (e.g. YRI or CEU), multiple 
#' allowed, default = "CEU". Passed to `LDlinkR::LDexpress()`.
#' @param r2d Either "r2" for LD r^2 or "d" for LD D', default = "r2". Passed 
#' to `LDlinkR::LDexpress()`.
#' @param token Personal access token for accessing 1000 Genomes LD data via 
#' LDlink API. See `LDlinkR` package documentation.
#' @param ... Optional arguments such as `genome_build` which are passed on to
#'   `LDlinkR::LDexpress()`
#' @return Returns an object of class 'locus' with an extra list element 'LDexp'
#'   containing a dataframe of information obtained via `LDexpress()`.
#' @details
#' The additional eQTL information obtained from LDlink web server can be
#' displayed using [eqtl_plot()] which generates a scatter plot with gene tracks
#' similar to a locus plot, or with [overlay_plot()] which tries to overlay the
#' EQTL analysis over the original locus results (e.g. GWAS).
#' 
#' @seealso [locus()] [eqtl_plot()] [overlay_plot()]
#' @export

link_eqtl <- function(loc,
                      pop = "CEU",
                      r2d = "r2",
                      token = "", ...) {
  if (!inherits(loc, "locus")) stop("Not a locus object")
  if (!requireNamespace("LDlinkR", quietly = TRUE)) {
    stop("Package 'LDlinkR' must be installed to use this feature",
         call. = FALSE)
  }
  labs <- loc$labs
  index_snp <- loc$index_snp
  
  if (token == "") stop("token is missing")
  LDexp <- mem_LDexpress(snps = index_snp, pop = pop, r2d = r2d, 
                         token = token, ...)
  for (i in c("R2", "D'", "Effect_Size", "P_value")) {
    LDexp[, i] <- as.numeric(LDexp[, i])
  }
  LDexp$Effect_Allele <- gsub("=.*", "", LDexp$Effect_Allele_Freq)
  loc$LDexp <- LDexp
  
  loc
}

mem_LDexpress <- memoise(LDlinkR::LDexpress)