Given harmonised data, this function conducts a two-sample MR analysis.
Arguments
- TwoSampleMRinput
Harmonised data.
- mr_plot
one of "None", "TwoSampleMR", "pQTLtools" for no, the original and the revised plots, respectively.
- prefix
a prefix for output files.
Details
As TwoSampleMR faces seemingly perplexing options, this function intends to simplify various steps in a two-sample MR as in Dimou and Tsilidis (2018) . It is particularly useful when a large numbher of MRs are necessary, e.g., multiple proteins and their cis/trans regions need to be examined, in which case prefix could direct the output to various directories.
References
Dimou NL, Tsilidis KK (2018). “A Primer in Mendelian Randomization Methodology with a Focus on Utilizing Published Summary Association Data.” In Evangelou E (ed.), Genetic Epidemiology: Methods and Protocols, chapter 13, 211-230. Springer New York, New York, NY. ISBN 978-1-4939-7868-7, doi:10.1007/978-1-4939-7868-7_13 .
Examples
library(TwoSampleMR)
library(pQTLtools)
outcomes <- "ebi-a-GCST007432"
prot <- "MMP.10"
type <- "cis"
f <- paste0(prot,"-",type,".mrx")
d <- read.table(file.path(find.package("pQTLtools",lib.loc=.libPaths()),"tests",f),
header=TRUE)
exposure <- TwoSampleMR::format_data(within(d,{P=10^logP}), phenotype_col="prot", snp_col="rsid",
chr_col="Chromosome", pos_col="Posistion",
effect_allele_col="Allele1", other_allele_col="Allele2",
eaf_col="Freq1", beta_col="Effect", se_col="StdErr",
pval_col="P", log_pval=FALSE,
samplesize_col="N")
clump <- TwoSampleMR::clump_data(exposure)
#> Please look at vignettes for options on running this locally if you need to run many instances of this command.
#> Clumping vmQoON, 1106 variants, using EUR population reference
#> Server code: 502; Server is possibly experiencing traffic, trying again...
#> Retry succeeded!
#> Removing 1102 of 1106 variants due to LD with other variants or absence from LD reference panel
outcome <- TwoSampleMR::extract_outcome_data(snps=exposure$SNP,outcomes=outcomes)
#> Extracting data for 1106 SNP(s) from 1 GWAS(s)
#> Finding proxies for 155 SNPs in outcome ebi-a-GCST007432
#> Extracting data for 155 SNP(s) from 1 GWAS(s)
#> Server code: 502; Server is possibly experiencing traffic, trying again...
#> Retry succeeded!
harmonise <- TwoSampleMR::harmonise_data(clump,outcome)
#> Harmonising MMP.10 (vmQoON) and FEV1 || id:ebi-a-GCST007432 (ebi-a-GCST007432)
prefix <- paste(outcomes,prot,type,sep="-")
run_TwoSampleMR(harmonise, mr_plot="pQTLtools", prefix=prefix)
#> Analysing 'vmQoON' on 'ebi-a-GCST007432'
#> Warning: Removed 1 rows containing missing values (`geom_errorbarh()`).
#> Warning: Removed 1 rows containing missing values (`geom_point()`).
#> Warning: Removed 1 rows containing missing values (`geom_errorbarh()`).
#> Warning: Removed 1 rows containing missing values (`geom_point()`).
caption <- "Table. MMP.10 variants and FEV1"
knitr::kable(read.delim(paste0(prefix,"-result.txt"),header=TRUE),
caption=paste(caption, "(result)"))
#>
#>
#> Table: Table. MMP.10 variants and FEV1 (result)
#>
#> |id.exposure |id.outcome |outcome |exposure |method | nsnp| b| se| pval|
#> |:-----------|:----------------|:-------------------------------------|:--------|:-------------------------|----:|---------:|-------:|-----:|
#> |vmQoON |ebi-a-GCST007432 |FEV1 || id:ebi-a-GCST007432 |MMP.10 |MR Egger | 4| -0.000577| 0.01001| 0.959|
#> |vmQoON |ebi-a-GCST007432 |FEV1 || id:ebi-a-GCST007432 |MMP.10 |Weighted median | 4| 0.002002| 0.00625| 0.749|
#> |vmQoON |ebi-a-GCST007432 |FEV1 || id:ebi-a-GCST007432 |MMP.10 |Inverse variance weighted | 4| 0.002067| 0.00594| 0.728|
#> |vmQoON |ebi-a-GCST007432 |FEV1 || id:ebi-a-GCST007432 |MMP.10 |Simple mode | 4| 0.004152| 0.00945| 0.690|
#> |vmQoON |ebi-a-GCST007432 |FEV1 || id:ebi-a-GCST007432 |MMP.10 |Weighted mode | 4| 0.001134| 0.00666| 0.876|
knitr::kable(read.delim(paste0(prefix,"-heterogeneity.txt"),header=TRUE),
caption=paste(caption,"(heterogeneity)"))
#>
#>
#> Table: Table. MMP.10 variants and FEV1 (heterogeneity)
#>
#> |id.exposure |id.outcome |outcome |exposure |method | Q| Q_df| Q_pval|
#> |:-----------|:----------------|:-------------------------------------|:--------|:-------------------------|-----:|----:|------:|
#> |vmQoON |ebi-a-GCST007432 |FEV1 || id:ebi-a-GCST007432 |MMP.10 |MR Egger | 0.114| 2| 0.944|
#> |vmQoON |ebi-a-GCST007432 |FEV1 || id:ebi-a-GCST007432 |MMP.10 |Inverse variance weighted | 0.222| 3| 0.974|
knitr::kable(read.delim(paste0(prefix,"-pleiotropy.txt"),header=TRUE),
caption=paste(caption,"(pleiotropy)"))
#>
#>
#> Table: Table. MMP.10 variants and FEV1 (pleiotropy)
#>
#> |id.exposure |id.outcome |outcome |exposure | egger_intercept| se| pval|
#> |:-----------|:----------------|:-------------------------------------|:--------|---------------:|-------:|-----:|
#> |vmQoON |ebi-a-GCST007432 |FEV1 || id:ebi-a-GCST007432 |MMP.10 | 0.00143| 0.00434| 0.774|
knitr::kable(read.delim(paste0(prefix,"-single.txt"),header=TRUE),
caption=paste(caption,"(single)"))
#>
#>
#> Table: Table. MMP.10 variants and FEV1 (single)
#>
#> |exposure |outcome |id.exposure |id.outcome | samplesize|SNP | b| se| p|
#> |:--------|:-------------------------------------|:-----------|:----------------|----------:|:-------------------------------|---------:|-------:|-----:|
#> |MMP.10 |FEV1 || id:ebi-a-GCST007432 |vmQoON |ebi-a-GCST007432 | 321047|rs11225354 | 0.006425| 0.01217| 0.598|
#> |MMP.10 |FEV1 || id:ebi-a-GCST007432 |vmQoON |ebi-a-GCST007432 | 321047|rs142915220 | -0.003993| 0.02856| 0.889|
#> |MMP.10 |FEV1 || id:ebi-a-GCST007432 |vmQoON |ebi-a-GCST007432 | 321047|rs17099622 | 0.004360| 0.02229| 0.845|
#> |MMP.10 |FEV1 || id:ebi-a-GCST007432 |vmQoON |ebi-a-GCST007432 | 321047|rs17860955 | 0.000616| 0.00739| 0.934|
#> |MMP.10 |FEV1 || id:ebi-a-GCST007432 |vmQoON |ebi-a-GCST007432 | 321047|All - Inverse variance weighted | 0.002067| 0.00594| 0.728|
#> |MMP.10 |FEV1 || id:ebi-a-GCST007432 |vmQoON |ebi-a-GCST007432 | 321047|All - MR Egger | -0.000577| 0.01001| 0.959|
knitr::kable(read.delim(paste0(prefix,"-loo.txt"),header=TRUE),
caption=paste(caption,"(loo)"))
#>
#>
#> Table: Table. MMP.10 variants and FEV1 (loo)
#>
#> |exposure |outcome |id.exposure |id.outcome | samplesize|SNP | b| se| p|
#> |:--------|:-------------------------------------|:-----------|:----------------|----------:|:-----------|--------:|-------:|-----:|
#> |MMP.10 |FEV1 || id:ebi-a-GCST007432 |vmQoON |ebi-a-GCST007432 | 321047|rs11225354 | 0.000703| 0.00681| 0.918|
#> |MMP.10 |FEV1 || id:ebi-a-GCST007432 |vmQoON |ebi-a-GCST007432 | 321047|rs142915220 | 0.002341| 0.00608| 0.700|
#> |MMP.10 |FEV1 || id:ebi-a-GCST007432 |vmQoON |ebi-a-GCST007432 | 321047|rs17099622 | 0.001891| 0.00617| 0.759|
#> |MMP.10 |FEV1 || id:ebi-a-GCST007432 |vmQoON |ebi-a-GCST007432 | 321047|rs17860955 | 0.004730| 0.01001| 0.636|
#> |MMP.10 |FEV1 || id:ebi-a-GCST007432 |vmQoON |ebi-a-GCST007432 | 321047|All | 0.002067| 0.00594| 0.728|
for (x in c("result","heterogeneity","pleiotropy","single","loo"))
unlink(paste0(prefix,"-",x,".txt"))