孟德爾隨機化小試從軟件安裝數據下載到多種檢驗

孟德爾隨機化（Mendelian Randomization, MR）作為一種利用基因數據推斷因果關系的強大工具，在流行病學研究中應用廣泛。本文將詳細講解MR的核心原理、完整分析流程，并附上關鍵代碼實現，幫助你從零開始完成一次MR分析。
在這里插入圖片描述

1. 什么是孟德爾隨機化？

孟德爾隨機化是一種基于全基因組關聯研究（GWAS）數據，利用單核苷酸多態性（SNP）作為工具變量（IV） 來揭示暴露與結局間因果關系的方法。

其核心邏輯類似“自然隨機對照試驗”：

隨機對照試驗（RCT）：人為將研究對象隨機分配到實驗組/對照組
孟德爾隨機化：通過基因變異的自然隨機分配（減數分裂時的隨機分離），將攜帶特定等位基因的個體視為“暴露組”，非攜帶者視為“對照組”

2. 為什么選擇MR？

相比隊列研究等觀察性研究，MR的優勢在于：

避免反向因果：基因在出生前已確定，暴露狀態（由基因影響）早于結局發生
減少混雜偏倚：基因與后天環境、行為等混雜因素通常無關
低成本驗證：可直接免費利用公開GWAS數據

MR的三大核心假設（必須滿足！）

MR結果的可靠性完全依賴以下三個假設，需嚴格驗證：

假設名稱	核心內容	驗證方法
關聯性假設	工具變量（SNP）與暴露因素顯著相關	計算p值（通常p<5e-8）、F統計量（建議F>10，避免弱工具變量偏倚）、R2
獨立性假設	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,  # 剔除r2>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){# 計算R2和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