主题
03 · 第 1 步 · 工具选择
执行日期：2026-08-07 状态：✅ 完成 —— 全部命中预期，16/16 产物与旧结果 MD5 相同
本页上半部分是在运行之前写的
「数据 / 版本 / 方法 / 预期结果」在跑之前落文，跑完只填「实际结果 / 解读 / 去向」。 ★ 跑前写下的预期与实际不符，一律先按缺陷处理，不许现场解释掉。
第 1 步与第 2 步在同一个脚本里跑
analysis/01_screen.R 一次运行同时完成第 1 步（工具选择：F/R²、Steiger filtering） 与第 2 步（MR + FDR）。代码结构与流程步骤不是一一对应，处理办法是 不动代码，一次运行产出两份记录，分两次检查——记录粒度可以细于运行粒度。 第 2 步的记录见 04 页。
结局命名
本页表格用内部 ID（T1D_gcst / T2D_mahajan 等）——它们是文件名与代码里的字符串，实录必须能对得上。论文显示名分别是 Type 1 diabetes / Type 2 diabetes，对照表见总览页。★ 内部 ID 绝不允许出现在图上或论文表格里。
一、这一步做什么
按 v4：① 统一坐标 build（必须最先）→ ② cis 窗口 ±500 kb → ③ F ≥ 10 · clumping · Steiger filtering。 作用单位是变异，全部发生在 MR 之前。
★ 不设 MHC 排除步骤（v4 §3.3）。
二、用什么数据、什么版本
| 项 | 值 |
|---|---|
| 暴露文件 | data/exposure/protein_info.csv → 软链至 /mnt/d/mrdata/exposure/protein_info.csv |
| 来源 | UKB-PPP（Sun BB et al. Nature 2023;622:329–338，PMID 37794186） |
| 取哪一列效应量 | BETA (discovery, wrt. A1) / SE (discovery) / log10(p) (discovery) —— discovery 集 |
| 对应样本量 | 34,557（欧裔 discovery；必须与所取效应量列一致） |
| 频率列 | A1FREQ (discovery) |
| 位置列 | GENPOS (hg38) ⚠️ 同表另有 hg37 坐标藏在变异 ID 里，取错列不报错 |
| 变异 ID 列 | Variant ID (CHROM:GENPOS (hg37):A0:A1:imp:v1) —— 只取 A0/A1，前两段 hg37 坐标弃用 |
| 筛选 | cis/trans == "cis" |
| build | GRCh38（与全部结局一致；LD 面板是 hg19 但按 rsID 匹配，无需 liftover） |
代码基线：F:\project\mr-pipeline-r9v3 提交 efb6584；config.yaml resume: false、缓存指纹 v3。
三、什么检验方法
3.1 ① 统一 build —— 已在数据层完成
全部数据源统一到 GRCh38，唯一例外 Mahajan T2D（GRCh37 且无 rsID）另由 tools/lift_mahajan_hg38.R 两跳桥接处理（共定位用）；在本步它走 position 匹配。
★ position 匹配键用「排序后的等位」而非原始 A0:A1 顺序： 固定顺序键要求外部数据集的 NEA:EA 与 UKB-PPP 的 A0:A1 恰好同序，否则整条漏掉—— 实测 Mahajan 因此漏掉 700/1656 个工具（42%）。排序等位使键与效应方向无关， 方向留给 harmonise 对齐。
3.2 ② cis 窗口
config.yaml iv.cis_window_kb: 500（±500 kb）[Zheng 2020]。
⚠️ 本步实际不重新划窗：工具直接采用 UKB-PPP 已判定为 cis 的哨兵 （原文 cis 定义为哨兵位于编码基因 ±1 Mb 内）。±500 kb 用于共定位取区域。 实测 cis 哨兵距基因均 < 500 kb，故两者不冲突；若需与原文 cis 定义完全对齐可放宽至 ±1 Mb。
3.3 ③ F ≥ 10 · clumping · Steiger filtering
F / R²（每工具单独算）[Pierce 2011]：
R² = 2β²·maf(1−maf) / [ 2β²·maf(1−maf) + 2N·maf(1−maf)·se² ]
F = R²(N−2) / (1−R²)clumping —— 本流水线不做二次 clumping，因为 UKB-PPP 上游已完成 （±1 Mb PLINK clumping，HLA 区整块当一个 locus，取 P 最小者为哨兵）。 这一点已按原文重写进 config.yaml 的 instrument_definition。
Steiger filtering [Hemani 2017]：config.yaml iv.steiger_filtering: true， 由 apply_steiger() 真正执行。 ⚠️ 二分类结局的 R² 走 get_r_from_lor()，需 units="log odds" + ncase/ncontrol/prevalence （2026-08-01 修；此前对数几率被当连续量算 R²，结局 R² 被低估约 11 倍）。
3.4 共用哨兵 SNP 的校正
protein_info.csv 里有 11 个 rsID 被 2 个蛋白共用（如 rs1859788 = PILRA + PILRB）， 但两蛋白的 cis 效应量并不相同。standardize_dataset 的 !duplicated(SNP) 只留第一个蛋白的 β/se，随后 ann 合并又把它扇给同组另一个蛋白 → 后者的 Wald 比用了别人的暴露效应量 （实测 11 组的 β 全不同，其中 2 组符号相反）。
fix_shared_instrument() 按真实 βexp 等比例改正。 ★ z = b/se = βout/seout 与 βexp 无关 → p 值、FDR、显著性判定不受影响； 受影响的只有 OR/CI 数值（符号相反那 2 组连方向也会翻转）。
四、★ 预期结果（跑之前写下）
4.1 工具池
| 项 | 预期 | 依据 |
|---|---|---|
| 工具表 cis 行数 | 1,955 | 与 Sun 2023 正文逐位吻合 |
| cis 唯一蛋白 | 1,954 | 同上 |
其中 MHC（作者 MHC 列 = 25.5–34.0 Mb） | 27 | _mhc.R 断言 |
| 标准化去重后 cis 工具 SNP | 1,875 | 旧运行日志 [R9_dm] cis 工具 SNP: 1875 |
4.2 每结局的工具与 F
应与旧结果逐位相同（同数据、同参数）：
| outcome | 行数 | 唯一 SNP | 唯一蛋白 | F 最小 | F < 10 |
|---|---|---|---|---|---|
| Retinopathy | 1587 | 1576 | 1587 | 43.8 | 0 |
| Maculopathy | 1587 | 1576 | 1587 | 43.8 | 0 |
| NeovascGlaucoma | 1587 | 1576 | 1587 | 43.8 | 0 |
| Neuropathy | 1587 | 1576 | 1587 | 43.8 | 0 |
| Nephropathy | 1587 | 1576 | 1587 | 43.8 | 0 |
| T1D_gcst | 1722 | 1712 | 1722 | 43.8 | 0 |
| T2D_mahajan | 1614 | 1603 | 1614 | 43.9 | 0 |
★ F 最小 43.8，无一低于 10 —— 即 F ≥ 10 这道门槛在本数据上不淘汰任何工具。 写 Methods 时必须如实这样说，不能写成"我们用 F≥10 筛掉了弱工具"。
★ Steiger filtering 预期剔除 0 个 SNP：cis-pQTL 暴露 R²（中位 2.2e-2） 远高于结局 R²（约 7e-6）。保留该步是为方法学描述可核，不是因为它起了作用。 实测（tools/steiger_prevalence_sens.R）prevalence 在 0.001–0.2 扫描下判定 0 处变化。
4.3 容差
| 项 | 容差 |
|---|---|
| 工具数、每结局行数/唯一 SNP/唯一蛋白 | 完全一致，差 1 即缺陷 |
| F / R² | 浮点精度内一致 |
| F < 10 的数量 | 必须为 0 |
⚠️ 若不一致，按此顺序排查：① 输入文件被动过 ② 包版本变了 ③ 代码改动的副作用。 不许先假设是正常波动。
4.4 预期产物
results/screen_R9_dm/
├ <outcome>_all.csv 7 个（含 SNP/protein/R2/Fstat/两侧等位/b/se/p/FDR/OR）
├ <outcome>_all.csv.fp 7 个缓存指纹（本轮 resume=false，仍会写）
├ <outcome>_significant.csv
└ SUMMARY_publication.csv五、实际结果
运行：Rscript analysis/01_screen.R R9_dm，退出码 0，10:13:27 → 10:22 左右。 日志：results/logs/screen_R9_dm_20260807.log
5.1 工具池 —— 全部命中预期
| 项 | 预期 | 实际 | |
|---|---|---|---|
| 工具表 cis 行数 | 1,955 | 1,955 | ✅ |
| standardize 丢弃 | — | 80 | — |
| 暴露 cis 工具 SNP | 1,875 | 1,875 | ✅ |
5.2 ★ 「丢弃 80 行」的完整拆解（跑前记录里没拆，这里补上）
复刻 standardize_dataset 的过滤链，逐步实测：
| 过滤步 | 剩余 | 丢弃 |
|---|---|---|
起点（cis/trans == "cis"） | 1955 | — |
① 缺失 / se > 0 | 1955 | 0 |
② 等位合法 ^[ACGT]+$ | 1955 | 0 |
③ EA != OA | 1955 | 0 |
④ !duplicated(SNP) | 1875 | 80 |
→ 80 行全部来自去重，其中：
- 69 行：
rsID是占位符"-"（工具表里共 70 行如此），去重只留 1 行 - 11 行：11 个 rsID 被 2 个蛋白共用，去重各留 1 行
共用哨兵 11 组明细（两蛋白的 cis 效应量并不相同，★ 其中 2 组符号相反）：
| rsID | 蛋白 A (β) | 蛋白 B (β) |
|---|---|---|
| rs11126696 | REG1A (0.251) | REG1B (0.345) |
| rs12459454 | CEACAM6 (−0.120) | CEACAM8 (−0.131) |
| rs1801689 | APOH (−1.462) | ERN1 (0.526) ← 符号相反 |
| rs1859788 | PILRA (1.000) | PILRB (1.147) |
| rs2236352 | PSME1 (0.094) | PSME2 (0.104) |
| rs28403575 | AP1G2 (−0.192) | THTPA (−0.750) |
| rs367070 | LILRA3 (−1.219) | LILRB1 (−0.671) |
| rs4244437 | IL12A_IL12B (0.543) | IL12B (0.560) |
| rs4804774 | CD209 (−0.414) | CLEC4M (−0.364) |
| rs4821544 | PVALB (−0.813) | TST (−0.108) |
| rs5030377 | ICAM1 (0.415) | ICAM4 (−0.115) ← 符号相反 |
fix_shared_instrument() 每个结局校正 11 行 Wald 比（日志实测），与预期一致。
5.3 ⚠️ 查出并已判定的一个隐患：占位符 "-" 被当成工具带进池子
事实：去重后有 1 行 rsID = "-" 留在 1,875 个工具里。
调查（不是推断，是实测）：
| 检查 | 结果 |
|---|---|
结局侧白名单 grep 用的是 grep -Fw -f <白名单> | 构造属实 |
grep -Fw 对 - 是否会匹配整字段为 - 的行 | 会（造样本实测确认） |
FinnGen R9 的 rsids 列是否用 - 作缺失占位 | 不用 —— 前 300 万行里 rsids == "-" 的有 0 行，缺失一律是空字符串 |
| 实际预过滤读入行数 | 1,663 行（若 - 命中会是几十万行） |
结果表里 SNP == "-" 的行数 | 0（01_screen.R:201 有 [grepl("^rs", SNP)] 过滤） |
判定：对本轮结果影响为零。 但有两点必须记下：
- 保护是偶然的，不是设计的。 结果之所以干净，靠的是 ① FinnGen 恰好不用
-作占位符 ② 下游有个grepl("^rs", SNP)过滤。 换一个用-作占位的结局数据源，白名单就会命中大量无关行， 经去重后与一个任意变异配成对。 - 报告的工具数 1,875 里含 1 个非工具。 写 Methods 时的诚实数字是 1,874 个可用工具，或写 1,875 并说明其中 1 个占位记录在下游被剔除。
★ 这条不影响任何数字，但属于换数据即爆的类型，记入待办。
5.4 逐结局的工具与 F —— 全部命中预期
| outcome | 行数 | 唯一 SNP | 唯一蛋白 | F 最小 | F 最大 | F<10 | 预期 |
|---|---|---|---|---|---|---|---|
| Retinopathy | 1587 | 1576 | 1587 | 43.8 | 23802.7 | 0 | ✅ |
| Maculopathy | 1587 | 1576 | 1587 | 43.8 | 23802.7 | 0 | ✅ |
| NeovascGlaucoma | 1587 | 1576 | 1587 | 43.8 | 23802.7 | 0 | ✅ |
| Neuropathy | 1587 | 1576 | 1587 | 43.8 | 23802.7 | 0 | ✅ |
| Nephropathy | 1587 | 1576 | 1587 | 43.8 | 23802.7 | 0 | ✅ |
| T1D_gcst | 1722 | 1712 | 1722 | 43.8 | 23802.7 | 0 | ✅ |
| T2D_mahajan | 1614 | 1603 | 1614 | 43.9 | 23802.7 | 0 | ✅ |
5.5 完整 SNP 流（以 Retinopathy 为例，日志实测）
工具表 cis 行 1955
↓ standardize 丢弃 80（全部是去重）
暴露 cis 工具 SNP 1875
↓ 结局侧 rsID 预过滤 (grep -Fw) 1663
↓ 结局 standardize 丢弃 10 1653
↓ 谐化：剔除 48 个「模糊回文」SNP 1576 ★跑前记录漏了这一道
↓ Steiger filtering 1576 ★剔除 0，与预期一致
↓ 共用哨兵校正 11 行 Wald 比被校正★ 两处补记：
- 模糊回文剔除 48 个：
harmonise_action: 2（按 MAF 推断链方向）下， A/T、C/G 且 eaf 落在 0.42–0.58 的位点靠频率也判不了链方向，必须剔除。 这一道我在跑前预期里漏写了 —— 它不属于 Steiger，属于谐化。 - Steiger 单位声明正确：日志
exposure=SD outcome=log odds，rsq.outcome中位 8.41e-06（Retinopathy），与 2026-08-01 修复后的预期量级一致 （此前对数几率被当连续量算，结局 R² 被低估约 11 倍）。
5.6 ★★ 与旧结果的比对：16/16 产物 CSV 的 MD5 完全相同
| 项 | 结果 |
|---|---|
results/screen_R9_dm/ 下 csv 文件数 | 16 |
| MD5 与旧结果相同 | 16 |
| MD5 不同 | 0 |
| 旧目录缺失 | 0 |
包括 7 个 _all.csv、7 个 _significant.csv、SUMMARY_publication.csv、overlap_matrix.csv。
六、如何解读
6.1 这一步的结论
工具选择完全复现：resume: false 从零重算，与三个月前的旧结果MD5 逐字节一致。 加上第 0 步 LDSC 的同样结果，可复现性证明已覆盖数据验收 + 工具选择 + MR 三层。
6.2 ★ 两个必须如实写进 Methods 的"门槛没起作用"
| 门槛 | 实测 | 该怎么写 |
|---|---|---|
| F ≥ 10 | F 最小 43.8，淘汰 0 个 | 「所有工具 F 值介于 43.8–23802.7，均远超常用的 F ≥ 10 门槛」——不能写成"我们筛掉了弱工具" |
| Steiger filtering | 剔除 0 个 | 「已执行 Steiger 定向过滤，未剔除任何工具；这与 cis-pQTL 暴露 R²（中位 2.2e-2）远高于结局 R²（中位 8.4e-06）一致」 |
真正起了作用的是模糊回文剔除（48 个），它反而没在流程图上出现。这条要补进 Methods。
6.3 一处方法学上值得说的地方
standardize 的 4 道过滤里，①②③ 各丢 0 行——即缺失、非法等位、EA==OA 在 UKB-PPP 哨兵表里一个都没有。这是数据源质量的旁证，值得在 Methods 用一句话带过。
七、进正文还是补充（暂定）
| 内容 | 暂定去向 | 依据 |
|---|---|---|
| 工具选择流程与参数 | 正文 Methods | 必写：数据源、cis 定义、F 计算式、Steiger、不做二次 clumping 的理由 |
| "F 最小 43.8、无工具被 F≥10 淘汰" | 正文 Methods 一句 | 如实交代门槛未起作用，比含糊写"应用了 F≥10"诚实 |
| "Steiger 剔除 0 个" | 正文 Methods 一句 | 同上 |
| 每工具 F/R² 全表 | 补充材料 | 审稿人常索要 |
| 共用哨兵校正（11 组，2 组符号相反） | 补充方法一段 | 不写会被质疑 OR 与 β 对不上 |
| MHC 区哨兵选择规则不同质 | Limitations | 见 00 页 §6.3 |
⚠️ 暂定，第 7 步统一复核。