Skip to contents

Given harmonised data, this function conducts a two-sample MR analysis.

Usage

run_TwoSampleMR(TwoSampleMRinput, mr_plot = "None", prefix = "")

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.

Value

No value is returned but several files.

Details

As TwoSampleMR faces seemingly perplexing options, this function intends to simplify various steps in a two-sample MR. 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.

Note

Adapted from script by Dimou NL, Tsilidis KK.

References

Dimou NL, Tsilidis KK. A Primer in Mendelian Randomization Methodology with a Focus on Utilizing Published Summary Association Data. In Evangelos Evangelou (ed.), Genetic Epidemiology: Methods and Protocols, Methods in Molecular Biology, vol. 1793, https://doi.org/10.1007/978-1-4939-7868-7_13, Springer Science+Business Media, LLC, part of Springer Nature 2018

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 oojDTr, 1106 variants, using EUR population reference
#> 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)
harmonise <- TwoSampleMR::harmonise_data(clump,outcome)
#> Harmonising MMP.10 (oojDTr) and FEV1 || id:ebi-a-GCST007432 (ebi-a-GCST007432)
prefix <- paste(outcomes,prot,type,sep="-")
run_TwoSampleMR(harmonise, mr_plot="pQTLtools", prefix=prefix)
#> Analysing 'oojDTr' 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|
#> |:-----------|:----------------|:-------------------------------------|:--------|:-------------------------|----:|---------:|-------:|-----:|
#> |oojDTr      |ebi-a-GCST007432 |FEV1 &#124;&#124; id:ebi-a-GCST007432 |MMP.10   |MR Egger                  |    4| -0.000577| 0.01001| 0.959|
#> |oojDTr      |ebi-a-GCST007432 |FEV1 &#124;&#124; id:ebi-a-GCST007432 |MMP.10   |Weighted median           |    4|  0.002002| 0.00634| 0.752|
#> |oojDTr      |ebi-a-GCST007432 |FEV1 &#124;&#124; id:ebi-a-GCST007432 |MMP.10   |Inverse variance weighted |    4|  0.002067| 0.00594| 0.728|
#> |oojDTr      |ebi-a-GCST007432 |FEV1 &#124;&#124; id:ebi-a-GCST007432 |MMP.10   |Simple mode               |    4|  0.004152| 0.00893| 0.674|
#> |oojDTr      |ebi-a-GCST007432 |FEV1 &#124;&#124; id:ebi-a-GCST007432 |MMP.10   |Weighted mode             |    4|  0.001134| 0.00747| 0.889|
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|
#> |:-----------|:----------------|:-------------------------------------|:--------|:-------------------------|-----:|----:|------:|
#> |oojDTr      |ebi-a-GCST007432 |FEV1 &#124;&#124; id:ebi-a-GCST007432 |MMP.10   |MR Egger                  | 0.114|    2|  0.944|
#> |oojDTr      |ebi-a-GCST007432 |FEV1 &#124;&#124; 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|
#> |:-----------|:----------------|:-------------------------------------|:--------|---------------:|-------:|-----:|
#> |oojDTr      |ebi-a-GCST007432 |FEV1 &#124;&#124; 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 &#124;&#124; id:ebi-a-GCST007432 |oojDTr      |ebi-a-GCST007432 |     321047|rs11225354                      |  0.006425| 0.01217| 0.598|
#> |MMP.10   |FEV1 &#124;&#124; id:ebi-a-GCST007432 |oojDTr      |ebi-a-GCST007432 |     321047|rs142915220                     | -0.003993| 0.02856| 0.889|
#> |MMP.10   |FEV1 &#124;&#124; id:ebi-a-GCST007432 |oojDTr      |ebi-a-GCST007432 |     321047|rs17099622                      |  0.004360| 0.02229| 0.845|
#> |MMP.10   |FEV1 &#124;&#124; id:ebi-a-GCST007432 |oojDTr      |ebi-a-GCST007432 |     321047|rs17860955                      |  0.000616| 0.00739| 0.934|
#> |MMP.10   |FEV1 &#124;&#124; id:ebi-a-GCST007432 |oojDTr      |ebi-a-GCST007432 |     321047|All - Inverse variance weighted |  0.002067| 0.00594| 0.728|
#> |MMP.10   |FEV1 &#124;&#124; id:ebi-a-GCST007432 |oojDTr      |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 &#124;&#124; id:ebi-a-GCST007432 |oojDTr      |ebi-a-GCST007432 |     321047|rs11225354  | 0.000703| 0.00681| 0.918|
#> |MMP.10   |FEV1 &#124;&#124; id:ebi-a-GCST007432 |oojDTr      |ebi-a-GCST007432 |     321047|rs142915220 | 0.002341| 0.00608| 0.700|
#> |MMP.10   |FEV1 &#124;&#124; id:ebi-a-GCST007432 |oojDTr      |ebi-a-GCST007432 |     321047|rs17099622  | 0.001891| 0.00617| 0.759|
#> |MMP.10   |FEV1 &#124;&#124; id:ebi-a-GCST007432 |oojDTr      |ebi-a-GCST007432 |     321047|rs17860955  | 0.004730| 0.01001| 0.636|
#> |MMP.10   |FEV1 &#124;&#124; id:ebi-a-GCST007432 |oojDTr      |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"))