CN114171110B - 一种基于联合似然的孟德尔随机化分析方法 - Google Patents

Abstract

本发明涉及疾病分析领域，具体提供了一种基于联合似然的孟德尔随机化分析方法，其可以同时对一组相互之间处于高LD(linkage disequilibrium，连锁不平衡)的初始候选SNPs(Single Nucleotide Polymorphism,单核苷酸多态性)工具变量进行建模，从中自动筛选合适的工具变量进行孟德尔随机化(Mendelian randomization,MR)分析,对在生物遗传中普遍存在的两种类型的水平多效性进行判别，并进行良好的控制，另外结合可扩展的抽样算法可以给出校准的p值，发明人将该方法称为工具变量自动筛选的两样本孟德尔随机化(MRAID)。发明人通过全面的模拟来证明MRAID的有效性，并且可以用于生物库大规模数据的实证分析，计算效率高，因此具有良好的实际应用价值。

Description

一种基于联合似然的孟德尔随机化分析方法

技术领域

本发明涉及疾病分析领域，具体而言，涉及一种基于联合似然的孟德尔随机化分析方法。

背景技术

研究复杂性状之间的因果关系并确定因果风险因素是理解复杂疾病的重要一步。孟德尔随机化(Mendelian randomization,MR)则是一种在观察性研究中进行因果推断常用的统计工具。MR是一种使用SNP作为工具变量推断暴露对结局因果效应的工具变量分析方法。MR只需要来自全基因组关联研究(GWASs)的汇总统计量，其通常在两样本研究中进行，暴露变量和结局变量是在两项独立研究中进行测量的。由于GWAS汇总数据的丰富可用性，大量用以确定各种常见疾病因果风险因素的MR分析不断涌现。许多最近发展的MR方法促进了这些MR研究，包括逆方差加权(IVW)方法，MR-Egger，基于中位数的MR，BWMR，RAPS，MRMix，CAUSE等等。不同的MR方法在建模假设和推断算法方面有所不同，但大多数方法都会遇到如下两个重要的建模和算法挑战，限制了其分析的有效性。

首先，几乎所有现有的MR方法都依赖于预先选择的一组独立SNP作为工具变量。选择的这些工具变量要求彼此独立，以确保在许多常见MR方法(如IVW)中使用的统计推断是有效的。独立的SNP通常是通过连锁不平衡聚集(LD clumping)来选择的，该过程首先根据SNP与暴露变量的边际关联证据对SNP进行排名，然后保留排名列表中低LD的SNP。但是，使用LD-clumping选择SNP可能不是最佳选择，因为所选SNP可能仅仅代表的是与因果SNP有LD关系的标记SNP(Tag SNP)，而不是因果SNP本身。使用标记SNP而不是因果SNP作为工具变量会降低MR分析的功效。另外，也许更重要的是，选择独立的SNP进行MR分析可能也不是理想的方法，因为复杂的性状可能会受到驻留在同一局部区域中的多个彼此存在LD的因果SNP的影响。因此，选择独立的SNP可能只捕获暴露变量中一小部分表型变异，再次导致随后的MR分析损失功效。事实上，在平行的研究领域全转录组关联分析中，已有文献证明，与仅仅使用独立SNP相比，结合关联SNP可以显着提高检验效能。因此，结合关联SNP来开发工具变量选择的有效方法对于充分发挥MR的潜力非常重要。

其次，只有少数MR方法对水平多效性进行建模，且难以有效控制水平多效性。当工具变量SNP通过暴露以外的途径对结局产生影响时就会发生水平多效性。在复杂性状分析中可以观察到大量水平多效性，并且通常有两种不同的类型。第一种水平多效性是通过独立于暴露的路径产生的，由此产生的水平多效性效应与SNP对暴露的效应无关。第二种类型的水平多效性是通过未观察到的暴露-结局混杂因素产生的，并导致水平多效性效应与SNP对暴露效应之间的关联。两种水平多效性都违反了标准的MR建模假设，并可能导致因果效应估计产生偏倚和假阳性的增加。早期MR分析通过简单地去除可能与结局变量相关的SNP工具变量来控制水平多效性。去除与结局相关的SNP将倾向选择一组保守的工具变量并造成随后MR分析中的功效损失。最近的MR方法通过指定不同水平多效性的建模假设清晰地将水平多效性纳入了模型中。例如，Egger假设SNP工具变量具有相同的水平多效性效应，而PMR-VC和BWMR假设水平多效性遵循正态分布，所有这些方法都是对第一类水平多效性进行了建模。相比之下，MRMix和CAUSE采用正态混合模型来控制两种类型的水平多效性。不幸的是，对这两种类型的水平多效性进行建模在技术上一直具有挑战性，因为由此产生的MR模型的似然函数通常包含一个无法解析求解的积分。因此，MRMix和CAUSE都基于非似然的方法来进行MR推断。具体地说，MRMix在因果效应候选网格上进行搜索，以确定在没有水平多效性的情况下最大限度地提高GWAS汇总统计量在预期子模型中的比例。CAUSE则是通过计算两个不同模型之间的期望对数后验密度，比较两个不同模型(一个有因果效应，另一个没有因果效应)的样本外预测精度，来进行因果推断。然而，基于非似然性的因果推断可能导致功效损失和(或)未经校准的检验统计量，而该统计量对于大规模筛选潜在疾病的因果风险因素至关重要。事实上，正如发明人将展示的那样，MRMix对于工具变量效应所服从的分布并不稳健，容易出现估计偏差，而CAUSE则会产生过于保守的p值。

发明内容

本发明的主要目的在于提供一种基于联合似然的孟德尔随机化分析方法，以解决相关技术中的问题。

为了实现上述目的，根据本发明的一个方面，提供了一种基于联合似然的孟德尔随机化分析方法，包括：

对一组互相之间处于高LD的初始候选SNP工具变量进行建模，并从中自动筛选合适的工具变量进行MR分析；

同时有效控制与工具变量和暴露之间的效应相关的或不相关的两种水平多效性效应；

采用可扩展的抽样算法来计算校准的p值。

进一步地，还包括在两样本MR框架下估计和检验暴露变量对结局变量的因果效应，其中暴露变量和结局变量是在两个独立的没有样本重叠GWAS中测量。

进一步地，所述重叠GWAS分别称为暴露GWAS和结局GWAS，在暴露GWAS中，发明人按照边际p值低于全基因组显著性阈值的标准进行初步筛选，选出与暴露变量相关的SNPs。

进一步地，对一组互相之间处于高LD的初始候选SNP工具变量进行建模包括通过以下三个线性回归对暴露，结局和基因型之间的关系进行建模：

x＝Z_xβ+ε_x,

其中，Z_x是一个n₁×p的基因型矩阵，表示的是在暴露GWAS中n₁个个体p个候选工具变量信息，Z_y是一个n₂×p的基因型矩阵，表示的是在结局GWAS中n₂个个体p个候选工具变量信息，y为结局变量，

为未观察到的暴露，β是p维向量，表示的是SNPs对暴露效应，η₀和η₁是表示水平多效性效应的p维向量，α代表的是暴露对结局的因果效应标量，ε_x是表示残差的n₁维向量，其每个元素独立同分布于正态分布

是表示残差的n₂维向量，其每个元素独立且服从相同的正态分布

ε_y是表示残差的n₂维向量，其每个元素独立且服从正态分布

进一步地，所述自动筛选合适的工具变量包括对SNP对暴露的影响效应(β)进行了稀疏性建模假设来进行工具变量自动筛选，其具体为：提出假设

其中δ₀是表示在零处的质点的狄拉克函数，第j个SNP对暴露的效应有1―π_β的概率为零，有π_β的概率效应不为零，若存在非零效应，其效应服从均值为0，方差为

的正态分布，其中方差参数

决定效应大小。

进一步地，还包括在η₀和η₁上指定了单独的模型假设，其具体为：假设每个选定的SNP工具变量都有π_c的概率通过混杂因素导致相关水平多效性，且SNP对结局的混杂效应是其对暴露效应的ρ倍，如果SNP对暴露的效应是β_j，选定的SNP通过混杂对结局的效应变成ρβ_j，则对

的假设，是

假设每个选定的SNP工具都有π₁的概率不通过暴露直接对结局表现出水平多效性，则对

的假设为

故对于被选定的工具变量其水平多效性为

而对于未被选定的非工具变量，假设

使用相同的方差参数

来对工具变量和非工具变量SNP的不相关水平多效性效应建模。

进一步地，所述采用可扩展的抽样算法来计算校准的p值包括：利用Gibbs抽样获得因果关系α的后验样本，利用该样本均值和样本标准差去归纳后验分布。

根据权利要求1所述的基于联合似然的孟德尔随机化分析方法，其特征在于，所述采用可扩展的抽样算法来计算校准的p值包括：给因果关系α指定一个正态的先验分布，获得其先验均值和标准差，利用矩估计的方法来获得近似最大似然估计

和标准误

构建一个近似的Wald检验统计量，获得假设检验的p值。

与现有技术相比，本发明具有以下有益效果：本发明表明MRAID在存在水平多效性(包括相关的和不相关的)的情况下为因果效应检验提供了校准的I类错误控制，降低了假阳性，比现有的MR方法更具效力。另外，作为一份副产品，MRAID还可以估计产生不相关或相关的水平多效性的SNP的比例。通过对英国生物银行(UK Biobank)中的645个性状对进行MR筛查分析，说明了MRAID的优点，确定了可能对心血管疾病相关性状有因果效应的多种生活方式风险因素。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其它的附图。

本说明书附图所绘示的结构、比例、大小等，均仅用以配合说明书所揭示的内容，以供熟悉此技术的人士了解与阅读，并非用以限定本发明可实施的限定条件，故不具技术上的实质意义，任何结构的修饰、比例关系的改变或大小的调整，在不影响本发明所能产生的功效及所能达成的目的下，均应仍落在本发明所揭示的技术内容得能涵盖的范围内。

图1为本发明实施例MRAID的原理图；

图2为本发明实施例不同MR方法在模拟中的一类错误控制率；

图3为本发明实施例不同MR方法在模拟中的功效；

图4为本发明实施例实际数据中不同MR方法在自身性状对分析中的点估计和95％置信区间；

图5为本发明实施例不同MR方法在检验生活方式风险因素和英国生物银行心血管疾病相关特征之间的因果关系的―log₁₀P值图；

具体实施方式

为使得本发明的发明目的、特征、优点能够更加的明显和易懂，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，下面所描述的实施例仅仅是本发明一部分实施例，而非全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

下面结合附图并通过具体实施方式来进一步说明本发明的技术方案。

个体数据的MRAID

发明人在此概述了MRAID，其推论和技术细节在图1的说明性图表中提供。发明人的目标是在两样本MR框架下估计和检验暴露变量对结局变量的因果效应，其中暴露变量和结局变量是在两个独立的没有样本重叠GWAS中测量的。发明人把这两个独立的GWAS分别称为暴露GWAS和结局GWAS。x是一个n₁维向量，表示在暴露GWAS中n₁个个体的暴露。y是一个n₂维向量，表示在结局GWAS中n₂个个体的结局。发明人将x和y标准化，使其均值为0，标准差为1。在暴露GWAS中，发明人按照边际p值低于全基因组显著性阈值(5×10^―8)的标准进行初步筛选，选出与暴露变量相关的SNPs。这些SNPs可能由于LD存在相互关联，会被选进候选工具变量的初始集合。Z_x是一个n₁×p的基因型矩阵，表示的是在暴露GWAS中n₁个个体p个候选工具变量信息，Z_y是一个n₂×p的基因型矩阵，表示的是在结局GWAS中n₂个个体p个候选工具变量信息。发明人对两个基因型矩阵的每一列进行标准化，使其均值为0，标准差为1。发明人通过以下三个线性回归对暴露，结局和基因型之间的关系进行建模：

x＝Z_xβ+ε_x, (1)

以上，公式(1)描述了暴露GWAS中基因型Z_x和暴露变量x之间的关系；公式(2)描述了结局GWAS中基因型Z_y与未观察到的暴露

之间的关系；公式(3)描述了结局GWAS中基因型Z_y、结局变量y和未观察到的暴露

之间的关系；β是p维向量，表示的是SNPs对暴露效应，η₀和η₁是表示水平多效性效应的p维向量，α代表的是暴露对结局的因果效应标量，ε_x是表示残差的n₁维向量，其每个元素独立同分布于正态分布

是表示残差的n₂维向量，其每个元素独立且服从于相同的正态分布

ε_y是表示残差的n₂维向量，其每个元素独立且服从正态分布

发明人注意到，虽然上述三个方程是基于两个独立的GWAS指定的，但它们通过共同参数β联系在一起。发明人下面仔细考虑了关于SNP对暴露变量的影响效应β以及水平多效性η₀和η₁的建模假设。

上述模型中包含的p个SNPs代表一组初始候选工具变量。虽然所有候选工具变量都与暴露有一定的关联，但其中大部分不太可能是暴露变量的因果SNP。相反，大多数候选工具变量很可能代表的是因为与暴露真正的因果SNP有LD关联而与暴露相关的标记SNPs。因此，最好是对候选工具变量进行额外筛选以识别暴露的因果SNP，并将其作为工具变量以最大限度提高MR分析的功效。为此，发明人借鉴了GWAS研究领域的精细定位思想，对SNPs对暴露的影响效应(β)进行了稀疏性建模假设来进行工具变量自动筛选。特别地，发明人提出假设

其中δ₀是表示在零处的质点的狄拉克函数。也就是说，第j个SNP对暴露的效应有1―π_β的概率为零，有π_β的概率效应不为零，若存在非零效应，其效应服从均值为0，方差为

的正态分布，其中方差参数

决定效应大小。对β的稀疏假设允许发明人筛选对暴露具有非零效应的SNP作为MR模型中的工具变量。

此外，上述模型中包含的p个SNPs也可能存在水平多效性，并通过暴露以外的路径影响结局变量。为了控制潜在的水平多效性效应并改善因果效应推断，发明人引入了两组参数η₀和η₁来对水平多效性效应进行建模。两组参数分别放置在两个SNP组(选定的工具变量SNP组和未选定的非工具变量SNP组)中，这些参数通过对β的稀疏建模假设进行分类。特别地，η₁表示β_j≠0的被选定SNP工具变量所表现出的水平多效性，而η₀表示具有β_j＝0的未被选择的非工具变量SNP所表现出的水平多效性。控制η₁可以帮助降低由工具变量SNP的水平多效性引起的因果效应估计的偏倚,而控制η₀可以减少公式(3)中的残差方差，从而有助于提高因果效应估计的统计效率。

为了有效控制两个SNP组的水平多效性，发明人在η₀和η₁上指定了单独的模型假设。具体来说，对于选定的SNP，发明人假设它们可以以两种不同的方式产生水平多效性：它们可以通过与暴露和结局相关的常见混杂影响结局，并且它们可以通过独立于暴露的路径影响结局。对于第一种水平多效性，发明人假设每个选定的SNP工具变量都有π_c的概率通过混杂因素导致水平多效性。参考先前的研究，发明人假设SNP对结局的混杂效应是其对暴露效应的ρ倍。因此，如果SNP对暴露的效应是β_j，则选定的SNP通过混杂对结局的效应变成ρβ_j，发明人对

的假设，是

对于第j个SNP，

作为η₁的一部分代表了第一种类型的水平多效性，其中I(·)是一个将水平多效性设置为ρβ_j指示函数。对于第二种类型的水平多效性，发明人假设每个选定的SNP工具都有π₁的概率不通过暴露直接对结局表现出水平多效性。对于第j个SNP，发明人用

代表η₁中代表第二类水平多效性效应的部分，针对其的假设为

其中方差

决定了水平多效性效应的大小。注意，第一类水平多效性因为经过混杂因素，所以其与SNP对暴露的效应相关，而第二类水平多效性与SNP对暴露的效应无关。整个水平多效性是两者的总和，也就是

当然，η₁是被选择的SNP具备的水平多效性效应，所以如果β_j＝0，则η_1j＝0。对于未被选择的β_j＝0的非工具变量SNP，发明人假设

也就是，非工具变量SNP有π₀的概率产生水平多效性，其效应服从与上述水平多效性相同方差

的正态分布。因为发明人通常没有足够数量的SNP来准确估计两个单独的参数，所以发明人使用相同的方差参数

来对工具变量和非工具变量SNP的不相关水平多效性效应建模。因为η₀是非工具变量的水平多效性效应，则当β_j≠0，也会有η_0j＝0。

上述水平多效性的参数是基于SNP工具变量选择建模的。水平多效性的等价和替代建模思路是将其分割成相关的水平多效性部分η_c和不相关的水平多效性部分η_u。具体来说，相关水平多效性仅在选定的SNP工具变量中发生，其效应服从η_cj|β_j≠0～π_cI(η_1j＝ρβ_j)+(1―π_c)δ₀，如果β_j＝0，那么η_cj＝0。另一方面，不相关的水平多效性在工具变量和非工具变量SNP均可存在，其分别服从如下分布

结合第一种建模思路，也就是，

和

总体而言，SNP对结局的影响通过三种不同的路径实现：通过暴露对结局的因果效应α，通过由未观测混杂介导的相关水平多效性效应和通过不相关的水平多效性效应。模型中的SNP可以不通过任何一种途径，也可以通过一种或多种路径对结局产生作用。值得注意的是，除非发明人做出进一步的建模假设，否则SNP对结局通过因果效应路和通过相关水平多效性效应路的影响是无法区分的。依据先前的研究，发明人假设π_c非常小。因此，在选定的对暴露具有非零效应的工具变量SNP中，只有一小部分对结果表现出相关的水平多效性。

发明人感兴趣的关键参数是因果关系α，在标准MR模型中对α的因果解释要求所选的SNP工具变量满足三个条件：(i)工具变量只与暴露相关(相关条件)；(ii)不与和暴露、结局相关的任何其他混杂相关(独立条件)；(iii)工具变量只通过暴露影响结局(排除限制条件)。发明人对β的建模假设使发明人能够选择满足相关性条件的SNPs。发明人对η₀和η₁的建模假设允许在违反独立条件和排除限制条件下建模。因此，发明人的模型用对β，η₀和η₁的具体假设有效地取代了一般条件的(ii)和(iii)。此外，通过利用ρ对工具变量-暴露效应和工具变量-结局效应间相关性的清晰建模，发明人的模型不再需要InSIDE假设，有时也被称作弱排除限制条件。

发明人感兴趣的是估计因果效应α和对零假设H₀:α＝0进行检验。但是，由于基于上述建模假设定义的似然函数非常复杂，而且包含了很多利用解析方法无法求解的积分，所以对α的推断具有很大的计算难度。发明人在极大似然框架下开发一种近似推理算法来实现似然函数的积分运算，并且可以获得用于检验α的近似p值。发明人的算法是基于α的似然函数可以被表示成其后验分布与先验分布的比值进行的。因为α的后验分布是近似正态分布，发明人可以利用Gibbs抽样获得α的后验样本，并且可以利用该样本均值和样本标准差去归纳后验分布。此外，发明人还可以给α指定一个正态的先验分布，获得其先验均值和标准差。因为α的似然可以表示成后验分布除以先验分布的形式，并且其自身是服从渐进正态分布的，则可以基于从后验分布和先验分布中获得的均值和标准差，利用矩估计的方法来获得近似最大似然估计

和标准误

接着，发明人可以构建一个近似的Wald检验统计量，获得假设检验的p值。需要说明的是，该算法虽然依赖于Gibbs抽样，但是发明人没有进行贝叶斯分析，而是把Gibbs抽样作为一种实用准确的数值近似工具，来获得在频率方法下难以获得的α的边际似然。

发明人把该模型和算法称为两样本工具变量自动筛选的MR分析(MRAID)。“工具变量自动筛选”突出了发明人的模型从一组相互间具有潜在高LD的候选变量中选择工具变量的思想。与现有的两样本MR方法相比，MRAID依赖于极大似然的推理框架，能够对相关的工具变量进行建模，进行工具变量的自动选择，控制相关和不相关的水平多效性，并且具备计算的可拓展性。

汇总数据的MRAID

尽管发明人利用个体数据对MRAID进行了描述，但其也可以拓展至仅使用汇总数据的情况。MRAID汇总数据版本的具体细节可以在详细信息中查看。简单来说，MRAID的汇总数据版本需要输入两种信息：SNP对暴露和结局边际效应大小的估计和暴露GWAS、结局GWAS中SNPs的相关矩阵。这两种输入的信息都是基于标准化的基因型数据获得的，其中每个SNP的基因型数据已标准化为零均值和单位标准差。发明人用p维向量

和

分别表示SNPs对暴露的边际效应和SNPs对结局的边际效应，用p×p维的矩阵Σ₁和Σ₂分别表示暴露GWAS和结局GWAS种的SNPs相关矩阵。Σ₁和Σ₂都是半正定的，并且可以从相同的LD参考面板中进行估计(例如，千人基因组中具有相同祖先的个体)。汇总数据的MRAID模型可基于以下两个方程构建

其中e_x，e_y是代表残差的p维向量，分别服从多元正态分布

和

针对汇总数据版本的MRAID也在最大似然框架下开发了类似的近似推断算法。

模拟

发明人进行了基于实际情况的模拟，来评估MRAID的性能，并将其与六种现有的MR方法进行比较。在模拟中，发明人从英国生物银行(UK Biobank)中随机挑选了60000个个体，随机地将这些个体分成等大小的两组样本：一组30000个人作为暴露GWAS，另一组其余30000个人作为结局GWAS。对于这些个体，发明人从1号染色体上的649695个SNPs中获取了他们的基因型，将每个SNP进行标准化使其均值为0和标准差为1，并使用标准化的基因型来模拟暴露和结局。具体来说，在暴露GWAS中，发明人随机选择对暴露具有非零效应的K个SNPs(K＝100或1000)。

表示K个SNPs构成的基因型矩阵。发明人从正态分布

中产生K个SNPs对暴露的效应β，其中标量

表示这些遗传变异解释的暴露变量中的方差比例。将K个SNPs的遗传效应相加，即为

从正态分布

中产生了的残差ε_x，然后，将总遗传效应和残差相加，产生了模拟的暴露变量x。在结局GWAS中，发明人获得了K个相同SNPs的基因型，作为

并使用与暴露GWAS相同的β来计算结局背后的遗传成分，即

发明人将因果效应α设置为

使得因果项

可以解释的结局变量中的方差比例恰为PVE_α。发明人从K个SNPs中随机获得π_cK个SNPs(四舍五入取整数)，使其具有相关水平多效性。相关水平多效性效应大小在模拟中被设置为ρβ，通过设置ρ以便相关水平多效性解释的结局变量中的方差比例为PVE_c。另外，发明人从K个SNPs中随机获得π₁K个SNPs(四舍五入为整数)，随机从剩余的非因果SNP中获取100―π₁K个SNPs，产生不相关的水平多效性，从而总共有100个SNPs具有不相关的水平多效性。发明人从正态分布中产生了这100个SNPs的不相关的水平多效性效应，并将其标准化，使不相关的水平多效性解释的结局表型方差的比例为PVE_u。残差ε_y从正态分布N(0,1―PVE_α-PVE_c―PVE_u)中产生。发明人将因果项、相关和不相关的水平多效性、残差相加，产生了模拟的结局变量y。按照事先无法确定因果SNPs的思路，发明人遵循标准MR步骤，选择合适的SNP作为工具变量。为此，发明人在GEMMA软件中利用线性回归模型在暴露GWAS中进行关联分析，并选定了p值低于5×10^―8的SNPs作为候选工具变量。对于被选定的SNPs，发明人获得了它们的效应估计、标准误和Z得分统计量，来用作汇总统计量输入。Z_x和Z_y分别表示暴露和结局GWAS中被选定的SNPs的标准化基因型矩阵。基于基因型矩阵，发明人获得了SNPs的相关矩阵，如

和

作为MR模型拟合的输入。

在模拟中，发明人首先进行了基线模拟设置，其参数设定为

PVE_α＝0,K＝100,PVE_u＝0,PVE_c＝0。在此基础上，发明人一次改变一个参数，以检查各种参数对方法性能的影响。对于

发明人将其设置为5％或10％，对于β，除了让其服从正态分布外，发明人还从贝叶斯稀疏线性混合模型(BSLMM)分布中产生。具体来说，发明人在K个SNPs中随机选择了1％或10％的SNPs有较大效应，这些大效应的SNP解释

的20％，剩余的SNP具有较小的效应，来解释其余的

对于K，发明人将其设置为100或1000。对于PVE_α，发明人将它在零模拟中设置为零，并在非零模拟中测试了不同值，当K＝100时，发明人设置了PVE_α为0.05％、0.15％和0.25％，当K＝1000时为0.5％、1.5％和2.5％，以确保充足的功效。对于不相关的水平多效性效应，发明人设置PVE_u为0、2.5％和5％。在没有不相关的水平多效性(PVE_u＝0)的情况下，进行零模拟(PVE_α＝0)时，发明人将K设置为100或1000。当存在不相关的水平多效性时，发明人设置K为100，并设置π₁为0，10％，20％和30％。发明人还模拟了相关的水平多效性，并设置了π_c为5％和10％，ρ为0.02和0.05。

对于零模拟，发明人在每个场景中重复了1000次模拟，来检验方法对I类错误的控制。在功效评价方面，发明人进行了100次非零模拟以及900次零模拟，通过这些模拟，发明人根据0.05的错误发现率(FDR)计算了功效。然后，发明人重复了五次这样的分析，并报告了这些重复分析的平均功效。需要说明的是，因为不同方法的相同p值可能对应于不同的I类错误，为了保证方法间比较的公平性，发明人根据FDR计算了功效而不是名义的p值阈值。

实际数据分析应用

发明人利用MRAID和其他MR方法识别了英国生物银行(UK Biobank)中38种生活方式风险因素与11种与心血管疾病相关性状之间的因果关联。英国生物银行由487298人和92693895个被填补的SNPs组成。

发明人在Nealelab遵循了相同的质量控制(QC)程序，共保留337129名欧洲血统的个体进行分析。发明人还筛选出HWE p值<10^―7，基因型调用率<95％或MAF<0.001的SNPs，以获得总计13876958个SNPs进行分析。对于被纳入分析的个体，发明人获得了其所有与生活方式相关的定量性状以及心血管疾病相关性状，删除了样本含量小于10000的特征，并专注于剩余的38个生活方式性状和11个心血管疾病相关性状进行分析。这38种生活方式性状包括8种身体活动性状、12种酒精摄入性状、10种饮食性状(如咖啡和水果)和8种与吸烟相关的性状。11个与心血管疾病相关性状包括四个脉冲波性状、两个血压性状(SBP和DBP)、四个脂质性状(LDL、HDL、TC和TG)和BMI。在观察性研究中，这些生活方式风险因素许多被发现与心血管疾病相关性状有关，但尚不清楚这些关联是否为因果关系。对于每个性状，依次去除其性别和前10个基因型主成分(PCs)的影响，得到性状残差，将残差标准化为均值为0，标准差为1，并将这些标准化后的残差用于MR分析。

为模拟两样本MR设计，发明人将337129个个体分成两个不重叠的集合：一个是包含168564个个体的暴露GWAS集合，另一个是包含168565个个体的结局GWAS集合。随机数据分割策略在之前的MR文献中被广泛应用，其确保了每个研究中样本的同质性和研究之间的独立性。发明人研究了在暴露GWAS中的38个生活方式性状和结局GWAS中11个与心血管疾病相关的性状。在这两个GWAS中，发明人在GEMMA软件中通过线性回归获得了每个性状的汇总统计量。当暴露GWAS中的生活方式性状作为暴露时，发明人选择了p值低于5×10^―8的SNPs作为每个暴露性状的候选工具。因为大多数MR方法要求至少两个工具SNPs，且当工具SNPs数量太大时，一些方法同样会变得不稳定，所以发明人去除了候选工具变量数量低于2或超过10000的暴露性状。这样，发明人通过少于2个候选工具变量的条件去除了三个性状，超过10000个候选工具变量的条件去除了四个性状。发明人将剩下的31个生活方式性状暴露和11个心血管疾病相关性状结局配对成341个性状对。31个性状间显著相关的SNPs平均数为286个。当结局GWAS中与心血管疾病相关的性状为暴露时，发明人去除了3个候选工具变量数少于2个的性状，发明人发现剩余8个性状的候选工具变量数均超过10000个。因此，对于这些剩余的性状，发明人使用一个更严格的1×10^―15的p值来选择工具SNPs，并且分析得到的304个性状对。8个与心血管疾病相关的性状中相关SNPs的平均数量是2,318。总计，发明人分析了645个性状对。

比较方法

发明人将MRAID的性能与以下现有的六种方法进行了比较：(1)IVW-R，IVW的随机效应版本，它通过加权和整合单个工具的效应估计来获得因果效应估计。它依赖随机效应来解释多效性和工具变量效应估计的异质性。(2)Weighted mode，IVW的众数版本。它将从单个工具获得的效应估计的众数，而不是平均值作为因果效应估计。(3)Robust，这是IVW的稳健版本。它使用MM-估计程序，包括初始S-估计，和进一步的与Tukey的双权重损失函数相结合的M-估计。发明人使用R包“MendelianRandomization”和默认设置来拟合方法(1)-(3)。(4)RAPS，这是一种调整轮廓得分的MR方法。它将随机效应和稳健损失函数结合到轮廓评分中，以解释系统的和特殊的基因多效性。发明人用R包“mr.raps”拟合RAPS。(5)MRMix，是一种依赖于混合模型来解释水平多效性效应和其与工具变量效应相关性的方法。发明人使用R包“MRMix”拟合MRMix。(6)CAUSE，其在解释相关的基因多效性时，识别了与因果效应模式一致的工具变量效应。发明人用R包“cause”拟合了CAUSE。发明人将MRAID和上述的六种方法进行了比较。在满足InSIDE假设的情况下IVW-R,Robust和RAPS都展现了良好的性能，而即使违反InSIDE假设，MRMix和Weighted mode都表现的很好。在模拟和实际数据应用中，发明人首先获得了达到全基因组显著性水平(p<5×10^―8)的SNPs，作为工具变量的候选集。发明人直接利用这个候选工具变量集进行MRAID分析。由于所有其他MR方法都需要独立的SNPs作为工具变量，因此发明人对候选工具变量集的SNPs进行了LD-clumping来选择独立的SNPs进行分析。LD-clumping是使用PLINK进行的，其中发明人将LD的r²参数设置为0.001。CAUSE还需要通过使用基因组中一组随机的SNPs来估计模型中一些冗余参数，发明人通过随机选择100000个SNPs来做到这一点。最后，发明人在功效模拟中探索了oracle方法，在该方法中发明人知道真实影响暴露的因果工具变量，获得了真实工具SNPs集合，通过剪枝从中选出独立的SNPs，并使用选定的SNPs集作为工具变量，利用IVW-R方法进行分析。

结果

方法综述

方法部分中描述了MRAID，其技术细节在详细信息中，方法示意图如图1所示。简单的说，MRAID是一个两样本MR方法，旨在使用GWAS汇总统计量推断暴露变量对结局变量的因果效应。MRAID对可能具有潜在LD的全基因组显著的SNPs联合建模，同时可以进行自动的工具变量筛选识别出适合进行MR分析的工具变量。此外，MRAID通过基于极大似然推理框架清晰地说明了两类水平多效性效应，并可扩展到生物库数据集进行运算(表1)。

表1各种两样本MR方法的平均计算时间(分钟)

计算是在Xeon Gold 6138 CPU的单个线程上执行的。计算了20次重复实验所耗时间的平均值，括号内的值表示标准差。#SNPs表示模型中包含的工具变量的数量。MRAID的计算时间是基于1000次Gibbs迭代采样进行的，其中丢掉了前200次。

模拟：I类错误控制

发明人在模拟中评估了MRAID的性能，并将其与六种现有的MR方法进行了比较(详见详细信息部分)。首先，发明人研究了不同情形下不同方法对I类错误的控制。在相关和不相关水平多效性效应都不存在的情况下，大多数方法，包括MRAID、IVW-R、Robust、RAPS和MRMix，都有合理的I类错误控制率(图2a)。另一方面，Weighted mode和CAUSE则对I类错误控制过于保守，这与前期研究结果是一致的。不论影响暴露的SNPs的数量还是它们对暴露的总效应怎样，在零模拟中不同方法获得的p值分布十分相似。发明人进一步研究了不同方法在不同情形下的稳健性，此时，设定SNPs对暴露的效应不再服从简单的正态分布，而是其中有一些SNPs具有比其他SNPs更大的效应。在这种设置下，MRAID、IVW-R和RAPS仍旧表现良好，而MRMix和Robust方法则显示出膨胀的I类错误，这可能是由于它们要求SNPs效应大小服从正态分布所致(图2b)。值得注意的是，发明人直接使用相关的SNPs进行MRAID分析，其他方法则使用的是通过clumping选择的独立SNPs。在没有clumping的情况下，其他所有MR方法都会产生过度膨胀的I类错误(图2)。

发明人研究了水平多效性效应对不同方法控制Ⅰ类错误的影响。当水平多效性效应存在，但与工具对暴露的效应不相关时，MRAID可以控制I类错误(图2c)。相反，Weightedmode和CAUSE都过于保守，而MRMix、Robust、IVW-R和RAPS产生了膨胀的p值(图2c)。无论不相关水平多效性的效应大小或具有不相关水平多效性效应的SNPs的比例如何，类似的结论都成立。当具有不相关水平多效性效应的SNPs比例减少时，MRMix的p值膨胀问题得到了缓解。当不相关的水平多效性效应和相关的水平多效性效应都存在时，MRAID保持有效的I类错误控制率(图2d)。相比之下，Weighted mode和CAUSE都过于保守，而MRMix、Robust、IVW-R和RAPS产生了膨胀的p值。无论相关水平多效性的效应多大，具有不相关水平多效性效应的SNPs的比例是多少，显示相关水平多效性效应的SNPs的比例是多少，相关水平多效性效应是如何产生的，类似的结论都成立。

模拟：功效比较

发明人检验了不同MR方法识别非零因果效应的功效。因为来自不同方法的相同p值可能对应于不同的I类错误，发明人基于0.05的FDR而不是名义p值阈值来计算功效，以便在方法之间进行公平比较。在不相关和相关的水平多效性效应都不存在的情况下，MRAID、IVW-R和RAPS在不同的情形下都比其他方法具有更高的功效。在这三种方法中，当工具对暴露的效应很小或因果效应很小时(图3a、3b)，MRAID比其他两种方法的功效略高一些，这可能是因为MRAID中实行了工具变量自动选择。当工具对暴露效应较大且因果效应较大时，MRAID的功效略小于其他两种方法，因为在其他方法中使用的简单工具选择方法可能在这些不太具有挑战性的设置中更有效。一般除了这三种方法，表现较好的是Robust。而Weighted mode，MRMix和CAUSE有轻微的低功效表现。

发明人研究了水平多效性对不同方法功效的影响。当水平多效性效应存在但与工具变量对暴露的效应不相关时，MRAID比其他MR方法更有效(图3C)。MRAID在功效方面的良好表现，随着水平多效性的增加变得更加明显，水平多效性的增加主要是体现在水平多效性效应的增大或和具有水平多效性效应的SNPs的比例增加。RAPS、Robust、CAUSE和IVW-R方法的表现次之，而MRMix和Weighted mode具有较低的功效。在这些方法中，IVW-R的表现对水平多效性效应的大小或具有水平多效性效应SNPs的比例特别敏感。当不相关的水平多效性效应和相关的水平多效性效应共同存在时，MRAID的功效仍然高于其他方法。无论相关水平多效性效应大小、具有相关水平多效性效应的工具SNPs的比例(图3d)、相关水平多效性效应的产生方式，MRAID的高功效表现都保持不变。随着具有不相关水平多效性效应的工具SNPs比例的增加，MRAID带来的功效优势尤其明显(图3d)。重要的是，MRAID的功效接近于oracle MR方法，特别是因果效应较大时，oracle方法使用真正的工具变量SNPs集进行MR推断，这种表现也是对MRAID中工具变量自动选择有效性的佐证。

接下来，发明人通过反向因果分析检验了不同MR方法辨别因果效应方向的能力。特别是，发明人在非零模拟中检验了结局对暴露的因果效应，其中暴露对结局有因果效应，但反之不成立。在存在水平多效性时，在反向MR分析中获得的与结局相关的SNPs工具变量将包含两种：一种是通过暴露间接与结局相关的SNPs工具变量，另一种则是直接与结局相关的SNPs，这种与结局直接相关的SNPs通常被认为它们对结局具有水平多效性效应。因为这两组SNPs表现出对暴露的异质性效应，发明人无法检测出结局对暴露的非零因果效应。因此，在水平多效性存在下的反向因果关系分析有效地用于零模拟的分析。事实上，发明人发现在一系列模拟情形中，MRAID提供了反向因果分析中有效的I类错误控制和校正的p值。相反，其他方法的I类错误控制高度依赖于水平多效性的程度。具体而言，当少量暴露工具变量SNPs表现出对结局的水平多效性时，在反向因果分析中，大多数结局的候选工具SNPs不会表现出对暴露的水平多效性效应。在这种情况下，CAUSE和Weighted mode都过于保守，而IVW-R、MRMix、RAPS和Robust产生了略为膨胀的p值。相反，当暴露的大多数工具SNPs对结局有水平多效性效应时，反向因果关系分析中结局的候选工具SNPs也将表现出水平多效性效应。在这种情况下，正如发明人在相应的零模拟情形中所示，MRMix、Robust、IVW-R和RAPS都产生了膨胀的p值。这些方法的p值在水平多效性效应更小时变得十分显著，在这种情况下，选择第二种SNPs作为工具变量就会更加困难。

最后，MRAID在零假设和各种备择假设下都产生了近似无偏的因果效应估计。

实际数据分析应用

发明人应用MRAID和其他MR方法分析了英国生物银行(UK Biobank)中38种生活方式危险因素和11个与心血管相关的性状(详细信息在具体实施方式部分)。具体来说，发明人将UK Biobank数据分为两个独立的、大小相等的子集，分别代表暴露GWAS和结局GWAS。发明人进行了两组分析。首先，发明人聚焦于8个与心血管疾病相关的性状，并检验了暴露GWAS中每个性状对结局GWAS中相同性状的因果效应，有效地检验了该性状对自身的因果效应。在这种分析中，真正的因果效应正好等于1。发明人发现所有的方法都能检测到8个与心血管疾病相关的性状对自身的非零因果效应(图4)。然而，只有MRAID和CAUSE能够产生涵盖所有8个性状对真正因果效应的95％置信区间，CAUSE产生的置信区间比MRAID长2.39-5.69倍(图4)。例如，在HDL-HDL分析中，只有MRAID(估计值＝0.98；95％CI:0.96-1.01)，CAUSE(0.95；0.82-1.09)和MRMix(0.96；0.90-1.02)能够正确推断因果关系，但MRAID可以提供最小的置信区间(图4h)。相反，其他四种方法获得的置信区间并不能涵盖真正的因果效应。在LDL-LDL分析中，MRAID(0.97；0.94-1.01)和CAUSE(0.96；0.84-1.08)正确推断了因果效应，MRAID提供了较小的置信区间(图4f)。而其他五种方法的置信区间也没有涵盖真正的因果关系。结果表明，MRAID和CAUSE都能对性状本身的分析产生精确的因果效应估计和校正的置信区间，其中MRAID比CAUSE的统计效能更好。

接下来，发明人研究了38个生活方式风险因素与11个心血管疾病相关性状之间的因果关系。生活方式风险因素与心血管疾病相关性状的关联性已被广泛证明。然而，由于一些关联效应在不同研究中有不同的估计，因此对于检测到的关联关系是否是因果关系仍存在争议。发明人既进行了正向因果分析，检验生活方式风险因素对心血管疾病相关性状的因果效应，又进行反向因果分析，检验心血管疾病相关性状对生活方式风险因素的因果效应。不同方法分析的性状对的p值分布如图5a所示。与模拟一致，发明人发现MRAID(膨胀系数GIF＝0.90)和MRMix(GIF＝0.78)的p值在分析中表现良好，略微保守，优于其他方法(图5a)。同样与模拟一致的是，发明人发现CAUSE的p值过于保守(GIF＝0.12)，而RAPS(GIF＝1.96)、Weighted mode(GIF＝1.70)、IVW-R(GIF＝2.12)和Robust(GIF＝2.00)的p值都显示出明显的膨胀(图5a)。事实上，在置换分析中，只有MRAID产生了校准的p值(图5b)。

基于Bonferroni校正的p值阈值(8×10^―5)，MRAID检测到8种具有很强的生物学证据的因果关系。例如，MRAID检测到吸烟对BMI具有负因果效应。吸烟和肥胖之间的负相关已经在观察性研究和MR研究中得到了充分的证明。具体来说，吸烟期间摄入尼古丁会降低静息代谢率，并抑制脂蛋白脂肪酶活性和其他激酶途径来减少脂多肽，所有这些都会导致脂肪组织中净能量存储的减少，从而导致体重减轻。尼古丁还会激活下丘脑中的乙酰胆碱受体，进而产生厌食神经元，从而抑制食欲和食物摄入量。还有，MRAID检测到有吸烟史的人开始吸烟的年龄对HDL的影响，表明吸烟行为对HDL的具有负向因果效应。众所周知，吸烟行为与HDL有因果关系。特别是吸烟可以改变脂质运输的关键酶的活性，降低卵磷脂-胆固醇酰基转移酶(LCAT)的活性，改变胆固醇酯转移蛋白(CETP)和肝脂肪酶的活性，所有这些都可能降低HDL代谢。此外，吸烟诱发的氧化修饰会使HDL功能失调，并剥夺其动脉粥样硬化保护特性。

重要的是，与其他方法相比，MRAID并没有错误地检测出许多错误的因果关联。吸烟对血压的影响便是一个众所周知的具有潜在假因果关联的例子。在观察性研究中发现，吸烟和血压之间存在负相关。然而，在大型数据集上的多个后续MR研究并不支持这两个性状之间的因果关系。事实上，在观察性研究中，吸烟和血压之间的联系很可能受到年龄、BMI、社会阶层、盐摄入量、饮酒习惯以及其他未测量的混杂因素的干扰。与这些先前的MR研究一致，MRAID没有检测到与吸烟相关的八种特征中的任何一种对收缩压或舒张压有显著的因果效应。相比之下，几乎所有其他的方法都错误地发现了一些与吸烟有关的特征对血压的因果效应。例如，不成功戒烟次数对SBP的因果效应并没有被MRAID(p＝0.44)、CAUSE(p＝0.01)或Weighted mode(p＝1.3×10^―4)检测到，但可以被IVW-R(p＝1.4×10^―6)、Robust(p＝4.1×10^―34)、RAPS(p＝5.5×10^―6)和MRMix(p＝2.8×10^―6)错误地识别出。同样地，有吸烟史的人开始吸烟年龄对SBP的因果影响也没有被MRAID(p＝0.06)和CAUSE(p＝1.3×10^―3)检测到，但也被IVW-R(p＝7.8×10^―5)、Robust(p＝1.5×10^―30)、RAPS(p＝8.9×10^―7)、Weighted mode(p＝1.7×10^―6)和MRMix(p＝4.1×10^―7)错误地识别出。BMI对驾驶时间的影响是另一个过度识别的例子，BMI不太可能对驾驶时间产生因果关系，至少不是正向的。事实上，MRAID(p＝0.01)和MRMix(p＝0.12)、CAUSE(p＝0.02)和Weighted mode(p＝0.04)都没有检测到BMI对驾驶时间有任何因果效应，但是，IVW-R(p＝2.5×10^―6)和RAPS(p＝3.1×10^―5)均检测到了BMI对驾驶时间的假阳性影响。

最后,发明人注意到MRAID一个重要的特性是它能够有效地将SNP对结局的影响分解为三个不同的路径:一个是SNP直接对结局进行作用，一个通过中介对结局进行作用，还有通过影响暴露和结局的混杂因素进行作用。因此，MRAID可以用来估计不同类别中SNPs的比例，包括与暴露相关的SNPs在全基因组显著的SNPs中所占的比例(π_β)，表现出相关水平多效性的SNPs的比例(π_c)，在选定的工具变量中具有不相关水平多效性的SNPs的比例(π₁)，以及在剩余的候选工具变量中具有不相关水平多效性的SNPs的比例(π₀)。在实际数据应用中，发明人估计了在645对性状对中，π_β,π_c,π₁和π₀的均值分别为14.6％、6.4％、16.4％和5％。此外，发明人估计了在8对有意义的性状对中它们平均值分别为6.2％，5.7％，11.4％和0.1％。具有相关水平多效性的SNPs所占比例与具有不相关水平多效性的SNPs所占比例也高度相关，后者一般大于前者。这些比例的估计支持了先前在复杂性状分析中识别出的广泛存在的水平多效性，并为这两种水平多效性对MR分析的影响程度提供了详细的量化数据。

讨论

MRAID是一种两样本孟德尔随机化的方法，该方法可以实现在相关的候选SNPs集合中自动筛选合适的工具变量，并且可以在基于似然的推理框架中控制相关和不相关的水平多效性。总的来说，通过工具SNPs的自动选择和似然框架下的统计推断，MRAID可以在各种各样的场景下给出校准的p值，并且具有比现存方法更高的统计效能。MRAID的优点在模拟和复杂性状的实际数据分析中均有体现。

在本研究中，发明人主要侧重于利用MRAID对数量性状进行建模。对于二值暴露和结局，可以把它们当作连续变量，直接应用MRAID进行MR分析。将二值暴露和结局视为连续变量可以通过将线性模型视为广义线性模型(如logistic回归)的一阶泰勒近似来证明。然而，只有当SNP对暴露和结局的效应相对较小时，这种近似才准确。虽然类似的方法已经在许多以前的MR研究中被应用，但发明人要提醒的是，当线性模型被用于拟合二值暴露和结局时，因果效应的解释可能具有挑战性，特别是当两阶段推断方法被用于MR分析时。例如，当二值暴露是基于连续风险因素的二分化产生时，不考虑潜在的连续风险因素，直接通过对二值暴露建模获得因果效应估计可能需要额外的建模假设。此外，在不考虑潜在的连续风险因素的情况下对二值暴露进行建模，可能会导致违反排除限制假设，因为即使二进制暴露不发生变化，这些工具变量也会通过连续风险因素影响结局。因此，对MRAID进行拓展，使其能够明确地对数量性状以外的数据类型进行建模，对于确保其被广泛应用非常重要。因为MRAID是数据驱动的建模方式，是通过极大似然框架对SNP-暴露模型和SNP-结局模型进行联合推断，所以它可以自然地被扩展到对其他类型的暴露或结局数据的建模，例如，推广到广义线性模型框架下。目前，发明人所知道的唯一能同时考虑二值风险因素和结局的基于概率的MR方法是IV-MVB。但是，IV-MVB需要个体水平的数据，且仅适用于单一样本的分析，其并不能高效地同时解决多个工具变量问题，尤其是当这些工具变量间存在相关关系的时候。因此，探索MRAID拓展至通用数据类型是未来研究的一个重要方向。

MRAID并非没有局限性。首先，当MRAID进行SNP工具变量的自动选择时，这种选择是建立在稀疏建模假设上的，该假设指定了SNP效应大小。稀疏建模假设包含多个超参数，而这些超参数的推断则是依赖于一种基于抽样的算法。对超参数的准确和可靠估计可能需要一定数量的候选工具变量。虽然在实际数据分析中，通过MRAID选择出的性状对的意义似乎并不依赖于被选择的候选工具数量(图5)，但仍需说明的是，当GWAS的样本量较小，暴露呈现非多基因结构时，MRAID的统计功效可能会降低。其次，MRAID主要遵循CAUSE的方法，通过引入一个潜在变量作为暴露和结局的混杂因子来校正相关水平多效性。由于其在仅建模单一未观察混杂因素方面的局限性，MRAID可能在暴露与结局之间存在多个或其他类型的混杂因素下并不完全有效。最后，MRAID的汇总数据版本需要两个LD矩阵作为输入，一个来自暴露GWAS，另一个来自结局GWAS。在目前的研究中，发明人利用个体数据估计了两个LD矩阵。在缺乏个体数据的情况下，可以从具有相同遗传祖先的参考面板(例如，千人基因组项目)估计这两种LD矩阵。然而，当选择的暴露GWAS和结局GWAS具有不同的遗传祖先时需要特别注意。

详细信息

1.个体水平数据的MRAID模型

发明人遵循发明内容中的符号，并以下列形式重写MRAID模型：

x＝Z_xβ+ε_x, (1)

其中v是二分类的p维列向量，指示工具变量SNPs是否具有相关多效性(v_j＝1)或者(v_j＝0)；

表示向量β和v的哈达玛积，也被称为元素对应乘积。对于j＝1,…,p，进一步有：

若β_j≠0,

若β_j＝0,

若β_j≠0，p(v_j＝1)＝π_cβ_j≠0

假设

服从一个逆伽马分布，

设置

和b_β＝0.2，以确保其先验均值为2/p。类似地，假设

其中

和b_η＝0.2来确保先验均值为1/p。假设比例参数服从贝塔分布：π_β～beta(λ_β1,λ_β2)，其中λ_β1＝0.5和λ_β2＝4.5，以确保先验均值为0.1；π_c～beta(λ_c1,λ_c2)，其中λ_c1＝0.5和λ_c2＝9.5，确保先验均值为0.05；π₁～beta(λ₂₁,λ₂₂)，其中λ₂₁＝0.5和λ₂₂＝1.5，确保先验均值为0.25；以及π₀～beta(λ₃₁,λ₃₂)，其中λ₃₁＝0.05和λ₃₂＝9.95，确保先验均值为0.005。这些先验均值被用来代表先验信息，即相对较小比例的SNPs对暴露具有非零的效应，相对较小比例的SNPs具有水平多效性，被选择的工具变量SNPs比未被选择的SNPs更容易表现出水平多效性，并且具有水平多效性的SNPs相较于相关多效性更易表现出不相关的多效性。对于其他超参数，发明人设置了

p(ρ)∝1。对于因果效应参数α，发明人设置了

其中

为了便于计算，发明人引进了一个p维的二分类指标γ＝(γ₁,…,γ_p)^T用来表示每一个SNP对暴露有非零效应(γ_j＝1)或者没有非零效应(γ_j＝0)。发明人也引入了另一个p维的二分类指标τ＝(τ₁,…,τ_p)^T来表示每一个SNP具有水平多效性效应(τ_j＝1)或者没有(τ_j＝0)。因此，有p(γ_j＝1)＝π_β，p(τ_j＝1|γ_j＝1)＝π₁，p(τ_j＝1|γ_j＝0)＝π₀和p(v_j＝1|γ_j＝1)＝π_c。

2.详细的抽样步骤和高效的计算方法

首先，发明人的目标是获得因果效应α的后验样本。为此，发明人依靠Gibbs抽样来获得参数的条件分布。发明人的后验似然的形式是

f(α,β,γ,η_u,τ,v,ρ|x,y)∝f(x,y|α,β,γ,η_u,τ,v,ρ)f(β|γ)f(η_u|τ)f(v|γ)f(τ|γ)f(γ|π_β)f(π_β)f(π₁)f(π₀)f(π_c)∝f(x|β)f(y|α,β,η_u,v,ρ)f(β|γ)f(η_u|τ)f(v|γ)f(τ|γ)f(γ|π_β)f(π_β)f(π₁)f(π₀)f(π_c)∝

其中

K是一个标量，它的表达式是

β_―j是一个p―1维向量，表示向量β不包含第j个元素β_j，Z_x,―j是一个n₁×(p―1)的矩阵，表示Z_x没有第j列向量Z_x,j；Z_y,―j是一个n₂×(p―1)矩阵，表示Z_y没有第j列向量Z_y,j。需要说明的是，由于Z_x，Z_y被标准化，其均值为0、标准差为1，则对于所有j＝1,…,p，有

基于上述似然，发明人对所有参数

进行Gibbs抽样。

具体地说，给定γ_j＝1和其他参数，β_j的后验条件分布是均值为

方差为

的正态分布，其中

给定τ_j＝1等参数，η_uj的后验条件分布是均值为

和方差为

的正态分布，其中

在将参数β_j积分后，可以获得γ_j的后验条件分布

在将参数η_j积分后，可以获得到τ_j的后验条件分布

同时，发明人还得到了v_j的后验条件分布

ρ的后验条件分布是一个方差为

和均值为

的正态分布。

α的后验条件分布是一个方差为

和均值为

的正态分布。

π_β的后验条件分布是

π₁的后验条件分布是

π₀的后验条件分布是

π_c的后验条件分布是

的后验条件分布是一个形状参数为

比例参数为

的逆伽马分布。

的后验条件分布是一个形状参数为

比例参数为

的逆伽马分布。

的后验条件分布是一个形状参数为

和比例参数为

的逆伽马分布。

的后验条件分布是一个形状参数为

和比例参数为

的逆伽马分布。

对于本研究中的所有分析，发明人运行了1000次Gibbs抽样迭代，前200次作为测试，并获得了用于假设检验的α的后验样本，这将在下一节中描述。

3.参数估计和p值计算

发明人将μ_α和

作为因果效应参数α的后验均值和后验方差，由于似然函数和后验分布都渐近地服从正态分布，同时发明人使用了正态分布

作为先验分布，发明人可以很容易地用矩方法得到近似的极大似然估计及其标准误

由此得到的

渐近地服从标准正态分布，这允许发明人获得一个相应的p值来检验零假设H₀:α＝0。在目前的研究中，由于发明人通过设置

对α给定了一个非信息性先验，因此发明人有

和

请注意，因为发明人仅对假设检验使用后验和先验之比表示作为似然^1,2，所以在α上选择的先验并不影响推理结果，。

4.汇总统计量的MRAID模型

汇总统计量的MRAID模型可以以边际效应大小估计和LD矩阵的形式给出。发明人在暴露GWAS数据中将工具变量SNPs的LD结构表示为Σ₁，在结局GWAS数据中为Σ₂，两者都是p×p的对称正定矩阵。通常，Σ₁和Σ₂来自相同的LD参考面板(例如，来自千人基因组项目的具有欧洲血统的个体)。汇总统计量的MRAID模型可以由以下两个方程式构造

其中，

是一个p维向量，表示SNPs对暴露的边际效应大小；

是一个p维向量，表示SNPs对结局的边际效应大小。e_x和e_y都是表示残差的p维向量，e_x服从多元正态分布

e_y服从多元正态分布

5.汇总统计量MRAID模型的详细抽样步骤

对于MRAID的汇总统计量版本，发明人可以得到所有参数的后验分布

其中

K是一个标量，它表达式是

β_―j是一个p―1维向量，表示向量β不包含第j个元素β_j；

基于上述似然，发明人对所有参数

进行Gibbs抽样。

给定γ_j＝1和其他参数，β_j的后验分布是均值为

方差为

的正态分布，其中

给定τ_j＝1和其它参数，η_uj的后验条件分布是均值为

方差为

的正态分布，其中

在将参数β_j积分后，可以获得γ_j的后验条件分布

在将参数η_uj积分后，可以获得τ_j的后验条件分布

同时，发明人还得到了v_j的后验条件分布

ρ的后验条件分布是一个方差为

均值为

的正态分布。

α的后验条件分布是一个方差为

均值为

的正态分布。

π_β的后验条件分布是

π₁的后验条件分布是

π₀的后验条件分布是

π_c的后验条件分布是

的后验条件分布是形状参数为

比例参数为

的逆伽马分布。

的后验条件分布是形状参数为

比例参数为

的逆伽马分布。

的后验条件分布是形状参数为

比例参数为

的逆伽马分布。

的后验条件分布是形状参数为

比例参数为

的一个逆伽马分布。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于联合似然的孟德尔随机化分析方法，其特征在于，包括：

对一组互相之间处于高LD的初始候选SNP工具变量进行建模，并从中自动筛选合适的工具变量进行MR分析；

同时有效控制与工具变量和暴露之间的效应相关的或不相关的两种水平多效性效应；

采用可扩展的抽样算法来计算校准的p值；

对一组互相之间处于高LD的初始候选SNP工具变量进行建模包括通过以下三个线性回归对暴露，结局和基因型之间的关系进行建模：

x＝Z_xβ+ε_x，

其中，Z_x是一个n₁×p的基因型矩阵，表示的是在暴露GWAS中n₁个个体p个候选工具变量信息，Z_y是一个n₂×p的基因型矩阵，表示的是在结局GWAS中n₂个个体p个候选工具变量信息，y为结局变量，

为未观察到的暴露，β是p维向量，表示的是SNPs对暴露效应，η₀和η₁是表示水平多效性效应的p维向量，α代表的是暴露对结局的因果效应标量，ε_x是表示残差的n₁维向量，其每个元素独立同分布于正态分布

是表示残差的n₂维向量，其每个元素服从相同的正态分布

ε_y是表示残差的n₂维向量，其每个元素独立且服从正态分布

2.根据权利要求1所述的基于联合似然的孟德尔随机化分析方法，其特征在于，还包括在两样本MR框架下估计和检验暴露变量对结局变量的因果效应，其中暴露变量和结局变量是在两个独立的没有样本重叠GWAS中测量。

3.根据权利要求1所述的基于联合似然的孟德尔随机化分析方法，其特征在于，所述重叠GWAS分别称为暴露GWAS和结局GWAS，在暴露GWAS中，按照边际p值低于全基因组显著性阈值的标准进行初步筛选，选出与暴露变量相关的SNPs。

4.根据权利要求1所述的基于联合似然的孟德尔随机化分析方法，其特征在于，所述自动筛选合适的工具变量包括对SNP对暴露的影响效应β进行了稀疏性建模假设来进行工具变量自动筛选，其具体为：提出假设

其中δ₀是表示在零处的质点的狄拉克函数，第j个SNP对暴露的效应有1-π_β的概率为零，有π_β的概率效应不为零，若存在非零效应，其效应服从均值为0，方差为

的正态分布，其中方差参数

决定效应大小。

5.根据权利要求1所述的基于联合似然的孟德尔随机化分析方法，其特征在于，还包括在η₀和η₁上指定了单独的模型假设，具体为：假设每个选定的SNP工具变量都有π_c的概率通过混杂因素导致相关水平多效性，其SNP对结局的混杂效应是其对暴露效应的ρ倍，如果SNP对暴露的效应是β_j，选定的SNP通过混杂对结局的效应变成ρβ_j，则对

的假设，是

假设每个选定的SNP工具都有π₁的概率不通过暴露直接对结局表现出不相关水平多效性，则对

的假设为

故对于被选定的工具变量，其水平多效性为

对于未被选则的非工具变量，假设

使用相同的方差参数

来对工具变量和非工具变量SNP的不相关水平多效性效应建模。

6.根据权利要求1所述的基于联合似然的孟德尔随机化分析方法，其特征在于，所述采用可扩展的抽样算法来计算校准的p值包括：利用Gibbs抽样获得因果效应α的后验样本，利用该样本均值和样本标准差去归纳后验分布。

7.根据权利要求1所述的基于联合似然的孟德尔随机化分析方法，其特征在于，所述采用可扩展的抽样算法来计算校准的p值包括：给因果关系α指定一个正态的先验分布，获得其先验均值和标准差，利用矩估计的方法来获得近似最大似然估计

和标准误

构建一个近似的Wald检验统计量，获得假设检验的p值。

Images