主题
02 · 工具变量选择
对应模块:R/04_iv_select.R参考:Zheng 2020;Staiger & Stock 1997;Pierce 2011
先说清楚:本课题的 cis-pQTL 筛查用的是「现成哨兵」,不走本页的 clumping
本页描述的是通用全量引擎(run_all.R)从原始 sumstats 里选工具的完整四步。但本课题的 主筛查(analysis/01_screen.R)不走这条——工具变量直接取自 UKB-PPP (Sun et al., Nature 2023) 的显著 cis-pQTL 表(Supplementary Table 9),每蛋白一个 LD 独立的哨兵变异。
- 数据实测:99.9% 的蛋白仅含 1 个 cis 工具,每「蛋白×基因组区域」恰好 1 个变异,无蛋白含 ≥3 个。
- 据论文方法学描述(经检索,尚未直接读到 Nature 原文全文,以下为原文摘要转述):这些 primary 关联系对全基因组显著变异做 PLINK 区域化 clumping(±1Mb,排除 HLA) 后的代表变异—— 也就是说,LD 独立性已由原文上游完成,本流水线直接采用其哨兵,不再二次 clumping。
- 论文另报告约 87% 的 cis 位点存在条件独立的 secondary 信号;本课题为保守起见仅用 primary 哨兵。
因此:审稿人若问"工具有没有做 LD 处理",答案是 做了,在数据源 UKB-PPP 上游用 PLINK 完成, 不是我们遗漏。本页的 clumping 步骤只在通用引擎(换非 pQTL 暴露、从原始 sumstats 起跑)时才触发。
作用
从标准化后的暴露数据里,选出满足 MR 相关性假设的工具变量(IV)。编排函数 select_instruments() 依次做四步,并全程记录 SNP 流水账(STROBE-MR item 10a 流程图):
raw → 显著性 → cis 窗口 → LD clumping → F≥10 剩余1. 显著性阈值(apply_significance)
按分析模式取阈值(iv: 段可调):
cis模式:pval_cis = 5e-6(cis 区域先验强,阈值放宽)genome_wide模式:pval_genome_wide = 5e-8(全基因组显著)
2. cis 窗口过滤(apply_cis_window,仅 cis 模式)
以基因的 TSS/TES 为中心,取 ±cis_window_kb(默认 500kb,可改 1000=±1Mb)内的 SNP。窗口由 resolve_cis_window()(01_config.R)解析:优先用暴露条目里显式的 chr/tss/tes,否则查 gene_annotation 按基因名定位。
window = [ min(tss,tes) − 窗口 , max(tss,tes) + 窗口 ]3. LD clumping(clump_iv)
去掉连锁不平衡(LD)冗余,保留近独立的工具。两种后端:
local:本地 PLINK 二进制 + 1000G LD 面板(/mnt/d/mrdata/ldpanel/EUR),走ieugwasr::ld_clump;api:OpenGWAS 在线 clumping。
参数 clump_r2 = 0.001、clump_kb = 10000。
LD 面板祖先必须匹配暴露人群(STROBE-MR)
clumping 用的 LD 结构必须与暴露 pQTL/GWAS 的人群祖先一致。用欧洲人面板去 clump 东亚/非洲人群的 pQTL,会因连锁结构错配选出有偏的工具变量。
为此每个暴露须声明 ancestry,ld 段须声明面板祖先(local 用 ld$ancestry,api 用 iv$clump_pop)。clump_iv() 在 clumping 前调用 assert_ld_ancestry_match():
- 暴露祖先 ≠ 面板祖先 → 直接报错拦截(
01_config.R加载配置时也会 fail-fast 提前拦下); - 暴露未声明
ancestry→ 仅告警放行(向后兼容)。
换非欧暴露:把 bfile 换成对应人群面板(如 /mnt/d/mrdata/ldpanel/EAS)、ld$ancestry 改成 EAS,并把暴露 ancestry 设为 EAS。单一全局 LD 配置只能服务同祖先的暴露;多祖先请分 config 分别跑。
失败守卫
clumping 失败时不中断:回退用未 clump 的 SNP,并打上 unclumped 标记 —— 报告会明确标注"结果可能因 LD 冗余而膨胀",而不是悄悄给出误导结果。SNP ≤1 个时直接跳过 clumping。
4. F 统计量与 R²,剔除弱工具(compute_f_stats)
衡量工具强度(相关性假设)。有 EAF 时按标准公式逐 SNP 计算:
每个 SNP 的方差解释比例 R²(maf = min(eaf, 1−eaf),n = 样本量):
每个 SNP 的 F 统计量:
整体 F(k 个强工具,r2_sum = ΣR²,n_mean 平均样本量):
弱工具剔除:去掉 F < f_stat_min(默认 10)的 SNP —— F<10 是弱工具偏倚的经典阈值(Staiger & Stock)。
EAF 缺失时的代理
若整列 EAF 缺失,无法算 R²,则用 (beta/se)² ≈ z² 作每-SNP F 的代理(此时不做弱工具剔除,仅提示)。
对应 config 字段
yaml
analysis_mode: "cis" # cis | genome_wide
iv:
pval_cis: 5.0e-6
pval_genome_wide: 5.0e-8
cis_window_kb: 500 # ±500kb
clump_r2: 0.001
clump_kb: 10000
clump_pop: "EUR" # api 模式的 LD 面板人群;须与暴露 ancestry 一致
f_stat_min: 10 # 弱工具阈值
ld:
method: "local" # local(PLINK) | api(OpenGWAS)
plink_bin: "plink1.9"
bfile: "/mnt/d/mrdata/ldpanel/EUR" # LD 面板 bfile 前缀
ancestry: "EUR" # local 模式必填:该面板人群祖先,须与暴露 ancestry 一致每个暴露/结局条目也需声明 ancestry(如 EUR/EAS/AFR),既供 STROBE-MR 报告样本祖先,也用于 clumping 前的面板匹配校验。
输出
一张工具变量表(含每 SNP 的 f_stat、r2),并在属性里带上 overall_f、overall_r2、unclumped 标记,以及 flow 计数(raw → significant → cis_window → clumped → strong_F)供报告画流程图。
下一步:工具变量与结局数据一起进入 谐化与定向过滤。