library(ieugwasr)
library(tidyverse)
library(TwoSampleMR)
library(rsnps)

Summary

Analysis

Depression GWAS:

gwasinfo("ieu-b-102") %>% str()
Classes ‘GwasInfo’ and 'data.frame':    1 obs. of  19 variables:
 $ id         : chr "ieu-b-102"
 $ trait      : chr "Major depression"
 $ year       : int 2019
 $ category   : chr "Disease"
 $ ncase      : int 170756
 $ priority   : int 0
 $ note       : chr "UKBiobank and PGC (excluding 23andme)"
 $ sex        : chr "Males and Females"
 $ consortium : chr "PGC"
 $ mr         : int 1
 $ build      : chr "HG19/GRCh37"
 $ subcategory: logi NA
 $ population : chr "European"
 $ author     : chr "Howard DM"
 $ ontology   : chr "EFO:0003761;MONDO:0002009"
 $ group_name : chr "public"
 $ unit       : logi NA
 $ ncontrol   : int 329443
 $ pmid       : int 30718901

C-reactive protein GWAS:

gwasinfo("ieu-b-35") %>% str()
Classes ‘GwasInfo’ and 'data.frame':    1 obs. of  19 variables:
 $ id         : chr "ieu-b-35"
 $ trait      : chr "C-Reactive protein level"
 $ year       : int 2018
 $ category   : chr "Continuous"
 $ priority   : int 0
 $ note       : logi NA
 $ sex        : chr "Males and Females"
 $ consortium : logi NA
 $ mr         : int 1
 $ build      : chr "HG19/GRCh37"
 $ subcategory: chr "Immune system"
 $ population : chr "European"
 $ author     : chr "Ligthart, S"
 $ ontology   : logi NA
 $ group_name : chr "public"
 $ unit       : logi NA
 $ sample_size: int 204402
 $ pmid       : int 30388399
 $ nsnp       : int 2414379

Perform MR

dat <- make_dat("ieu-b-35", "ieu-b-102")
Extracting data for 57 SNP(s) from 1 GWAS(s)
Finding proxies for 2 SNPs in outcome ieu-b-102
Extracting data for 2 SNP(s) from 1 GWAS(s)
Harmonising C-Reactive protein level || id:ieu-b-35 (ieu-b-35) and Major depression || id:ieu-b-102 (ieu-b-102)
res %>% select(method, nsnp, b, se, pval)

Use only CRP variant

ress <- mr_singlesnp(dat)
mr_forest_plot(ress)
$`ieu-b-35.ieu-b-102`

attr(,"split_type")
[1] "data.frame"
attr(,"split_labels")

Find CRP SNP

crp_snp <- subset(dat, chr == 1) %>%
  mutate(dist=abs(pos.exposure - 159682079)) %>%
  arrange(dist) %>%
  select(dist, SNP, pval.exposure) %>%
  {.$SNP[1]}

Find IL6R SNP

dbsnp_info <- ncbi_snp_query(dat$SNP)
Getting info about the following rsIDs: rs10240168, rs10512597, rs1051338, rs10521222, rs10778215, rs10832027, rs10925027, rs11108056, rs12202641, rs12587622, rs1260326, rs12960928, rs12995480, rs13233571, rs13409371, rs1441169, rs1490384, rs1509394, rs1558902, rs1582763, rs17658229, rs178810, rs1800961, rs1805096, rs1880241, rs2064009, rs2239222, rs2293476, rs2315008, rs2352975, rs2710804, rs2794520, rs2836878, rs2852151, rs2891677, rs3122633, rs340005, rs387976, rs4092465, rs4129267, rs4246598, rs4420638, rs4655802, rs4656849, rs469772, rs4767920, rs4841132, rs6001193, rs644234, rs6485751, rs6601302, rs6672627, rs7121935, rs7310409, rs9271608, rs9284725

MR using just CRP

subset(ress, SNP == crp_snp) %>% str()
'data.frame':   1 obs. of  9 variables:
 $ exposure   : chr "C-Reactive protein level || id:ieu-b-35"
 $ outcome    : chr "Major depression || id:ieu-b-102"
 $ id.exposure: chr "ieu-b-35"
 $ id.outcome : chr "ieu-b-102"
 $ samplesize : num NA
 $ SNP        : chr "rs2794520"
 $ b          : num -0.0258
 $ se         : num 0.0252
 $ p          : num 0.307

MR using just IL6R

subset(ress, SNP == il6r_snp) %>% str()
'data.frame':   1 obs. of  9 variables:
 $ exposure   : chr "C-Reactive protein level || id:ieu-b-35"
 $ outcome    : chr "Major depression || id:ieu-b-102"
 $ id.exposure: chr "ieu-b-35"
 $ id.outcome : chr "ieu-b-102"
 $ samplesize : num NA
 $ SNP        : chr "rs4129267"
 $ b          : num -0.0548
 $ se         : num 0.0503
 $ p          : num 0.275

Multivariable MR of CRP and BMI on MDD

mvexp <- mv_extract_exposures(c("ieu-b-40", "ieu-b-35"))
Please look at vignettes for options on running this locally if you need to run many instances of this command.
Clumping 1, 563 variants, using EUR population reference
Removing 81 of 563 variants due to LD with other variants or absence from LD reference panel
Extracting data for 482 SNP(s) from 2 GWAS(s)
Finding proxies for 1 SNPs in outcome ieu-b-40
Extracting data for 1 SNP(s) from 1 GWAS(s)
Finding proxies for 2 SNPs in outcome ieu-b-35
Extracting data for 2 SNP(s) from 1 GWAS(s)
Harmonising body mass index || id:ieu-b-40 (ieu-b-40) and C-Reactive protein level || id:ieu-b-35 (ieu-b-35)
Removing the following SNPs for being palindromic with intermediate allele frequencies:
rs10778215, rs10887578, rs1144387, rs138289, rs1454687, rs1521527, rs189843, rs2163188, rs2284746, rs2643452, rs486359, rs6011457, rs6595205, rs6815910
mvout <- extract_outcome_data(mvexp$SNP, "ieu-b-102")
Extracting data for 481 SNP(s) from 1 GWAS(s)
mvdat <- mv_harmonise_data(mvexp, mvout)
Harmonising body mass index || id:ieu-b-40 (ieu-b-40) and Major depression || id:ieu-b-102 (ieu-b-102)
Removing the following SNPs for being palindromic with intermediate allele frequencies:
rs10778215, rs10887578, rs1144387, rs138289, rs1454687, rs1521527, rs189843, rs2163188, rs2284746, rs2643452, rs486359, rs6011457, rs6595205, rs6815910
mvres <- mv_multiple(mvdat) 
mvres$result %>% select(exposure, nsnp, b, se, pval)

Reverse MR

datr <- make_dat("ieu-b-102", "ieu-b-35")
Extracting data for 50 SNP(s) from 1 GWAS(s)
Finding proxies for 27 SNPs in outcome ieu-b-35
Extracting data for 27 SNP(s) from 1 GWAS(s)
Harmonising Major depression || id:ieu-b-102 (ieu-b-102) and C-Reactive protein level || id:ieu-b-35 (ieu-b-35)
resr <- mr(datr)
Analysing 'ieu-b-102' on 'ieu-b-35'
mr_scatter_plot(resr, datr)
$`ieu-b-102.ieu-b-35`

attr(,"split_type")
[1] "data.frame"
attr(,"split_labels")

resr %>% select(method, nsnp,b, se, pval)

Human immune system traits from Roederer et al 2015

str(hi_exp)
'data.frame':   479 obs. of  15 variables:
 $ chr.exposure          : chr  "1" "3" "2" "5" ...
 $ se.exposure           : num  0.0226 0.0209 0.0341 0.0254 0.0318 ...
 $ pval.exposure         : num  1.30e-14 3.13e-09 1.11e-08 1.13e-08 1.44e-08 ...
 $ beta.exposure         : num  -0.174 -0.124 0.195 -0.145 0.18 ...
 $ samplesize.exposure   : num  460 448 419 456 429 415 420 415 414 450 ...
 $ pos.exposure          : int  161662484 154835722 68273851 75395725 13470196 20307219 30026413 18046200 72583375 161479745 ...
 $ id.exposure           : chr  "met-b-268" "met-b-128" "met-b-159" "met-b-159" ...
 $ SNP                   : chr  "rs4657090" "rs10513469" "rs10496160" "rs7702044" ...
 $ effect_allele.exposure: chr  "A" "A" "C" "C" ...
 $ other_allele.exposure : chr  "G" "C" "G" "T" ...
 $ eaf.exposure          : num  0.2946 0.3817 0.0621 0.1721 0.0734 ...
 $ exposure              : chr  "1c+mDC:%32+ || id:met-b-268" "B:%Mature || id:met-b-128" "CD123 on 11c+123+DC || id:met-b-159" "CD123 on 11c+123+DC || id:met-b-159" ...
 $ mr_keep.exposure      : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...
 $ pval_origin.exposure  : chr  "reported" "reported" "reported" "reported" ...
 $ data_source.exposure  : chr  "igd" "igd" "igd" "igd" ...
hi_out <- extract_outcome_data(hi_exp$SNP, "ieu-b-102")
Extracting data for 438 SNP(s) from 1 GWAS(s)
Finding proxies for 20 SNPs in outcome ieu-b-102
Extracting data for 20 SNP(s) from 1 GWAS(s)
hi_dat <- harmonise_data(hi_exp, hi_out)
Harmonising MDC:%32+ || id:met-b-266 (met-b-266) and Major depression || id:ieu-b-102 (ieu-b-102)
Removing the following SNPs for being palindromic with intermediate allele frequencies:
rs6997616
Harmonising NKeff:%314-158a+ || id:met-b-241 (met-b-241) and Major depression || id:ieu-b-102 (ieu-b-102)
Removing the following SNPs for being palindromic with intermediate allele frequencies:
rs6469144, rs937352
Harmonising CD4:%RTE || id:met-b-178 (met-b-178) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD123 on 11c+123+DC || id:met-b-159 (met-b-159) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising B:%Mature || id:met-b-128 (met-b-128) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD4mem:%preTh17 (2) || id:met-b-204 (met-b-204) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD4:%Treg(39+73+) || id:met-b-194 (met-b-194) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD161 on CD4mem || id:met-b-147 (met-b-147) and Major depression || id:ieu-b-102 (ieu-b-102)
Removing the following SNPs for being palindromic with intermediate allele frequencies:
rs1385912
Harmonising CD4mem:%cFTH (1) || id:met-b-211 (met-b-211) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD27 on CD8 T || id:met-b-157 (met-b-157) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD27 on IgG+B || id:met-b-146 (met-b-146) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD27 on CD4 T || id:met-b-155 (met-b-155) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising PD1 on CD4mem || id:met-b-148 (met-b-148) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD8:%TM-like || id:met-b-165 (met-b-165) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD27 on IgA+B || id:met-b-144 (met-b-144) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD32 on pDC || id:met-b-143 (met-b-143) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising NK:%Eff || id:met-b-121 (met-b-121) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising NKearly:%337+158b+ || id:met-b-249 (met-b-249) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD4mem:%preTh17 (3) || id:met-b-234 (met-b-234) and Major depression || id:ieu-b-102 (ieu-b-102)
Removing the following SNPs for being palindromic with intermediate allele frequencies:
rs2241006
Harmonising NKeff:%314-R7- || id:met-b-247 (met-b-247) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising NKearly:%335+314- || id:met-b-250 (met-b-250) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD34 on Stem || id:met-b-151 (met-b-151) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD8:%39+ || id:met-b-185 (met-b-185) and Major depression || id:ieu-b-102 (ieu-b-102)
Removing the following SNPs for being palindromic with intermediate allele frequencies:
rs4074424
Harmonising 1c+mDC:%32+ || id:met-b-268 (met-b-268) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD4:%Treg(39+73-) || id:met-b-189 (met-b-189) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD4:%Treg(39+) || id:met-b-191 (met-b-191) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising NKeff:%337+ || id:met-b-239 (met-b-239) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD39 on CD4 T || id:met-b-156 (met-b-156) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD8:%SCM || id:met-b-179 (met-b-179) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD32 on 11c+123+DC || id:met-b-162 (met-b-162) and Major depression || id:ieu-b-102 (ieu-b-102)
Harmonising CD32 on mDC || id:met-b-150 (met-b-150) and Major depression || id:ieu-b-102 (ieu-b-102)
hi_res <- mr(hi_dat)
Analysing 'met-b-121' on 'ieu-b-102'
Analysing 'met-b-128' on 'ieu-b-102'
Analysing 'met-b-143' on 'ieu-b-102'
Analysing 'met-b-144' on 'ieu-b-102'
Analysing 'met-b-146' on 'ieu-b-102'
Analysing 'met-b-147' on 'ieu-b-102'
Analysing 'met-b-148' on 'ieu-b-102'
Analysing 'met-b-150' on 'ieu-b-102'
Analysing 'met-b-151' on 'ieu-b-102'
Analysing 'met-b-155' on 'ieu-b-102'
Analysing 'met-b-156' on 'ieu-b-102'
Analysing 'met-b-157' on 'ieu-b-102'
Analysing 'met-b-159' on 'ieu-b-102'
Analysing 'met-b-162' on 'ieu-b-102'
Analysing 'met-b-165' on 'ieu-b-102'
Analysing 'met-b-178' on 'ieu-b-102'
Analysing 'met-b-179' on 'ieu-b-102'
No SNPs available for MR analysis of 'met-b-185' on 'ieu-b-102'
Analysing 'met-b-189' on 'ieu-b-102'
Analysing 'met-b-191' on 'ieu-b-102'
Analysing 'met-b-194' on 'ieu-b-102'
Analysing 'met-b-204' on 'ieu-b-102'
Analysing 'met-b-211' on 'ieu-b-102'
No SNPs available for MR analysis of 'met-b-234' on 'ieu-b-102'
Analysing 'met-b-239' on 'ieu-b-102'
Analysing 'met-b-241' on 'ieu-b-102'
Analysing 'met-b-247' on 'ieu-b-102'
Analysing 'met-b-249' on 'ieu-b-102'
Analysing 'met-b-250' on 'ieu-b-102'
Analysing 'met-b-266' on 'ieu-b-102'
Analysing 'met-b-268' on 'ieu-b-102'

Does anything surpass 5% FDR threshold?

subset(hi_res, method %in% c("Inverse variance weighted", "Wald ratio")) %>% mutate(fdr=p.adjust(pval, "fdr")) %>%
  arrange(fdr) %>%
  filter(fdr < 0.05) %>%
  select(exposure, nsnp, b, se, pval)
NA
hi_res[which.min(hi_res$pval), ] %>% str()
'data.frame':   1 obs. of  9 variables:
 $ id.exposure: chr "met-b-249"
 $ id.outcome : chr "ieu-b-102"
 $ outcome    : chr "Major depression || id:ieu-b-102"
 $ exposure   : chr "NKearly:%337+158b+ || id:met-b-249"
 $ method     : chr "Wald ratio"
 $ nsnp       : num 1
 $ b          : num -0.11
 $ se         : num 0.0371
 $ pval       : num 0.00302
LS0tCnRpdGxlOiAiSW5mbGFtbWF0aW9uIGFuZCBkZXByZXNzaW9uIgphdXRob3I6ICJHaWJyYW4gSGVtYW5pIgpkYXRlOiAiMjAyMi0wMi0xMCIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3IsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CmxpYnJhcnkoaWV1Z3dhc3IpCmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KFR3b1NhbXBsZU1SKQpsaWJyYXJ5KHJzbnBzKQpgYGAKCgojIyBTdW1tYXJ5CgotIE5vdCBtdWNoIHN1cHBvcnQgZm9yIENSUCBvbiBjYXVzYWwgcGF0aHdheSB0byBNREQsIG1vZGVzdCBwcm90ZWN0aXZlIGVmZmVjdCBvZiBoaWdoZXIgaW5mbGFtbWF0aW9uCi0gRWl0aGVyIHVzaW5nIGFsbCB2YXJpYW50cyBvciBqdXN0IENSUCBvciBJTDZSIHZhcmlhbnRzCi0gTVZNUiBvZiBDUlAgKyBCTUkgb24gTUREIHNob3dzIEJNSSBhc3NvY2lhdGVzIGluZGVwZW5kZW50IG9mIENSUAotIFJldmVyc2UgTVIgb2YgTUREIG9uIENSUCBpcyBudWxsLCBtb2Rlc3QgcG9zaXRpdmUgZWZmZWN0IG9mIGhpZ2hlciBNREQgbGlhYmlsaXR5IHdpdGggaGlnaGVyIENSUAoKIyMgQW5hbHlzaXMKCkRlcHJlc3Npb24gR1dBUzoKCmBgYHtyfQpnd2FzaW5mbygiaWV1LWItMTAyIikgJT4lIHN0cigpCmBgYAoKQy1yZWFjdGl2ZSBwcm90ZWluIEdXQVM6CgpgYGB7cn0KZ3dhc2luZm8oImlldS1iLTM1IikgJT4lIHN0cigpCmBgYAoKUGVyZm9ybSBNUgoKYGBge3J9CmRhdCA8LSBtYWtlX2RhdCgiaWV1LWItMzUiLCAiaWV1LWItMTAyIikKcmVzIDwtIG1yKGRhdCkKbXJfc2NhdHRlcl9wbG90KHJlcywgZGF0KQpgYGAKCmBgYHtyfQpyZXMgJT4lIHNlbGVjdChtZXRob2QsIG5zbnAsIGIsIHNlLCBwdmFsKQpgYGAKClVzZSBvbmx5IENSUCB2YXJpYW50CgpgYGB7cn0KcmVzcyA8LSBtcl9zaW5nbGVzbnAoZGF0KQptcl9mb3Jlc3RfcGxvdChyZXNzKQpgYGAKCkZpbmQgQ1JQIFNOUAoKYGBge3J9CmNycF9zbnAgPC0gc3Vic2V0KGRhdCwgY2hyID09IDEpICU+JQogIG11dGF0ZShkaXN0PWFicyhwb3MuZXhwb3N1cmUgLSAxNTk2ODIwNzkpKSAlPiUKICBhcnJhbmdlKGRpc3QpICU+JQogIHNlbGVjdChkaXN0LCBTTlAsIHB2YWwuZXhwb3N1cmUpICU+JQogIHsuJFNOUFsxXX0KYGBgCgpGaW5kIElMNlIgU05QCgpgYGB7cn0KZGJzbnBfaW5mbyA8LSBuY2JpX3NucF9xdWVyeShkYXQkU05QKQppbDZyX3NucCA8LSBzdWJzZXQoZGJzbnBfaW5mbywgZ2VuZSA9PSAiSUw2UiIpJHJzaWQKYGBgCk1SIHVzaW5nIGp1c3QgQ1JQCgpgYGB7cn0Kc3Vic2V0KHJlc3MsIFNOUCA9PSBjcnBfc25wKSAlPiUgc3RyKCkKYGBgCk1SIHVzaW5nIGp1c3QgSUw2UgoKYGBge3J9CnN1YnNldChyZXNzLCBTTlAgPT0gaWw2cl9zbnApICU+JSBzdHIoKQpgYGAKCk11bHRpdmFyaWFibGUgTVIgb2YgQ1JQIGFuZCBCTUkgb24gTURECgpgYGB7cn0KbXZleHAgPC0gbXZfZXh0cmFjdF9leHBvc3VyZXMoYygiaWV1LWItNDAiLCAiaWV1LWItMzUiKSkKbXZvdXQgPC0gZXh0cmFjdF9vdXRjb21lX2RhdGEobXZleHAkU05QLCAiaWV1LWItMTAyIikKbXZkYXQgPC0gbXZfaGFybW9uaXNlX2RhdGEobXZleHAsIG12b3V0KQptdnJlcyA8LSBtdl9tdWx0aXBsZShtdmRhdCkgCm12cmVzJHJlc3VsdCAlPiUgc2VsZWN0KGV4cG9zdXJlLCBuc25wLCBiLCBzZSwgcHZhbCkKYGBgCgpSZXZlcnNlIE1SCgpgYGB7cn0KZGF0ciA8LSBtYWtlX2RhdCgiaWV1LWItMTAyIiwgImlldS1iLTM1IikKcmVzciA8LSBtcihkYXRyKQptcl9zY2F0dGVyX3Bsb3QocmVzciwgZGF0cikKYGBgCgpgYGB7cn0KcmVzciAlPiUgc2VsZWN0KG1ldGhvZCwgbnNucCxiLCBzZSwgcHZhbCkKYGBgCgojIyBIdW1hbiBpbW11bmUgc3lzdGVtIHRyYWl0cyBmcm9tIFJvZWRlcmVyIGV0IGFsIDIwMTUKCgpgYGB7cn0KZ2kgPC0gZ3dhc2luZm8oKQpoaV9leHAgPC0gZXh0cmFjdF9pbnN0cnVtZW50cyhnaSRpZFtncmVwbCgibWV0LWItIiwgZ2kkaWQpXSkKc3RyKGhpX2V4cCkKYGBgCgpgYGB7cn0KaGlfb3V0IDwtIGV4dHJhY3Rfb3V0Y29tZV9kYXRhKGhpX2V4cCRTTlAsICJpZXUtYi0xMDIiKQpoaV9kYXQgPC0gaGFybW9uaXNlX2RhdGEoaGlfZXhwLCBoaV9vdXQpCmhpX3JlcyA8LSBtcihoaV9kYXQpCmBgYAoKRG9lcyBhbnl0aGluZyBzdXJwYXNzIDUlIEZEUiB0aHJlc2hvbGQ/CgpgYGB7cn0Kc3Vic2V0KGhpX3JlcywgbWV0aG9kICVpbiUgYygiSW52ZXJzZSB2YXJpYW5jZSB3ZWlnaHRlZCIsICJXYWxkIHJhdGlvIikpICU+JSBtdXRhdGUoZmRyPXAuYWRqdXN0KHB2YWwsICJmZHIiKSkgJT4lCiAgYXJyYW5nZShmZHIpICU+JQogIGZpbHRlcihmZHIgPCAwLjA1KSAlPiUKICBzZWxlY3QoZXhwb3N1cmUsIG5zbnAsIGIsIHNlLCBwdmFsKQpgYGAKCmBgYHtyfQpoaV9yZXNbd2hpY2gubWluKGhpX3JlcyRwdmFsKSwgXSAlPiUgc3RyKCkKYGBgCgo=