孟德尔随机化（文献复现）

两样本孟德尔随机化-文献复现

---

前言

关于孟德尔随机化的复现，使用R语言
文献原文：PMID: 36816044

一、写作思路

1.先前的研究表明吸烟对衰弱的发生起促进作用，然而在饮酒方面有着矛盾的发现，一些研究表明与从不饮酒的人群相比，重度饮酒可能会降低衰弱风险，其他研究则没有这一结论。
2.目前对于吸烟和饮酒对于衰弱影响的研究多是观察性研究，易受混杂因素干扰。
3.可能存在反向因果关系

二、数据来源

作者采取两样本孟德尔随机化（ two-sample Mendelian randomization ）来探究吸烟和饮酒与衰弱之间的因果关系，本文只复现吸烟与衰弱关系。
1.暴露数据
GWAS原文：PMID: 36477530
在suppl.中，表单ST2meta分析结果，选择暴露为SmkInit（开始吸烟 文献中有具体说明）人种为EUR（欧洲世系），粘贴到文本 txt格式（方便下文分析）

2.结局数据
IEU-GWAS库，在线搜索frailty index，复制GWAS ID

三、MR代码

1.读取暴露数据

代码如下：

library('TwoSampleMR')

#读取暴露数据
exp_dat<-read_exposure_data(filename = "smkinit.txt",
                            sep='\t',
                            snp_col = "rsID",
                            beta_col = "Variant Effect Size",
                            se_col = "Standard Error",
                            effect_allele_col = "Alternate Allele",
                            other_allele_col = "Reference Allele",
                            eaf_col = "Alternate Allele Frequency",
                            pval_col = "P-value",
                            samplesize_col = "N")
write.csv(exp_dat, file="exp_dat.csv")

2.读取结局数据

代码如下：

out_dat<-extract_outcome_data(snps=exp_dat$SNP,
                               outcomes = 'ebi-a-GCST90020053')#在线抓取 来自网络 默认找代理

3.协调数据

dat <- harmonise_data(exposure_dat =  exp_dat, outcome_dat = out_dat,action=1) #默认action=2
write.csv(dat, file="dat.csv")

这里action取值我也不是很确定，只是取2和文献作者的结果更接近

4.进行mr分析

res <- mr(dat,method_list = c("mr_ivw","mr_egger_regression","mr_ivw_radial","mr_weighted_median","mr_weighted_mode"))
generate_odds_ratios(res) #置信区间

置信区间的图就不放了
可以看到，本文与文献结果仍有些许差异，可能是由于代理寻找不同

5.异质性检验

mr_heterogeneity(dat)

异质性检验Q统计量 P<0.05，存在异质性，文中处理方式：
去掉潜在离群值，重新进行分析

run_mr_presso(dat,NbDistribution = 1000)

在跑离群值时，
报错：Error in MRPRESSO::mr_presso(BetaOutcome = “beta.outcome”, BetaExposure = “beta.exposure”, :
Not enough elements to compute empirical P-values, increase NbDistribution
因此离群值这部分没做处理。

6.水平多效性检验

mr_pleiotropy_test(dat)

可以看到，截距项 P>0.05，不存在水平多效性，满足工具变量核心假设

7.单个snp效应值+leave one out

#获取单独的每个SNP的效应值
res_single <- mr_singlesnp(dat) #方便森林图
#去一法分析
res_loo <- mr_leaveoneout(dat)
#散点图 
p1 <- mr_scatter_plot(res, dat)
p1[[1]] 
#森林图
p2 <- mr_forest_plot(res_single)
p2[[1]]
#去一法图 
p3 <- mr_leaveoneout_plot(res_loo)
p3[[1]]

森林图和leave one out 作图实在太丑了，但是没有相应的参数可以调整

总结

1.工具变量的选择时即处理暴露数据没有做独立性与相关性检验，是因为暴露数据来源于GWAS原文，已经满足独立性与相关性要求；
2.敏感性分析：单个snp效应驱动（leave-one-out）以及去掉离群值重新进行mr分析，讨论效应大小变化。