616_孟德尔随机化分析¶

# === 安装TwoSampleMR包 ===
# 方法1：从r-universe安装（推荐）
install.packages(
    "TwoSampleMR",
    repos = c("https://mrcieu.r-universe.dev", 
              "https://cloud.r-project.org")
)

# 方法2：从GitHub安装
# install.packages("remotes")
remotes::install_github("MRCIEU/TwoSampleMR")

# 安装其他MR相关包
install.packages("MendelianRandomization")  # 另一个MR包
remotes::install_github("rondolab/MR-PRESSO")  # MR-PRESSO

# 加载包
library(TwoSampleMR)  # 加载TwoSampleMR

# === 2024年重要更新：OpenGWAS需要API认证 ===
# 访问 https://api.opengwas.io 获取API key
# 在R中设置：
Sys.setenv(OPENGWAS_JWT = "your_api_key_here")  # 设置API密钥

# === 方法1：从OpenGWAS数据库获取 ===
# 搜索可用的GWAS数据集
ao <- available_outcomes()  # 查看所有可用GWAS
bmi_ids <- ao[grep("body mass index", ao$trait, ignore.case = TRUE), ]
head(bmi_ids[, c("id", "trait", "nsnp", "sample_size")])  # 查看BMI相关

# 提取暴露的工具变量SNP
exposure_dat <- extract_instruments(
    outcomes = "ieu-b-40",  # BMI的GWAS数据集ID
    p1 = 5e-8,              # 全基因组显著性阈值
    clump = TRUE,           # LD clumping去除冗余SNP
    r2 = 0.001,             # clumping的LD阈值
    kb = 10000              # clumping的距离窗口（10Mb）
)
cat("提取到", nrow(exposure_dat), "个工具变量SNP\n")

# === 方法2：从本地GWAS结果读取 ===
exposure_dat <- read_exposure_data(
    filename = "bmi_gwas_summary.txt",  # 本地文件
    sep = "\t",                          # 分隔符
    snp_col = "SNP",                     # SNP列
    beta_col = "BETA",                   # 效应值列
    se_col = "SE",                       # 标准误列
    effect_allele_col = "A1",            # 效应等位基因
    other_allele_col = "A2",             # 其他等位基因
    pval_col = "P",                      # p值列
    eaf_col = "EAF"                      # 效应等位基因频率
)

# === 从OpenGWAS提取结局数据 ===
outcome_dat <- extract_outcome_data(
    snps = exposure_dat$SNP,  # 使用暴露的SNP
    outcomes = "ieu-b-7"      # 2型糖尿病的GWAS数据集ID
)

# === 从本地文件读取 ===
outcome_dat <- read_outcome_data(
    filename = "t2d_gwas_summary.txt",
    snps = exposure_dat$SNP,  # 只提取暴露的SNP
    sep = "\t",
    snp_col = "SNP",
    beta_col = "BETA",
    se_col = "SE",
    effect_allele_col = "A1",
    other_allele_col = "A2",
    pval_col = "P",
    eaf_col = "EAF"
)

# 确保暴露和结局的等位基因编码一致
dat <- harmonise_data(
    exposure_dat = exposure_dat,  # 暴露数据
    outcome_dat = outcome_dat,    # 结局数据
    action = 2                     # 自动协调（推荐）
)
# action=1: 假设所有等位基因都在正链
# action=2: 尝试推断正确的链方向（推荐）
# action=3: 严格模式，无法确定的SNP全部去除

cat("协调后保留", nrow(dat), "个SNP\n")

# === 主要MR分析 ===
res <- mr(dat)  # 运行所有默认MR方法
print(res)  # 查看结果

# 输出包含5种方法的结果：
# MR Egger        — 能检测多效性
# Weighted median  — 允许50%的SNP无效
# IVW             — 主要分析方法（逆方差加权）
# Simple mode     — 简单模式
# Weighted mode   — 加权模式

# === 指定特定方法 ===
res_specific <- mr(
    dat,
    method_list = c(
        "mr_ivw",                    # IVW（主要方法）
        "mr_egger_regression",       # MR-Egger
        "mr_weighted_median",        # 加权中位数
        "mr_penalised_weighted_median"  # 惩罚加权中位数
    )
)

# === 1. 异质性检验 ===
het <- mr_heterogeneity(dat)  # Cochran's Q检验
print(het)
# p < 0.05 说明SNP间存在异质性，需要关注

# === 2. 水平多效性检验（MR-Egger截距）===
pleio <- mr_pleiotropy_test(dat)  # 多效性检验
print(pleio)
# p < 0.05 说明存在水平多效性，IVW结果可能有偏

# === 3. 留一法分析（Leave-one-out）===
loo <- mr_leaveoneout(dat)  # 每次去掉一个SNP重新分析
# 如果去掉某个SNP后结果变化很大，说明该SNP可能有问题

# === 4. MR-PRESSO（检测和去除离群SNP）===
library(MRPRESSO)
presso_res <- mr_presso(
    BetaOutcome = "beta.outcome",
    BetaExposure = "beta.exposure",
    SdOutcome = "se.outcome",
    SdExposure = "se.exposure",
    OUTLIERtest = TRUE,     # 检测离群SNP
    DISTORTIONtest = TRUE,  # 检测校正前后差异
    data = dat,
    NbDistribution = 1000,  # 分布模拟次数
    SignifThreshold = 0.05  # 显著性阈值
)

# === 1. 散点图（Scatter plot）===
p1 <- mr_scatter_plot(res, dat)  # 暴露-结局效应的散点图
ggsave("mr_scatter.pdf", p1[[1]], width = 8, height = 6)

# === 2. 森林图（Forest plot）===
res_single <- mr_singlesnp(dat)  # 单SNP分析
p2 <- mr_forest_plot(res_single)  # 森林图
ggsave("mr_forest.pdf", p2[[1]], width = 8, height = 10)

# === 3. 留一法图（Leave-one-out plot）===
p3 <- mr_leaveoneout_plot(loo)  # 留一法图
ggsave("mr_leaveoneout.pdf", p3[[1]], width = 8, height = 10)

# === 4. 漏斗图（Funnel plot）===
p4 <- mr_funnel_plot(res_single)  # 漏斗图
ggsave("mr_funnel.pdf", p4[[1]], width = 8, height = 6)

# === 计算F-statistic（工具变量强度）===
# F > 10 说明是强工具变量，可以信赖MR结果

dat$F_stat <- (dat$beta.exposure / dat$se.exposure)^2  # 计算F值
cat("F-statistic范围:", range(dat$F_stat), "\n")
cat("平均F-statistic:", mean(dat$F_stat), "\n")
cat("F < 10的弱工具变量数:", sum(dat$F_stat < 10), "\n")

# 如果有弱工具变量，考虑去除
dat_strong <- dat[dat$F_stat >= 10, ]  # 只保留强IV
cat("去除弱IV后剩余SNP数:", nrow(dat_strong), "\n")

# === 计算总R²（暴露的方差解释比例）===
r2 <- sum(2 * dat$eaf.exposure * (1 - dat$eaf.exposure) * 
          dat$beta.exposure^2)
cat("总R²:", r2, "\n")
# R²太小（<1%）说明工具变量解释暴露的能力弱

# === 整理最终结果 ===
cat("\n========== MR分析结果报告 ==========\n")
cat("暴露: BMI\n")
cat("结局: 2型糖尿病\n")
cat("工具变量SNP数:", nrow(dat), "\n")
cat("平均F-statistic:", round(mean(dat$F_stat), 1), "\n\n")

# 主要结果
cat("--- 主要分析 ---\n")
for(i in 1:nrow(res)) {
    or <- exp(res$b[i])  # 转为OR
    ci_low <- exp(res$b[i] - 1.96 * res$se[i])  # 95%CI下限
    ci_high <- exp(res$b[i] + 1.96 * res$se[i])  # 95%CI上限
    cat(sprintf("%s: OR=%.2f (95%%CI: %.2f-%.2f), p=%.2e\n",
        res$method[i], or, ci_low, ci_high, res$pval[i]))
}

# 敏感性分析
cat("\n--- 敏感性分析 ---\n")
cat("异质性检验 p:", het$Q_pval, "\n")
cat("MR-Egger截距 p:", pleio$pval, "\n")