Given harmonised data, this function conducts a two-sample MR analysis.
Details
As TwoSampleMR faces seemingly perplexing options, this function intends to simplify various steps in a two-sample MR as in dt18;textualpQTLtools. 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.
Check your authentication token if the example below fails to run.
Examples
suppressMessages(require(dplyr))
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 <- exposure[sample(1:nrow(exposure),nrow(exposure)/80),] # TwoSampleMR::clump_data(exposure)
# outcome <- TwoSampleMR::extract_outcome_data(snps=exposure$SNP,outcomes="ebi-a-GCST007432")
outcome <- pQTLtools::import_OpenGWAS("ebi-a-GCST007432","11:102090035-103364929","gwasvcf") %>%
as.data.frame() %>%
dplyr::mutate(outcome="FEV1",LP=10^-LP) %>%
dplyr::select(ID,outcome,REF,ALT,AF,ES,SE,LP,SS,id) %>%
setNames(c("SNP","outcome",paste0(c("other_allele","effect_allele","eaf","beta","se",
"pval","samplesize","id"),".outcome")))
#> [1] "https://gwas.mrcieu.ac.uk/files/ebi-a-GCST007432/ebi-a-GCST007432.vcf.gz"
unlink("ebi-a-GCST007432.vcf.gz.tbi")
harmonise <- TwoSampleMR::harmonise_data(clump,outcome)
#> Harmonising MMP.10 (ImbABK) and FEV1 (ebi-a-GCST007432)
prefix <- paste(prot,type,sep="-")
pQTLtools::run_TwoSampleMR(harmonise, mr_plot="pQTLtools", prefix=prefix)
#> Analysing 'ImbABK' on 'ebi-a-GCST007432'
#> Warning: Removed 1 row containing missing values or values outside the scale range (`geom_errorbarh()`).
#> Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).
#> Warning: Removed 1 row containing missing values or values outside the scale range (`geom_errorbarh()`).
#> Warning: Removed 1 row containing missing values or values outside the scale range (`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|
#> |:-----------|:----------------|:-------|:--------|:-------------------------|----:|--------:|-------:|-----:|
#> |ImbABK |ebi-a-GCST007432 |FEV1 |MMP.10 |MR Egger | 11| 0.00883| 0.02618| 0.744|
#> |ImbABK |ebi-a-GCST007432 |FEV1 |MMP.10 |Weighted median | 11| 0.00156| 0.01176| 0.895|
#> |ImbABK |ebi-a-GCST007432 |FEV1 |MMP.10 |Inverse variance weighted | 11| -0.00365| 0.00836| 0.662|
#> |ImbABK |ebi-a-GCST007432 |FEV1 |MMP.10 |Simple mode | 11| 0.00866| 0.01998| 0.674|
#> |ImbABK |ebi-a-GCST007432 |FEV1 |MMP.10 |Weighted mode | 11| 0.01204| 0.01875| 0.535|
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|
#> |:-----------|:----------------|:-------|:--------|:-------------------------|----:|----:|------:|
#> |ImbABK |ebi-a-GCST007432 |FEV1 |MMP.10 |MR Egger | 7.97| 9| 0.537|
#> |ImbABK |ebi-a-GCST007432 |FEV1 |MMP.10 |Inverse variance weighted | 8.22| 10| 0.607|
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|
#> |:-----------|:----------------|:-------|:--------|---------------:|------:|-----:|
#> |ImbABK |ebi-a-GCST007432 |FEV1 |MMP.10 | -0.00156| 0.0031| 0.627|
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 |ImbABK |ebi-a-GCST007432 | 321047|rs12418070 | -0.032669| 0.02311| 0.157|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs12807063 | 0.011429| 0.03238| 0.724|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs17098913 | 0.018028| 0.01972| 0.361|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs17099562 | 0.014941| 0.03349| 0.655|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs17860942 | 0.000563| 0.03151| 0.986|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs17860975 | 0.009259| 0.03056| 0.762|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs2186602 | 0.034356| 0.02822| 0.223|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs2671375 | -0.052478| 0.03499| 0.134|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs4559658 | -0.009709| 0.03689| 0.792|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs76344628 | -0.019860| 0.02103| 0.345|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs7926065 | -0.023011| 0.03356| 0.493|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|All - Inverse variance weighted | -0.003653| 0.00836| 0.662|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|All - MR Egger | 0.008831| 0.02618| 0.744|
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 |ImbABK |ebi-a-GCST007432 | 321047|rs12418070 | 0.000711| 0.00896| 0.937|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs12807063 | -0.004729| 0.00865| 0.585|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs17098913 | -0.008397| 0.00922| 0.363|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs17099562 | -0.004887| 0.00863| 0.571|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs17860942 | -0.003972| 0.00867| 0.647|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs17860975 | -0.004696| 0.00869| 0.589|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs2186602 | -0.007304| 0.00875| 0.404|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs2671375 | -0.000700| 0.00860| 0.935|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs4559658 | -0.003325| 0.00858| 0.698|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs76344628 | -0.000615| 0.00910| 0.946|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|rs7926065 | -0.002374| 0.00863| 0.783|
#> |MMP.10 |FEV1 |ImbABK |ebi-a-GCST007432 | 321047|All | -0.003653| 0.00836| 0.662|
for (x in c("result","heterogeneity","pleiotropy","single","loo"))
unlink(paste0(prefix,"-",x,".txt"))