A robust cis-Mendelian randomization method with application to drug target discovery
Nature Communications
摘要
这篇文章提出了一种新的稳健的顺式孟德尔随机化（cis-MR）方法，称为cisMR-cML，用于药物靶点发现。孟德尔随机化（MR）使用基因变异作为工具变量来推断特征之间的因果关系，而cis-MR则专注于单一基因组区域，只使用顺式SNPs（cis-SNPs）。这种方法的一个重要应用是药物靶点的发现，例如使用蛋白质的cis-pQTL作为暴露变量来评估疾病的因果效应，从而为药物靶点发现提供一种经济有效的途径。
cisMR-cML方法基于约束最大似然估计，能够处理顺式SNPs的基因多效性和连锁不平衡（LD）问题，相较于现有方法在工具变量假设的违反情况下表现出更强的鲁棒性。实例分析表明，cisMR-cML在模拟研究和冠状动脉疾病（CAD）药物靶点分析中的表现优于其他现有的方法，发现了三个潜在的CAD药物靶点：PCSK9、COLEC11和FGFR1。
引言
现有研究：
文章提到了MR方法，这是一种利用遗传变异作为工具变量（instrumental variables, IVs）来研究特征之间因果关系的统计方法。这种方法的优势在于遗传变异在概念上是随机分配的，从而减少了混杂因素和逆向因果关系的偏倚。 研究存在的问题：
传统的MR方法主要关注全基因组显著的SNPs与暴露的关联，但这种方法可能无法有效处理顺式SNPs（cis-SNPs）的多效性（pleiotropy）和连锁不平衡（linkage disequilibrium, LD）。 文章指出，现有的cis-MR方法在处理顺式SNPs时存在局限性，因为这些SNPs通常位于基因组的同一区域，存在较强的LD，这可能导致IV假设的违反。 为什么要做这项研究：
为了解决现有方法在处理顺式SNPs时的不足，作者提出了一种新的方法——cisMR-cML，这是一种基于约束最大似然估计的方法，它对IV假设的违反具有鲁棒性，并有强大的理论支持。 该方法旨在更准确地识别和利用顺式SNPs作为IVs，从而提高MR分析的准确性和可靠性。 能改变什么：
cisMR-cML方法的提出，可以改变我们对药物靶点发现的理解和应用。通过更精确地估计基因表达与疾病之间的因果关系，这种方法有助于识别潜在的药物靶点，特别是在冠心病（CAD）等疾病的研究中。 该方法的应用可以提高我们对遗传学和疾病之间复杂关系的理解，为个性化医疗和精准治疗提供科学依据。
方法
cisMR-cML 方法的核心是如何基于约束最大似然估计（cMLE）来估计暴露变量（Exposure）和结局变量（Outcome）之间的因果效应。cisMR-cML 特别适用于存在无效工具变量（Instrumental Variables, IVs）的情况，即一些 SNP 对结局的直接效应未通过暴露变量中介。这一方法包括了从基本模型到对数似然函数的推导，并通过约束最大化方法来估计因果效应。
以下是方法部分的详细步骤解析和推导。
1.基本模型设定
假设我们有 个 SNP 作为 IV，每个 SNP 与暴露 和结局 都有不同的效应。模型的基本设定为：
暴露的模型：其中 表示第 个 SNP 的基因型， 是第 个 SNP 对暴露 的效应， 是独立于 SNP 的随机误差项。
结局的模型：其中 是因果效应， 是第 个 SNP 对结局的直接效应，未通过暴露变量中介， 是独立于 SNP 的随机误差项。
2.联合模型的构建
通过将暴露模型（1）代入结局模型（2），我们可以得到一个联合模型：
记，则上式可以简化为：
其中 是新的误差项。
3.估计暴露和结局的联合效应
“Transformation between Marginal and Conditional SNP-Effect Estimates” 这部分的内容讨论了如何在孟德尔随机化（MR）分析中将 SNP 对暴露和结局的边际效应估计转化为条件效应估计。这个转换是因为边际效应估计无法准确捕捉 SNPs 的共线性结构或连锁不平衡（LD）信息，而条件效应则考虑了不同 SNPs 之间的相关性，从而提供了更精确的估计。这种转换在 cisMR-cML 方法中尤为关键，因为该方法依赖于 SNPs 的条件效应估计来计算因果效应。 以下是这部分内容的逐步详细讲解和推导： 在 GWAS（基因组关联研究）中，边际效应估计是每个 SNP 单独与暴露或结局变量之间的关联，忽略了其他 SNP 的影响。这意味着每个 SNP 的边际效应估计只考虑了它对暴露或结局的独立贡献。然而，由于 SNPs 之间可能存在 LD，边际效应并不能准确反映 SNP 对暴露或结局的真实影响。 相反，条件效应估计通过调整其他 SNPs 的影响来得到每个 SNP 的条件效应，从而更好地考虑了 SNPs 之间的相关性（如 LD 结构）。这种条件效应估计可以通过 LD 矩阵进行计算。 假设我们有暴露 和结局 的 GWAS 边际效应估计为：其中： 在本节中，目标是将这些边际效应估计转换为条件效应估计，以便准确反映每个 SNP 对暴露和结局的联合效应。 从暴露的模型可以得出：其中： 在引入 LD 矩阵后，条件效应估计可以用如下关系表示：其中： 通过这种转换，可以得到条件效应估计，即暴露 中每个 SNP 对结局的影响，同时考虑 SNPs 之间的相关性。 为了进一步准确估计条件效应，还需要计算条件效应估计的协方差矩阵。暴露效应的条件估计协方差矩阵 可以通过边际估计的协方差矩阵进行转换： 其中： 在实际应用中，如果无法获得 LD 矩阵，则可以通过公开的参考样本（如 1000 Genomes 项目）中的 LD 矩阵 近似替代。 当 GWAS 中的边际效应不是在标准化的基因型和表型上计算时，可以通过以下方法近似得到标准化的效应估计： 通过以上近似，cisMR-cML 可以在没有标准化的情况下依然获得合理的条件效应估计。 同样的步骤可以应用于结局 的条件效应估计。我们可以类似地表示为： 其中 是结局的条件效应估计， 是结局 GWAS 样本的 LD 矩阵，。 在 cisMR-cML 方法中，选择条件效应估计是因为边际效应在存在 LD 的情况下可能不满足有效 IV 的要求。例如： cisMR-cML 方法中的“边际效应与条件效应估计之间的转换”步骤，通过 LD 矩阵将 SNP 对暴露和结局的边际效应估计转换为条件效应估计。这样可以确保 SNP 的估计值在存在 LD 时更为准确，同时提高因果推断的可靠性。这种转换是 cisMR-cML 方法中实现有效 IV 筛选和因果效应估计的关键步骤。3.1 边际效应与条件效应
3.2 边际效应与条件效应之间的转换
3.3 条件效应估计的推导
3.4 协方差矩阵的转换
3.5 近似条件效应的估计
3.6 对结局的条件效应估计
3.7 重要性：条件效应 vs. 边际效应
总结
在 cisMR-cML 方法中，我们将 SNP 对暴露 和结局 的效应分别看作是多变量正态分布的估计值：
其中：
和 分别是 和 的协方差矩阵。
根据模型的设定，结局的效应可以表示为暴露的效应和直接效应的线性组合，即：
其中 表示 SNP 对结局的直接效应向量。
4.对数似然函数的构建
基于上述模型的设定，给定观测到的 GWAS 估计 和，可以构建以下对数似然函数：
其中：
第一项 表示 SNP 对暴露的估计误差。 第二项 表示 SNP 对结局的估计误差。
5.约束条件
在 cisMR-cML 中，假设有 个无效工具变量，即在模型中存在 个 SNP 满足。因此，我们对 的稀疏性进行约束：
其中 是指示函数，用于标记哪些 SNP 对结局有直接效应。
6.约束最大似然估计的优化问题
通过最大化对数似然函数，在给定约束条件下估计、 和，即：
这个优化问题本质上是一个最佳子集选择问题（Best Subset Selection Problem），目标是在约束条件下找到最优的 和。
在求解一个子集选择问题时，为了估计无效工具变量（invalid IVs）的集合 和因果参数，所使用的几种不同算法的选择及其效果。这些算法包括一种启发式算法（heuristic algorithm）、拼接算法（splicing algorithm）和带 Lasso 正则化的最小化算法。每种算法在理论上或实践中各有优缺点，以下是对每种方法的详细解释和比较。 在此上下文中，问题是一个最佳子集选择问题（Best-Subset Selection Problem），其目的是从一组候选的工具变量（IVs）中识别出无效 IVs。无效 IVs 是具有水平多效性（pleiotropy）的 SNPs，即那些通过其他路径直接影响结局而不是通过暴露的 SNPs。无效 IVs 可以表示为非零的，因此问题可以转化为一个求解最优子集以满足条件的约束优化问题。 由于工具变量的数量 较多时计算代价很高，因此直接获得全局解（global solution）是计算上不现实的。因此提出了几种近似或替代的算法来解决该问题。 算法原理： 优缺点： 适用场景：启发式算法适合快速分析，特别是在时间和计算资源有限的情况下非常实用。 算法原理： 优缺点： 适用场景：拼接算法适合需要高精度结果且计算资源充足的场景，特别是在 IV 数量较少或计算效率要求不高的情况下。 算法原理： 优缺点： 适用场景：Lasso 方法适合处理较大规模的数据集，尤其在目标是筛选大量候选变量中的少量有效变量时较为有效。 在本文的数值实验中，启发式算法表现良好，并且具有较高的计算效率。因此，研究者最终选择了启发式算法来估计无效 IVs 集合和因果参数 。尽管拼接算法能更接近全局最优解，但其计算代价过高，不适合较大规模的 IV 集合。同样，尽管 Lasso 正则化能够进行自动变量选择，但在本实验中效果不如启发式算法，表现出次优的结果。因此，为了在计算效率和精度之间取得平衡，最终选择了启发式算法。算法讲解
6.1 问题背景：最佳子集选择问题
6.2 启发式算法（Heuristic Algorithm）
6.3 拼接算法（Splicing Algorithm）
6.4 带 Lasso 正则化的最小化算法（Lasso Penalty Method）
6.5 算法的对比总结
6.6 算法选择的原因
7.使用 BIC 选择最优 值
为了确定合适的无效工具变量数量，使用贝叶斯信息准则（BIC） 作为模型选择标准：
其中 是样本大小的最小值。最终的 选择为使得 BIC 最小化的值：
贝叶斯信息准则（BIC）是一种广泛使用的模型选择标准，用于评估和选择统计模型。本文中BIC 被用于确定合适的无效工具变量数量 。下面将详细讲解这个公式及其组成部分。 似然函数 ： 惩罚项： 整体结构： BIC 在这个上下文中的作用是通过评估不同数量的无效工具变量对模型拟合的影响，从而选择出最佳的模型。通过结合拟合优度和复杂度的考虑，BIC 提供了一种有效的方式来平衡模型的准确性与简单性，进而提高模型在新数据上的泛化能力。这种方法在许多统计和计量经济学研究中都有广泛应用。BIC
BIC公式
组成部分解析
最优化选择
总结
8.标准误差的估计
一旦确定了，可以基于Fisher 信息矩阵 估计其标准误差。最终通过数据扰动技术（如 Bootstrap）获得更稳定的估计结果和统计推断。
总结
cisMR-cML 方法的核心在于构建对数似然函数，并在无效工具变量数量 的约束下进行优化，最终通过选择合适的 值来获得因果效应的约束最大似然估计。这一方法的提出解决了孟德尔随机化分析中因无效工具变量而引起的偏差问题，并通过多步计算得到了稳健的因果效应估计。
Box 1 通常用于列出算法的主要步骤，以帮助读者更清晰地理解实现的过程。根据你的需求，以下是 Box 1 的一个示例，详细描述了在估计无效工具变量（IVs）集合和因果参数 时的算法步骤。这个算法涉及到逐步选择和排除无效 IVs 的过程，以最终确定无效工具变量集合。 Box 1：Heuristic Algorithm for Estimating Invalid IVs and Causal Parameter Input: GWAS summary statistics,, LD matrices,, maximum number of invalid IVs, initial values Initialize:[初始化] Iterative Update:[迭代更新] Update. [更新。] Update residuals:. Update invalid IV set:. [重新计算残差，并更新无效 IV 集合。] Set,,. [将最新的参数更新为最终输出。] Update parameters. [对给定的无效 IV 集合，在约束条件下最大化对数似然，更新。] Set. Set. [将 更新为。] Maximize under the constraint for: If: [如果新的对数似然值 高于之前的，则进行以下更新：] Return:[输出结果] 该算法通过逐步迭代更新无效 IV 集合来接近最优解。初始选择具有最大残差的 IV 作为无效 IVs，并在每次迭代中重新计算残差，更新无效 IV 集合。这样，算法逐渐逼近对数似然最大化的解，从而得到因果参数的稳定估计。这种方法是一种启发式算法，在计算效率和结果准确性之间取得平衡，适合处理具有多效性的无效工具变量筛选问题。 本文选择使用启发式算法是因为最佳子集选择问题本身计算代价高，特别是在存在大量候选工具变量（IVs）的情况下。我们可以分析几种可能的替代算法，包括拼接算法（Splicing Algorithm）和带 Lasso 正则化的最小化算法（Lasso Penalty Method），并讨论它们各自的优缺点和实现步骤。 启发式算法的主要优势在于计算效率高，适用于大型数据集。但由于它并不系统性地探索所有可能的变量组合，因此可能会收敛到局部最优解，而不是全局最优解。这在模型复杂或噪声较大的情况下可能影响准确性。 以下是两个替代算法的选择，它们可以解决局部最优解的问题，但各自也有不同的挑战。 拼接算法是一种全局优化算法，它通过系统地搜索变量的组合来找到一个全局最优解。这种方法能更准确地识别无效 IVs，但由于计算代价高，通常适合小规模的数据集。 步骤： 优缺点： Lasso 正则化方法通过在优化过程中引入惩罚项，将效应较小或无效的变量系数（即无效 IVs）压缩到接近零，从而达到选择有效变量的目的。这种方法适合高维数据，但在存在强相关性时，可能会导致变量选择不准确。 步骤： 设置目标函数：定义包含 Lasso 惩罚项的目标函数：其中 是惩罚系数，控制对 的压缩力度。 初始化：选择合适的 值，根据经验或通过交叉验证确定。 迭代优化：使用梯度下降或坐标下降等优化算法，迭代更新参数 和，在每一步更新中将较小的 压缩至零。 变量选择：在收敛时，取所有 的 IVs 作为无效 IV 集合，(r_i = 0$ 的 IVs 则视为有效工具变量。 输出：返回最终的无效 IV 集合及估计的因果参数。 优缺点： 启发式算法在这一问题上被优先选择的原因在于： 在最佳子集选择问题上，不同算法各有优缺点。拼接算法虽然能够找到全局最优解，但计算代价太高；Lasso 正则化方法计算效率高，但存在选择偏差。相较之下，启发式算法在大规模数据下更具优势，能够较快逼近合理的解。因此，该问题最终选择了启发式算法，以在准确性和效率之间取得平衡。box1
算法的直观解释
算法选择问题
① 启发式算法的局限性和原因
② 替代算法
A. 拼接算法（Splicing Algorithm）
B. 带 Lasso 正则化的最小化算法（Lasso Penalty Method）
③ 为什么选择启发式算法
④ 总结
模拟
在这篇文章中，作者进行了两组模拟研究，目的是验证所提出的cisMR-cML方法在各种不同情境下的性能。以下是关于模拟研究的详细讲解，包括数据生成、多种情境的考虑和具体设置。
第一组模拟研究
第一组模拟研究中，作者模拟了具有不同强度LD结构的遗传数据，目的是比较cisMR-cML与其他现有方法在不同情境下的表现。以下是具体的步骤和情境设置：
数据生成：
SNPs的生成：作者生成了10个SNPs，SNP之间的LD模式遵循自回归（autoregressive）LD结构，相关系数分别为弱（ρ = 0.2）、中等（ρ = 0.6）和强（ρ = 0.8）三种水平。 两种场景： 场景1：所有10个SNP对暴露均有作用，即 |IX| = 10。 场景2：只有一半的SNP对暴露有作用，即 |IX| = 5，且IY ∩ IX 的大小为2。 工具变量的选择：
LD Clumping：在某些分析中，使用LD Clumping的策略获得部分独立的工具变量，再应用于现有的稳健MR方法。但作者发现这种方法会导致工具变量数量的减少，进而影响方法的统计功效。 方法比较：
使用的cis-MR方法： 独立IV版本的方法： 多基因型MR方法： cisMR-cML：在条件效应的基础上，计算所有10个SNP的条件估计，并进行稳健分析。 LEgger：条件估计的MR-Egger方法。 GIVW和GEgger：使用边际效应计算。 IVW-IND、Egger-IND和cML-IND：使用独立工具变量的IVW、Egger和MR-cML方法。 包括MR.LDP、MR.Corr2、MR.CUE和MRAID。 模拟结果：
在场景1中（所有10个SNP均对暴露有影响），当所有10个工具变量均为有效工具变量时（即没有无效工具变量），所有方法都能很好地控制I类错误率。cisMR-cML相对其他方法更加保守。 当存在4个无效工具变量时，只有cisMR-cML能够有效控制I类错误率并保持较高的功效。 GIVW和GEgger随着SNP之间相关性的增加，I类错误率显著增加。而其他使用独立工具变量的方法也出现了I类错误率膨胀的情况。 使用边际效应与条件效应的对比：
文章还进行了cisMR-cML与使用边际效应估计的cisMR-cML-Marg的对比，发现后者由于违背了多数条件而导致I类错误率膨胀。
第二组模拟研究
在第二组模拟研究中，作者使用了更具现实性的LD结构来生成数据，目的是评估方法在接近真实数据情况下的表现。
数据生成：
真实LD模式的应用：通过从UK Biobank的基因型数据中选取cis区间的SNPs，生成暴露和结局数据。这些数据分别用于暴露和结局的GWAS，之后使用COJO分析来选择与暴露和结局相关的SNPs。 两种场景： 场景1：从蛋白质组分析中随机选择50种蛋白质，这些蛋白质至少有五个SNP与暴露相关，并且至少有一个SNP与结局相关。 场景2：随机生成50种蛋白质及其对应的基因效应。 方法比较：
使用的方法： cisMR-cML和LEgger：使用COJO选择的SNPs进行条件估计。 cisMR-cML-X、GIVW-X和GEgger-X：仅使用与暴露相关的SNPs进行分析。 IVW-IND和其他多基因型MR方法（MR.LDP、MR.Corr2、MR.CUE和MRAID）。 模拟结果：
在场景1和场景2中，cisMR-cML是唯一有效控制I类错误且具有最小均方根误差（RMSE）的方法。 LEgger的功效较低，使用仅与暴露相关的SNPs的其他方法则表现出不良的表现，包括I类错误率膨胀。
总结
数据生成：使用了自回归LD结构和基于真实LD模式的模拟。 情境设置：考虑了不同强度的LD相关性和不同数量的无效工具变量，以考察方法在各种情形下的表现。 方法比较：cisMR-cML相比其他方法表现出更强的鲁棒性，特别是在处理水平多效性和连锁不平衡的情况下。
实例
在这篇文章的实例研究中，作者使用了cisMR-cML方法进行了两个实例分析，旨在展示其在真实数据中的应用，并比较现有方法的性能。
实例研究 1：使用下游生物标志物评估CAD的药物靶点
1. 数据来源
暴露和结局数据： 暴露：低密度脂蛋白胆固醇（LDL），其与PCSK9蛋白的浓度和活性有关。 结局：冠状动脉疾病（CAD）的发病风险。 使用GCTA-COJO方法在PCSK9基因区域选择了9个SNPs，其中8个与LDL相关，1个与CAD相关。
2. 研究目的
使用PCSK9基因区域的遗传变异，评估LDL水平与CAD风险之间的因果关系。 确认已有研究发现，即PCSK9抑制剂可以通过降低LDL水平来降低CAD风险。
3. 现有研究的缺陷
现有的MR分析通常只使用与暴露相关的工具变量（例如与LDL水平显著相关的SNPs），而忽略了可能与结局（例如CAD风险）相关的工具变量，这可能会导致多效性偏倚。 在许多现有的cis-MR方法中，对SNPs之间存在的连锁不平衡（LD）未能很好地处理，这会导致分析中I类错误率膨胀，尤其是在SNPs之间存在较强相关性时。
4. 结果
因果估计结果：cisMR-cML、LEgger、GIVW-X和GEgger-X均提示LDL对CAD风险具有显著的正因果效应。各自的p值分别为 7.3 × 10^-4、9.2 × 10^-3、8.6 × 10^-8 和 0.02。 对比结果： cisMR-cML：通过使用与暴露和结局相关的SNPs并建模条件效应，能更好地避免多效性偏倚，从而得出稳健的因果推断。 其他方法：仅使用与暴露相关的SNPs进行分析的方法在结果中存在一定程度的偏倚，可能因为未包含与结局相关的SNPs而导致多效性效应通过连锁不平衡传播到暴露。
5. 新发现
PCSK9蛋白：cisMR-cML分析结果表明LDL对CAD的显著因果效应，验证了现有关于PCSK9抑制剂降低CAD风险的研究。该方法有效确认了LDL作为潜在药物靶点的因果影响，并支持PCSK9抑制剂的治疗价值。
6. 讨论
研究显示，使用cisMR-cML方法，在使用下游生物标志物作为暴露的情况下，可以更有效地识别与疾病相关的潜在因果关系，从而为药物靶点的验证和再定位提供支持。
实例研究 2：CAD的蛋白质组广泛分析
1. 数据来源
暴露和结局数据： 暴露：从ARIC欧洲血统（EA）队列中提取的蛋白质水平数据，样本量为7213。 结局：冠状动脉疾病（CAD）风险。 使用GCTA-COJO对蛋白质组中的773种蛋白进行了pQTL（蛋白质定量性状位点）分析。
2. 研究目的
进行蛋白质组的全范围扫描，以评估多个蛋白质的表达水平对CAD风险的因果作用。 目标是识别潜在的CAD药物靶点，并评估其与CAD的因果关系及共定位证据。
3. 现有研究的缺陷
现有的cis-MR方法在处理蛋白质组数据时通常只关注与暴露相关的SNPs，可能忽略了与结局相关的遗传变异，从而导致多效性偏倚的存在。 由于蛋白质与SNPs之间存在复杂的相关性，传统的方法在处理这些数据时往往功效不足或者存在高I类错误率的问题。
4. 结果
蛋白质-CAD因果关系：
显著蛋白质靶点：在FDR（假发现率）小于0.05的情况下，cisMR-cML识别出三个蛋白质与CAD风险之间存在因果关系：PCSK9、COLEC11和FGFR1。其中PCSK9和COLEC11的因果变异具有共定位证据（H4-PP≥0.7），表明它们可能是CAD的潜在靶点。 FGFR1蛋白：虽然共定位分析未能找到与CAD有关的显著因果变异（H1-PP约为96%），但该基因与心脏病的潜在联系仍然值得进一步研究。 对比结果：
GIVW-X识别了18种蛋白，其中只有4种具有共定位证据，且有5种蛋白质可能因LD与CAD有关，从而违反了MR假设。 其他方法如GEgger-X和Wald-ratio也存在I类错误率膨胀的问题，而LEgger并未识别出任何显著结果。
5. 新发现
COLEC11蛋白：与CAD相关，参与了凝集素补体激活途径，对先天免疫系统有重要作用，可能成为新的CAD治疗靶点。 FGFR1蛋白：虽然在共定位分析中没有找到明确的因果变异，但其在细胞增殖和血管生成中的作用使其成为未来潜在的研究对象。
6. 讨论
多效性与共定位：cisMR-cML在使用条件效应模型、包含与结局相关的SNPs的情况下，表现出稳健的因果推断能力，并能够有效控制I类错误率。这一点在蛋白质组广泛扫描中得到了体现。 方法的有效性：cisMR-cML相较于其他方法表现出更好的稳健性，尤其是在处理复杂的SNP相关结构和多效性时，其因果推断结果更具可信度。
总结
通过这两组实例分析，cisMR-cML方法在处理复杂基因与表型之间的因果关系时表现出显著的优势，特别是在解决多效性和遗传相关性问题方面。
Lin, Z., Pan, W. A robust cis-Mendelian randomization method with application to drug target discovery. Nat Commun 15, 6072 (2024). https://doi.org/10.1038/s41467-024-50385-y参考文献
