孟德尔随机化完整教程

🧬

孟德尔随机化（Mendelian Randomization, MR） 是一种利用遗传变异作为工具变量来推断暴露与结局之间因果关系的方法。由于基因型在受精时随机分配且不受环境混淆因素影响，MR被视为"自然的随机对照试验"。

📑 目录

MR基本原理

🎯 核心思想

利用基因型作为工具变量推断暴露对结局的因果效应：

基因型 (G) → 暴露 (X) → 结局 (Y)

MR三大假设

1 关联性假设

工具变量（SNP）与暴露（X）强相关

2 独立性假设

SNP与结局的混淆因素无关

3 排他性假设

SNP仅通过暴露影响结局

⚠️ 违反假设的后果：
• 违反假设1 → 弱工具变量偏倚（需F统计量检验）
• 违反假设2 → 混杂偏倚
• 违反假设3 → 水平多效性偏倚（需MR-Egger和MR-PRESSO检验）

前置准备

R包安装

# 安装核心MR包
install.packages("TwoSampleMR")
install.packages("MRPRESSO")

# 安装数据处理包
install.packages("data.table")
install.packages("dplyr")
install.packages("openxlsx")

# 加载包
library(TwoSampleMR)
library(MRPRESSO)
library(data.table)
library(dplyr)

GWAS数据来源

数据库	规模	特点
IEU OpenGWAS	50,000+ 暴露	免费，API访问，最常用
GWAS Catalog	100,000+ 研究	文献发表数据，质量参差
UK Biobank	500,000 样本	大队列，表型丰富
FinnGen	500,000 芬兰人	疾病终点，更新快

TwoSampleMR标准工作流

# ============================================================
# TwoSampleMR 标准分析流程
# ============================================================

library(TwoSampleMR)
library(data.table)

# ------------------
# 1. 提取暴露数据
# ------------------

# 从IEU OpenGWAS数据库提取
exp_dat <- read_exposure_data(
    filename = "ieu-a-2",  # BMI的GWAS ID
    snp_col = "SNP",
    beta_col = "beta",
    se_col = "se",
    effect_allele_col = "effect_allele",
    other_allele_col = "other_allele",
    eaf_col = "eaf",
    pval_col = "pval"
)

# 或从本地文件读取
exp_dat <- fread("path/to/exposure_gwas.txt")
exp_obj <- format_data(
    exp_dat,
    type = "exposure",
    snp_col = "SNP",
    beta_col = "beta",
    se_col = "se",
    effect_allele_col = "EA",
    other_allele_col = "OA",
    eaf_col = "EAF",
    pval_col = "P"
)

# ------------------
# 2. 提取结局数据
# ------------------

# 从IEU OpenGWAS提取（基于暴露的SNPs）
out_dat <- extract_outcome_data(
    snps = exp_dat$SNP,
    outcomes = "ieu-a-7"  # 冠心病的GWAS ID
)

# 或从本地文件读取
out_dat <- fread("path/to/outcome_gwas.txt")
out_obj <- format_data(
    out_dat,
    type = "outcome",
    # ... 同上参数 ...
)

# ------------------
# 3. 数据协调
# ------------------

dat <- harmonise_data(
    exposure_dat = exp_obj,
    outcome_dat = out_obj
)

# 查看协调后的数据
head(dat)

# ------------------
# 4. MR分析
# ------------------

res <- mr(dat, method_list = c(
    "mr_ivw",           # 逆方差加权法（主要方法）
    "mr_weighted_median", # 加权中位数法
    "mr_egger_regression" # MR-Egger回归
))

# ------------------
# 5. 敏感性分析
# ------------------

# 异质性检验
het <- mr_heterogeneity(dat)

# 多效性检验（MR-Egger截距）
pleio <- mr_pleiotropy_test(dat)

# Leave-one-out分析
loo <- mr_leaveoneout(dat)

# ------------------
# 6. 结果汇总
# ------------------

# 合并所有结果
res_all <- rbind(res, het, pleio)

# 查看主要结果
res <- subset(res, mr_method == "Inverse variance weighted")

💡 常用MR方法对比：
• IVW：主要方法，假设所有SNP均有效（无多效性）
• Weighted Median：允许≤50%的SNP无效
• MR-Egger：允许所有SNP存在多效性，但统计效能较低

F统计量计算（最重要）

🔴 为什么重要？

F统计量衡量工具变量强度。F > 10是避免弱工具变量偏倚的经验法则。F < 10会导致因果效应估计向观察性关联偏倚。

计算公式

F = (beta²) / (se²)

R代码实现

# ============================================================
# F统计量计算脚本
# ============================================================

library(data.table)
library(dplyr)
library(openxlsx)

# 读取暴露GWAS数据
exp_dat <- fread("path/to/exposure_gwas.txt")

# 计算F统计量（单SNP情况）
# F = (beta^2) / (se^2)
if ("beta.exposure" %in% names(exp_dat)) {
    exp_dat$F_stat <- (exp_dat$beta.exposure^2) / (exp_dat$se.exposure^2)
} else if ("beta" %in% names(exp_dat)) {
    exp_dat$F_stat <- (exp_dat$beta^2) / (exp_dat$se^2)
}

# 汇总统计
f_summary <- data.frame(
    N_SNPs = nrow(exp_dat),
    F_mean = mean(exp_dat$F_stat, na.rm = TRUE),
    F_median = median(exp_dat$F_stat, na.rm = TRUE),
    F_min = min(exp_dat$F_stat, na.rm = TRUE),
    F_max = max(exp_dat$F_stat, na.rm = TRUE),
    N_F_above_10 = sum(exp_dat$F_stat > 10, na.rm = TRUE),
    N_F_above_20 = sum(exp_dat$F_stat > 20, na.rm = TRUE),
    Pct_F_above_10 = round(sum(exp_dat$F_stat > 10, na.rm = TRUE) / nrow(exp_dat) * 100, 1),
    Pct_F_above_20 = round(sum(exp_dat$F_stat > 20, na.rm = TRUE) / nrow(exp_dat) * 100, 1)
)

print(f_summary)

# 保存结果
write.xlsx(f_summary, "F_statistic_report.xlsx")

解读标准

指标	通过标准	说明
Mean F	> 10	平均F统计量应>10
F>10比例	≥ 80%	至少80%的SNP满足F>10
F>20比例	越高越好	用于敏感性分析，工具变量更强

⚠️ F统计量不足怎么办？
• 降低P值阈值（如P < 5×10⁻⁸ → P < 1×10⁻⁵）增加SNP数量
• 考虑使用暴露的汇总数据而非单SNP数据
• 报告时注明弱工具变量限制

MR-PRESSO离群值检测

🎯 目的

检测并校正水平多效性离群值（某些SNP通过非暴露途径影响结局）。MR-PRESSO通过模拟识别离群值并重新计算MR效应值。

R代码实现

# ============================================================
# MR-PRESSO 分析
# ============================================================

library(MRPRESSO)
library(TwoSampleMR)

# 准备数据（需为MRInputObject格式）
presso_res <- mr_presso(
    BetaOutcome = "beta.outcome",
    BetaExposure = "beta.exposure",
    SdOutcome = "se.outcome",
    SdExposure = "se.exposure",
    OUTLIERtest = TRUE,
    DISTORTIONtest = TRUE,
    data = dat,
    NbDistribution = 1000,  # 模拟次数，增加提高准确性
    SignifThreshold = 0.05
)

# 查看离群值检测结果
print(presso_res)

# Global test：检测是否存在水平多效性
# P < 0.05 表示存在显著多效性

# 离群值SNP列表
outliers <- presso_res$`MR-PRESSO results`$`Outliers test`$RSNP

# 校正后的MR结果（移除离群值）
presso_corrected <- presso_res$`MR-PRESSO results`$`Distortion test`

📌 注意事项：
• MR-PRESSO需要至少3个SNP才能运行
• NbDistribution设为10000可提高准确性，但计算时间更长
• 校正后的效应值更可靠，应报告校正前后对比

Steiger方向性检验

🎯 目的

验证因果方向是 暴露→结局，而非 结局→暴露（反向因果）。基于SNP对暴露和结局的解释方差（R²）比较。

R代码实现

# ============================================================
# Steiger 方向性检验
# ============================================================

library(TwoSampleMR)

# 需要协调后的数据
# dat <- harmonise_data(exp_obj, out_obj)

# 执行Steiger检验
steiger_res <- directionality_test(dat)

# 查看结果
print(steiger_res)

# 关键输出指标：
# - snp_r2.exposure: SNP对暴露的解释方差
# - snp_r2.outcome: SNP对结局的解释方差
# - correct_causal_direction: TRUE表示方向正确
# - steiger_pval: 方向性的显著性P值

# 批量分析多个暴露-结局对时
steiger_results <- data.frame()

for (i in 1:nrow(matched_paths)) {
    # ... 数据准备 ...
    dat <- harmonise_data(exp_dat, out_dat)

    if (nrow(dat) > 0) {
        steiger <- directionality_test(dat)

        steiger_row <- data.frame(
            id.exposure = exp_id,
            id.outcome = out_id,
            snp_r2.exposure = steiger$snp_r2.exposure,
            snp_r2.outcome = steiger$snp_r2.outcome,
            correct_causal_direction = steiger$correct_causal_direction,
            steiger_pval = steiger$steiger_pval
        )

        steiger_results <- rbind(steiger_results, steiger_row)
    }
}

解读标准

指标	通过标准	说明
correct_causal_direction	TRUE	因果方向正确（暴露→结局）
steiger_pval	< 0.05	方向性统计显著
snp_r2.exposure > snp_r2.outcome	是	SNP对暴露的解释方差更大

⚠️ 方向错误怎么办？
• 检查暴露和结局是否弄反
• 可能存在反向因果（结局影响暴露）
• 考虑双向MR分析验证

异质性检验

🎯 目的

评估不同SNP的因果效应估计是否一致。异质性可能提示存在多效性或效应修饰。

R代码实现

# ============================================================
# 异质性检验（Cochran's Q）
# ============================================================

library(TwoSampleMR)

# 执行异质性检验
het <- mr_heterogeneity(dat)

# 查看结果
print(het)

# 关键输出指标：
# - Q: Cochran's Q统计量
# - Q_df: 自由度
# - Q_pval: 异质性显著性P值

# 解读：
# Q_pval < 0.05 → 存在显著异质性
# 此时建议使用随机效应IVW模型

💡 异质性处理：
• P < 0.05：使用随机效应模型（IVW默认为固定效应）
• 使用MR-PRESSO移除离群值
• 亚组分析寻找异质性来源

多效性检验（MR-Egger截距）

🎯 目的

检测有向水平多效性（SNP通过多条途径影响结局的平均效应）。MR-Egger截距显著不为零表明存在多效性偏倚。

R代码实现

# ============================================================
# MR-Egger 截距检验
# ============================================================

library(TwoSampleMR)

# 执行多效性检验
pleio <- mr_pleiotropy_test(dat)

# 查看结果
print(pleio)

# 关键输出指标：
# - egger_intercept: Egger回归截距
# - se: 截距标准误
# - pval: 截距显著性P值

# 解读：
# pval > 0.05 → 截距与0无差异 → 无显著多效性
# pval < 0.05 → 截距显著 ≠ 0 → 存在多效性偏倚

📌 MR-Egger vs IVW：
• MR-Egger允许所有SNP存在多效性（通过截距校正）
• 但统计效能较低，需要更多SNP
• 当截距P>0.05时，IVW结果是主要报告

Leave-one-out敏感性分析

🎯 目的

检验单个SNP是否驱动整体MR结果。逐一移除每个SNP重新计算MR效应，观察结果是否稳定。

R代码实现

# ============================================================
# Leave-one-out 分析
# ============================================================

library(TwoSampleMR)

# 执行LOO分析
loo <- mr_leaveoneout(dat)

# 查看结果
print(loo)

# 可视化（推荐）
mr_leaveoneout_plot(loo)

# 解读：
# 如果移除某个SNP后效应值大幅变化/翻转，
# 说明该SNP驱动结果，需谨慎报告

💡 稳定性判断：
• 所有LOO点估计值置信区间重叠 → 结果稳定
• 移除某SNP后方向改变 → 需报告该SNP影响
• 结合单SNP森林图综合判断

MR结果可视化

常用图表

# ============================================================
# MR结果可视化
# ============================================================

library(TwoSampleMR)

# 1. 森林图（不同MR方法对比）
mr_forest_plot(res)

# 2. 散点图（SNP水平效应）
mr_scatter_plot(res, dat)

# 3. 漏斗图（对称性检验）
mr_funnel_plot(res)

# 4. Leave-one-out图
mr_leaveoneout_plot(loo)

# 5. 单SNP森林图
mr_singlesnp_plot(dat)

📊 森林图

展示不同MR方法的效应值和置信区间，用于方法间对比

📍 散点图

展示SNP对暴露和结局的效应关系，直观显示因果斜率

🔺 漏斗图

SNP精确度vs效应值，对称性提示无多效性

📉 LOO图

逐一移除SNP后的效应值变化，检验稳定性

📝 论文Methods报告模板

Sensitivity Analyses and Quality Control

Instrument strength was assessed using F-statistics, with F > 10 indicating adequate strength to avoid weak instrument bias. The mean F-statistic was XX.X (range: X.X-X.X), with XX.X% of instruments having F > 10.

Horizontal pleiotropy was evaluated using MR-PRESSO global test and MR-Egger intercept. The MR-Egger intercept was not significantly different from zero for XX.X% of associations (P > 0.05), suggesting no directional pleiotropy.

Heterogeneity was assessed using Cochran's Q statistic. Significant heterogeneity (P < 0.05) was observed in XX (X.X%) associations, for which random-effects IVW estimates were reported.

Steiger filtering confirmed the correct causal direction (exposure → outcome) for XX.X% of associations.

Leave-one-out analysis confirmed that no single SNP was driving the observed associations.

✅ MR质控检查清单

步骤	质控项目	目的	必要性
1	F统计量筛选	排除弱工具变量偏倚	⭐⭐⭐⭐⭐ 必须
2	MR-PRESSO	检测水平多效性离群值	⭐⭐⭐⭐⭐ 必须
3	Steiger检验	验证因果方向	⭐⭐⭐⭐ 强烈推荐
4	异质性检验	Cochran's Q检验	⭐⭐⭐⭐ 推荐
5	MR-Egger截距	检测有向多效性	⭐⭐⭐⭐ 推荐
6	Leave-one-out	检测单SNP驱动结果	⭐⭐⭐ 推荐

常见错误

❌ 错误1：忽略F统计量检验

问题：F < 10会导致因果效应向观察性关联偏倚

✅ 正确做法：报告每个分析的F统计量均值和范围

❌ 错误2：只报告IVW结果

问题：单一方法结果可能不可靠

✅ 正确做法：同时报告IVW、Weighted Median、MR-Egger结果

❌ 错误3：样本重叠未说明

问题：暴露和结局样本重叠会导致偏倚

✅ 正确做法：在Methods中说明样本重叠情况

❌ 错误4：未进行Steiger方向检验

问题：可能存在反向因果

✅ 正确做法：进行双向MR或Steiger检验验证方向