作者：麦子
转载请注明：解螺旋·临床医生科研成长平台
“相关不等于因果”是一句如雷贯耳的统计学训诫了，每当我们观察到某风险因素与某疾病的相关性时，这句话一定会敲着我们的脑袋，提醒我们注意寻找那些相关性的内在逻辑，再排查一下其他可能的影响因素等等，但这又谈何容易呢？
要排除其他因素的干扰，我们通常采用严格控制变量的RCT设计，但其实RCT也有一些防不胜防的干扰因素，比如社会经济条件啦、生活环境等等，更别说那些由于伦理原因无法设计RCT的课题。那还有没有别的办法来帮我们做因果推断呢？
1986年，Katan提出了一种思想，他想到了大名鼎鼎的豌豆射手孟德尔~
孟德尔发现的遗传规律可概括为：“亲代等位基因随机分配给子代。”也就是说，每个个体出生之前，已经经历了一次生命的大随机~这就成为一个可以利用的“先天分组”。
如果我们因为一些前期观察数据或经验怀疑X表型（比如HDL-C水平降低）会导致Y结局（比如冠心病），那么通常可以设计一个RCT，用一种药物来升高HDL-C水平，如果治疗组比对照组的冠心病（CHD）发生率显著降低了，那么HDL-C和CHD之间有因果关系的证据就大大加强了。
如果很不巧，两个课题组分别做出了互相矛盾的结果又该怎么办？Meta分析是一条路，但现在还有越来越多的人在尝试走另一条路，嗯，就是豌豆大法。
因为现在基因组学的发展也欣欣向荣，十年来GWAS已鉴定出了许多基因变异与表型之间的关系，方便我们利用。能否选择一些与HDL-C相关的基因变异为工具变量（Instrumental Variable，IV），来衔接HDL-C与CHD的关系呢？
一群人携带某些变异，使他们的HDL-C水平较高，另一群人没有这些变异，HDL-C较低，这就分成了两组，可以观察他们发生CHD的风险。基因变异不易受外界环境、社会行为等混杂因素的干扰，事件发生的时序上也比较合理，用来辅助因果推断有很大的优势。
不过要注意的是，IV在其中起了非常关键的作用。在一份去年的综述中，Sekula P等人总结了MR的概念模型（下图A），同时也提醒我们，一个成功的孟德尔随机化分析的IV要符合以下三个前提：
1） IV与目标表型（暴露因素）必须有强关联，并且可重复性良好（图B）；
2） IV与混杂因素没有关联（图C）
3） IV与结局事件没有直接关联，只能通过暴露因素发生关系（图D）。
JASN. 2016 Nov; 27(11): 3253–3265.
做孟德尔随机化分析在数学上也还有好几种方法，这些方法计算出来的结果会有些许差异，也有些方法相对更受推崇（比如Egger法）。但我们暂时不用管太多，因为有个比较简便的操作，即今年8月刚发布的R语言函数包MendelianRandomization，把这些方法打包在一起。
除去前边做回归分析的步骤，它的流程就分为三步：整合数据、选择MR方法、出图。我们具体来看一下。
整合数据
MendelianRandomization包中提供了两份示例数据，一份是28个遗传变异（SNP）与LDL-C、HDL-C、甘油三酯、冠心病（CHD）风险的关系，数据来自Arterioscler Thromb Vasc Biol. 2010 Nov;30(11):2264-76.
另一份是CASR基因上的7个变异与血钙、快测血糖的关系，数据来自Eur J Epidemiol. 2015 Jul;30(7):543-52. 。两份数据的区别在于，前者的基因变异被认为是互相独立的，而后面的7个变异是互相有关联的，它们的关系还需要提供一个相关系数的矩阵，本例中保存在名为calc.rho中的数据框里。
大家可以找这两篇文章来读一读，了解它们的研究背景，这里就不再多说。两种类型的数据在计算时只有一个参数的差别。下边的演示以第一份数据为主，有差别的地方再简单提醒。
我们需要准备的数据主要有四种：基因变异和暴露因素（干预方法）的相关系数（bx），以及bx的标准误（bxse），它们是在前期的单变量回归分析中获得；基因变异和结局事件的相关性（by），以及by的标准误（byse），同样也是由回归分析中获得。回归分析就不在本文中演示了。
其他有用的数据就是基因变异的名称，如rs12785878等。运算时不需要，画图的时候如果有名字就可以标出来了。另外就是相关性系数的矩阵，如前所述的calc.rho，如果没有提供，就直接认为那些变异是互相独立的。
在示例中，暴露因素是各种脂代谢的指标，我们先用LDL-C来演示，其与SNP的相关系保存在ldlc变量中，相应的标准误为ldlcse；结局事件为CHD风险，SNP与其相关性用比值比对数表示，保存在chdlodds变量中，相应的标准误为chdloddsse。
下面将它们用mr_input()函数整合到一个MRInput对象中：
library('MendelianRandomization')
MRInput <- mr_input(bx = ldlc,
bxse = ldlcse,
by = chdlodds,
byse = chdloddsse)
看看这个Input长啥样：
SNP的名字和编号都是默认的，如果你准备好了自己的SNP名，可以在上述语句中加上snps参数来指定。
当几个变异之间有相关性，需要提供一个矩阵时，就加一个corr参数，比如第二份数据的情况就写成这样：
MRInput.cor <- mr_input(bx = calcium,
bxse = calciumse,
by = fastgluc,
byse = fastglucse,
corr = calc.rho)
区别就这点了。后面的步骤咱们采用第一个，MRInput。
选择方法
MendelianRandomization包提供了四种方法，分别为逆方差加权法（Inverse-variance weighted，IVW），中位数法、Egger法、最大似然法（Maximun likelihood），分别用对应的四个函数，每个函数中的method参数又可指定不同的亚类。这四种方法比较常用，不过开发者说MR包后续还会更新，增加其他方法和功能。
比如逆方差加权法，可使用mr_ivw()函数最简单的一句就在括号中填上刚才做好的MRInput：
IVW <- mr_ivw(MRInput)
这样就返回一个表，有估计值、标准误、95%CI、P值等信息：
如果需要更详细的设置，都可加上相应的参数，如加权方式（weighted）、是否采用惩罚加权（penalized）、robust回归（robust）等。如果基因变异之间有关联，别忘了加上correl = T。
其他方法的函数，如中位数法的mr_median()、Egger法的mr_egger()、最大似然法的mr_maxlik()，使用起来都大同小异，把MRInput放括号里就行了，更详细的参数大家可通过help查询。
不过最惊艳的是，本包还有一个全方法一起上的霸道函数mr_allmethods()。使用方法也非常简单：
MRAll <- mr_allmethods(MRInput)
这样便可得到所有方法生成的结果汇总，包括估计值、标准误、95% CI、P值。
作图
作图仿佛是使用R的过程中最有成就感的事，因为效果直观且赏心悦目嘛。那么来吧。该包中只有一个作图函数mr_plot()，能做两种图，区别在于输入的数据。
如果输入咱们第一步做好的MRInput，就得到各个基因变异分别与暴露因素和结局指标关系的交互散点图：
mr_plot(MRInput)
交互就是可以跟用户互动，图上本没有标签，当鼠标悬停于某个点时，就可在弹窗中看到该点的相关信息。这样看起来很清爽，但可能发表不方便，所以你也可以加一个interactive参数将它变成静态。
mr_plot(MRInput, interactive = F, labels = T)
加上labels = T之后每个标签都会显示，就是看着有点犯密恐。自己选一个吧。
此外还可以添加line参数来指定拟合线，是选择IVW的估计值还是Egger的估计值。
另一种图，是输入刚才第二步做好的MRALL数据，用于各种方法的比较。
mr_plot(MRAll)
上图是用LDL-C作为暴露因素做出来的图。如果我们把第一步中的LDL-C改成HDL-C，再跑一遍试试？整个流程如下：
## 整合数据
HDLInput <- mr_input(bx = hdlc,
bxse = hdlcse,
by = chdlodds,
byse = chdloddsse)
## 全方法分析
HDLAll <- mr_allmethods(HDLInput)
## 作图比较
mr_plot(HDLAll)
可见二者有很明显的区别，LDL-C与CHD风险的关系，用各种方法计算时差别都不大，但HDL-C则见到Egger系列与其他方法有着很大的区别。这也是它作为CHD风险因素有争议的来源之一。
这个包要说有缺点，那就是需要安装的东西比较多，因为它后台需要调用几十个其他的函数包。可以理解，既好用、表面又简单的东西，里面一定是个迷宫。如果你之前没有安装过那些包，那么这次会几十个一起安装，所以请耐心、耐心、耐心等待。
参考资料：
1. Emdin, C. A., Khera, A. V. & Kathiresan, S. Mendelian Randomization. JAMA 318, 1925–1926 (2017).
2. 王莉娜 & Zhang Zuofeng. 孟德尔随机化法在因果推断中的应用. 中华流行病学杂志 38, 547–552 (2017).
3. Sekula, P., Del Greco M, F., Pattaro, C. & Köttgen, A. Mendelian Randomization as an Approach to Assess Causality Using Observational Data. J. Am. Soc. Nephrol. 27, 3253–3265 (2016).
4. https://cran.r-project.org/web/packages/MendelianRandomization/vignettes/Vignette_MR.pdf
联系客服
