Helper/wrapper functions for executing coloc & associated analyses.

• Author: David Zhang
• Maintainer: Regina H. Reynolds
• Contributor(s): Alejandro Martinez Carrasco

If you use any of the code from this repository, please cite the original coloc R package and our doi.

The code in this repository is released under an MIT license. This repository is distributed in the hope that it will be useful to the wider community, but without any warranty of any kind. Please see the LICENSE file for more details.

There is no plan to ever submit this code to CRAN or Bioconductor. This code was developed for use by the Ryten Lab and collaborators thereof. While most of the code has been separately tested by David and Regina, it is highly likely bugs exist.

If you would like to install the development version from GitHub you can use the following command:

if(!require("remotes"))
   install.packages("remotes") # if necessary
library(remotes)
install_github("RHReynolds/colochelpR")

The functions herein are simply wrapper functions for executing coloc. Please refer to the coloc vignette for more background on use of coloc. Please remember to cite coloc, if you use this helper package.

For an example of how the wrapper functions in this package have been used in analyses, please refer to the following repository.

For these helper functions to work, it is important that GWASs/QTL datasets have been "tidied" and the necessary columns are present. So far, we have only used coloc for GWAS and eQTL datasets, the necessary columns for which are presented in tables below.

• According to the original coloc paper, it is possible to compute posterior probabilities from single-variant association p.value and maf, but estimated single SNP regression co-efficients (beta) and their variances (varbeta) or standard errors (se) are preferred.
• If you do not have the columns beta and se, but you do have the odds ratio (often OR) and the Z-score, you can generate beta and se via the following formulas:

GWAS %>% 
  dplyr::mutate(beta = log(OR, base = exp(1)),
                se = log(OR, base = exp(1))/Z_SCORE)

• You will also need the number of individuals within your GWAS.
  • This does not need to be a column in your GWAS dataframe. A vector with the number will suffice.
  • If the GWAS is a case-control dataset, you will additionally need the number of cases and controls, in order to calculate the proportion of individuals that are cases (i.e. n_cases/n_total).

• As with the GWAS dataset, if columns beta and se are not available, p.values and maf can be used.
• Many eQTL datasets do not provide MAFs. It is, however, possible to use a reference database to look up MAFs for SNPs within the eQTL datasets. For example, the package MafDb.1Kgenomes.phase3.hs37d5 contains MAFs for a number of populations (remember to select the population that matches your eQTL dataset). See example code below:

mafdb <- MafDb.1Kgenomes.phase3.hs37d5

example_df <- data.frame(rs_id = c("rs12921634", "rs1476958", "rs56189750"))

mafs <- GenomicScores::gscores(x = mafdb, ranges = example_df$rs_id %>% as.character(), pop = "EUR_AF")
mafs <- mafs %>%
  as.data.frame() %>%
  tibble::rownames_to_column(var = "SNP") %>%
  dplyr::rename(maf = EUR_AF) %>% 
  dplyr::select(SNP, maf)

example_df <- example_df %>% 
  inner_join(mafs, by = c("rs_id" = "SNP"))

• The format that SNPs are provided in often varies between location-based (CHR:BP) or RS-ID-based formats.

• As a result, there may sometimes be a need to convert between formats.

• The most common conversion we have run into is conversion of RS IDs to chromosome locations, and have therefore written a function to do this (colochelpR::convert_rs_to_loc()).

  • This functions uses the BSgenome::snpsById() function, and therefore requires a SNPlocs object, which is a container for storing known SNP locations.
  • Use BSgenome::available.SNPs() to return a vector of the available SNPlocs packages.
  • Remember to use a SNPlocs package with the same genome build as the genome build you are intending to match these RS IDs to.

• We have also written a wrapper function for conversion of chromosome locations to RS ids (colochelpR::convert_loc_to_rs()).

  • This functions uses the BSgenome::snpsByOverlaps() function, and therefore requires a SNPlocs object, which is a container for storing known SNP locations (as described above).
  • Users should note that some chromosome-bp locations may have more than one associated RS-ID, thus the user may have to check for duplicates after conversion. We have chosen not to do this, as some users may wish to simply remove duplicated SNPs, or keep one of the duplicates.
  • Users should also note that some chromosome-bp locations may not have an associated RS ID. These will be returned as NA in the RS ID column (SNP) after the conversion.