作者:麦子
转载请注明:解螺旋·临床医生科研成长平台
“相关不等于因果”是一句如雷贯耳的统计学训诫了,每当我们观察到某风险因素与某疾病的相关性时,这句话一定会敲着我们的脑袋,提醒我们注意寻找那些相关性的内在逻辑,再排查一下其他可能的影响因素等等,但这又谈何容易呢?
要排除其他因素的干扰,我们通常采用严格控制变量的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
联系客服