---
title: "Getting started with gwasrapidd"
output: rmarkdown::html_vignette
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).

There are three ways to access the Catalog database:

- Via the Web interface at
[https://www.ebi.ac.uk/gwas](https://www.ebi.ac.uk/gwas)
- By downloading a
database dump in TSV or RDF/OWL formats:
[https://www.ebi.ac.uk/gwas/docs/file-downloads](https://www.ebi.ac.uk/gwas/docs/file-downloads)
- Or through the GWAS Catalog REST API:
[https://www.ebi.ac.uk/gwas/rest/docs/api](https://www.ebi.ac.uk/gwas/rest/docs/api)

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, and 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 1).


<img src="../man/figures/get_fns.png" title="Figure 1 | gwasrapidd retrieval functions." alt="Figure 1 | gwasrapidd retrieval functions." width="500" style="display: block; margin: auto auto auto 0;" />

You can use a combination of several search criteria with each retrieval function as shown in Figure 2. For example, if you want to get studies using either one of these two criteria: *study accession identifier* and *variant identifier*, you could run the following code:


```r
library(gwasrapidd)
my_studies <- get_studies(study_id = 'GCST000858', variant_id = 'rs12752552')
```

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.

<img src="../man/figures/get_criteria.png" title="Figure 2 | 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)." alt="Figure 2 | 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)." width="500" style="display: block; margin: auto auto auto 0;" />

## Example 1 | 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] 8
my_studies@studies$study_id
#> [1] "GCST003097"   "GCST011008"   "GCST007071"   "GCST009873"   "GCST011005"   "GCST011009"   "GCST90029015" "GCST90029016"
```

Apparently only 8 studies: GCST003097, GCST011008, GCST007071, GCST009873, GCST011005, GCST011009, GCST90029015, GCST90029016. 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."
```

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
