孟德尔随机化小试从软件安装数据下载到多种检验

孟德尔随机化（Mendelian Randomization, MR）作为一种利用基因数据推断因果关系的强大工具，在流行病学研究中应用广泛。本文将详细讲解MR的核心原理、完整分析流程，并附上关键代码实现，帮助你从零开始完成一次MR分析。
在这里插入图片描述

1. 什么是孟德尔随机化？

孟德尔随机化是一种基于全基因组关联研究（GWAS）数据，利用单核苷酸多态性（SNP）作为工具变量（IV） 来揭示暴露与结局间因果关系的方法。

其核心逻辑类似“自然随机对照试验”：

随机对照试验（RCT）：人为将研究对象随机分配到实验组/对照组
孟德尔随机化：通过基因变异的自然随机分配（减数分裂时的随机分离），将携带特定等位基因的个体视为“暴露组”，非携带者视为“对照组”

2. 为什么选择MR？

相比队列研究等观察性研究，MR的优势在于：

避免反向因果：基因在出生前已确定，暴露状态（由基因影响）早于结局发生
减少混杂偏倚：基因与后天环境、行为等混杂因素通常无关
低成本验证：可直接免费利用公开GWAS数据

MR的三大核心假设（必须满足！）

MR结果的可靠性完全依赖以下三个假设，需严格验证：

假设名称	核心内容	验证方法
关联性假设	工具变量（SNP）与暴露因素显著相关	计算p值（通常p<5e-8）、F统计量（建议F>10，避免弱工具变量偏倚）、R²
独立性假设	SNP与结局之间不存在通过混杂因素的关联（SNP仅通过暴露影响结局）	确保SNP与已知混杂因素无关联，通过敏感性分析验证
排他性假设	SNP不直接影响结局，仅通过暴露间接影响	计算SNP与结局的独立关联性（p>0.05），MR-Egger回归可弱化此假设要求

关键提示：若假设不成立，因果推断结果可能完全错误。例如，若某SNP既影响BMI（暴露），又直接影响心脏病（结局），则违反排他性假设，不能用于推断BMI与心脏病的因果关系。

MR分析完整流程与代码实现

1. 环境配置（R语言）

MR分析涉及到的R包可以一次性安装

MRCIEU/TwoSampleMR
VariantAnnotation
mrcieu/gwasglue
VariantAnnotation
mrcieu/gwasglue 
CMplot
MendelianRandomization
LDlinkR
ggplot2
ggfunnel
cowplot
plinkbinr
FastDownloader
FastTraitR
MRPRESSO
SNPlocs.Hsapiens.dbSNP155.GRCh37
SNPlocs.Hsapiens.dbSNP155.GRCh38
tidyverse
MungeSumstats
GenomicFiles
explodecomputer/plinkbinr
MRPRESSO
pedropark99/ggfunnel
xiechengyong123/friendly2MR
NightingaleHealth/ggforestplot
BSgenome.Hsapiens.1000genomes.hs37d5
BSgenome.Hsapiens.NCBI.GRCh38
glmnet
meta
psych

2. 数据来源

MR分析依赖GWAS汇总数据，常用公开数据库包括：

IEU汇总数据库（MRCIEU）：全球最大的GWAS数据平台（https://gwas.mrcieu.ac.uk/）
精神病学基因组联盟（PGC）：专注精神疾病GWAS数据
社会科学遗传学联盟（SSGAC）：包含教育、收入等社会表型数据
UK Biobank：包含超50万样本的多表型数据
自建数据：自己开展的GWAS分析结果

3. 暴露数据处理（关键步骤）

暴露数据即与“暴露因素”相关的GWAS数据，需筛选符合“关联性假设”的SNP。

3.1 从VCF文件获取暴露数据

# 读取VCF格式的GWAS数据
VCF_dat <- VariantAnnotation::readVcf('ieu-a-2.vcf.gz')
# 转换为TwoSampleMR所需格式
exp_dat <- gwasglue::gwasvcf_to_TwoSampleMR(vcf = VCF_dat)# 筛选符合关联性假设的SNP（p<5e-8）
exp_dat <- subset(exp_dat, pval.exposure < 5e-08)

3.2 处理连锁不平衡（LD）

SNP间的连锁不平衡会导致工具变量相关性过高，需通过“聚类”去除：

# 自定义本地LD聚类函数（避免调用云端API）
fix_ld_clump_local <- function (dat, tempfile, clump_kb, clump_r2, clump_p, bfile, plink_bin) {shell <- ifelse(Sys.info()["sysname"] == "Windows", "cmd", "sh")# 输出SNP和p值到临时文件write.table(data.frame(SNP = dat[["rsid"]], P = dat[["pval"]]),file = tempfile, row.names = F, col.names = T, quote = F)# 调用PLINK进行LD聚类fun2 <- paste0(shQuote(plink_bin, type = shell), " --bfile ",shQuote(bfile, type = shell), " --clump ", shQuote(tempfile,type = shell), " --clump-p1 ", clump_p, " --clump-r2 ",clump_r2, " --clump-kb ", clump_kb, " --out ", shQuote(tempfile,type = shell))print(fun2)system(fun2)# 读取聚类结果res <- read.table(paste(tempfile, ".clumped", sep = ""), header = T)unlink(paste(tempfile, "*", sep = ""))# 输出去LD后的SNPreturn(subset(dat, dat[["rsid"]] %in% res[["SNP"]]))}# 运行LD聚类（以欧洲人群为例）
fuck <- fix_ld_clump_local(dat = dplyr::tibble(rsid=exp_dat$SNP, pval=exp_dat$pval.exposure),tempfile = file.path(getwd(),'tmp.ld_clump.exp_dat'),clump_kb = 10000,  # 10kb内的SNP进行聚类clump_r2 = 0.001,  # 剔除r²>0.001的SNPclump_p = 1,plink_bin = 'path/to/plink',  # PLINK路径bfile = "1kg/EUR"  # 1000G欧洲人群参考面板
)# 提取去LD后的暴露数据
exp_dat_clumped <- exp_dat[exp_dat$SNP %in% fuck$rsid,]
saveRDS(file = 'ieu-a-2.exp_gwas', exp_dat_clumped)

3.3 从其他来源获取暴露数据

# 从GWAS目录获取（例如BMI数据）
df_gwas <- subset(MRInstruments::gwas_catalog,grepl("Speliotes", Author) & Phenotype == "Body mass index")exp_dat <- TwoSampleMR::format_data(df_gwas)# 从GTEx获取eQTL数据（例如特定基因在脂肪组织的表达）
df_gwas <- subset(MRInstruments::gtex_eqtl, gene_name == "IRAK1BP1" & tissue == "Adipose Subcutaneous")exp_dat <- TwoSampleMR::format_gtex_eqtl(df_gwas)# 从IEU数据库直接提取（通过ID）
exp_gwas <- TwoSampleMR::extract_instruments(outcomes = 'ieu-a-2')  # 'ieu-a-2'为某表型ID

3.4 评估工具变量强度（计算F统计量）

F统计量用于判断是否为“弱工具变量”（F<10提示可能存在偏倚）：

MR_calc_r2_F <- function(beta, eaf, N, se){# 计算R²和F统计量r2 <- (2 * (beta^2) * eaf * (1 - eaf)) /(2 * (beta^2) * eaf * (1 - eaf) + 2 * N * eaf * (1 - eaf) * se^2)F <- r2 * (N - 2) / (1 - r2)print(paste("平均F值：", mean(F)))return(dplyr::tibble(r2=r2, F=F))
}# 计算并筛选F>10的SNPf_stats <- MR_calc_r2_F(beta = exp_dat_clumped$beta.exposure,  # 暴露的beta值（log(OR)）eaf = exp_dat_clumped$eaf.exposure,    # 等位基因频率N = exp_dat_clumped$samplesize.exposure,  # 样本量se = exp_dat_clumped$se.exposure       # 标准误)exp_dat_final <- exp_dat_clumped[f_stats$F > 10, ]  # 保留强工具变量

4. 结局数据处理

结局数据即与“结局变量”相关的GWAS数据，需满足“排他性假设”。

4.1 从PGC数据库获取结局数据（以ADHD为例）


# 读取ADHD的GWAS汇总数据（meta分析结果）df_gwas =read.table(gzfile('ADHD2022_iPSYCH_deCODE_PGC.meta.gz'), header = T)# 仅保留与暴露SNP重叠的数据
df_gwas <- df_gwas[df_gwas$SNP %in% exp_gwas$SNP,]# 转换为TwoSampleMR格式
out_gwas <- data.frame(SNP = df_gwas$SNP,chr = as.character(df_gwas$CHR),pos = df_gwas$BP,beta.outcome = log(df_gwas$OR),  # 将OR转换为log(OR)se.outcome = df_gwas$SE,samplesize.outcome = df_gwas$Nca + df_gwas$Nco,  # 总样本量（病例+对照）pval.outcome = df_gwas$P,eaf.outcome = with(df_gwas, (FRQ_A_38691*Nca + FRQ_U_186843*Nco)/(Nca + Nco)),  # 计算整体等位基因频率effect_allele.outcome = df_gwas$A1,other_allele.outcome = df_gwas$A2,outcome = 'ADHD',id.outcome = 'ADHD2022_iPSYCH_deCODE_PGC'
)# 筛选符合排他性假设的SNP（p>0.05，与结局无直接关联）
out_gwas <- subset(out_gwas, pval.outcome > 0.05)

4.2 从其他来源获取结局数据

opengwas_jwt需要注册ieu并申请token后才可以从网站下载数据

# 从IEU数据库提取
out_gwas <- TwoSampleMR::extract_outcome_data(snps = exp_gwas$SNP, outcomes = 'ieu-a-7'，opengwas_jwt = opengwas_jwt)  # 'ieu-a-7'为结局ID# 从UK Biobank提取
anxiety_outcome <- TwoSampleMR::extract_outcome_data(snps = exp_dat_clumped$SNP, outcomes = "ukb-b-11311")

5. 数据协调（Harmonization）

确保暴露和结局数据中SNP的等位基因方向一致（关键步骤，否则结果完全错误）：

# 协调暴露和结局数据
dat <- TwoSampleMR::harmonise_data(exposure_dat = exp_gwas,  # 处理后的暴露数据outcome_dat = out_gwas    # 处理后的结局数据
)
# 查看协调结果（检查是否有方向冲突的SNP被剔除）
head(dat)

关键说明：协调过程会自动处理：

等位基因方向不一致的SNP（如暴露中A为效应等位基因，结局中a为效应等位基因，会自动转换）
回文SNP（如A/T和T/A，无法判断方向时会被剔除）
完全不匹配的SNP（直接剔除）

6. 因果效应估计

使用多种方法计算暴露对结局的因果效应，结果一致更可信：

# 查看所有可用的MR方法
TwoSampleMR::mr_method_list()# 选择常用方法进行分析（IVW、Egger、加权中位数）
mr_regression <- TwoSampleMR::mr(dat,method_list = c('mr_ivw', 'mr_egger_regression', 'mr_weighted_median')
)
print(mr_regression)# 若结局为分类变量（如疾病），转换为OR值（方便解读）
mr_regression_or <- TwoSampleMR::generate_odds_ratios(mr_res = mr_regression)
print(mr_regression_or)# 绘制散点图（直观展示各SNP的效应及整体趋势）
pdf(file = 'MR散点图.pdf', width = 6, height = 6)
print(TwoSampleMR::mr_scatter_plot(mr_results = mr_regression, dat = dat))
dev.off()

7. 假设验证与敏感性分析

7.1 异质性检验（判断工具变量效应是否一致）

# IVW方法的Q统计量检验
heterogeneity_ivw <- TwoSampleMR::mr_heterogeneity(dat)
print(heterogeneity_ivw)  # Q_pval < 0.05提示存在异质性# MR-PRESSO检验（检测离群值导致的异质性）
heterogeneity_presso <- TwoSampleMR::run_mr_presso(dat, NbDistribution = 3000)  # 3000次模拟# 全局检验（是否存在异质性）
print(heterogeneity_presso[[1]]$`MR-PRESSO results`$`Global Test`$Pvalue)# 离群值检测（若存在，需剔除后重新分析）
print(heterogeneity_presso[[1]]$`MR-PRESSO results`$`Distortion Test`$`Outliers Indices`)

若存在异质性：使用随机效应IVW模型，或剔除离群值后重新分析
若无异质性：固定效应IVW模型更可靠

7.2 水平多效性检验（验证独立性假设）

# MR-Egger截距检验
pleiotropy_test <- TwoSampleMR::mr_pleiotropy_test(dat)
print(pleiotropy_test)
# 若p < 0.05：拒绝截距为0的假设，提示存在水平多效性（结果不可靠）

7.3 Leave-one-out分析（敏感性检验）

# 逐一剔除每个SNP后重新分析，判断结果是否依赖某个SNP
mr_loo <- TwoSampleMR::mr_leaveoneout(dat)# 绘制leave-one-out图
TwoSampleMR::mr_leaveoneout_plot(leaveoneout_results = mr_loo)
# 若剔除某SNP后结果显著变化（如效应值跨过0），提示结果不稳定

7.4 单SNP分析（查看个体工具变量效应）

# 每个SNP单独进行MR分析
mr_res_single <- TwoSampleMR::mr_singlesnp(dat)# 绘制漏斗图（判断是否存在小研究效应）
TwoSampleMR::mr_funnel_plot(mr_res_single)# 绘制森林图（展示每个SNP的效应）
TwoSampleMR::mr_forest_plot(mr_res_single)

7.5 方向性检验（判断因果方向）

TwoSampleMR::directionality_test(dat) # TRUE表示暴露→结局的方向更可能成立

8. 一键生成报告（可选）

# 生成Markdown格式的完整报告
TwoSampleMR::mr_report(dat, output_type = "md")

四、MR的局限性与注意事项

依赖GWAS数据质量：样本量过小、人群分层会导致结果偏倚
工具变量选择至关重要：弱工具变量（F<10）会导致效应估计偏差
无法完全避免多效性：即使通过检验，仍可能存在未识别的水平多效性
仅反映平均效应：无法推断个体水平的因果关系
数据重叠问题：暴露和结局数据若来自同一人群，需控制样本重叠比例（建议<10%）

参考

https://mrcieu.github.io/TwoSampleMR/articles/perform_mr.html
https://hexo.limour.top/Mendelian-Randomization