Perform fast queries against a massive database of complete GWAS summary data

Authentication

From 1st May 2024, most queries to the OpenGWAS API will require user authentication. For more information on why this is necessary, see this blog post.

To authenticate, you need to generate a token from the OpenGWAS website. The token behaves like a password, and it will be used to authorise the requests you make to the OpenGWAS API. Here are the steps to generate the token and then have ieugwasr automatically use it for your queries:

1. Login to https://api.opengwas.io/profile/
2. Generate a new token
3. Add OPENGWAS_JWT=<token> to your .Renviron file. This file could be either in your home directory or in the working directory of your R session. You can check the location of your .Renviron file by running Sys.getenv("R_ENVIRON_USER") in R.
4. Restart your R session
5. To check that your token is being recognised, run ieugwasr::get_opengwas_jwt(). If it returns a long random string then you are authenticated.
6. To check that your token is working, run user(). It will make a request to the API for your user information using your token. It should return a list with your user information. If it returns an error, then your token is not working.

Now any query to OpenGWAS will automatically include your token to authorise the request.

IMPORTANT NOTE: Do not share this token with others as it is equivalent to a password. If you think your token has been compromised, you can generate a new one from the OpenGWAS website.

General API queries

The API has a number of endpoints documented here. A general way to access them in R is using the api_query function. There are two types of endpoints - GET and POST.

• GET - you provide a single URL which includes the endpoint and query. For example, for the association endpoint you can obtain some rsids in some studies, e.g.

  api_query("associations/ieu-a-2,ieu-a-7/rs234,rs123")

• POST - Here you send a “payload” to the endpoint. So, the path specifies the endpoint and you add a list of query specifications. This is useful for long lists of rsids being queried, for example

  api_query("associations", query = list(rsid = c("rs234", "rs123"), id = c("ieu-a-2", "ieu-a-7")))

The api_query function returns a response object from the httr package. See below for a list of functions that make the input and output to api_query more convenient.

Get API status

library(ieugwasr)
api_status()

Get list of all available studies

gwasinfo()

Get list of a specific study

gwasinfo("ieu-a-2")

Extract particular associations from particular studies

Provide a list of variants to be obtained from a list of studies:

associations(variants = c("rs123", "7:105561135"), id = c("ieu-a-2", "ieu-a-7"))

By default this will look for LD proxies using 1000 genomes reference data (Europeans only, the reference panel has INDELs removed and only retains SNPs with MAF > 0.01). This behaviour can be turned off using proxies=0 as an argument.

Note that the queries are performed on rsids, but chromosome:position values will be automatically converted. A range query can be done using e.g.

associations(variants = "7:105561135-105563135", id = c("ieu-a-2"), proxies = 0)

Get the tophits from a study

The tophits can be obtained using

tophits(id = "ieu-a-2")

Note that it will perform strict clumping by default (r2 = 0.001 and radius = 10000kb). This can be turned off with clump=0.

Perform PheWAS

Lookup association of specified variants across every study, returning at a particular threshold. Note that no LD proxy lookups are made here.

phewas(variants = "rs1205", pval = 1e-5)

PheWAS can also be performed in only specific subsets of the data. The datasets in the OpenGWAS database are organised by batch, you can see info about it here: https://gwas.mrcieu.ac.uk/datasets/ or get a list of batches and their descriptions using:

batches()

You can perform PheWAS in only specified batches using:

phewas(variants = "rs1205", pval = 1e-5, batch = c('ieu-a', 'ukb-b'))

By default PheWAS is performed in all batches (which is of course somewhat slower).

LD clumping

The API has a wrapper around plink version 1.90 and can use it to perform clumping with an LD reference panel from 1000 genomes reference data.

a <- tophits(id = "ieu-a-2", clump = 0)
b <- ld_clump(
    dplyr::tibble(rsid = a$rsid, pval = a$p, id = a$id)
)

There are 5 super-populations that can be requested via the pop argument. By default this will use the Europeans subset (EUR super-population). The reference panel has INDELs removed and only retains SNPs with MAF > 0.01 in the selected population.

Note that you can perform the same operation locally if you provide a path to plink and a bed/bim/fam LD reference dataset. e.g.

ld_clump(
    dplyr::tibble(rsid = a$rsid, pval = a$p, id = a$id),
    plink_bin = "/path/to/plink",
    bfile = "/path/to/reference_data"
)

See the following vignette for more information: Running local LD operations

LD matrix

Similarly, a matrix of LD r values can be generated using

ld_matrix(b$variant)

This uses the API by default but is limited to only 500 variants. You can use, instead, local plink and LD reference data in the same manner as in the ld_clump function, e.g.

ld_matrix(b$variant, plink_bin = "/path/to/plink", bfile = "/path/to/reference_data")

There are 5 super-populations that can be requested via the pop argument. By default this will use the Europeans subset (EUR super-population). The reference panel has INDELs removed and only retains SNPs with MAF > 0.01 in the selected population.

Super-populations:

• EUR = European
• EAS = East Asian
• AMR = Admixed American
• SAS = South Asian
• AFR = African

See the following vignette for more information: Running local LD operations

Variant information

Translating between rsids and chromosome:position, while also getting other information, can be achieved.

The chrpos argument can accept the following

• <chr>:<position>
• <chr>:<start>-<end>

For example

a <- variants_chrpos(c("7:105561135-105563135", "10:44865737"))

This provides a table with dbSNP variant IDs, gene info, and various other metadata. Similar data can be obtained from searching by rsid

b <- variants_rsid(c("rs234", "rs333"))

And a list of variants within a particular gene region can also be found. Provide a ensembl or entrez gene ID (e.g. ENSG00000123374 or 1017) to the following:

c <- variants_gene("ENSG00000123374")

Extracting GWAS summary data based on gene region

Here is an example of how to obtain summary data for some datasets for a gene region. As an example, we’ll extract CDK2 (HGNC number 1017) from a BMI dataset (ieu-a-2)

Use the mygene bioconductor package to query the mygene.info API.

library(mygene)
a <- mygene::getGene("1017", fields = "genomic_pos_hg19")
r <- paste0(a[[1]]$genomic_pos_hg19$chr, ":", a[[1]]$genomic_pos_hg19$start, "-", a[[1]]$genomic_pos_hg19$end)
b <- ieugwasr::associations(r, "ieu-a-2")

1000 genomes annotations

The OpenGWAS database contains a database of population annotations from the 1000 genomes project - the alternative allele frequencies and the LD scores for each variant, calculated for each super population separately. Only variants are present if they are MAF > 1% in at least one super population. You can access this info in different ways

1. Look up a particular set of rsids

   afl2_rsid(c("rs234", "rs123"))

2. Look up a set of positions or regions

   afl2_chrpos("1:100000-900000")

3. Extract annotations for a list of 20k variants that are common in all super populations, and evenly spaced across the genome

   afl2_list()

4. Extract annotations for a 1.3 million HapMap3 variants

   afl2_list("hapmap3")

5. Infer the ancestry of a particular study by comparing the allele frequencies with different super population reference frequencies

   snplist <- afl2_list()
   eur_example <- associations(snplist$rsid, "ieu-a-2")
   infer_ancestry(eur_example, snplist)
   eas_example <- associations(snplist$rsid, "bbj-a-10")
   infer_ancestry(eur_example, snplist)