---
title: "Genetic colocalisation"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Genetic colocalisation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval=FALSE
)
```

Here we'll perform colocalisation analysis for a particular region, and plot the regions as well. We'll do the same analysis two ways:

- querying the association data from the [IEU GWAS database](https://gwas.mrcieu.ac.uk/), and 
- downloading the [GWAS VCF](https://github.com/MRCIEU/gwas_vcf_spec) files and querying those.

We'll use the example of LDL cholesterol [ieu-a-300](https://gwas.mrcieu.ac.uk/datasets/ieu-a-300/) and coronary heart disease [ieu-a-7](https://gwas.mrcieu.ac.uk/datasets/ieu-a-7/). 


Load libraries:

```{r}
suppressPackageStartupMessages(suppressWarnings({
	library(gwasglue)
	library(dplyr)
	library(gassocplot)
	library(coloc)
}))
```

## ieugwasr

First find a region that we know to be associated with LDL cholesterol.

```{r}
top <- ieugwasr::tophits('ieu-a-300') %>% arrange(p)
top
```

Choose the best signal and create a range

```{r}
chrpos <- paste0(top$chr[1], ":", top$position[1] - 90000, "-", top$position[1] + 90000)
chrpos
```

Extract, harmonise and format the data

```{r}
out <- ieugwasr_to_coloc(id1='ieu-a-300', id2='ieu-a-7', chrompos=chrpos)
```

Run colocalisation analysis

```{r}
res <- coloc::coloc.abf(out[[1]], out[[2]])
```

Plot

```{r}
temp <- coloc_to_gassocplot(out)
gassocplot::stack_assoc_plot(temp$markers, temp$z, temp$corr, traits=temp$traits)
```


## gwasvcf

Let's do the same with the vcf files (and the indexes). Download from here:

- [https://gwas.mrcieu.ac.uk/files/ieu-a-300/ieu-a-300.vcf.gz](https://gwas.mrcieu.ac.uk/files/ieu-a-300/ieu-a-300.vcf.gz)
- [https://gwas.mrcieu.ac.uk/files/ieu-a-300/ieu-a-300.vcf.gz.tbi](https://gwas.mrcieu.ac.uk/files/ieu-a-300/ieu-a-300.vcf.gz.tbi)
- [https://gwas.mrcieu.ac.uk/files/ieu-a-7/ieu-a-7.vcf.gz](https://gwas.mrcieu.ac.uk/files/ieu-a-7/ieu-a-7.vcf.gz)
- [https://gwas.mrcieu.ac.uk/files/ieu-a-7/ieu-a-7.vcf.gz.tbi](https://gwas.mrcieu.ac.uk/files/ieu-a-7/ieu-a-7.vcf.gz.tbi)


Set a region to plot:

```{r}
chrpos <- "19:11112306-11292306"
```

Extract region from vcf files and convert to coloc object

```{r}
vout <- gwasvcf_to_coloc("ieu-a-300.vcf.gz", "ieu-a-7.vcf.gz", chrpos)
```

Run colocalisation analysis

```{r}
vres <- coloc::coloc.abf(vout[[1]], vout[[2]])
```

Plot

```{r}
library(gassocplot)
temp <- coloc_to_gassocplot(vout)
gassocplot::stack_assoc_plot(temp$markers, temp$z, temp$corr, traits=temp$traits)
```