前面给大家介绍了孟德尔随机化的概念和它所包含的内容以及一些注意事项(详见上篇文章:一文带你读懂孟德尔随机化),那么今天就带着各位同学实际操作一下,我们如何进行孟德尔随机化分析。
01
#第一步 安装加载包
install.packages("remotes")
remotes::install_github("MRCIEU/TwoSampleMR")
#加载R包
library(TwoSampleMR)
02
#第二步 读取在线数据
exposure_dat <- extract_instruments("ukb-b-18336")
#参数默认为P=5*10-8,r方=0.001,kb=10000
注:(在线读取的数据来自于IEU数据库(https://gwas.mrcieu.ac.uk/),
其他数据库数据需先下载本地读入)
03
#第三步 计算F值
#建立计算公式
MR_F <- function(sample_size,num_IVs,r_square){
numberator <- r_square * (sample_size - 1 - num_IVs)
denominator <- (1 - r_square) * num_IVs
f <- numberator / denominator return(f) }
my_PVE <- function(beta, se_beta, maf, N){
numberator <- 2 * beta^2 * maf * (1 - maf)
denominator <- 2 * beta^2 * maf * (1 - maf) + se_beta^2 * 2 * N*maf * (1 - maf) return(numberator/denominator) }
#计算R^2
exposure_dat$pve<- my_PVE(beta = exposure_dat$beta.exposure,
se_beta = exposure_dat$se.exposure,
maf = exposure_dat$eaf.exposure,
N = exposure_dat$samplesize.exposure)
#得到F值
IV_infor <- exposure_dat %>% group_by(exposure) %>% summarise(samplesize=mean(samplesize.exposure),nIV=n(),R2=sum(pve))
IV_infor$F_stat<-MR_F(sample_size=IV_infor$samplesize, num_IVs = IV_infor$nIV, r_square = IV_infor$R2)
04
#第四步 提取结局数据
outcome_dat<-extract_outcome_data(snps=exposure_dat$SNP,outcomes="ukb-b-13348")
05
#第五步 将暴露数据和结局数据合并
dat <- harmonise_data(exposure_dat, outcome_dat)
06
#第六步 孟德尔随机化分析
mrResult=mr(dat)
mrTab=generate_odds_ratios(mrResult)
07
#第七步 异质性分析
heterTab=mr_heterogeneity(dat)
08
#第八步 多效性检验
pleioTab=mr_pleiotropy_test(dat)
09
#第九步 绘制散点图
mr_scatter_plot(mrResult, dat)
#mr的结果可视化
10
#第十步 森林图
res_single=mr_singlesnp(dat)
mr_forest_plot(res_single)
#得到每个工具变量对结局的影响
11
#第十一步 漏斗图
mr_funnel_plot(singlesnp_results = res_single)
12
#第十二步 留一法敏感性分析
mr_leaveoneout_plot(leaveoneout_results = mr_leaveoneout(dat))
小结:我们在读取孟德尔随机化数据的时候,一般有两种读取方式,包括在线读取和本地读取方式。通常建议初学者可以先使用在线读取方式进行练习,然后再进行本地读取。好了,今天的分享就到这里,后面我们会继续给大家讲解本地读取的方式,感兴趣的小伙伴赶紧练习起来吧。