Once instruments for the exposure trait have been specified, those variants need to be extracted from the outcome trait.

Available studies in IEU GWAS database

The IEU GWAS database (IGD) contains complete GWAS summary statistics from a large number of studies. You can browse them here:

https://gwas.mrcieu.ac.uk/

To obtain details about the available GWASs programmatically do the following:

ao <- available_outcomes()
#> API: public: http://gwas-api.mrcieu.ac.uk/
head(ao)
#> # A tibble: 6 x 22
#>   id      trait     group_name  year consortium author sex   population    unit 
#>   <chr>   <chr>     <chr>      <int> <chr>      <chr>  <chr> <chr>         <chr>
#> 1 ieu-b-~ Number o~ public      2021 MRC-IEU    Clare~ Male~ European      SD   
#> 2 prot-a~ Leptin r~ public      2018 NA         Sun BB Male~ European      NA   
#> 3 prot-a~ Adapter ~ public      2018 NA         Sun BB Male~ European      NA   
#> 4 ukb-e-~ Impedanc~ public      2020 NA         Pan-U~ Male~ Greater Midd~ NA   
#> 5 prot-a~ Dual spe~ public      2018 NA         Sun BB Male~ European      NA   
#> 6 eqtl-a~ ENSG0000~ public      2018 NA         Vosa U Male~ European      NA   
#> # ... with 13 more variables: sample_size <int>, nsnp <int>, build <chr>,
#> #   category <chr>, subcategory <chr>, ontology <chr>, note <chr>, mr <int>,
#> #   pmid <int>, priority <int>, ncase <int>, ncontrol <int>, sd <dbl>

For information about authentication see https://github.com/MRCIEU/ieugwasr/blob/master/README.md#authentication.

The available_outcomes function returns a table of all the available studies in the database. Each study has a unique ID. e.g.

head(subset(ao, select=c(trait, id)))
#> # A tibble: 6 x 2
#>   trait                                                    id                   
#>   <chr>                                                    <chr>                
#> 1 Number of children                                       ieu-b-4760           
#> 2 Leptin receptor                                          prot-a-1725          
#> 3 Adapter molecule crk                                     prot-a-664           
#> 4 Impedance of arm (right)                                 ukb-e-23109_MID      
#> 5 Dual specificity mitogen-activated protein kinase kinas~ prot-a-1841          
#> 6 ENSG00000167699                                          eqtl-a-ENSG000001676~

Extracting particular SNPs from particular studies

If we want to perform MR of BMI against coronary heart disease, we need to identify the SNPs that influence the BMI, and then extract those SNPs from a GWAS on coronary heart disease.

Let’s get the Locke et al 2014 instruments for BMI as an example:

bmi_exp_dat <- extract_instruments(outcomes='ieu-a-2')
head(bmi_exp_dat)
#>   pval.exposure samplesize.exposure chr.exposure se.exposure beta.exposure
#> 1   2.18198e-08              339152            1      0.0030       -0.0168
#> 2   4.56773e-11              339065            1      0.0031        0.0201
#> 3   5.05941e-14              313621            1      0.0087        0.0659
#> 4   5.45205e-10              338768            1      0.0029        0.0181
#> 5   1.88018e-28              338123            1      0.0030        0.0331
#> 6   2.28718e-40              339078            1      0.0037        0.0497
#>   pos.exposure id.exposure        SNP effect_allele.exposure
#> 1     47684677     ieu-a-2   rs977747                      G
#> 2     78048331     ieu-a-2 rs17381664                      C
#> 3    110082886     ieu-a-2  rs7550711                      T
#> 4    201784287     ieu-a-2  rs2820292                      C
#> 5     72837239     ieu-a-2  rs7531118                      C
#> 6    177889480     ieu-a-2   rs543874                      G
#>   other_allele.exposure eaf.exposure                      exposure
#> 1                     T       0.5333 Body mass index || id:ieu-a-2
#> 2                     T       0.4250 Body mass index || id:ieu-a-2
#> 3                     C       0.0339 Body mass index || id:ieu-a-2
#> 4                     A       0.5083 Body mass index || id:ieu-a-2
#> 5                     T       0.6083 Body mass index || id:ieu-a-2
#> 6                     A       0.2667 Body mass index || id:ieu-a-2
#>   mr_keep.exposure pval_origin.exposure data_source.exposure
#> 1             TRUE             reported                  igd
#> 2             TRUE             reported                  igd
#> 3             TRUE             reported                  igd
#> 4             TRUE             reported                  igd
#> 5             TRUE             reported                  igd
#> 6             TRUE             reported                  igd

We now need to find a suitable GWAS for coronary heart disease. We can search the available studies:

ao[grepl("heart disease", ao$trait), ]
#> # A tibble: 28 x 22
#>    id      trait       group_name  year consortium author sex   population unit 
#>    <chr>   <chr>       <chr>      <int> <chr>      <chr>  <chr> <chr>      <chr>
#>  1 ieu-a-8 Coronary h~ public      2011 CARDIoGRAM Schun~ Male~ European   log ~
#>  2 finn-a~ Valvular h~ public      2020 NA         NA     Male~ European   NA   
#>  3 finn-a~ Other or i~ public      2020 NA         NA     Male~ European   NA   
#>  4 ukb-b-~ Diagnoses ~ public      2018 MRC-IEU    Ben E~ Male~ European   SD   
#>  5 ukb-a-~ Diagnoses ~ public      2017 Neale Lab  Neale  Male~ European   SD   
#>  6 ukb-e-~ I25 Chroni~ public      2020 NA         Pan-U~ Male~ South Asi~ NA   
#>  7 ukb-b-~ Diagnoses ~ public      2018 MRC-IEU    Ben E~ Male~ European   SD   
#>  8 finn-a~ Ischemic h~ public      2020 NA         NA     Male~ European   NA   
#>  9 finn-a~ Major coro~ public      2020 NA         NA     Male~ European   NA   
#> 10 finn-a~ Other hear~ public      2020 NA         NA     Male~ European   NA   
#> # ... with 18 more rows, and 13 more variables: sample_size <int>, nsnp <int>,
#> #   build <chr>, category <chr>, subcategory <chr>, ontology <chr>, note <chr>,
#> #   mr <int>, pmid <int>, priority <int>, ncase <int>, ncontrol <int>, sd <dbl>

The most recent CARDIOGRAM GWAS is ID number ieu-a-7. We can extract the BMI SNPs from this GWAS as follows:

chd_out_dat <- extract_outcome_data(
    snps = bmi_exp_dat$SNP,
    outcomes = 'ieu-a-7'
)
#> Extracting data for 79 SNP(s) from 1 GWAS(s)

The extract_outcome_data() function is flexible. The snps argument only requires an array of rsIDs, and the outcomes argument can be a vector of outcomes, e.g.

chd_out_dat  <- extract_outcome_data(c("rs234", "rs17097147"), c('ieu-a-2', 'ieu-a-7'))

will extract the two SNPs from each of the outcomes ieu-a-2 and ieu-a-7.

LD proxies

By default if a particular requested SNP is not present in the outcome GWAS then a SNP (proxy) that is in LD with the requested SNP (target) will be searched for instead. LD proxies are defined using 1000 genomes European sample data. The effect of the proxy SNP on the outcome is returned, along with the proxy SNP, the effect allele of the proxy SNP, and the corresponding allele (in phase) for the target SNP.

The parameters for handling LD proxies are as follows:

  • proxies = TRUE or FALSE (TRUE by default)
  • rsq = numeric value of minimum rsq to find a proxy. Default is 0.8, minimum is 0.6
  • palindromes = Allow palindromic SNPs? Default is 1 (yes)
  • maf_threshold = If palindromes allowed then what is the maximum minor allele frequency of palindromes allowed? Default is 0.3.

Using local GWAS summary data

If you have GWAS summary data that is not present in IEU GWAS database, this can still be used to perform analysis.

Supposing there is a GWAS summary file called “gwas_summary.csv” with e.g. 2 million rows and it looks like this:

rsid,effect,SE,a1,a2,a1_freq,p-value,Units,Gene,n
rs10767664,0.19,0.030612245,A,T,0.78,5.00E-26,kg/m2,BDNF,225238
rs13078807,0.1,0.020408163,G,A,0.2,4.00E-11,kg/m2,CADM2,221431
rs1514175,0.07,0.020408163,A,G,0.43,8.00E-14,kg/m2,TNNI3K,207641
rs1558902,0.39,0.020408163,A,T,0.42,5.00E-120,kg/m2,FTO,222476
...
...

To extract the exposure SNPs from this data, we would use the following command:

outcome_dat <- read_outcome_data(
    snps = bmi_exp_dat$SNP,
    filename = "gwas_summary.csv",
    sep = ",",
    snp_col = "rsid",
    beta_col = "effect",
    se_col = "SE",
    effect_allele_col = "a1",
    other_allele_col = "a2",
    eaf_col = "a1_freq",
    pval_col = "p-value",
    units_col = "Units",
    gene_col = "Gene",
    samplesize_col = "n"
)

This returns an outcome data frame with only the SNPs that were requested (if those SNPs were present in the “gwas_summary.csv” file).

Outcome data format

The extract_outcome_data function returns a table of SNP effects for the requested SNPs on the requested outcomes. The format of the data is similar to the exposure data format, except the main columns are as follows:

  • SNP
  • beta.outcome
  • se.outcome
  • samplesize.outcome
  • ncase.outcome
  • ncontrol.outcome
  • pval.outcome
  • eaf.outcome
  • effect_allele.outcom
  • other_allele.outcome
  • units.outcome
  • outcome
  • consortium.outcome
  • year.outcome
  • pmid.outcome
  • id.outcome
  • originalname.outcome
  • proxy.outcome
  • target_snp.outcome
  • proxy_snp.outcome
  • target_a1.outcome
  • target_a2.outcome
  • proxy_a1.outcome
  • proxy_a2.outcome
  • mr_keep.outcome
  • data_source.outcome

More advanced use of local data

We have developed a summary data format called “GWAS VCF,” which is designed to store GWAS results in a strict and performant way. It is possible to use this format with the TwoSampleMR package. Going down this avenue also allows you to use LD proxy functionality using your own LD reference files (or ones that we provide). For more details, see this package that explains the format and how to query it in R:

https://github.com/mrcieu/gwasvcf

and this package for how to connect the data to other packages including TwoSampleMR

https://github.com/MRCIEU/gwasglue