多变量MR分析 - 医研AI工具库

📊 新手导读：MVMR解决什么问题？

当多个暴露彼此相关（如血脂组分）时，单变量MR难区分独立效应，MVMR用于估计“条件因果效应”。

背景介绍

MVMR可在同一模型中同时纳入多个暴露，降低暴露间相关导致的混杂，适合做机制拆分与药物靶点优先级评估。

备注（统一三段）

适用人群：暴露变量高度相关（如血脂组分）并需估计独立效应的人群。

常见错误：忽略条件F统计量与共线性诊断，会让MVMR系数不稳定。

论文写法：结果报告建议优先给条件效应，并补充与单变量MR方向一致性比较。

📖 什么是多变量MR？

多变量孟德尔随机化（Multivariate Mendelian Randomization, MVMR）是MR的扩展方法，允许同时估计多个暴露对同一结局的因果效应。

🎯

识别独立效应

区分多个暴露的独立作用

🛡️

控制混杂

调整暴露之间的相关性

💊

药物靶向

评估药物靶点效应

🔄 单变量 vs 多变量MR

特征	单变量MR	多变量MR
暴露数量	1个	≥2个
混杂控制	手动调整	模型内同时调整
工具变量	独立筛选	联合筛选
适用场景	初步探索	精细分析

💻 方法一：TwoSampleMR包

使用 mv_extract_exposures 在线提取

library(TwoSampleMR)

# 定义多个暴露的GWAS ID
# 示例：HDL、LDL、TG对冠心病的效应
id_exposure <- c(
  "ieu-a-299", # HDL cholesterol
  "ieu-a-300", # LDL cholesterol 
  "ieu-a-302"  # Triglycerides
)

# 定义结局ID
id_outcome <- "ieu-a-7" # Coronary heart disease

# 提取多个暴露的SNP
exposure_dat <- mv_extract_exposures(id_exposure)

# 提取结局数据
outcome_dat <- extract_outcome_data(
  snps = exposure_dat$SNP,
  outcomes = id_outcome
)

# 数据协调
mvdat <- mv_harmonise_data(exposure_dat, outcome_dat)

# 多变量MR分析
res <- mv_multiple(mvdat)

# 生成OR值
res_OR <- generate_odds_ratios(res$result)
print(res_OR)

SNP选择规则

# 分别提取各暴露的SNP，然后合并去重
exp1 <- extract_instruments(outcomes = "ieu-a-299") # HDL
exp2 <- extract_instruments(outcomes = "ieu-a-300") # LDL
exp3 <- extract_instruments(outcomes = "ieu-a-302") # TG

# 合并所有SNP
exp_all <- rbind(exp1, exp2, exp3)
exp_all <- exp_all["SNP"] # 只保留SNP列
exp_all <- unique(exp_all) # 去重

# 查看最终SNP数量
cat("最终SNP数量:", nrow(exp_all), "\n")

LASSO特征选择

当暴露高度相关时，使用LASSO去除冗余SNP

# LASSO特征选择
lasso_features <- mv_lasso_feature_selection(mvdat)
cat("保留的特征:", lasso_features, "\n")

# 基于LASSO结果进行MR
mv_lasso <- mv_subset(
  mvdat,
  features = lasso_features,
  intercept = FALSE,
  instrument_specific = FALSE,
  pval_threshold = 5e-08,
  plots = FALSE
)

# 结果
mv_lasso_OR <- generate_odds_ratios(mv_lasso$result)
print(mv_lasso_OR)

💻 方法二：MendelianRandomization包

数据准备

library(MendelianRandomization)
library(data.table)

# 假设已从GWAS提取以下数据（示例数据）
# 暴露1: LDL
ldlc <- c(-0.15, -0.22, 0.08, -0.31, 0.12, -0.19)
ldlcse <- c(0.03, 0.04, 0.05, 0.03, 0.04, 0.02)

# 暴露2: HDL
hdlc <- c(0.05, -0.08, 0.11, 0.03, -0.06, 0.09)
hdlcse <- c(0.02, 0.03, 0.04, 0.02, 0.03, 0.02)

# 暴露3: TG
trig <- c(0.18, -0.12, 0.25, 0.07, -0.15, 0.21)
trigse <- c(0.04, 0.05, 0.06, 0.03, 0.04, 0.05)

# 结局: CHD
chdlodds <- c(0.22, 0.15, -0.08, 0.31, -0.12, 0.18)
chdloddsse <- c(0.05, 0.06, 0.07, 0.04, 0.05, 0.06)

构建MVMR输入对象

# 构建多变量MR输入对象
MRMVInputObject <- mr_mvinput(
  bx = cbind(ldlc, hdlc, trig),      # 暴露效应值
  bxse = cbind(ldlcse, hdlcse, trigse),  # 暴露标准误
  by = chdlodds,              # 结局效应值
  byse = chdloddsse            # 结局标准误
)

MRMVInputObject

IVW方法

# 逆方差加权 (IVW)
MRMVObject_ivw <- mr_mvivw(
  MRMVInputObject, 
  model = "default",
  correl = FALSE,
  distribution = "normal",
  alpha = 0.05
)

print(MRMVObject_ivw)

MR-Egger方法

# MR-Egger
MRMVObject_egger <- mr_mvegger(
  MRMVInputObject,
  orientate = 1, # 第一个暴露作为方向参考
  correl = FALSE,
  distribution = "normal",
  alpha = 0.05
)

print(MRMVObject_egger)

加权中位数

# 加权中位数 (Weighted Median)
MRMVObject_median <- mr_mvmedian(
  MRMVInputObject,
  distribution = "normal",
  alpha = 0.05,
  iterations = 10000,
  seed = 314159265
)

print(MRMVObject_median)

LASSO方法

# LASSO (自动去除水平多效性)
MRMVObject_lasso <- mr_mvlasso(
  MRMVInputObject,
  orientate = 1,
  distribution = "normal",
  alpha = 0.05,
  lambda = numeric(0) # 自动选择lambda
)

print(MRMVObject_lasso)

💻 方法三：MVMR包（工具变量强度评估）

MVMR包专注于评估多变量MR中工具变量的强度，特别适合检测弱工具变量问题。

# 安装MVMR包
# remotes::install_github("WSpiller/RMVMR", build_opts = c("--no-resave-data", "--no-manual"))

library(MVMR)

# 准备数据格式
# 需要: bx (暴露效应), by (结局效应), seBX (暴露SE), seBY (结局SE), RSID (SNP ID)
F.data <- format_mvmr(
  BXGs = cbind(exp_a$beta, exp_b$beta), # 多个暴露的beta
  BYG = outcome$Beta,          # 结局的beta
  seBXGs = cbind(exp_a$SE, exp_b$SE),  # 多个暴露的SE
  seBYG = outcome$SE,           # 结局的SE
  RSID = exp_a$SNP            # SNP ID
)

head(F.data)

弱工具变量检验

# 检验工具变量强度
# Q统计量: 评估异质性
# F统计量: 评估弱工具变量 (F > 10 为佳)

sres <- strength_mvmr(r_input = F.data, gencov = 0)
print(sres)

# 结果解读:
# F > 10: 工具变量足够强
# F < 10: 存在弱工具变量偏倚

水平多效性检验

# 检验水平多效性 (Sargan检验)
pres <- pleiotropy_mvmr(r_input = F.data, gencov = 0)
print(pres)

# 结果解读:
# p > 0.05: 不存在显著的水平多效性
# p < 0.05: 存在水平多效性，需谨慎解读

IVW效应估计

# IVW方法估计因果效应
res <- ivw_mvmr(r_input = F.data)
print(res)

# 输出: 各暴露的beta、SE、P值

💻 方法四：MR-PRESSO多变量

MR-PRESSO不仅可以做单变量分析，也支持多变量离群值检测。

library(MRPRESSO)

# 准备数据 (从mv_harmonise_data获取)
SummaryStats <- cbind(
  mvdat[["outcome_beta"]],    # 结局beta
  mvdat[["exposure_beta"]][,1], # 暴露1 beta
  mvdat[["exposure_beta"]][,2], # 暴露2 beta
  mvdat[["exposure_beta"]][,3], # 暴露3 beta
  mvdat[["outcome_se"]],     # 结局SE
  mvdat[["exposure_se"]][,1],  # 暴露1 SE
  mvdat[["exposure_se"]][,2],  # 暴露2 SE
  mvdat[["exposure_se"]][,3]  # 暴露3 SE
)
SummaryStats <- data.frame(SummaryStats)

# MR-PRESSO多变量分析
mr_presso(
  BetaOutcome = "X1",
  BetaExposure = c("X2", "X3", "X4"),
  SdOutcome = "X5",
  SdExposure = c("X6", "X7", "X8"),
  OUTLIERtest = TRUE,
  DISTORTIONtest = TRUE,
  data = SummaryStats,
  NbDistribution = 1000, # Bootstrap次数
  SignifThreshold = 0.05
)

📈 结果解读

方法	暴露1 (LDL)	暴露2 (HDL)	暴露3 (TG)
IVW	0.42 (0.31-0.53)	-0.18 (-0.27--0.09)		MR-Egger	0.38 (0.22-0.54)	-0.15 (-0.29--0.01)	Weighted Median	0.40 (0.27-0.53)	-0.16 (-0.28--0.04)	0.33 (0.20-0.46)

结论：LDL升高CHD风险(OR=1.52)，HDL降低CHD风险(OR=0.84)，TG升高CHD风险(OR=1.42)。

⚠️ 常见问题与解决方案

Q1: 暴露之间高度相关怎么办？

使用LASSO进行特征选择，或选择生物学上独立的暴露

Q2: F统计量<10怎么处理？

使用更宽松的P值阈值(如1e-5)获取更多SNP，或使用MVMR包评估

Q3: 样本重叠问题

优先使用独立的GWAS数据，或使用Egger类方法调整

Q4: 水平多效性检验显著

使用MR-PRESSO去除离群SNP，或使用LASSO/Egger方法

📚 推荐阅读

原始文献

• Sanderson E et al. Multivariable MR. IJE 2022
• Burgess S et al. MVMR. AJHG 2020
• Yang J et al. MVMR with GWAS. Nat Genet 2020

配套代码

• MVMR.R (玮瑜科研课程)
• 多变量MR.R (实操教程)