In RajLabMSSM/catalogueR: QTL data extraction and colocalization

library(catalogueR)

Introduction

Here, we aim to more robustly test whether the genetic signals underlying two dataset (e.g. GWAS vs. eQTL in the same locus)are indeed the same, using a methodology called colocalization.
Specifically, we will use coloc, which infers the probability that each SNP is causal in a given locus in each of the datasets. It then tests the hypothesis that those signals show substantially similar association distributions.

List datasets

We have previously queried eQTL Catalogue using several Parkinson's disease GWAS loci (with eQTLcatalogue_query()).
Let's first gather those files saved as csv files.

gwas.qtl_paths <- catalogueR::eQTLcatalogue_example_queries(fnames = c(
  "BST1__Alasoo_2018.macrophage_IFNg+Salmonella.tsv.gz",
  "BST1__Alasoo_2018.macrophage_naive.tsv.gz",
  "BST1__Alasoo_2018.macrophage_Salmonella.tsv.gz"
))

Run coloc

Since its original release, coloc has been updated so that it can now model multiple causal variants within a given dataset (see the paper), whereas previously it could assume one causal variants. Thus, it may be able to better estimate the colocalization probability between two datasets.

coloc_QTLs <- catalogueR::COLOC_run(gwas.qtl_paths = gwas.qtl_paths,
                                    top_snp_only = TRUE,
                                    split_by_group = FALSE,
                                    method = "abf")

Plot results

First, let's plot only the results with >80% colocalization probability. This is a colocalization threshold commonly used in the field.

coloc_plot <- catalogueR::COLOC_heatmap(coloc_QTLs = coloc_QTLs, 
                                             coloc_thresh = .8)

Uh oh! Looks like there were no colocalizations at >=80% probability. Let's explore the data a bit and lower the threshold to 0 (not advisable except for demonstration purposes).

coloc_plot <- catalogueR::COLOC_heatmap(coloc_QTLs = coloc_QTLs,
                                             coloc_thresh = 0)

Too many results:

In other cases, you may have too many colocalized results to show in a plot all at once.
As a general rule, the more tests you do to higher your probability threshold should be, since Bayesian methods don't calculate P-values (and this cannot used multiple-testing correction methods like FDR).
For example, if you ran many thousands of coloc test across many GWAS-eQTL pairs, it's advisable to raise your PP threshold to .90 or even .99. This also has the benefit of narrowing down your (sometimes many!) results so that you may focus on only the strongest colocalizations.