基于贝叶斯网络的孟德尔随机化用于变异优先级排序和表型因果推断

摘要

背景

一般来说，孟德尔随机化（MR）分析依赖于三个核心假设（图1a）：（i）相关性：IV与暴露之间存在相关性;（ii）外生性或可交换性：IV与暴露和结局之间的任何混杂因素无关（ ）;（iii）排他性：IV仅通过暴露影响结局（ ）。

然而严格的假设经常被违反（图1b），这使得确定合适的IV变得具有挑战性。首先，尽管GWAS研究已经确定了许多风险位点但对多基因复合性状的影响通常很小，导致弱工具偏倚。多重检验负担、“赢家诅咒”、连锁不平衡（LD） 和人口分层增加了GWAS 中假阳性信号的风险，这可以通过应用多种工具来改进，然而相关的工具也会导致不稳定的估计，并在纳入非因果变异时引入额外的遗传混杂。已提出的策略如 LD 逐步修剪、主成分分析（PCA）和惩罚等，旨在从大量相关的弱工具中提取适当数量的独立工具，但面临对稳健性的考虑。

另一个问题是，由于水平多效性（变异通过感兴趣的暴露以外的其他途径影响结局），许多IV实际上是无效的。基因相互作用如LD和异位显性，也可能违反类似于多效性的排他性。对于个体层面的数据，lasso类型的方法，如sisVIVE和post-adaptive Lasso有助于控制多效性效应的影响。最近的一些方法，如 TSHT 和 CIIV通过从候选集合中识别有效IV来减轻多效性。

上述方法依赖于许多难以置信的假设，很难对复杂的真实遗传模式进行建模。特别是，由于复杂的基因相互作用（GxG），将非因果变异作为 IV 不仅会导致不稳定的估计并影响统计功效，还可能引入 GxG 诱导的相关多效性，违反InSIDE假设（图 1c）。因果图模型提供了一种表示潜在因果关系的替代方法。通过因果图，机器学习技术如因果贝叶斯网络（BN）等目前应用于识别遗传相互作用和因果变异。它们也将成为传统 MR 的有效的补充。

在本文中，我们通过整合因果发现和推理提出了一种基于贝叶斯网络的孟德尔随机化（BNMR）方法（图1d）。我们的目标是在学习阶段解决相关的弱工具的问题，在推理阶段应对多效性问题。随机图森林 （RGF）是一种由一系列BN结构学习过程组成的集成方法。我们应用 RGF，优先考虑效应较小且相互作用的变异，并将对暴露有直接影响的变异识别为工具。然后，我们通过贝叶斯MR框架估计因果效应，并使用收缩先验应对潜在的水平多效性。我们通过模拟证明了BNMR在工具变量选择和效应估计方面都优于传统方法。通过应用到UK Biobank数据，我们研究了从血液学特征到血压和精神疾病的因果途径，带来了新的生物学见解。

方法

BNMR模型概览

BNMR是一个使用个体级数据的两阶段MR框架。在学习阶段，我们应用RGF从大量相关的弱工具中选择具有可靠相关性的变异。我们利用BN来表征复杂的条件概率关系，并根据它们与感兴趣暴露的关系将变量集划分为三个子集（DIE分割），

中的变异直接影响暴露，中的变异通过基因相互作用或者连锁间接影响暴露，即中的变异与被中的变异d分离（），中的变异不影响暴露（）。这三个子集在因果马尔科夫、忠实性和充分性假设下通过 BN 进行区分，且只能是因果图中的父节点集合。

在推理阶段，我们显式地模拟了潜在的水平多效性。由于数量性状是由遗传和环境因素决定的，假设线性且没有相互作用，我们有

其中和表示对和的效应。变异通过两种不同的途径影响结局：感兴趣的因果路径，即通过暴露的路径；水平多效性路径，即直接影响或通过以外的其他路径。

其中是  对  的因果效应,而 表示多效性效应。通过引入  和 之间的相关矩阵 得到

重写方程(2)(3)得到

其中，，。和相关但都与独立。事实上只有一部分变异需要作为 IV 被纳入，因此需要进行高精度的识别。我们对于冗余参数施加一个收缩先验以使得可识别。

随机图森林 RGF中的BN结构学习

为了降低结构学习的计算复杂性并评估每条边的置信度，我们提出了受随机森林启发的RGF。在RGF 中，使用 bootstrapping 或下采样创建个子图，在每个子图s中，个个体中的个个体和个变异中的个变异被采样。因此，我们将DIE分割的过程归结为一系列因果BN的结构学习。因此，我们将DIE分割的过程归结为一系列因果BN的结构学习。

由于变异是自然性状的原因，我们可以简化结构学习以绘制骨架确定图。我们通过扫描每个图中与暴露直接相邻的变异并计算每个变异的邻接分数（所有子图中出现边的频率）来识别，这是平均因果图中变异-暴露相关性的置信度。邻接分数较高的变异具有较高的置信水平，可以被划分为。我们可以选择指定数量的主导变异或分数高于给定阈值的变异作为IV。

目前已经提出了很多用于BN结构学习的算法。基于得分的方法通过详尽或启发式地探索候选图并最大化网络得分来确定最佳网络，而基于约束的方法利用一系列条件独立性测试来建立节点之间的边缘约束，并随后细化方向。我们选择了基于得分的方法，包括爬山算法（hc） 和禁忌搜索 （tabu），基于约束的方法，包括稳定PC （pc.stable）、iamb和gs，以及混合学习方法，包括 Max-Min 爬山 （mmhc） 和限制最大化 （r smax2）。所有这些方法都是使用 R 包 bnlearn 实现的。

具有收缩先验的贝叶斯MR估计

我们在贝叶斯的框架下指定模型（6）。总误差项和可以分解为一个混杂相关项和完全随机项，即，，。假设两个完全随机项是不相关的，我们有。我们只需要选择的子集作为IV。贝叶斯 MR 的模型为

BNMR可以将方程(8)修改为 probit或者 logistic 回归来拓展到二元结局，

其中连接函数可以是inverse-probit或者inverse-logit。

我们将 BNMR 的估计值与其他 IV 选择和 MR 估计方法进行比较。我们使用 R 包统计信息实现 PCA，并使用 R 包 glmnet 实现惩罚回归，其中 10折交叉验证用于确定的最佳 值。比较方法是使用R 包 AER、ivmodel、MendelianRandomization、cause、R2BGLiMS 和 CIIV 实现的。我们将 BNMR 实现为 R 包，源代码可在 https://github.com/sjl-sjtu/bnmr2 上获得。

模拟

我们在模拟中使用来自UK Biobank的模拟和真实基因组学 数据。对于模拟基因组学，从多项分布中采样的个独立位点，其基因型频率满足Hardy-Weinberg 平衡 （HWE），效应等位基因频率。根据从采样的LD平方相关系数模拟每个位点的个相关位点，并合成具有位点的基因组学。用于模拟表型的真实基因组数据来自UK Biobank欧洲血统人群中10、17 和 22 号染色体上的变异。

表型由线性模型生成，

其中因果效应为，是的子集。混杂且。从模拟基因组中随机挑选出影响的遗传变异，其效应大小为。对于多效性变异，（平衡多效性）或者（方向多效性）。我们对每个情景重复100次并报告平均值。

结果

BNMR可以在学习阶段以高精度从众多弱的相互作用的变异中有效地识别效应变异

首先，我们比较了RGF在模拟数据集中不同超参数和结构学习算法的表现（图2）。RGF 表现出较低的错误发现率（FDR）和较高的 AUC，随着子样本个数和子样本的样本量的增加，尽管这种改进伴随着时间消耗的增加。基于约束的方法产生较低的 FDR，而基于分数的方法在速度上更胜一筹。随着选择阈值的增加，已识别变异的数量随着精度的提高而减少。

    图 3 模拟中不同 IV 选择策略的错误发现率（FDR）。a：RGF的FDR 。b：LD逐步修剪的FDR。c：惩罚回归的FDR。d：GWAS的主导 SNP 的 FDR。

与传统的关联检验（线性回归）相比，RGF、LD逐步修剪和惩罚回归（尤其是套索和弹性网）都可以降低FDR，而RGF实现了最高的精度，是确定候选效应变异优先级和识别真实效应变异的有效工具（  ）（图 3）。在采用这些变量选择策略之前，我们进行预过滤以减少候选变异的数量。在 RGF 之前，更严格的阈值可提高top变异的精度，但可能会减少召回率。大约至的阈值可能是合适的。此外，邻接分数可以看作是变体所属的置信度，因此，它也是对IV强度的评估。邻接分数与常用的F统计量之间的相关性表明，RGF能够选择可靠的工具来减少弱工具偏倚。

图 4 BNMR 在基于英国生物样本库真实基因型的模拟表型上的表现。a：不同样本量（n）的真实基因型数据的ROC曲线。这里的 p 固定为 20,000。b： 具有不同基因组大小的真实基因型数据的ROC曲线（p）。这里的 n 固定为 5000

由于基因型是难以模拟的，我们使用方程(10)生成表型，基于来自英国生物银行的真实基因型数据。使用合成数据集，我们展示了RGF对不同规模基因组的适应性（图4）。RGF能够处理真实遗传数据的复杂结构，对不同尺度的数据集表现出良好的适应性。

BNMR可以有效降低估计值的均方误差，增强统计功效，并且在推理阶段对水平多效性具有鲁棒性

我们首先将BNMR的性能与其他现有的MR方法进行比较，包括两阶段最小二乘法（TSLS），有限信息最大似然比（LIML），逆方差加权（IVW），MR-Egger，加权中位数，加权众数，JAM-MR，CAUSE和 CIIV（图 5a）。我们在模拟中包括两种类型的多效性位点：独立影响暴露和结果的多效性位点，或由基因相互作用（如连锁）导致的对暴露和结果的影响。对于每种情况，我们检查了所有基因座的预期平均多效性效应为 0（平衡多效性）或非零（定向多效性）的设置下各种方法的性能。大多数流行的方法在平衡多效性方面表现相对较好，但无法应对具有复杂方向性和相关多效性的情况。由于MR-Egger对InSIDE假设的敏感性，当相关多效性变异的数量增加时，其表现明显更差。一些双样本方法，如 CAUSE，在应用于单样本研究时，会面临夸大的估计方差和样本重叠的偏差。广义上基于多元的方法，如加权中位数和众数估计，表现出稳定性。一般来说，BNMR在均方误差（MSE）方面优于现有方法，特别是由于其估计方差较小，从而增强了统计能力。 尽管依赖于 InSIDE 假设，但使用 RGF 进行 IV 选择的过程增强了 InSIDE 假设的稳健性，使其更适用于对基因-基因相互作用产生的相关多效性。

为了显示 BN 为传统贝叶斯 MR 带来的加成，我们随后评估了使用从 BN 获得的 IV （BNMR） 与直接使用 GWAS 先导 SNP 作为 IV （BMR） 相比，使用贝叶斯 MR 估计的改进。当存在基因相互作用的相关多效性时，可以观察到MSE的显着降低（图5b）。这种增强主要表现在平衡情景下估计方差的衰减，而当多效性效应是方向性的时，偏差和方差都会减弱。为了更好地了解BN的作用，我们检查了RGF选择的IV在其他传统MR方法上的性能（图5c）。由于单样本设计中的样本重叠，以及 IVW 和 MR-Egger 使用独立 IV 的要求，而 RGF 旨在识别对暴露有直接影响的 IV（可能彼此不独立），而使用 RGF 选择的 IV 在使用 IVW 估计器时效果不佳，甚至会增加偏差。即便如此，另一方面，当使用 MR-Egger 估计时，与通过 LD 修剪获得的 IV 相比，使用 RGF 选择的 IV 减少了估计的偏差和方差。我们认为这是因为 MR-Egger 和贝叶斯 MR 基于类似的假设，即 InSIDE 假设，要求  对  的影响和从  到  的多效性通路是独立的。当  影响 和  的混杂（即相关的水平多效性）时违反这一假设。如果这种相关的水平多效性是由基因-基因相互作用引起的，其中是另一个遗传位点（图1c），RGF将倾向于选择 作为IV，以确保仍然满足InSIDE假设。因此，RGF增强了InSIDE假设的鲁棒性。

我们还对工具数量和迭代次数以及不同的收缩先验进行灵敏度分析（图 5）。当工具过多或过少时，偏倚会增加。BNMR估计通常对先验不敏感，尽管基于误差线的均匀尖峰和平板先验比其他先验效率低一些（图5d）。贝叶斯套索先验显示出最快的采样速度和放缩的标准误差，但偏差略高。马蹄先验虽然效率略低，但由于Rhat最低，收敛性能更胜一筹。

BNMR应用于大规模生物样本库级别的数据，证明了从红细胞相关性状到血压的因果关系

我们利用来自英国生物银行的数据，提供了两项具有连续和二元结果的真实世界研究的说明性示例。血液学指标通常在各种生理过程中有所不同，是相关疾病的潜在指标。

相关研究表明，红细胞相关特征，包括红细胞计数 （RBC）、血红蛋白浓度 （HGB）、血细胞比容 （HCT）、红细胞与血浆的比例和平均红细胞体积 （MCV），与收缩压和舒张压（SBP 和 DBP）密切相关，红细胞异常可能是某些心脑血管疾病（如高血压）的指标。为了检查红细胞参数对血压的因果影响，我们涉及 246,659 名白种人血统的参与者，他们自我报告没有高血压或其他心血管疾病（英国生物银行非癌症疾病编码 1065-1094），并在入组时进行了可用的血液常规测量。使用PLINK 2.0进行基因组质量控制，SNP缺失率、次要等位基因频率（MAF）和HWE检验的相应阈值分别为0.05、0.01和1e−6。使用 Python 包 PyMC 和 JAX 对大型数据集进行快速后验采样。为了提高功效，我们使用来自同一种族群体但具有不同样本的不同数据集的汇总统计进行了初步的GWAS过滤。

我们采用两种预过滤策略来减少候选变体的数量，然后进行 BNMR 分析。第一种策略采用更严格的 GWAS阈值1e−20，而第二种策略采用更宽松的阈值 5e−8，然后是 LD clumping。结果（图6e）一致表明，RBC、HGB和HCT对DBP和SBP均表现出显著的正向影响，并且对SBP的影响幅度大于对DBP的影响幅度。而 MCV 对血压的负面影响不显着。替代方法使用LD聚集后GWAS的顶级显著SNP作为工具，不同甚至相互矛盾的结果（表S5）提醒我们MR方法的重要性。MR-Egger检验表明，所有因果关系均不显著。另一方面，TSLS 表明 RBC、HGB 和 HCT 对 DBP 有显着影响，但对 SBP 没有显着影响，而 CIIV 估计对 DBP 有正向影响，对 SBP 有负向影响。

比较 RGF 和概率精细映射方法（如 Susie）确定的因果变异非常有趣。一般来说，RGF侧重于基因组全球景观，而精细定位方法则更侧重于局部特征。以 RBC 为例，如果我们在初步筛选后使用这两种方法对所有候选基因座进行粗略的精细映射，我们会发现大多数信号在两种方法之间是共享的（图 6f）。当我们专门检查GWAS峰附近的局部结构时（例如，6号染色体上从41,600,000到42,200,000的区域），尽管Susie倾向于识别更多的因果位点，但两种方法的最显着信号（rs112233623）是相同的（图6g）。考虑到在IV选择中，我们更关注GWAS中由遗传相关性和“赢家诅咒”引起的假阳性信号，RGF的相对保守性并不是劣势。

潜在的机制可能与血液粘度有关。较高的 RBC、HGB 和 HCT 意味着血液粘度和外周血流阻力增加，导致高血压。此外，红细胞和血红蛋白还影响一氧化氮生物利用度，这是血管舒张、血栓形成抑制和血管形成等血管心理调节的关键信号。尽管分子机制仍有待发现，但研究结果表明，这些血液学指标可能不仅是高血压的指标，而且是高血压的潜在治疗靶点。

BNMR表明，白细胞增加会导致患抑郁症的风险

近年来，神经免疫相互作用一直是一个吸引人的话题。两个系统之间存在广泛的双向电路。神经系统通过交感神经和副交感神经的直接连接调节免疫活性和细胞因子平衡，一些神经递质和神经内分泌激素也可以作为免疫调节剂。同时，免疫系统在发育过程中参与突触的消除和可塑性，并调节大脑活动。免疫相关血液生物标志物为许多精神疾病的病理机制提供了新的见解。例如，免疫失调长期以来一直被认为与包括抑郁症在内的心理障碍有关。最近的研究报告了白细胞计数与抑郁症之间的相关性。

我们利用英国生物样本库的疾病记录，通过随机选择相同数量的相同种族的健康个体来评估白细胞计数及其两种亚型淋巴细胞和单核细胞计数是否会因果影响抑郁症，从而构建病例对照研究。排除极值超过的受试者，分析包括 22,324 例病例和 22,861 例对照。所有参与者都是高加索血统。