---
title: "Getting started with gwasrapidd"
output: 
  bookdown::html_document2:
    base_format: rmarkdown::html_vignette
    fig_caption: yes
    number_sections: false
# as_is: for figure captions in pkgdown output
pkgdown:
  as_is: true
vignette: >
  %\VignetteIndexEntry{Getting started with gwasrapidd}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: gwasrapidd.bib
---

## The GWAS Catalog

The GWAS Catalog is a service provided by the EMBL-EBI and NHGRI that offers a
manually curated and freely available database of published genome-wide
association studies (GWAS). The Catalog website and infrastructure is hosted by
the [EMBL-EBI](https://www.ebi.ac.uk).

`{gwasrapidd}` facilitates the access to the Catalog via the REST API,
allowing you to programmatically retrieve data directly into R.

## GWAS Catalog Entities

The Catalog REST API is organized around four core entities:

- studies
- associations
- variants
- traits

`{gwasrapidd}` provides four corresponding functions to get each of the
entities: `get_studies()`, `get_associations()`, `get_variants()`, and
`get_traits()`.

Each function maps to an appropriately named S4 classed object: [studies](https://rmagno.eu/gwasrapidd/reference/studies-class.html), [associations](https://rmagno.eu/gwasrapidd/reference/associations-class.html), [variants](https://rmagno.eu/gwasrapidd/reference/variants-class.html), and [traits](https://rmagno.eu/gwasrapidd/reference/traits-class.html) (see Figure \@ref(fig:fns)).


```{r fns, out.width = "70%", echo=FALSE, fig.cap='`{gwasrapidd}` retrieval functions.'}
knitr::include_graphics('../man/figures/get_fns.png') 
```

You can use a combination of several search criteria with each retrieval
function as shown in Figure \@ref(fig:criteria). For example, if you want to get studies using
either one of these two criteria:

- study accession identifier (`study_id`)
- variant identifier (`variant_id`)

You could run the following code:

```r
library(gwasrapidd)

my_studies <- get_studies(study_id = 'GCST000858', variant_id = 'rs12752552')
my_studies@studies
```

This command returns all studies that match either `'GCST000858'` or
`'rs12752552'`. This is equivalent to running `get_studies` separately on each
criteria, and combining the results afterwards:


```r
s1 <- get_studies(study_id = 'GCST000858')
s2 <- get_studies(variant_id = 'rs12752552')
my_studies <- gwasrapidd::union(s1, s2)
```

All four retrieval functions accept the `set_operation` parameter which defines the way the results obtained with each criterion are combined. The two options for this parameter are `'union'` (default) or `'intersection'`, resulting, respectively, in an OR or AND operation.

```{r criteria, out.width = "70%", echo=FALSE, fig.cap='`{gwasrapidd}` arguments for retrieval functions. Colors indicate the criteria that can be used for retrieving GWAS Catalog entities: studies (green), associations (red), variants (purple), and traits (orange).'}
knitr::include_graphics('../man/figures/get_criteria.png') 
```

## Finding Risk Alleles Associated with Autoimmune Disease

As a first example, take the work by @Light2014. In this work the authors focused on variants that had been previously reported in genome-wide association studies (GWAS) for autoimmune disease.

With **gwasrapidd** we can interrogate the GWAS Catalog for the study/studies by searching by *autoimmune disease* (an EFO trait). To do that let's load gwasrapidd first:


```r
library(gwasrapidd)
```

Then query the GWAS Catalog by EFO trait:


```r
my_studies <- get_studies(efo_trait = 'autoimmune disease')
```

We can now check how many GWAS studies we got back:


```r
gwasrapidd::n(my_studies)
#> [1] 9
my_studies@studies$study_id
#> [1] "GCST003097"   "GCST011008"   "GCST007071"   "GCST009873"   "GCST011005"   "GCST011009"   "GCST90029015"
#> [8] "GCST90029016" "GCST90012738"
```

Apparently only 9 studies: GCST003097, GCST011008, GCST007071, GCST009873, GCST011005, GCST011009, GCST90029015, GCST90029016, GCST90012738. Let's see the associated publication titles:


```r
my_studies@publications$title
#> [1] "Meta-analysis of shared genetic architecture across ten pediatric autoimmune diseases."                                                            
#> [2] "Genetic factors underlying the bidirectional relationship between autoimmune and mental disorders - findings from a Danish population-based study."
#> [3] "Leveraging Polygenic Functional Enrichment to Improve GWAS Power."                                                                                 
#> [4] "Meta-analysis of Immunochip data of four autoimmune diseases reveals novel single-disease and cross-phenotype associations."                       
#> [5] "Genetic factors underlying the bidirectional relationship between autoimmune and mental disorders - findings from a Danish population-based study."
#> [6] "Genetic factors underlying the bidirectional relationship between autoimmune and mental disorders - findings from a Danish population-based study."
#> [7] "Mixed-model association for biobank-scale datasets."                                                                                               
#> [8] "Mixed-model association for biobank-scale datasets."                                                                                               
#> [9] "Clinical laboratory test-wide association scan of polygenic scores identifies biomarkers of complex disease."
```

If you want to further inspect these publications, you can quickly browse the respective PubMed entries:

``` r
# This launches your web browser at https://www.ncbi.nlm.nih.gov/pubmed/26301688
open_in_pubmed(my_studies@publications$pubmed_id)
```

Now if we want to know the variants previously associated with autoimmune disease, as used by @Light2014, we need to retrieve statistical association information on these variants, and then filter them based on the same level of significance $P < 1\times 10^{-6}$ [@Light2014].

So let's start by getting the associations by `study_id`:


```r
# You could have also used get_associations(efo_trait = 'autoimmune disease')
my_associations <- get_associations(study_id = my_studies@studies$study_id)
```

Seemingly, there are 182 associations.


```r
gwasrapidd::n(my_associations)
#> [1] 182
```

However, not all variants meet the level of significance, as required by @Light2014:


```r
# Get association ids for which pvalue is less than 1e-6.
dplyr::filter(my_associations@associations, pvalue < 1e-6) %>% # Filter by p-value
  tidyr::drop_na(pvalue) %>%
  dplyr::pull(association_id) -> association_ids # Extract column association_id
```

Here we subset the `my_associations` object by a vector of association identifiers (`association_ids`) into a smaller object, `my_associations2`:


```r
# Extract associations by association id
my_associations2 <- my_associations[association_ids]
gwasrapidd::n(my_associations2)
#> [1] 180
```

Of the 182 associations found in GWAS Catalog, 180 meet the p-value threshold of $1\times 10^{-6}$. Here are the variants, and their respective risk allele and risk frequency:


```r
my_associations2@risk_alleles[c('variant_id', 'risk_allele', 'risk_frequency')]
#> # A tibble: 180 × 3
#>    variant_id risk_allele risk_frequency
#>    <chr>      <chr>                <dbl>
#>  1 rs11580078 G                     0.43
#>  2 rs6679677  A                     0.09
#>  3 rs34884278 C                     0.3 
#>  4 rs6689858  C                     0.29
#>  5 rs2075184  T                     0.23
#>  6 rs36001488 C                     0.48
#>  7 rs4676410  A                     0.19
#>  8 rs4625     G                     0.31
#>  9 rs62324212 A                     0.42
#> 10 rs7725052  C                     0.43
#> # … with 170 more rows
```

## References
