01 · 共享函数与数据读取

对应脚本：common/mr_funcs.R被谁用：01/02/03/05 各步全部 source() 它。改这里会影响所有分析——动它要谨慎。

作用

集中放"每一步都要用"的东西：结局数据配置、按 rsID 从大文件里抠数、Z-score→β 换算、LD clump、单组 MR、共定位。

结局配置 `OUTCOMES`

一个 list，每个结局记 file / ncase / ncontrol / disease，IAMDGC 额外标 fmt="iamdgc"（格式不同要走另一套解析）。

OUTCOMES <- list(
  DR_strict = list(file=".../finngen_R13_DM_RETINOPATHY_STRICT.gz", ncase=6970,  ncontrol=88076,  disease="DR"),
  AMD       = list(file=".../finngen_R13_H7_AMD.gz",                ncase=13947, ncontrol=458803, disease="AMD"),
  WetAMD    = list(...), DryAMD = list(...),
  AMD_iamdgc= list(file=".../IAMDGC_lateAMD_related_EUR_annotated.txt", ncase=15616, ncontrol=16723, disease="AMD", fmt="iamdgc")
)

换数据只改这里（详见结果页"换数据怎么改"）。disease 字段是 FDR 分家的依据——DR 一族、AMD 一族。

按 rsID 提取结局 `extract_finngen` / `extract_iamdgc`

避免把 800MB 结局整表读进内存：用 awk 只捞工具变量那几百个 rsID 的行。

# FinnGen：rsids 在第 5 列
cmd <- sprintf("zcat %s | awk -F'\\t' 'NR==FNR{a[$1]=1;next} ($5 in a){print}' %s -", oc$file, tmp)

等位一致性（关键正确性）：

FinnGen effect_allele = alt、eaf = af_alt —— eaf 对应的正是 effect allele，一致 ✓
IAMDGC 是 REGENIE 输出，effect_allele = ALLELE1、A1FREQ 是 ALLELE1 的频率、BETA 是每 ALLELE1 效应 —— 全部对齐到 ALLELE1 ✓（这是 IAMDGC 最容易搞错的点，代码是对的）

已知限制（诚实记录）

FinnGen 多等位位点的 rsids 列可能是逗号连的多个 rsID（如 rs1,rs2），($5 in a) 精确匹配整字段 → 这类行会被漏掉。丢的多是罕见变异，量小，但发表方法里应提一句。

Z-score → β/se 换算 `z_to_beta`

eQTLGen 只给 Z 分数，没给 β。用等位频率 + 样本量近似（Zhu 2016 / Zheng 2020）：

\hat{β} = \frac{Z}{\sqrt{2 p (1 - p) (N + Z^{2})}}, se (\hat{β}) = \frac{1}{\sqrt{2 p (1 - p) (N + Z^{2})}}

z_to_beta <- function(z, p, n) {
  denom <- sqrt(2 * p * (1 - p) * (n + z^2))
  list(beta = z / denom, se = 1 / denom)
}

$p (1 - p)$ 对等位对称，所以用哪个等位的频率结果一样——不怕填反。

LD clump `clump_local`

用本地 plink1.9 + 1000G EUR 参考做独立工具筛选（不走 OpenGWAS API，避免跨墙/限流）。

system2(PLINK, c("--bfile", LDREF, "--clump", ...,
                 "--clump-r2", r2, "--clump-kb", kb, ...))

eQTL/视网膜：r2<0.1, 1Mb（cis 区独立工具）
代谢物：r2<0.001, 10Mb（全基因组，更严）

plink 失败时保留全部 SNP（if(!file.exists(cf)) return(df)）并在报告里标 caveat —— 这是"失败不静默"的设计。

单组 MR `run_mr_one`

按工具数自适应选方法（和主流水线一致）：

工具数	方法
1	Wald ratio
2	IVW
≥3	IVW + MR-Egger + 加权中位数 + 异质性 Q + 多效性截距 + Steiger

Steiger 定向、Egger 截距、Cochran Q 都算好挂在结果属性上，供下游判 robust。

共定位 `run_coloc`

coloc.abf（Giambartolomei 2014）判"暴露和疾病是不是同一个因果变异"。

三个实现要点：

两侧按 rsID 去重（coloc.abf 遇重复 SNP 直接报错）—— 暴露保留最显著行、结局保留第一行。这是踩过坑后加的。
等位翻转：effect/other 互换时把结局 β 取负对齐。
暴露当 type="quant"、疾病当 type="cc"（带病例比例 s）。

d1 <- list(beta=..., varbeta=se^2, type="quant", N=n_exp, MAF=pmin(eaf,1-eaf))
d2 <- list(beta=..., varbeta=se^2, type="cc", N=ncase+ncontrol, s=ncase/(ncase+ncontrol), MAF=...)
cr <- coloc.abf(d1, d2)

nsnp<5 直接返回 PP.H4=NA（不硬算假值）—— 视网膜层大量返回 NA 就是这里，属数据稀疏不是阴性。

金标准

下游判"因果稳健"= MR FDR<0.05 且 coloc PP.H4≥0.8 两条都满足（见 02）。

一个可讨论的假设

quant 暴露没传 sdY，coloc.abf 会用 MAF/N/β 自己估。标准做法，但属附加假设，方法里可提一句。

01 · 共享函数与数据读取 ​

作用 ​

结局配置 OUTCOMES ​

按 rsID 提取结局 extract_finngen / extract_iamdgc ​

Z-score → β/se 换算 z_to_beta ​

LD clump clump_local ​

单组 MR run_mr_one ​

共定位 run_coloc ​