多组学因果推断实操：孟德尔随机化与中介效应建模

实操技巧：通过 IEU GWAS 数据库的 extract_instruments() 函数可直接获取标准化工具变量，无需本地存储几十 G 原始数据，代码如下（R 语言）：

# 安装并加载 TwoSampleMR 包
install.packages("TwoSampleMR")
library(TwoSampleMR)
# 从 IEU 数据库提取免疫表型（CD39+Treg）的工具变量
exposure_dat <- extract_instruments(outcomes = "ieu-a-1234", p1 = 1e-5, clump = TRUE)
# 查看工具变量基本信息
head(exposure_dat[, c("snp", "beta", "se", "pval")])

（二）Step2：数据预处理与质量控制

1. 工具变量筛选

遵循'三步骤筛选法'：

1. 关联性筛选：保留与暴露 / 中介显著关联的 SNP（P＜1e-5），例如 CD39⁺Treg 表型筛选出 26 个有效 SNP。
2. 独立性筛选：采用 LD 剪枝（clump 参数设置 r²＜0.01，距离＞1000kb），排除连锁不平衡的 SNP。
3. 有效性筛选：计算 F 统计量，剔除 F＜10 的弱工具变量，确保 βₓₘ估计无偏。

2. 数据标准化与合并

• 效应量统一：将所有数据的效应量转换为'每增加 1 个标准差暴露 / 中介，结局的变化量'。
• 等位基因对齐：通过 harmonise_data() 函数统一 SNP 的等位基因方向，避免链翻转导致的效应量偏差：

import pandas as pd
import numpy as np
from TwoSampleMR.py import harmonise

# 读取暴露、中介、结局数据
exposure = pd.read_csv("exposure_snps.csv")
mediator = pd.read_csv("mediator_snps.csv")
outcome = pd.read_csv("outcome_snps.csv")

# 等位基因对齐
harmonised_exp_med = harmonise.harmonise_data(exposure, mediator)
harmonised_med_out = harmonise.harmonise_data(mediator, outcome)

（三）Step3：孟德尔随机化核心分析

采用'三步法'验证因果关系链：暴露→中介、中介→结局、暴露→结局。

1. 暴露对中介的 MR 分析（X→M）

以 CD39⁺Treg（暴露）对 GPI（中介）的效应分析为例，采用多种 MR 方法交叉验证：

# 1. 加载数据
med_dat <- extract_outcome_data(snps = exposure_dat$snp, outcomes = "ieu-b-5678")
dat_xm <- harmonise_data(exposure_dat, med_dat)

# 2. 多方法分析
mr_results_xm <- mr(dat_xm, method_list = c("ivw", "mr_egger_regression", "weighted_median"))

# 3. 结果输出
print(mr_results_xm[, c("method", "b", "se", "pval", "or")])

关键结果解读：若 IVW 法显示 βₓₘ=0.23（95% CI：0.11~0.35，P=0.002），表明 CD39⁺Treg 水平每增加 1SD，GPI 浓度显著升高 23%，满足中介效应的前提条件。

2. 中介对结局的 MR 分析（M→Y）

验证 GPI 对 AA 的因果效应，需控制暴露因素的混杂：

import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor

# 构建多变量回归模型（控制暴露效应）
X = harmonised_med_out[["beta_med", "beta_exp"]] # 中介效应量 + 暴露效应量
X = sm.add_constant(X)
y = harmonised_med_out["beta_out"] # 计算 VIF 值，排除共线性（VIF＜5 为合格）
vif = pd.DataFrame()
vif["变量"] = X.columns
vif["VIF"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]

# 拟合模型
model_my = sm.OLS(y, X).fit()
print(model_my.summary())

3. 暴露对结局的总效应分析（X→Y）

采用 IVW 法作为核心方法，MR-Egger 和加权中位数法作为敏感性验证：

# 提取结局数据
out_dat <- extract_outcome_data(snps = exposure_dat$snp, outcomes = "ieu-c-9012")
dat_xy <- harmonise_data(exposure_dat, out_dat)

# 总效应估计
mr_results_xy <- mr(dat_xy, method_list = c("ivw", "mr_egger_regression", "weighted_median"))

# 森林图可视化
mr_forest_plot(mr_results_xy, dat_xy)

（四）Step4：中介效应建模与量化

采用乘积法计算中介效应，公式为：中介效应 = βₓₘ × βₘᵧ；中介比例 =（βₓₘ × βₘᵧ）/ βₓᵧ × 100%。

1. 单中介效应计算

# 提取关键效应量
beta_xm <- mr_results_xm$b[mr_results_xm$method == "ivw"] # 暴露对中介的效应
beta_my <- model_my.params["beta_med"] # 中介对结局的效应
beta_xy <- mr_results_xy$b[mr_results_xy$method == "ivw"] # 暴露对结局的总效应

# 计算中介效应及比例
mediation_effect <- beta_xm * beta_my
mediation_proportion <- (mediation_effect / beta_xy) * 100

# 输出结果
cat(paste0("中介效应：", round(mediation_effect, 3), "\n"))
cat(paste0("中介比例：", round(mediation_proportion, 2), "%\n"))

案例结果：如 GPI 在 CD39⁺Treg 与 AA 间的中介效应为 - 0.003，中介比例为 - 8.80%（95% CI：-16.00%~-1.63%，P=0.016），表明 GPI 在该通路中发挥负向中介作用。

2. 多中介效应分析

当存在多个代谢物作为潜在中介时，采用多变量 MR（MVMR）模型控制中介间的相互干扰：

# 加载 MVMR 包
install.packages("MVMR")
library(MVMR)

# 构建多中介数据集
mv_data <- format_mvmr(exposure_dat = exposure_multi, mediator_dat = mediator_multi, outcome_dat = out_dat)

# 拟合多中介模型
mv_result <- mvmr(mv_data, method = "ivw")

# 提取各中介的效应量
mediation_effects <- mv_result$beta_mediator

（五）Step5：敏感性分析与结果稳健性验证

敏感性分析是 MR 研究的'生命线'，需从异质性、多效性、稳定性三方面展开：

1. 异质性检验

采用 Cochran's Q 检验评估工具变量间的效应异质性：

# 异质性检验
heterogeneity_xm <- mr_heterogeneity(dat_xm, method = "ivw")
cat(paste0("Cochran's Q 统计量：", round(heterogeneity_xm$q_stat, 2), "\n"))
cat(paste0("P 值：", round(heterogeneity_xm$q_pval, 3), "\n"))

判断标准：P＞0.05 表明无显著异质性，模型拟合良好。

2. 水平多效性评估

• MR-Egger 截距检验：若截距项 95% CI 覆盖 0 且 P＞0.05，提示无定向多效性。
• MR-PRESSO 检验：识别并剔除异常值 SNP，验证结果稳定性。

# MR-PRESSO 检验（Python 实现）
from pymrpress import mrp

# 执行检验
presso_result = mrp.mr_presso(dat_xy["beta_exp"], dat_xy["beta_out"], dat_xy["se_exp"], dat_xy["se_out"])

# 查看结果
print("全局异常值检验 P 值:", presso_result.global_pval)
print("校正后效应量:", presso_result.corrected_beta)

3. 稳定性检验

采用**留一法（Leave-One-Out）**分析单个 SNP 对总效应的影响：

# 留一法分析
loo_result <- mr_leaveoneout(dat_xy)

# 可视化留一法结果
mr_leaveoneout_plot(loo_result)

若剔除任一 SNP 后效应量仍稳定在 95% CI 内，表明结果不受单个工具变量影响。

四、进阶实战：多组学整合的典型案例解析

以'SGLT1 抑制→血浆蛋白 / 代谢物→老年人衰弱'研究为例，展示多组学中介分析的完整应用。

（一）研究设计亮点

1. 暴露变量创新：以 SLC5A1 基因变异模拟 SGLT1 抑制效应，避免药物干预的伦理限制；
2. 多组学中介整合：同时纳入 1558 种血浆蛋白和 1352 种代谢物作为中介候选；
3. 分层分析策略：按性别和衰弱诊断标准（FNIH/EWGSOP）分层验证，提升结果可靠性。

（二）核心发现解读

1. 总效应验证：SGLT1 抑制与老年人虚弱指数（FI）降低显著相关（IVW: OR=0.82, 95% CI:0.71~0.94, P=0.005）；
2. 蛋白中介通路：神经内分泌转化酶 1（PCSK1）中介比例达 13.49%，提示其在代谢调控中的关键作用；
3. 代谢物中介通路：棕榈酰肉碱同时介导 SGLT1 抑制对 FI 和低握力的效应，中介比例 8.35%。

（三）关键代码片段：多组学中介筛选

# 批量筛选蛋白中介
protein_mediators <- list()
for (i in 1:nrow(protein_list)) {
    # 提取单个蛋白的工具变量
    med_snps <- extract_instruments(outcomes = protein_list$id[i], p1 = 1e-5)
    
    # 验证 X→M 和 M→Y 效应
    if (nrow(med_snps) > 2) { # 至少 3 个工具变量
        dat_xm <- harmonise_data(exposure_dat, med_snps)
        res_xm <- mr(dat_xm, method = "ivw")
        dat_my <- harmonise_data(med_snps, out_dat)
        res_my <- mr(dat_my, method = "ivw")
        
        # 筛选显著中介
        if (res_xm$pval < 0.05/1558 & res_my$pval < 0.05/1558) {
            protein_mediators[[i]] <- data.frame(
                protein = protein_list$name[i],
                beta_xm = res_xm$b,
                beta_my = res_my$b
            )
        }
    }
}

# 合并结果
mediator_results <- do.call(rbind, protein_mediators)

五、常见问题与解决方案

（一）工具变量相关问题

1. 弱工具变量偏倚
   • 问题：F＜10 导致效应量估计偏倚；
   • 解决方案：降低 P 值阈值（如从 1e-5 降至 1e-6），或合并多个 GWAS 数据集提升样本量。
2. 工具变量数量不足
   • 问题：部分代谢物 / 蛋白对应的 pQTL 数量＜3；
   • 解决方案：采用基因水平工具变量（如以基因内所有 SNP 构建多基因风险评分），或参考 GTEx 等数据库补充 eQTL 数据。

（二）多组学数据整合问题

1. 数据异质性
   • 问题：不同数据库的人群 ancestry、检测平台差异导致结果不一致；
   • 解决方案：限定 ancestry 为欧洲人群（减少种群分层），采用 ComBat 方法校正批次效应。
2. 多重检验负担
   • 问题：上万种多组学分子分析导致假阳性升高；
   • 解决方案：采用分层校正策略，先对组学层面进行 FDR 校正（如蛋白组 FDR＜0.1），再对显著分子进行 Bonferroni 校正。

（三）中介效应解读问题

1. 负向中介效应
   • 问题：中介比例为负值（如 - 8.80%）难以理解；
   • 解决方案：结合生物学背景，可能提示中介因子具有'缓冲效应'，如 GPI 通过抑制炎症反应抵消部分 CD39⁺Treg 的促 AA 效应。
2. 中介效应与因果机制的区别
   • 问题：将统计中介等同于生物学机制；
   • 解决方案：需通过湿实验验证，如在细胞模型中干预 GPI 水平，观察 CD39⁺Treg 对造血细胞凋亡的影响。

六、前沿拓展：多组学因果推断的未来方向

1. 动态中介分析：结合纵向多组学数据，捕捉中介效应随时间的变化，例如青少年时期的代谢物如何介导肥胖对老年心血管疾病的风险。
2. 多模态整合：融合影像组学数据（如脑结构影像），解析'基因→蛋白→脑结构→认知下降'的跨尺度因果通路。
3. 机器学习赋能：采用自编码器对高维多组学数据降维，或用随机森林筛选关键中介因子，减少人工筛选偏差。
4. 药物靶点优先级排序：通过中介比例量化靶点价值，例如中介比例＞30% 的蛋白（如 PON3，中介比例 44.90%）可作为优先开发的药物靶点。

七、结语

多组学孟德尔随机化与中介效应建模的结合，打破了传统关联分析的局限性，为解析复杂疾病机制提供了'因果透镜'。从工具变量筛选到敏感性验证，从单中介到多组学整合，每一步都需要严格的方法学规范与生物学逻辑支撑。

本指南提供的分析框架与代码，可直接应用于免疫、代谢、肿瘤等多个研究领域。但需牢记：统计因果推断是生成假设的工具，而非验证机制的终点。只有将 MR 发现与湿实验验证、临床观察相结合，才能真正实现从'数据'到'机制'再到'干预'的转化，为精准医学研究注入新动能。