孟德尔随机化分析怎么做：从 GWAS 数据到敏感性分析

# 1. 在 https://api.opengwas.io/profile 登录并生成令牌（有效期 14 天）
# 2. 写入 ~/.Renviron（或项目目录下的 .Renviron），变量名必须是 OPENGWAS_JWT
usethis::edit_r_environ()
# 在打开的文件里加一行：OPENGWAS_JWT=eyJhbGciOi...（不加引号，末尾留一个换行）

# 3. 重启 R 后核对
ieugwasr::get_opengwas_jwt()   # 返回一长串字符说明已读到
ieugwasr::user()               # 返回账号信息说明令牌有效；401 说明过期或复制错

# 4. 当前版本的参数名是 opengwas_jwt；旧教程里的 access_token 已不存在
exp <- TwoSampleMR::extract_instruments("ieu-a-300", p1 = 5e-8, r2 = 0.001, kb = 10000)

library(TwoSampleMR); library(data.table)   # fread 读 .gz 还需要 R.utils

## A. 暴露：作者网站或 GWAS Catalog 下载的文本摘要统计（以 GLGC 2013 LDL-C 为例）
ldl <- fread("jointGwasMc_LDL.txt.gz")
setnames(ldl, c("P-value", "Freq.A1.1000G.EUR"), c("pval", "eaf"))
class(ldl$pval)              # 本页实测为 "character"：文件里有 1.24e-652 这类超出 double 范围的值
ldl[, pval := as.numeric(pval)]                 # 不转换时 pval < 5e-8 是字符比较，筛选失效
ldl[, `:=`(A1 = toupper(A1), A2 = toupper(A2))] # GLGC 的等位基因是小写，A1 为效应等位基因
exp_dat <- format_data(as.data.frame(ldl[pval < 5e-8]), type = "exposure",
  snp_col = "rsid", beta_col = "beta", se_col = "se", eaf_col = "eaf",
  effect_allele_col = "A1", other_allele_col = "A2", pval_col = "pval", samplesize_col = "N")

## B. 结局：GWAS Catalog 协调后文件（*.h.tsv.gz，GRCh38，hm_ 前缀列已统一方向）
# 先用 awk 只保留工具变量的行，370 MB 文件约 80 秒，R 里不用读 860 万行
#   gzip -dc GCST003116.h.tsv.gz | awk -F'\t' 'NR==FNR{a[$1];next} FNR==1 || ($2 in a)' snps.txt - > cad_subset.tsv
cad <- fread("cad_subset.tsv")
out_dat <- format_data(as.data.frame(cad), type = "outcome", snps = exp_dat$SNP,
  snp_col = "hm_rsid", beta_col = "hm_beta", se_col = "standard_error",
  eaf_col = "hm_effect_allele_frequency", effect_allele_col = "hm_effect_allele",
  other_allele_col = "hm_other_allele", pval_col = "p_value", chr_col = "hm_chrom", pos_col = "hm_pos")

## C. IEU OpenGWAS 的 VCF：FORMAT 为 ES:SE:LP:AF:SS:ID，LP 是 -log10(P)，效应等位基因是 ALT
#   bcftools query -f '%ID\t%CHROM\t%POS\t%ALT\t%REF\t[%ES]\t[%SE]\t[%LP]\t[%AF]\t[%SS]\n' ieu-a-300.vcf.gz > ieu-a-300.tsv
v <- fread("ieu-a-300.tsv", col.names = c("SNP","chr","pos","effect_allele","other_allele","beta","se","lp","eaf","samplesize"))
exp_vcf <- format_data(as.data.frame(v), type = "exposure", pval_col = "lp", log_pval = TRUE)

## D. Pan-UKB：P 值列是 neglog10_pval_EUR，同样用 log_pval = TRUE；效应等位基因是 alt，坐标 GRCh37

# FinnGen 文件带 tabix 索引，只取工具变量所在位置，不用下载整个 810 MB 文件
# regions_hg38.tsv 两列：染色体（1-23，不带 chr）与 GRCh38 位置；可从 GWAS Catalog 协调后文件的 hm_chrom、hm_pos 得到
tabix -h -R regions_hg38.tsv \
  https://storage.googleapis.com/finngen-public-data-r13/summary_stats/finngen_R13_I9_CHD.gz \
  > finngen_chd_subset.tsv
# 列：#chrom pos ref alt rsids nearest_genes pval mlogp beta sebeta af_alt ...（alt 为效应等位基因）
# 本页实测：77 个位置，4 分 26 秒，返回 79 行（含多等位位点）；当前目录会留下 .tbi 索引文件

# 参考面板：http://fileserve.mrcieu.ac.uk/ld/1kg.v3.tgz（1.57 GB，含 EUR/EAS/AFR/AMR/SAS）
# 只解压需要的人群：tar -xzf 1kg.v3.tgz EUR.bed EUR.bim EUR.fam
# plink 1.9 可用 conda 安装（bioconda::plink）或 genetics.binaRies::get_plink_binary()
library(ieugwasr)
# 先看有多少显著 SNP 不在参考面板里：这些 SNP 会被 ld_clump 直接删除
bim <- data.table::fread("ref/EUR.bim", select = 2)
miss <- exp_dat[!exp_dat$SNP %in% bim$V2, ]
nrow(miss); min(miss$pval.exposure)   # 本页实测 132 个，最小 P 为 4.3e-127

clumped <- ld_clump(
  dplyr::tibble(rsid = exp_dat$SNP, pval = exp_dat$pval.exposure, id = exp_dat$id.exposure),
  clump_kb = 10000, clump_r2 = 0.001, clump_p = 1,
  bfile = "ref/EUR",                  # 写到文件名前缀，不带 .bed
  plink_bin = Sys.which("plink"))
exp_dat <- exp_dat[exp_dat$SNP %in% clumped$rsid, ]

# 单 SNP F 统计量（近似式）
exp_dat$F <- (exp_dat$beta.exposure / exp_dat$se.exposure)^2

# 解释方差 R2（含 SE 与 N 的写法；beta 为 SD 单位时可直接用于连续性状）
b <- exp_dat$beta.exposure; s <- exp_dat$se.exposure
f <- exp_dat$eaf.exposure;  N <- exp_dat$samplesize.exposure
exp_dat$R2 <- 2*b^2*f*(1-f) / (2*b^2*f*(1-f) + 2*s^2*N*f*(1-f))

# 整体 F：k 个 SNP、总 R2
k <- nrow(exp_dat); R2 <- sum(exp_dat$R2, na.rm = TRUE)
F_all <- R2 * (median(N) - k - 1) / (k * (1 - R2))
# 本页实测（79 个 SNP）：单 SNP F 最小 27.8，中位 58.4；总 R2 0.085；整体 F 201

# 结局里缺某个工具变量时，在本地参考面板里找 r2 >= 0.8 的代理 SNP
plink --bfile ref/EUR --r2 in-phase with-freqs \
  --ld-snp rs646776 --ld-window-kb 500 --ld-window 99999 --ld-window-r2 0.8 \
  --out proxy_rs646776
# 输出 PHASE 列如 CG/TA：目标 SNP 的 C 等位基因与代理 SNP 的 G 等位基因在同一单倍型上
# 用代理 SNP 的结局效应时，按这个对应关系把效应等位基因换回目标 SNP 的等位基因
# 本页实测：单个 SNP 约 5 秒，rs646776 找到 10 个代理（含 r2=1 的 rs12740374）

dat <- harmonise_data(exp_dat, out_dat, action = 2)
attr(dat, "log")                      # 记录交换、剔除的 SNP 数，写进补充材料
dat <- dat[dat$mr_keep, ]             # 被剔除的行仍在数据框里，导出补充表前先过滤

res <- mr(dat, method_list = c("mr_ivw", "mr_ivw_fe", "mr_egger_regression",
                               "mr_weighted_median", "mr_weighted_mode"))
generate_odds_ratios(res)
mr_heterogeneity(dat)                 # Cochran Q
mr_pleiotropy_test(dat)               # MR-Egger 截距
Isq(dat$beta.exposure, dat$se.exposure)   # I2GX，低于 0.9 时 MR-Egger 有回归稀释
loo <- mr_leaveoneout(dat); ss <- mr_singlesnp(dat)

# Steiger：结局为二分类时先算 r.outcome，否则按连续性状近似
dat$r.outcome  <- get_r_from_lor(dat$beta.outcome, dat$eaf.outcome,
                                 ncase = 60801, ncontrol = 123504, prevalence = 0.05)
dat$r.exposure <- get_r_from_bsen(dat$beta.exposure, dat$se.exposure, dat$samplesize.exposure)
directionality_test(dat)

# MR-PRESSO：NbDistribution 至少为 SNP 数的 20 倍，否则离群检验精度不足
library(MRPRESSO)
nb <- max(1000, ceiling(nrow(dat) * 20 / 1000) * 1000)
pr <- mr_presso(BetaOutcome = "beta.outcome", BetaExposure = "beta.exposure",
  SdOutcome = "se.outcome", SdExposure = "se.exposure", data = dat,
  OUTLIERtest = TRUE, DISTORTIONtest = TRUE, NbDistribution = nb, seed = 2026)

# 作图
mr_scatter_plot(res, dat); mr_forest_plot(ss); mr_leaveoneout_plot(loo); mr_funnel_plot(ss)

# 同一暴露、两个独立结局来源的 IVW 估计合并（meta 包 8.5）
library(meta)
m <- metagen(TE = c(0.4123, 0.2741), seTE = c(0.0511, 0.0318),
             studlab = c("CARDIoGRAMplusC4D 2015", "FinnGen R13 I9_CHD"),
             sm = "OR", common = TRUE, random = TRUE)
forest(m)
# 本页实测：共同效应 OR 1.37 (1.30-1.44)；随机效应 1.40 (1.22-1.60)；Q = 5.28，P = 0.022，I2 = 81%