Hi，大家好，我是晨曦

今天这期推文，我们再次踏入孟德尔随机化分析的大门，前面有小伙伴反映希望这期的专栏可以多推出几期，那么自然是没有问题
那么今天这期推文，我们通过一篇6分左右的文献带领大家去掌握孟德尔随机化分析，然后晨曦会在每篇推文的最后再附上关于孟德尔随机化（MR）分析的知识点，希望可以通过这种方式帮助大家快速、准确的掌握MR分析
那么，我们开始吧
今天我们选择的文献题目为：Coffee Consumption and Cardiovascular Diseases: A Mendelian Randomization Study（咖啡摄入量和心血管疾病：孟德尔随机化研究）
注意：我们并不是要复现这篇文献，因为这里面的部分数据晨曦并没有弄到，我们是通过这篇文献的方法学，然后用自己的数据进行演练
我们想要做的是BMI指数和心血管疾病的孟德尔随机化分析
第一步：准备工作
#准备工作library(TwoSampleMR) #加载R包
注意一：TwoSampleMR包可以进行双样本的孟德尔随机化分析即一组样本为工具变量-暴露因素，另一组为暴露因素-结局变量
第二步：获取工具变量-暴露因素的GWAS数据
#读取工具变量和暴露因素数据exposure_dat <-extract_instruments("ieu-a-2")#提取工具变量-暴露因素
注意一：读取数据我们主要有两种方式，一种就是通过各种GWAS数据库或者GWAS文献，找到GWAS数据然后下载下来用TwoSampleMR包进行读取，或者因为TwoSampleMR包本身就相当于一个 端口（简单理解为GEOquery包）所以可以很方便的从IEU OpenGWAS project (mrcieu.ac.uk)这个GWAS数据库上获得数据
注意二：我们这里是通过extract_instruments函数进行读取的，需要了解一下参数：
总结成一句话就是，我们通过设置p1参数找到与暴露因素具有显著相关的工具变量；然后通过设置clump参数去掉连锁不平衡（LD）的工具变量（简单理解就是彼此工具变量太相关了，研究起来没啥意义）；然后我们通过设置r2和kb参数来制定去除LD的标准（默认即可）
注意三：我们这里可以看一下文献中是如何对其进行描述的：
通过设置p1参数为p < 5 × 10−8找到与暴露因素具有显著相关的工具变量，然后通过设置clump参数为TRUE确定去掉连锁不平衡（LDA）的工具变量，然后我们通过设置r2和kb参数分别为1000和0.01，则是说明去掉在1000kb范围内与最显著SNP的r2大于0.01的SNP，这里我们需要保证我们的SNP也就是工具变量大概有100以上进入下轮分析，所以可以发现随着r2的变小与kb的变大，被去除的存在连锁不平衡的SNP会越来越多，而最终剩下的SNP会越来越少
上一步我们获取了工具变量和暴露因素的GWAS数据，那么下面就是获取暴露因素和结局变量的GWAS数据
第三步：获取暴露因素-结局变量的GWAS数据
#获取暴露因素和结局变量数据outcome_dat <-extract_outcome_data(snps=exposure_dat$SNP, outcomes="ieu-a-7")
注意一：既然是获取数据，那么自然也有两种方法，外在以及内在，获取数据这里除了函数不同以外，其实和第二步骤差不多（结局变量的ID号可以通过IEU OpenGWAS project (mrcieu.ac.uk)这个数据库查询得到）
注意二：尽管默认的参数可以直接运行，但是我们必须要了解参数的作用以及意义，才可以让我们的MR分析更加准确，需要了解的参数如下：
总结成一句话就是，我们通过设置snps参数获得工具变量与暴露因素的SNP，然后我们需要与暴露因素和结局变量的SNP取交集，构建工具变量-暴露因素-结局变量这样的关系，所以我们需要提供工具变量与暴露因素的SNP来作为索引；通过设置outcomes参数确定我们的结局变量是什么疾病（下面的示例中我们设定的是冠心病）；通过设置proxies参数为TRUE或者FALSE，可以帮助我们决定是否使用代理SNP（简单理解就是提供的SNP并没有在结局变量中存在，然后设置为TRUE的时候就会通过其它标准选择一些SNP加入其中），这里通常我们可以选择FALSE；通过设置rsq参数我们明确的是LD决定系数（r2），这个参数以及后续的align_alleles和palindromes参数均是在proxies参数为TRUE的时候才会使用（简单理解，通过设置这些参数可以找到当输入的SNP找不到的时候，可以找到代理SNP也就是与输入SNP存在强连锁不平衡的SNP）；通过设置maf_threshold参数可以确定SNP在结局变量中的最小等位基因频率，默认值是0.3，不过大样本GWAS可以适当调低（比如0.01）；大陆用户access_token参数不需要设置或者设置为NULL
这一步骤在文献中就没有对应的解释，但是我们可以看到文献中的描述，其实重点就是挑选工具变量，然后调到了哪些工具变量应该是完成分析以后的事情了
第四步：对数据进行预处理，使其效应等位与效应量保持统一
dat <- harmonise_data( exposure_dat = exposure_dat, outcome_dat = outcome_dat)
注意一：必备步骤，不需要理解，执行即可
第五步：MR分析
res <- mr(dat)#> Analysing 'ieu-a-2' on 'ieu-a-7'res#> id.exposure id.outcome outcome#> 1 ieu-a-2 ieu-a-7 Coronary heart disease || id:ieu-a-7#> 2 ieu-a-2 ieu-a-7 Coronary heart disease || id:ieu-a-7#> 3 ieu-a-2 ieu-a-7 Coronary heart disease || id:ieu-a-7#> 4 ieu-a-2 ieu-a-7 Coronary heart disease || id:ieu-a-7#> 5 ieu-a-2 ieu-a-7 Coronary heart disease || id:ieu-a-7#> exposure method nsnp b#> 1 Body mass index || id:ieu-a-2 MR Egger 79 0.5024935#> 2 Body mass index || id:ieu-a-2 Weighted median 79 0.3870065#> 3 Body mass index || id:ieu-a-2 Inverse variance weighted 79 0.4459091#> 4 Body mass index || id:ieu-a-2 Simple mode 79 0.3401554#> 5 Body mass index || id:ieu-a-2 Weighted mode 79 0.3888249#> se pval#> 1 0.14396056 8.012590e-04#> 2 0.07347663 1.386190e-07#> 3 0.05898302 4.032020e-14#> 4 0.15326612 2.936644e-02#> 5 0.09714470 1.417399e-04#扩展：添加额外的MR方法mr_method_list()#查看支持的方法mr(dat, method_list=c("mr_egger_regression", "mr_ivw"))#添加
注意一：mr函数默认使用五种算法来进行MR分析（ MR Egger，Weighted median，Inverse variance weighted，Simple mode ，Weighted mode ），这五个算法其中最重要的就是Inverse variance weighted（IVW），当然其它算法也是有着各自的用途
注意二：这里我们需要重点解释一下结果：
b（效应值）：正值代表暴露因素的增加会导致结局变量增加（BMI增加会导致心血管疾病的风险增加）（OR是风险比，针对的是二分类变量，如果是连续型变量比如说是身高或者血红蛋白的数量，那么就会变成Beta（Beta和OR可以相互转化），也就是说如果针对的是二分类结局那么需要把B转化为OR，而OR值则是大于1增加风险，小于1则降低风险）
Pval：小于0.05则代表我们的暴露因素与结局变量是否具有统计学意义，那么很显然是具有的
注意三：我们接下来需要重点解释一下五种方法的区别：
简单理解：最重要的方法就是IVW方法，这里可能就会存在一种可能，假设IVW方法是p小于0.05，而其它四种方法都是p大于0.05，那么这个时候我们暴露因素是否与结局变量具有因果关系呢，答案为——是的，也就是IVW的结果就是这么具有决断力，但是需要满足其它四种方法的b（效应值）的方向要与IVW的b值方向相同（需要注意OR和b可以相互转化）
第六步：敏感性分析（有三部分组成|结果出来需要检验是否可靠）
第一部分：异质性检验 mr_heterogeneity
mr_heterogeneity(dat)
注意一：SNP位点是固定的，也就是说不管我测什么疾病/暴露，都是那些SNP，然后双样本的含义其实就是，我选择SNP-暴露和SNP-结局，这就是两组样本，那么自然就会存在人群不同、测序方法不同等很多情况，我们通过异质性检验如果P值小于0.05说明存在异质性，如果存在异质性，我们需要在文章中提及虽然存在异质性，但是不影响IVW的结果，我们的结论是可靠的这一句话，然后结果也是可以接受的
第二部分：水平多效性检验Horizontal pleiotropy
mr_pleiotropy_test(dat)
注意一：P值大于0.05则说明没有水平多效应（简单理解，水平多效应其实就是是否你的研究中存在混杂因素）
第三部分：Leave-one-out analysis
res_loo <- mr_leaveoneout(dat)mr_leaveoneout_plot(res_loo)
注意一：这个分析的本质其实类似于构建预测模型筛选预测变量的时候使用的经典统计迭代法，通过逐个剔除SNP来判断某个SNP是否对结果造成显著性改变
注意二：这块的方法默认使用的反方向加权（inverse variance weighted method），通常我们只需要看All是否大于0，如果大于0证明我们的结果可靠（重点看中心点所在的位置即可）
注意三：我们通过下面的文献可以知道，作者并没有使用常规的敏感性检测方法，这个是很正常的，就像是生信分析中差异分析就有三个R包可以执行，我们这里只需要知道，我们进行敏感性分析的目的就是评估我们的结果是否可靠（后期晨曦也会推送一些其它关于敏感性分析的推文）
#散点可视化p1 <- mr_scatter_plot(res, dat)
第七步：可视化结果（散点图、森林图、漏斗图）

注意一：图上的每一个点代表着一个SNP位点，横坐标是SNP对暴露（BMI）的效应，纵坐标是SNP对结局（心血管疾病）的效应，彩色的线表示的是MR拟合结果。从图中我们不难看出，随着BMI的升高，心血管疾病的发病风险也在升高，同时从这张图上，我们还能看到，当SNP对BMI的效应为0，也就是工具变量的效应为0的时候，我们的结局变量并不是0，说明了我们的数据内存在水平多效应（简单理解就是存在混杂因素）
注意二：我们重点关注的是趋势，截距只要不是偏差很多都是可以接受的
#森林图可视化res_single <- mr_singlesnp(dat)mr_forest_plot(res_single)
注意一：上述森林图中的每一条水平实线反映的是单个SNP利用Wald ratio方法估计出来的结果：有的实线完全在0的左边，说明由这个SNP估计出来的结果是BMI增加能降低心血管疾病的风险；有的实线完全在0的右边，说明由这个SNP估计出来的结果是BMI增加能升高心血管疾病的发病风险；而那些跨过0的说明结果不显著。所以我们单看某一个SNP的结果可能是有问题的，只有把结果综合起来看才能得到合理结果，这就是最底下红线，它反映出IVW方法下BMI的升高增加心血管疾病的发病风险
注意二：然后我们来看一下文献中是如何展示这种可视化结果的（文献中做了一个心血管疾病的大类，小伙伴们可以看红色框内的，其实我们把数据得到以后，就可以单独绘制这个小分栏）
#绘制漏斗图mr_funnel_plot(res_single)
注意一：这个可视化我们需要关注IVW那根线左右两边的点是否大致对称（这里面涉及到算法，简单理解就是MR分析之所以可以被认为是RCT就是因为孟德尔第二定律随机分组，所以按照道理来说左右的点应该是对称分布的，如果有特别离群的点，比如说最左边那个，我们可以去除，然后再次重复上述分析流程）
#如果要去除点，我们需要知道，可视化的本质是数据，所以只需要查看可视化数据即可p1 <- mr_funnel_plot(res_single)p1

至此，全部分析流程讲解结束
然后接下来，我们再来回顾一些关于MR分析的一些理论知识，方便各位小伙伴补充基本的知识基础（评论区好评加更哦~）
1.为了探索疾病的机制，我们通常将风险因素（暴露因素）与终点结局（疾病）之间的关系分类为因果关系与非因果关系，后者其实就是相关，这个通常是观察性研究告诉我们的
2.孟德尔随机化分析（MR分析）的基础是工具变量的创建
3.以往要想回答因果关系的问题，我们需要使用适当的研究设计（区别于传统的观察性研究），比如说使用前瞻性随机实验（RCT），但是RCT的弊端就是费时、费力、费钱，同时有些RCT并不被伦理所支持
4.目前来看，疾病的发展是多基因以及多因素（环境、饮食等等）共同造成的结果，比如说冠心病往往据集在有心脏病家族史的家族内，但是同时高盐高糖高脂肪的饮食也会造成疾病的发生，目前我们通过全基因组关联研究（GWAS）已经发现了数十万甚至百万的遗传变异与疾病结果的关联，这是MR分析的数据基础即一种利用遗传数据来评估可改变的非遗传暴露因素的因果效应的技术
5.例子：免疫炎症是心血管疾病的一个重要原因，目前来说存在炎症假说即免疫炎症机制导致心血管疾病的发生，而干预这一途径将会降低心血管疾病的发病风险，免疫炎症中的biomarker举一个出来——CRP，2010年的观察性研究告诉我们，CRP与冠心病的风险相关，但是我们并不知道这两者是否存在因果关系，是CRP升高导致冠心病发病风险增加，还是冠心病会导致CRP升高，那么这里我们首先可以做出一个假设：CRP的长期升高是否会导致冠心病的发病风险增加，那么我们需要找到一个与CRP水平相关的遗传变异，然后我们就可以通过MR分析得出因果关系结论，这就是MR分析
6.工具变量的限定条件（即满足以下条件才可以被认为是工具变量）
a.工具变量需要与暴露因素相关；
b.工具变量不可以与结局变量相关（工具变量只可以通过暴露因素影响结局）
c.工具变量本身不与其它暴露因素相关（比如说想要研究吸烟和冠心病的关系，但是BMI与冠心病的发生也有关系，那么我们的工具变量只能与吸烟相关，不能与BMI相关）
d.工具变量不能与其它已知的混杂因素相关（性别、年龄等等）
7.工具变量的选择十分重要，其实我们这里可以概括一下MR分析的流程，找寻工具变量和暴露因素的数据、找寻暴露因素和终点结局的数据、MR分析、敏感性分析、可视化，是不是特别像我们预测模型的构建即准备输入数据、构建模型、模型评估、模型可视化，其实本质上都是一样的
好啦~潦潦草草也快1W字了，感谢各位小伙伴的阅读，如果感兴趣MR分析专题，欢迎各位小伙伴在评论区评论哦~
那么，本期推文到这里就结束了
我是晨曦，我们下期再见~
参考文献：
1.TwoSampleMR官网：https://mrcieu.github.io/TwoSampleMR/
2.Hemani G, Tilling K, Davey Smith G. Orienting the causal relationship between imprecisely measured traits using GWAS summary data. PLoS Genet. 2017 Nov 17;13(11):e1007081. doi: 10.1371/journal.pgen.1007081. Erratum in: PLoS Genet. 2017 Dec 29;13(12 ):e1007149. PMID: 29149188; PMCID: PMC5711033.
3.Hemani G, Zheng J, Elsworth B, Wade KH, Haberland V, Baird D, Laurin C, Burgess S, Bowden J, Langdon R, Tan VY, Yarmolinsky J, Shihab HA, Timpson NJ, Evans DM, Relton C, Martin RM, Davey Smith G, Gaunt TR, Haycock PC. The MR-Base platform supports systematic causal inference across the human phenome. Elife. 2018 May 30;7:e34408. doi: 10.7554/eLife.34408. PMID: 29846171; PMCID: PMC5976434.
—
END
—

撰文丨晨   曦
排版丨三叶虫
编辑丨三叶虫
欢迎大家关注解螺旋生信频道-挑圈联靠公号~
继续阅读
阅读原文