--- title: "Outcome data" author: "Gibran Hemani" output: rmarkdown::html_vignette: toc: true toc_depth: 2 bibliography: refs.bib vignette: > %\VignetteIndexEntry{Outcome data} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", cache = TRUE ) ``` ```{r message=FALSE} library(TwoSampleMR) ``` ```{r, eval=FALSE, echo=FALSE} # Get data for the vignette ao <- available_outcomes() bmi_exp_dat <- extract_instruments(outcomes = 'ieu-a-2') chd_out_dat1 <- extract_outcome_data( snps = bmi_exp_dat$SNP, outcomes = 'ieu-a-7' ) chd_out_dat2 <- extract_outcome_data( snps = c("rs234", "rs17097147"), outcomes = c('ieu-a-2', 'ieu-a-7') ) save(ao, bmi_exp_dat, chd_out_dat1, chd_out_dat2, file = file.path("inst", "extdata", "vig_outcome.RData"), compress = "xz") ``` ```{r, echo=FALSE, eval=TRUE, warning=FALSE} load(system.file("extdata", "vig_outcome.RData", package = "TwoSampleMR")) ``` 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: ```{r, eval=FALSE} ao <- available_outcomes() ``` ```{r} head(ao) ``` For information about authentication see . The `available_outcomes` function returns a table of all the available studies in the database. Each study has a unique ID. e.g. ```{r} head(subset(ao, select = c(trait, id))) ``` ## 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: ```{r, eval=FALSE} bmi_exp_dat <- extract_instruments(outcomes = 'ieu-a-2') ``` ```{r} head(bmi_exp_dat) ``` We now need to find a suitable GWAS for coronary heart disease. We can search the available studies: ```{r} ao[grepl("heart disease", ao$trait), ] ``` The most recent CARDIOGRAM GWAS is ID number `ieu-a-7`. We can extract the BMI SNPs from this GWAS as follows: ```{r, eval=FALSE} chd_out_dat1 <- extract_outcome_data( snps = bmi_exp_dat$SNP, outcomes = 'ieu-a-7' ) ``` 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. ```{r eval=FALSE} chd_out_dat2 <- extract_outcome_data( snps = c("rs234", "rs17097147"), outcomes = 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: ```{r eval=FALSE} 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