孟德尔随机化（Mendelian randomization，MR）是借助遗传变异（genetic variation）作为工具变量（instrumental variable，Ⅳ）来推断暴露因素与结局之间因果关联的方法，有效避免了反向因果关联和潜在混杂因素导致的偏倚。随着全基因组关联研究（genome-wide association study，GWAS）的增多，孟德尔随机化分析广泛应用于观察性流行病学研究中。传统孟德尔随机化分析方法，如逆方差加权法（inverse-variance weighted，IVW），假定所有的遗传变异均满足工具变量的3个核心假设^[1]：（1）工具变量与暴露因素之间有强相关（关联性假设）；（2）工具变量与混杂因素无关联（独立性假设）；（3）工具变量只能通过暴露因素对结局产生作用，不能通过其他途径影响结局发生（排他性假设）。当工具变量存在多效性时，因果效应的估计则会有偏^[2]。MR-Egger回归是在IVW的基础上修正而来的一种基于汇总数据的多工具变量孟德尔随机化方法。与IVW不同的是，该方法仅需满足工具变量多效性效应独立于工具变量与暴露因素之间的关联（instrument strength independent of direct effect，InSIDE）假设和无测量误差（no measurement error，NOME）假设，不如工具变量的3个核心假设要求严格。同时，该方法既能检测多效性又能校正多效性偏倚，所以在以多个遗传变异作为工具变量的研究中，采用MR-Egger回归能够保持孟德尔随机化方法的有效性。本研究将对MR-Egger回归进行阐述，并用此方法探索高密度脂蛋白胆固醇（high-density lipoprotein cholesterol，HDL-C）与冠状动脉疾病（coronary artery disease，CAD）之间的因果关联。

资料和方法

基本原理  MR-Egger回归是由Bowden等^[3]于2015年提出，广泛运用于采用多个单核苷酸多态性（single nucleotide polymorphism，SNP）作为工具变量进行因果推断的孟德尔随机化分析中，尤其适用于遗传变异存在定向多效性的情况。该方法的核心在于加权线性回归时考虑截距项的存在，利用截距项来衡量工具变量间平均多效性的大小，斜率则是因果效应的无偏估计。MR-Egger回归的效用取决于两个前提：第一，InSIDE假设要求SNP对暴露因素的影响独立于它们对结局的多效性影响，弱化了排他性假设；第二，NOME假设要求工具变量与暴露因素的关联中不存在测量误差，$ {I}^{2} $统计量可以评估是否满足NOME假设，如果不满足，则会产生弱工具变量偏倚^[4]。

记J个工具变量为Gj（j=1，2，…，J），暴露因素为X，结局为Y，工具变量Gj对暴露因素X和结局Y的效应估计值分别为$ {\widehat{\beta }}_{{X}_{j}} $和$ {\widehat{\beta }}_{{Y}_{j}} $及其对应的标准误为$ {\sigma }_{{\widehat{\beta }}_{{X}_{j}}} $和$ {\sigma }_{{\widehat{\beta }}_{{Y}_{j}}} $，遗传变异与结局之间的关联可表示为：

$ {\beta }_{{Y}_{j}}={\alpha }_{j}+\theta {\beta }_{{X}_{j}} $

其中，αj表示工具变量Gj对结局Y的直接影响，当αj≠0时，工具变量具有多效性；θ表示暴露因素X对结局Y的影响；$ \theta {\beta }_{{X}_{j}} $则为工具变量Gj通过暴露因素X对结局Y的间接影响^[5]。

对于满足工具变量所有假设的单一遗传变异Gj（αj = 0），暴露因素X对结局Y的因果效应可被一致性地估计为关联估计的比值，比值及相应的标准误分别为：$ {\widehat{\theta }}_{j}={\widehat{\beta }}_{{Y}_{j}}/{\widehat{\beta }}_{{X}_{j}}$, ${\sigma }_{{\widehat{\theta }}_{j}}={\sigma }_{{\widehat{\beta }}_{{Y}_{j}}}/{\widehat{\beta }}_{{X}_{j}}$ 。

由于MR-Egger回归并没有限制截距项为0，回归方程可表示为：

$ {\widehat{\beta }}_{{Y}_{j}}={\theta }_{0E}+{\theta }_{E}{\widehat{\beta }}_{{X}_{j}}+{\epsilon }_{{E}_{j}}, {\epsilon }_{{E}_{j}} \sim N\left(0, {\sigma }^{2}{\sigma }_{{\widehat{\beta }}_{{Y}_{j}}}^{2}\right) $

其中，权重为$ {\omega }_{j} $，$ {\omega }_{j}=\sigma {\left({\widehat{\beta }}_{{Y}_{j}}\right)}^{-2} $；$ {\epsilon }_{{E}_{j}} $为残差项；$ {\theta }_{0E} $为截距项，$ {\widehat{\theta }}_{0E} $表示J个工具变量平均多效性效应的估计值；$ {\theta }_{E} $为斜率，$ {\widehat{\theta }}_{E} $表示经过偏倚校正后的真实因果效应$ \theta $的估计值，等价于X-G和Y-G关联效应值协方差与X-G关联效应值方差之比：

$ {\widehat{\theta }}_{E}=\mathrm{c}\mathrm{o}\mathrm{v}\left({\widehat{\beta }}_{{Y}_{j}}, {\widehat{\beta }}_{{X}_{j}}\right)/\mathrm{v}\mathrm{a}\mathrm{r}\left({\widehat{\beta }}_{{X}_{j}}\right) $

当InSIDE假设和NOME假设均满足时，$ \theta $是因果效应的无偏估计^[6]：

$ \mathrm{v}\mathrm{a}\mathrm{r}\left({\widehat{\beta }}_{{X}_{j}}\right)=\mathrm{v}\mathrm{a}\mathrm{r}\left({\beta }_{{X}_{j}}\right)={\sigma }_{{\beta }_{{X}_{j}}}^{2} $

$ E\left[\mathrm{c}\mathrm{o}\mathrm{v}\left({\widehat{\beta }}_{{Y}_{j}, }{\widehat{\beta }}_{{X}_{j}}\right)\right]=E\left[\mathrm{c}\mathrm{o}\mathrm{v}\left({\widehat{\beta }}_{{Y}_{j}, }{\beta }_{{X}_{j}}\right)\right]=\theta \mathrm{v}\mathrm{a}\mathrm{r}\left({\beta }_{{X}_{j}}\right)=\theta {\sigma }_{{\beta }_{{X}_{j}}}^{2} $

$ E\left[{\widehat{\theta }}_{E}\right]=\theta {\sigma }_{{\beta }_{{X}_{j}}}^{2}/{\sigma }_{{\beta }_{{X}_{j}}}^{2}=\theta $

若满足InSIDE假设但不满足NOME假设，因果效应估计值则会有偏：

$ {\widehat{\theta }}_{E}\approx \theta \frac{\mathrm{v}\mathrm{a}\mathrm{r}\left({\beta }_{{X}_{j}}\right)}{\mathrm{v}\mathrm{a}\mathrm{r}\left({\widehat{\beta }}_{{X}_{j}}\right)}=\theta \frac{{\sigma }_{{\beta }_{{X}_{j}}}^{2}}{{\sigma }_{{\beta }_{{X}_{j}}}^{2}+{s}^{2}} $

其中，$ {s}^{2} $为测量误差，$ {s}^{2}=0 $时才满足NOME假设；$ {I}^{2} $统计量是在Meta分析中评估研究间异质性的工具，$ {\sigma }_{{\beta }_{{X}_{j}}}^{2}/\left({\sigma }_{{\beta }_{{X}_{j}}}^{2}+{s}^{2}\right) $的大小可用$ {I}_{GX}^{2}=\left({Q}_{GX}-\left(J-1\right)\right)/{Q}_{GX} $来估计，$ {I}_{GX}^{2} $代表SNP-暴露关联的真实方差，取值范围为0~1，$ {I}_{GX}^{2} $越接近1，MR-Egger回归的估计越接近真实因果效应。其中，Cochran’s Q统计量是评估工具变量间异质性的工具，$ {Q}_{GX}=\frac{{\sum }_{j=1}^{J}{\left({\widehat{\beta }}_{{X}_{j}}-{\overline{\widehat{\beta }}}_{X}\right)}^{2}}{{\sigma }_{{\beta }_{{X}_{j}}}^{2}} $，$ {\overline{\widehat{\beta }}}_{X} $是SNP-暴露关联的平均值。

IVW法是由Burgess等^[2]于2013年提出的一种加权线性回归模型。与MR-Egger回归不同，该方法限制回归直线通过原点，将两个或多个SNP的效应值进行整合，因果效应估计值为：

$ {\widehat{\theta }}_{ivw}=\frac{{\sum }_{j=1}^{J}V{\left({\widehat{\beta }}_{{Y}_{j}}\right)}^{-1}{\widehat{\beta }}_{{Y}_{j}}{\widehat{\beta }}_{{X}_{j}}}{{\sum }_{j=1}^{J}V{\left({\widehat{\beta }}_{{Y}_{j}}\right)}^{-1}{\left({\widehat{\beta }}_{{X}_{j}}\right)}^{2}} $

实例分析  HDL-C资料来自全球脂质遗传学（Global Lipids Genetics Consortium，GLGC）数据库（http://csg.sph.umich.edu/willer/public/lipids2013/），该数据库中有关基因位点与HDL-C的GWAS研究结果发表于2013年，该研究样本量为94 595人，分析了2 418 527个与HDL-C相关的SNP位点^[7]。与HDL-C相关的SNP筛选标准：（1）达到全基因组统计显著性水平，即P < 5×10^-8；（2）去除连锁不平衡（linkage disequilibrium，LD）的SNPs，排除标准为r² > 0.01，$ kb=5\mathrm{ }000 $，筛选出125个与HDL-C相关且相互独立的SNP，提取每个SNP的rs编号、效应等位基因（effect allele，EA）和次要等位基因（non-effect allele，non-EA）、β系数、P值和标准误。CAD数据来源于2018年Van Der Harst等^[8]发表的GWAS统计结果，该研究将CARDIoGRAMplusC4D[Coronary Artery Disease Genome wide Replication and Meta-analysis（CARDIoGRAM）plus The Coronary Artery Disease（C4D）Genetics]的数据与英国生物银行（UK Biobank，UKBB）的CAD数据进行GWAS分析，共包括547 216人，其中病例组122 733人，对照组424 528人，分析了69 033个与CAD相关的SNP位点。设置最小等位基因频率MAF=0.01，从GLGC数据库筛选的SNP在CAD来源的GWAS研究中均存在，经等位基因对齐后，rs11637365、rs3790106、rs4986970、rs6589581和rs964184由于存在回文结构被剔除，最终有120个SNP作为工具变量分析HDL-C与CAD的因果关联。上述研究人群均为欧洲血统，具有相似的年龄和性别，符合两样本孟德尔随机化分析的要求。利用R4.0.3中的TwoSampleMR包进行两样本孟德尔随机化分析，检验水准α=0.05。因果效应分析使用IVW法和MR-Egger回归，对比两者差异；敏感性分析包括利用MR-Egger回归进行多效性检验（pleiotropy test）和逐个剔除检验（leave-one-out sensitivity test），“leave-one-out”是指逐一剔除SNP，分别计算剩下SNP的合并效应，如果剔除某个SNP后其他SNP估计出来的MR结果和总结果相差很大，说明MR结果对该SNP是敏感的。

结果

本研究针对HDL-C和CAD选取有效SNP作为Ⅳ，经筛选确定了120个SNP作为Ⅳ，表 1列出了前10个SNP的相关信息。采用MR-Egger回归和IVW法估计因果效应（表 2）。IVW结果支持HDL-C和CAD之间存在因果关系（OR=0.82，95%CI：0.75~0.89），且HDL-C每增加1个标准差，CAD风险降低18%。而MR-Egger法不支持（OR=0.96，95%CI：0.83~1.11），即HDL-C对CAD的发生风险没有影响。多效性检验结果显示，MR-Egger回归的截距为-0.01，P=0.008，说明因果分析结果会受到水平多效性的影响。在水平多效性存在的情况下，传统孟德尔随机化方法IVW的估计出现了偏差，MR-Egger回归则能识别并校正多效性，给出接近真实值的因果估计。敏感性分析结果见图 1，所有线条均在0的左侧，说明无论去除哪个SNP都不会对结果产生根本影响，即孟德尔随机化的结果稳健。

讨论

早期孟德尔随机化研究倾向于使用单一遗传变异，并集中在单一研究人群中的特定危险因素与疾病的关联研究。近年GWAS中大量增长的基因型-表型关联导致众多暴露因素的大量遗传变异被识别出来^[9]。许多遗传变异被认为具有多效性效应，要证明排他性假设的成立几乎不可能。MR-Egger回归区别于传统孟德尔随机化方法之处在于放宽了SNP之间不存在水平多效性的要求，它假设在基因-暴露关联和基因变异对结果的直接影响之间没有相关性，是比更严格的排他性标准更弱的要求，即使所有工具变量是无效的，MR-Egger回归也能得到准确的因果效应估计^[10]。因此，在工具变量存在多效性的情况下，MR-Egger回归是一个更好的选择。

本研究借助GLGC数据库和Van Der Harst的GWAS统计结果，提取与HDL-C和CAD均密切相关且相互独立的SNP作为工具变量，利用MR-Egger回归和IVW法进行两样本孟德尔随机化分析，估计HDL-C和CAD的因果关联并进行比较。IVW结果显示两者之间具有因果关系，HDL-C对于CAD有保护作用，但IVW法假设截距为0，未考虑所纳入的工具变量的基因多效性；而MR-Egger回归并不局限于截距为0，是在IVW基础上的修正，同时考虑基因多效性对结果的影响，因此给出了并不一致的估计结果，即HDL-C与CAD之间不存在因果关联，与以往孟德尔随机化研究结果相吻合^[11]。随后的逐一剔除检验显示MR-Egger结果更具有稳健性。因此，MR-Egger回归在孟德尔随机化分析中有较好的应用价值。

与传统孟德尔随机化方法相比，MR-Egger回归的偏倚更小，而且因果无效假设的排斥率更接近名义上的5%^[3]。MR-Egger回归的局限性是它在因果估计方面的统计效力较低。如果没有证据表明IVW和MR-Egger之间存在系统性差异，那么通常采用IVW，虽然IVW可能不那么稳健，但在此情况下，IVW估计值的精度更高。有研究者指出MR-Egger回归应被视为一种敏感性分析，用来检验是否违反工具变量假设，而不是作为上述更常见方法的替代^[12]。实际应用时应结合多种方法的结果进行综合考量，如果多种方法的结果一致，则结论更为可靠，从而更好地应用于传统观察性流行病学研究，为病因推断提供更强有力的证据，并有利于为未来疾病预防的干预目标设定优先级。

作者贡献声明  徐艺耘  查阅文献，论文构思、撰写和修订，数据统计分析。刘振球，樊虹，张欣，施婷婷，吴声  结果解释，论文修改。张铁军  论文构思、指导和修订。

利益冲突声明  所有作者均声明不存在利益冲突。