R数据分析：孟德尔随机化分析文献解析和实例操练

最近抽空研读了一篇探讨高血压和肾功能关系的文献，记录下来分享给大家，主要也是想看看孟德尔随机化的统计分析结果在论文中是如何呈现的，之后我会给大家写写孟德尔随机化的统计分析在R语言中的做法，希望可以帮助到大家。

文章的题目是A bidirectional Mendelian randomization study supports causal effects of kidney function on blood pressure，这篇文章用到的统计技巧叫做Two-sample MR----两样本孟德尔随机分析。我还查阅了别的孟德尔随机化的文献，这个Two-sample MR的分析其实是非常常用的。

Two-sample MR分析的一般步骤

第一步是找工具变量，我们要的是基因作为工具变量这些个基因都是从别人的研究中挑出来的，所有的基因研究有个专门的库叫做genome wide association studies (GWAS)。我们需要做的就是从这个库中挑出来我们自己需要的和我们暴露相关的基因变量SNPs。这是第一步。

第二步就是估计我们的工具变量对结局的作用，工具变量对结局的作用也是从所有的研究中估计出来的整体效应，这样可以拒绝单个研究的偏倚。

第三步就是合并多个SNP的效应量，这个效应量是我们得到暴露和结局因果效应的前提。

第四步就是用合并后的数据进行孟德尔随机化分析和相应的敏感性分析。

做分析的整个流程就在下面的图中啦：

总体来看就是在孟德尔随机化研究中我们的工具变量可以不需要你收集，工具变量的效应也不需要你计算，这些都只需要你在GWAS挖掘合并就行。就是说做孟德尔随机化研究是不需要你有原始数据的。

我们把上面的步骤具体在刚刚提到的文献中走一遍：

这篇文献是要研究高血压和肾功能的因果方向的，就是到底是高血压导致肾功能下降，还是肾功能下降导致的高血压，具体地就是研究eGFRcr和BP的因果方向。

首先作者从别人的基因研究中找自己研究变量的工具变量，别人的研究的情况如下表：注意下表是包含一个联盟的很多个研究的(肾功能的工具变量是从CKDGen Consortium找来的，血压的工具变量是UKB-ICBP中找来的)，是需要进行meta整合的：

通过meta分析作者就筛选出了两个变量可能的工具变量，因为每个变量的工具变量其实是比较多的，为了保证同一个变量工具变量间的独立性，作者有做一个叫LD clumping的操作：

To ensure independence among genetic instruments, we applied LD clumping60 with a clumping window of 10 MB and an r2 cutoff of 0.001 (default of the ld_clumpfunction)

作者有把筛出来的变量的暴露和结局的工具变量展示在文献中（但是放在补充材料中的，我并没有能找到，欸）

然后就到第二步和第三步，估计工具变量对暴露和结局的作用，这个时候要考虑工具变量一定不能直接影响结局（叫做pleiotropy），所以作者会用好几个算法来估计SNP的作用，并将多个SNP的效应合并，用到的是harmonise_data这个函数。

工具变量其实有很多的，所以就有上面提到的pleiotropic问题，作者是用不同的方法来估计参数（inverse variance weighted method, mendelian randomisation-Egger (MR-Egger) method, weighted median method, and weighted mode based estimation）来均衡pleiotropic问题的影响，最后得到一个总的合并后的效应，这个效应被认为是比较稳健的：

We applied four complementary methods of two sample mendelian randomisation (inverse variance weighted method, mendelian randomisation-Egger (MR-Egger) method, weighted median method, and weighted mode based estimation), which make different assumptions about horizontal pleiotropy. A consistent effect across the four methods is less likely to be a false positive.

有了合并后的工具变量对结局的效应，最后一步就是mr分析，估计暴露与结局的直接因果效应。

整个文献就是这么回事，正文中呈现出来的图表包括下面两个，一个是肾功能对血压的作用的森林图和效应的置信区间和p值，另一个是血压对肾功能的作用的森林图和效应的置信区间和p值：

最终作者得到肾功能是影响血压的原因，而非相反。整个文献就完啦，这个就是一个完整的孟德尔随机化研究的统计呈现过程。

Two-sample MR实操

接下来给大家写如何做出文献中的结果：

Two sample Mendelian randomisation (2SMR) is a method to estimate the causal effect of an exposure on an outcome using only summary statistics from genome wide association studies (GWAS).

具体地，作者使用的R包是‘TwoSampleMR’。

使用TwoSampleMR的基本流程包括4步：

选择工具变量

从GWAS数据库提取工具变量

合并效应量

做MR分析，敏感性分析，画图，出报告

在分析中我们要使用SNPs来作为工具变量：

Mendelian randomization is a method to assess the causal effect of an exposure on an outcome using an instrument, defined by one or more SNPs, as a proxy for the exposure.

这个工具变量怎么来呢，有人已经把所用的可能用到的SNPs总结好了放在了R包里，我们直接用就行，从GWAS中获取，我们需要在GWAS中找到暴露相关的工具变量和结局相关的工具变量从而形成两个样本，相应的分析跑两次，所以一般文献中的分析方法就叫做two-sample Mendelian randomization

Methodological advances mean that Mendelian randomization can be implemented using summary statistics from GWAS, without individual level data. This requires SNP-exposure associations and SNP-outcome associations obtained from separate datasets and is known as two-sample Mendelian randomization

这里又有许多术语需要给大家做做铺垫：

孟德尔随机化：Mendelian randomization is a method to assess the causal effect of an exposure on an outcome using an instrument, defined by one or more single nucleotide polymorphisms, as a proxy for exposure.（SNP就是工具变量）

Genome-wide association study (GWAS)：Genome-wide association studies identify the genetic variants that are associated with a given phenotype.和暴露（表型）相关的基因都是从GWAS中找的

Heterogeneity：Heterogeneity is defined as the variation in the causal estimate across SNPs.这个是不同SNP效应的异质性，所以多个SNP的效应是需要合并的

其实还有好多，感觉写不完。。。直接上例子吧

实例解析

我现在想通过孟德尔随机化研究血压和冠心病的关系，首先我们得接入孟德尔随机基础库MR-Base这个数据库，找到相应的结局（即血压和冠心病）的相关基因，代码如下：

ao<- available_outcomes()exposure_dat<- extract_instruments(c('ukb-a-360'))outcome_dat<- extract_outcome_data(exposure_dat$SNP, c('7'), proxies=1, rsq = 0.8, align_alleles = 1, palindromes=1, maf_threshold = 0.3)

运行上面的代码然后你去查看对象exposure_dat或者outcome_dat就会看到相关基因研究的基本特征，给大家截个图吧：

然后计算基因对结局的合并效应量，代码如下：

dat<-harmonise_data(exposure_dat,outcome_dat,action=2)

最后一步就是进行mr分析：

mr_results<- mr(dat)

运行所有的代码后我们就可以得到不同算法模型下我们的结果如下图两图：

文献中的结果

我跑出来的结果

可以看到上面的图是文献中的，底下的图是我自己复现出来的，一模一样哦。可以看出来血压是对冠心病有因果影响的，我们还可以得到相应系数的置信区间：

generate_odds_ratios(mr_results)

同时因为我们是选了很多个SNP进行的效应估计，我们可以出整个SNP效应的森林图和漏斗图，漏斗图可以帮我们看SNP的异质性，我们还可以用Cochran’s Q and I2统计量检验不同SNP的异质性，然后选择正确的方法对各个SNP的效应进行合并得到总效应：

从森林图就可以看到，整个血压的工具变量对肾功能的作用也是显著的，同时我们可以出散点图直接看出来血压对肾功能的系数，如下：

上图就是不同参数估计方法下面估计出来的暴露（血压）对结局（肾功能）的因果效应，这个图的纵轴是工具变量对结局的影响效应，横轴是工具变量对暴露的因果效应，两个效应做比值就是暴露对结局的效应，也就是图中线的斜率，也可以看出来这个线不同的算法估计出来是有些许差异的，但总体是斜向上的，就意味着暴露（血压）对结局（冠心病）是有着正向的因果关系的。

然后再用同样的方法做一遍冠心病对血压的因果关系的分析（流程完全一样就不重复了），做完之后发现冠心病对血压其实是没有因果作用的，这样我们就知道了其实是血压导致了冠心病而非冠心病导致了血压升高，这样一个完整的孟德尔随机化研究就完成了。

可以看到一个孟德尔随机化研究好像也不难哦。好啦，今天就到这儿吧，今天有粉丝在通话中鼓励我继续（赶紧）写公众号文章，有被感动到，所以特更一期，哈哈。

小结

今天从一篇文献出发给大家写了孟德尔随机化的流程和做法，希望可以帮助到大家，感谢大家耐心看完，自己的文章都写的很细，代码都在原文中，希望大家都可以自己做一做，请关注后私信回复“数据链接”获取所有数据和本人收集的学习资料。如果对您有用请先收藏，再点赞转发。

也欢迎大家的意见和建议，大家想了解什么统计方法都可以在文章下留言，说不定我看见了就会给你写教程哦，另欢迎私信。

本文参考文献：

Wootton RE, Lawn RB, Millard LA, et al. : Evaluation of the causal effects between subjective wellbeing and cardiometabolic health: mendelian randomisation study. BMJ. 2018;362:k3788. 10.1136/bmj.k3788

A bidirectional Mendelian randomization study supports causal effects of kidney function on blood pressureYu, Zhi et al.Kidney International, Volume 98, Issue 3, 708 - 716

Clinical effect of naturally random allocation to lower systolic blood pressure beginning before the development of hypertension.Ference BA, Julius S, Mahajan N, Levy PD, Williams KA Sr, Flack JMHypertension. 2014 Jun; 63(6):1182-8.

Bowden, Jack, George Davey Smith, and Stephen Burgess. 2015. “Mendelian randomization with invalid instruments: effect estimation and bias detection through Egger regression.” International Journal of Epidemiology In press.