前面给大家介绍了孟德尔随机化的概念和它所包含的内容以及一些注意事项,忘记的同学可以复习一下哦(一文带你读懂孟德尔随机化 - 简书 (jianshu.com))。
那么今天就要带着各位同学实际操作一下,我们如何进行孟德尔随机化分析。(这里我们以抑郁症为暴露因素,前列腺癌为结局带着大家进行操作)
#第一步,下载安装R包
install.packages("remotes")
remotes::install_github("MRCIEU/TwoSampleMR")
BiocManager::install("VariantAnnotation")
BiocManager::install("gwasglue")
#读取R包
library(TwoSampleMR)
library(VariantAnnotation)
library(gwasglue)
#第二步读取暴露数据
#在线读取
exposure_dat <- extract_instruments("ukb-b-18336")
#参数默认为P=5*10-8,r方=0.001,kb=10000
注:(在线读取的数据来自于IEU数据库(https://gwas.mrcieu.ac.uk/),
其他数据库数据需先下载本地读入)
#本地读取
bim_vcf <- readVcf("ukb-b-18336.vcf.gz")
a <- gwasvcf_to_TwoSampleMR(vcf=bim_vcf,type="exposure")
b<-subset(a,P<5e-8)
write.csv(b, file="exposure.csv")
(#此步骤是最不同于在线读取的地方。我们需将写出的文件,整理为下方代码读取的格式,然后放在TwoSampleMR包所在的位置处)
bmi<-system.file("exposure.csv",package="TwoSampleMR")
bmi_exp_dat<-read_exposure_data(filename = bmi,sep =",",snp_col ="SNP",beta_col="beta",se_col="standard_error",effect_allele_col="effect_allele",other_allele_col="other_allele",clump = TRUE)
#第三步提取结局数据
outcome_dat<-extract_outcome_data(snps=exposure_dat$SNP,outcomes="ukb-b-13348")
#第四步,将暴露数据和结局数据合并
dat <- harmonise_data(exposure_dat, outcome_dat)
#第五步,孟德尔随机化分析
mrResult=mr(dat)
#对结果进行OR值计算
mrTab=generate_odds_ratios(mrResult)
#异质性分析
heterTab=mr_heterogeneity(dat)
#多效性检验
pleioTab=mr_pleiotropy_test(dat)
#绘制散点图
mr_scatter_plot(mrResult, dat)
#漏斗图
mr_funnel_plot(singlesnp_results = res_single)
#留一法敏感性分析
mr_leaveoneout_plot(leaveoneout_results = mr_leaveoneout(dat))
小结:我们在读取孟德尔随机化数据的时候,一般有两种读取方式,包括在线读取和本地读取方式。通常建议初学者可以先使用在线读取方式进行练习,然后在进行本地读取。而这两种方式的区别就是我们在线读取的数据一般比较久远,而本地数据读取可能因为数据保存的格式不同,整理起来稍显麻烦。好了,今天的分享就到这里,感兴趣的小伙伴赶紧练习起来吧。