两样本孟德尔随机化：代码实操，一学就会

#install.packages("rlang")
library(rlang)

# Install packages
#install.packages("remotes")
#remotes::install_github("MRCIEU/TwoSampleMR")
#remotes::install_github("MRCIEU/ieugwasr")
#remotes::install_github("MRCIEU/MRInstruments")
#remotes::install_github("n-mounier/MRlap")
#install.packages("data.table")
#install.packages("devtools")
#install.packages("dplyr")
#install.packages("xlsx")

# Load Packages
library(TwoSampleMR)
library(ieugwasr)
library(MendelianRandomization)
library(MRInstruments)
library(MRlap)
library(dplyr)
library(xlsx)

> gwasinfo("ukb-b-12493")  ##原发性高血压
           id                                                             trait
1 ukb-b-12493 Diagnoses - secondary ICD10: I10 Essential (primary) hypertension
                                                                                 note
1 41204#I10: Output from GWAS pipeline using Phesant derived variables from UKBiobank
  sample_size       build year category ncase consortium ncontrol    nsnp
1      463010 HG19/GRCh37 2018   Binary 54358    MRC-IEU   408652 9851867
  group_name population priority       author ontology unit               sex mr
1     public   European        1 Ben Elsworth       NA   SD Males and Females  1
  subcategory
1          NA

hyperten_tophits <- tophits(id="ukb-b-12493", clump=0)
hyperten_gwas <- rename(hyperten_tophits, c(
  "SNP"="rsid",
  "effect_allele.exposure"="ea",
  "other_allele.exposure"="nea",
  "beta.exposure"="beta",
  "se.exposure"="se",
  "eaf.exposure"="eaf",
  "pval.exposure"="p",
  "N"="n"))

hyperten_liberal <- clump_data(hyperten_gwas, clump_kb = 10000, clump_r2 = 0.05,) ##参数可调
## Removing 2685 of 2773 variants due to LD with other variants or absence from LD reference panel
hyperten_liberal$id.exposure = "hyperten - liberal clump"

# Calculate R2 and F stat for exposure data
# Liberal hypertension F stat
hyperten_liberal$r2 <- (2 * (hyperten_liberal$beta.exposure^2) * hyperten_liberal$eaf.exposure * (1 - hyperten_liberal$eaf.exposure)) /
  (2 * (hyperten_liberal$beta.exposure^2) * hyperten_liberal$eaf.exposure * (1 - hyperten_liberal$eaf.exposure) +
     2 * hyperten_liberal$N * hyperten_liberal$eaf.exposure * (1 - hyperten_liberal$eaf.exposure) * hyperten_liberal$se.exposure^2)

hyperten_liberal$F <- hyperten_liberal$r2 * (hyperten_liberal$N - 2) / (1 - hyperten_liberal$r2)
hyperten_liberal_meanF <- mean(hyperten_liberal$F)

# Find SNPs in outcome data
ao <- available_outcomes()

anxiety_hyperten_liberal <- extract_outcome_data(snps = hyperten_liberal$SNP, outcomes = "ukb-b-11311")
##
Extracting data for 88 SNP(s) from 1 GWAS(s)
Finding proxies for 34 SNPs in outcome ukb-b-11311
Extracting data for 34 SNP(s) from 1 GWAS(s)
Server code: 502; Server is possibly experiencing traffic, trying again...
Server code: 502; Server is possibly experiencing traffic, trying again...
Server code: 502; Server is possibly experiencing traffic, trying again...
Retry succeeded!

gwasinfo("ukb-b-11311") ##焦虑
  id                                                            trait
1 ukb-b-11311 Diagnoses - secondary ICD10: F41.9 Anxiety disorder, unspecified
                                                                                  note
1 41204#F419: Output from GWAS pipeline using Phesant derived variables from UKBiobank
  sample_size       build year category ncase consortium ncontrol    nsnp
1      463010 HG19/GRCh37 2018   Binary  1523    MRC-IEU   461487 9851867
  group_name population priority       author ontology unit               sex mr
1     public   European        1 Ben Elsworth       NA   SD Males and Females  1
  subcategory
1          NA

# Harmonise datasets
harm_anxiety_hyperten_liberal <- harmonise_data(exposure_dat = hyperten_liberal, outcome_dat = anxiety_hyperten_liberal)
##
Harmonising hyperten (hyperten - liberal clump) and Diagnoses - secondary ICD10: F41.9 Anxiety disorder, unspecified || id:ukb-b-11311 (ukb-b-11311)
Removing the following SNPs for being palindromic with intermediate allele frequencies:
rs1870735, rs35184780

> dim(harm_anxiety_hyperten_liberal)
[1] 55 43

anxiety_hyperten_liberal_results <- mr(dat = harm_anxiety_hyperten_liberal)

# F stats for each analysis
anxiety_hyperten_liberal_meanF <- mean(harm_anxiety_hyperten_liberal$F)