主题
方法学审计与已知限制（2026-07-20）
这一页记录一次系统性方法学审计的全部发现：查了什么、发现了什么、改了什么、还剩什么限制。 投稿写 Methods 与 Limitations 时直接从这页取材。
审计边界要先说清楚：本轮覆盖了工具变量定义、谐化与等位方向、效应量与显著性口径、 多重检验、共定位候选选择、LDSC 参数、样本重叠、外部数据源核对。未覆盖：hyprcoloc/coloc 算法内部实现、PheWAS 外部库本身的质量、单细胞图谱的注释质量。
一、审计缘起
先前发现三处「静默失败被当成结论」的工程 bug（见阶段七的订正块）。 既然工程层出过这种问题，就有必要把方法学层也过一遍——毕竟前者只是让分析没跑成， 后者会让跑出来的数字本身不对。
二、已修复：影响已发表数值的问题
1. 共用哨兵 SNP 的蛋白，用错了暴露效应量 ⚠️ 影响候选靶点
protein_info.csv 里有 11 个 rsID 被 2 个蛋白共用（如 rs1859788 同时是 PILRA 和 PILRB 的 cis 哨兵），但两个蛋白的 cis 效应量并不相同。standardize_dataset 的 !duplicated(SNP) 只留下 文件里第一个蛋白的 β/se，随后按 SNP 合并注释时又把它扇给了同组另一个蛋白。
实测 11 组的 β 全都不同，其中 2 组符号相反：
| 受影响蛋白 | 真实 β | 误用的 β | 说明 |
|---|---|---|---|
| PILRB | 1.147 | 1.000（PILRA 的） | A 级候选靶点，OR 需按 ×0.872 修正 |
| ICAM4 | −0.115 | +0.415（ICAM1 的） | 符号相反 → 方向反转 |
| ERN1 | +0.526 | −1.462（APOH 的） | 符号相反 |
| REG1B / CEACAM8 / PSME2 / THTPA / LILRB1 / IL12B / CLEC4M / TST | 均不同 | — | 偏差 0.26–7.5 倍 |
为什么显著性没受影响
Wald 比 b = βout/βexp、se = seout/βexp，所以 z = b/se = βout/seout 与 βexp 无关。 ⇒ p 值、FDR、显著性判定、靶点入选全部不变；错的只有 OR 与置信区间的数值 （符号相反的那两组连方向也错）。
已修：按各蛋白真实 βexp 等比例校正 b/se，并用校正后的暴露量重算 R²/F。
2. 外部 T2D（Mahajan）按位置匹配时漏掉 42% 的工具
read_pos_outcome() 自己拼了个 chr:pos:NEA:EA 的顺序敏感键去匹配 UKB-PPP 的 A0:A1， 外部数据集只要等位顺序相反就整条漏掉。
| 匹配方式 | 命中工具数 |
|---|---|
固定顺序 A0:A1（旧实现） | 956 |
| 按排序等位（顺序无关） | 1656 |
R/02_standardize.R 里其实早就写了顺序无关的排序等位方案，还注明"比固定 NEA:EA 顺序更稳健"， 批量筛查这条线自己绕开了它。已改为与之同口径，方向交给后面的 harmonise 对齐。
3. LDSC 的 h² 用了总人数而非二分类有效样本量
munge 传的是 N = ncase + ncontrol。本课题各结局病例占比从 2.2% 到 19.7%， Neff/N 在 0.086–0.633 之间（相差 7.4 倍），等于每个性状的 h² 被一个各不相同的 因子压低 → h² 柱状图跨性状根本不可比。
已改为 Neff = 4/(1/ncase + 1/ncontrol)。改后 h² 不再被病例占比牵着走—— 病例占比最高的 Retinopathy（19.7%）h² 反而最低（0.051），而 T1D_wide（2.7%）最高（0.238）， 说明数值已由性状本身而非抽样结构决定。
| 性状 | 病例占比 | h²（Neff 口径） |
|---|---|---|
| T1D_wide | 0.027 | 0.238 |
| WetAMD | 0.022 | 0.178 |
| T2D | 0.185 | 0.162 |
| AMD | 0.030 | 0.127 |
| Retinopathy | 0.197 | 0.051 |
rg 是「近似不变」，不是严格不变——修正一处早先的判断
最初判断"rg 不受影响（分子分母同比例缩放）"。实测并非严格成立： AMD–WetAMD 0.925→0.931、WetAMD–DryAMD 0.851→0.865、AMD–DryAMD 0.966→0.972， 最大变动 0.014。原因是 LDSC 的回归权重本身依赖 N，不只是 h² 的等比例缩放。 Retinopathy–Nephropathy(0.884) 与 Retinopathy–Maculopathy(0.920) 则完全没变。 ⇒ 变动幅度 <2%，基于 rg 的定性结论不受影响，但论文里的 rg 数值须用新表。
4. Bonferroni/FDR 的检验数口径
FDR 原先在"把共用哨兵展开成各蛋白行"之前计算，检验数 m = 唯一变异数（1648）， 而汇总表报告的 n_instruments 是蛋白数（1659）。二者不一致，且 m 偏小 → FDR 偏松。 已把 FDR 移到展开之后，m = 该结局实际报告的蛋白数。
5. 外部 T1D 的病例/对照数写错了
outcome_manifest.csv 原写 18,942 / 501,638，而 GWAS Catalog 与数据文件本身都是 16,971 / 418,767（文件里是单一常数值，逐行核对无异）。这个数字会直接进投稿汇总表。已订正。
三、已修正：声明与实现不符
Steiger 定向过滤——声明了但从未执行
config.yaml 写着 steiger_filtering: true，但批量筛查 01_screen.R 直接调 harmonise_pair()， 绕过了 apply_steiger()；只有没用来出结果的 run_all.R 才走。
按 TwoSampleMR 同款口径复算：三个结局各 1659 个工具，方向错误 0 个 （暴露 R² 中位 2.2×10⁻² vs 结局 R² 中位 ~7×10⁻⁶，差四个数量级）。
⇒ 对结果零影响，但方法学描述不能写"做了 Steiger 过滤"。已让批量筛查真正执行该步， 使描述与实现一致。
四、设计层面的选择——影响候选靶点清单，值得重新考虑
只在 1 个疾病上显著的蛋白，从设计上进不了候选清单
05_hyprcoloc.R 只对跨 ≥2 个并发症显著的蛋白跑共定位（"至少 2 个疾病才值得多性状共定位"）， 而 A 级候选靶点必须有共定位支持 ⇒ 单病蛋白被结构性排除。
技术上并不需要这个限制——2 性状（1 蛋白 + 1 疾病）的 hyprcoloc 完全跑得通， 08_replicate_iamdgc.R 做的正是这个。
这会让候选清单系统性偏向"跨病共享"的蛋白，而一个只对 AMD 有强特异效应的蛋白反而出局。 被排除的蛋白里有若干生物学上很有说服力的 AMD 候选（如 MERTK——RPE 吞噬光感受器外节的 核心基因；CFP（备解素）——补体旁路通路，正是 AMD 的核心通路）。
这不是 bug，是选择；但要么改，要么在 Limitations 里写明
两条路：① 对单病蛋白补跑 2 性状共定位，把候选池扩起来；② 保留现设计，但在文中明确 "本研究的候选靶点定义要求跨 ≥2 个并发症共享，因此对疾病特异性蛋白不敏感"。
五、需在论文中如实交代的限定
1. LDSC 的遗传相关被共享对照结构性抬高
已用官方 FinnGen R13 manifest 逐个核对：Retinopathy 与 Nephropathy 的对照数 完全相同（62,519）——它们共用同一批对照。AMD / WetAMD / DryAMD 之间同样大量共用样本。
⇒ rg 0.85–0.97 不能当作"疾病间生物学同源"的独立证据，写作时必须限定为 "同一生物库、共享对照下的遗传相关"。
2. 工具变量的覆盖率
1,955 条 cis 记录 → 70 条 rsID = "-"（这 70 个蛋白从未进入分析）→ 1,874 个唯一 rsID → 谐化成功进入结果的 1,659 个。已核实被排除的 70 个蛋白不含任何候选靶点。
3. 阈值都是可辩护但需明示的选择
| 参数 | 取值 | 说明 |
|---|---|---|
| hyprcoloc 簇后验概率 | 0.7 | 低于常用的 0.8，config 里有声明 |
高多效标记 PLEIO_MAX | 100 | 硬编码在 06_network.R，不在 config，属任意选择 |
| FDR 分组 | by_outcome | 改成全局会让 R13 显著数从 333 降到 289 |
4. 样本重叠
本节结论已于 2026-07-29 更正，2026-07-31 补记
原文写的是「逐个核实后确认无重叠 ⇒ 全部结局与暴露样本独立」，这句话现在是错的， 已撤回。错因有两层：
- 它论证的是已退役的 T1D 源
GCST90475661（MVP）——该源经标志位点核验实为 T2D 主导，早已换成GCST90824163； - 换上去的新源本身含 UK Biobank（1,445 例 + 362,050 对照，占 44.5%）—— 我们此前写的「非 UKB、无重叠」对新源同样不成立。
更正后的口径（现行）：结局侧重叠 5–6%、病例侧 0.9%， 工具 F 值最低 40.6、中位 763，估计偏倚约 0.16%——可忽略，但必须在正文披露， 不能再写成「无重叠」。零重叠备选 Crouch 2025（GCST90013791）已下载并实测， 方向一致率 84.6%。详见外部 T1D 数据源与样本重叠。
现行逐结局核实：
- FinnGen（芬兰）、IAMDGC（Fritsche 2016 国际联盟）与 UK Biobank 无重叠 ✅
- 外部 T2D 用 Mahajan 2018 的 noUKBB 版本，原作者已剔除 UKBB ✅
- 外部 T1D
GCST90824163含 UKB，重叠 5–6% ⚠ 须披露
config.yaml 尚未同步
config.yaml 的 sample_overlap.note 仍是旧文本（还在写 GCST90475661 与 「全部结局与暴露样本独立」）。它会进报告模板，改稿前必须一并更新。
六、审计通过、确认没问题的部分
- 等位与方向处理：FinnGen AMD 与 IAMDGC AMD 的强信号（p<1e-5）14/14 方向一致， 谐化后效应等位逐个统一 ⇒ 等位对齐正确
- 效应量转换：
-log10(p) → p、OR → logOR、GCST 的se_from_ci重构均正确 - F/R² 公式、OR/CI 计算、
target_action方向逻辑（OR>1 → 抑制） mr_singlesnp的汇总行被^rs正确滤掉，未混入单 SNP 结果- 工具变量定义：实测 1953/1954 个蛋白恰好 1 个 cis 哨兵（99.95%），与 config 中 引用 UKB-PPP 方法学的描述一致
- PheWAS 的检验数：SNP 数 × OpenGWAS 数据集数，溯源字段完整，且有 "查询阈值比 Bonferroni 还严"的守卫
七、修复后全量重跑的结果变化（2026-07-20 实测）
01 筛查（R9+R13）与全部下游已用修复后代码重跑完毕。
变了什么
| 项目 | 修复前 | 修复后 | 说明 |
|---|---|---|---|
| 显著关联总数 | 616 | 662 | 新增 46 条，消失 0 条 |
| T2D_mahajan 工具数 | 955 | 1614 | 等位顺序修复 |
| T2D_mahajan 显著数 | 7 | 28 | 同上（R9/R13 同一份数据，各 +21） |
| AMD 显著数 | 30 | 33 | 检验计数口径；新增 SEPTIN8/FABP6/NPPB |
| WetAMD 显著数 | 25 | 26 | 同上 |
| 其余 26 个结局 | — | 完全未变 | — |
| PILRB OR（AMD/湿性/干性） | 1.0643 / 1.0705 / 1.0591 | 1.0558 / 1.0612 / 1.0514 | 用回自己的 βexp |
| ICAM4 OR（T2D） | 0.9578 | 1.1685 | 方向反转 |
| PheWAS 显著 SNP / 阈值 | 228 / 4.38e-9 | 236 / 4.23e-9 | 随显著集合变化 |
新增的 3 个 AMD 关联是边缘显著，写作时必须标注
SEPTIN8 / FABP6 / NPPB 的 FDR 全是 0.049（旧口径下为 0.054），是检验计数口径变化导致的 临界翻转，不是新证据。不能与 TNFRSF10A 那类强信号并列陈述。
没变什么——主要结论全部稳住
- 候选靶点清单不变：仍是 8 个 A 级（PILRA / PILRB / TGFB1 / WARS / CSF2 / CASP10 / TNFRSF10A / SPRY2）+ APOE（B 级高多效）
- 共定位判定不变：31 个蛋白参与，9 个获支持
- TNFRSF10A 仍是唯一跨队列共定位复制的靶点
- p 值逐条不变：616 条可比对关联最大相对差 9×10⁻¹³（浮点误差）， 实证了"Wald 比的 z 与 βexp 无关"这一推断
⇒ 这轮修复改的是效应量数值与外部 T2D 的功效，没有动摇任何一条主要结论。
八、出图脚本审阅（2026-07-20 补做）
图是投稿的门面，且画错了比数据错更难发现——图能生成、数字也对，错的是"图在声称什么"。 本轮逐个审了 02_visualize.R 与 12–15。
查出并已修的问题
① h² 图注与图的内容自相矛盾（14_figure_downstream.R）
图注写着 "Observed-scale estimates depend on sample case ratio and are not comparable across endpoints"， 但该图已按 Neff 重新生成、恰恰是可比的。这行字会原样印进论文。已改为说明 Neff 口径并注明误差棒含义。
② h² 与 rg 全程没有标准误
ldsc() 返回的 res$V 是 S 下三角拉直后的抽样协方差矩阵，里面有不确定度， 但脚本只取了 diag(S) 的点估计。补出来之后解读直接改变：
| 估计值 | SE | z | |
|---|---|---|---|
| AMD–DryAMD rg | 0.972 | 0.015 | 63.5 |
| Retinopathy–Maculopathy rg | 0.920 | 0.124 | 7.4 |
| Maculopathy–T1D_wide rg | 0.940 | 0.151 | 6.2 |
| T2D h² | 0.162 | 0.008 | 19.4 |
| Maculopathy h² | 0.074 | 0.029 | 2.5 |
| Neuropathy h² | 0.065 | 0.029 | 2.3 |
rg 的 SE 跨度 0.015–0.281，h² 的 z 从 2.3 到 19.4。此前文档反复引用的 "糖网↔黄斑 rg 极高"，置信区间其实很宽；Maculopathy 与 Neuropathy 的 h² 勉强离零。不带 SE 报这些数字会误导审稿人。
已修：04_ldsc.R 现输出 h2_observed.csv（含 se/z/95%CI）与新的 genetic_correlation_rg_long.csv（每对 rg/se/z/p，rg 的 SE 用 delta 法： rg = S_ij/√(S_ii·S_jj)，对 (S_ij, S_ii, S_jj) 求梯度后 JᵀΣJ）；h² 图加 95% CI 误差棒。
③ 图上排版：数值标签被误差棒横线穿过、图注右侧被裁掉——均已修。
审阅通过的部分
13_figure_locus.R的-log10(p)用-(pnorm(-|z|, log.p=TRUE) + log(2))/log(10)在对数域计算，避免下溢（做法正确）；facet_grid(trait ~ ., scales="free_y")共享 x 轴、各轨道独立 y 轴，符合 locus 图规范12_figure_screen.R的显著性口径与01_screen.R一致（FDR<0.05）；热图色标按 |log2(OR)| 的 98 分位截断且在图注里写明，没有偷偷压缩极端值02_visualize.R的标题明确标注 "MR only (cis-pQTL, FDR<0.05)"， 如实说明未经共定位过滤、仅供探索，没有夸大15_figure_singlecell.py用 Mann-Whitney U + AUC 效应量 + 全部"基因×细胞类型" 一起 BH 校正；脚本自己就写明"26 万细胞下 P 值必然极小，真正有解释力的是 AUC"
单细胞检验的一个方法学限定（需写进论文）
同一供体的细胞不是独立观测（pseudo-replication）。把 26 万个细胞当独立样本做 Mann-Whitney，p 值会系统性偏松。脚本靠"以 AUC 为准、不看 p"规避了这个问题， 这是正确的应对，但论文里必须写明，否则那些极小的 p 值会被当成强证据。
九、run_all.R 引擎路径审阅（2026-07-20 补做）
前八节覆盖的是 analysis/01→16 批量主线——本课题的全部结果都由它产出。 本节补审此前未查的 6 个文件，它们属于 run_all.R 单暴露引擎，本课题一次都没走过， 故以下问题不影响现有任何结果；但换课题若要走这条线，必须先修。
查出并已修
| 文件 | 问题 | 严重度 |
|---|---|---|
R/12_bidirectional.R + run_all.R | 反向 MR 结构性失效 | 高 |
R/04_iv_select.R | 弱工具过滤"失败即放行"；overall F 无守卫 | 中 |
R/07_sensitivity.R | E-value 把 OR 当 RR；MR-PRESSO 分布数写死 | 中 |
① 反向 MR 拿到的是被过滤过的结局数据
run_all.R 的 Phase 2 为省内存，按「暴露工具 SNP 并集」预过滤读取结局； 而 run_bidirectional() 又从这份数据里选「结局的反向工具」——等于只能在暴露自己的 cis-pQTL 里找结局的全基因组显著位点。结果要么是 0 个（静默返回 insufficient_data）， 要么更糟：碰巧通过几个，跑出一个被标成 status = "ok" 的无意义反向 MR。
已修：新增 out_std_full 参数，调用方必须显式提供未过滤的结局数据；拿不到就 明确拒绝出结果而不是给一个错误答案。run_all.R 调用点改为单独全量读一次。
② 弱工具过滤失败即放行
d[is.na(f_stat) | f_stat >= f_stat_min] —— F 算不出来的 SNP 会被保留。 弱工具偏倚正是要靠这一步挡住。已改为算不出即剔除并告警。
同时给 overall F 加守卫：逐 SNP R² 直接相加，若 Σ R² ≥ 1 则 (1 − ΣR²) 归零或变负， 公式会算出荒谬值甚至负数。现改为返回 NA 并告警"工具间可能非独立"。
③ E-value 把 OR 当 RR
E = RR + √(RR(RR−1)) 的原始定义针对风险比。结局不罕见时 OR 明显高估 RR， 直接代入会高估 E-value——把证据说得比实际更稳健。本课题结局患病率差异极大 （AMD 约 3%、糖网约 20%），不能一概而论。已改为 evalue_or(est, rare)： 罕见结局用 OR≈RR，常见结局用 √OR 近似 RR；未声明时告警。
MR-PRESSO 的 NbDistribution 原写死 1000，使全局检验 p 值下限锁死在 0.001。 已改为可配置（sensitivity$presso_nb），并在 p 触及分辨率下限时告警。
审阅通过
R/06_mr_core.R：按 IV 数自适应选方法（1→Wald、2→IVW 固定、≥3→IVW 随机+Egger+ 加权中位数/众数）、OR/CI 提取、主分析行标注均正确R/08_coloc_abf.R：setkey排序后match()对齐两侧 SNP ✅；maf = pmin(eaf, 1−eaf)使等位编码翻转不影响 ✅；s = ncase/(ncase+ncontrol)符合 coloc 的 case 比例定义 ✅R/11_report.R：无问题
coloc.abf 的一个固有性质，解读时别搞错
coloc.abf 基于 Z² 计算近似贝叶斯因子，不看效应方向——两个性状效应方向相反时， 它照样会判定共定位。这不是 bug，但不能把"共定位成立"读成"同向因果"。
smoke test 暴露的两个环境级阻塞（已修）
跑 tests/ 下的 smoke test 时，4 个里有 2 个直接报错：
| 阻塞 | 后果 |
|---|---|
config.yaml 的 ld.bfile 指向 reference/ld/EUR，本机不存在 | IV 选择阶段直接崩 |
mr.raps 包未安装，而 pick_methods() 在 IV≥3 时会带上它 | run_mr() 硬崩，主分析跑不完 |
已修：LD 面板路径改为实际位置 /mnt/d/mrdata/ldpanel/EUR、plink 改为 plink1.9； MR-RAPS 改为可用才启用，缺失时降级并告警（少一个稳健性方法，不该让主分析跑不动）。
修完 4 个 smoke test 全部通过：m1 读取标准化 ✅、m2 单 IV Wald ✅、 m2b 多 IV（IVW/Egger 截距 P=0.727/MR-PRESSO 全局 P=0.014 检出 1 个离群/LOO/方向性）✅、 position 位置匹配 ✅。
⇒ 这条引擎路径此前是完全跑不动的，这也解释了为什么本课题从未用它出过结果。
十、手工复算——脱离流水线的独立验证 ✅
前九节都是代码审阅：读代码找问题。但审阅有个根本局限——看不出"我以为它这么算、它其实那么算"。 唯一能真正验证整条链的办法，是完全不调用流水线任何一行代码，从原始文件手工算一遍再对。
脚本已归档：tools/manual_verify_tnfrsf10a.py、tools/manual_verify_flip.py （只用 Python 标准库，连 scipy 都不依赖，p 值用 math.erfc 直接算）。
复算 1：TNFRSF10A（主打靶点）
手工完成：解析 UKB-PPP 复合变异 ID 8:23082971:G:T:imp:v1 → A0=G(other) / A1=T(effect)； 从 FinnGen R13 AMD 原始 .gz 里 zgrep 出 rs13278062 那一行；按 alt=effect 约定取 β； 手工谐化；手工算 Wald 比与 R²/F。
| 量 | 手算 | 流水线 | 相对差 |
|---|---|---|---|
| beta.exposure | −0.462 | −0.462 | 0 |
| beta.outcome | +0.0689603 | +0.0689603 | 0 |
| b (Wald) | −0.1492647186 | 同 | 2.6×10⁻¹⁵ |
| se | 0.02815454545 | 同 | 1.7×10⁻¹⁵ |
| OR | 0.8613410717 | 同 | 3.9×10⁻¹⁶ |
| 95%CI | [0.8150974, 0.9102083] | 同 | ~1×10⁻¹⁶ |
| p | 1.147791753×10⁻⁷ | 同 | 5.4×10⁻¹⁵ |
| R² / F | 0.11194202 / 4355.75 | 同 | ~1×10⁻¹⁵ |
10 个量全部吻合到浮点精度。
复算 2：ABO / rs505922（补验等位翻转分支）
复算 1 有个缺口：TNFRSF10A 的等位在两侧恰好一致（G/T vs G/T），没走到翻转分支—— 而 ICAM4 那个 bug 说明方向错误正是最危险的一类。故另找一个两侧编码相反的位点：
- 暴露
A0=T, A1(effect)=C，β = +1.1710 - 结局
ref=C, alt(effect)=T，β = −0.048420 ← 效应等位与暴露相反 - 手工翻转后 β = +0.048420 → b=+0.04134910、se=0.01115636、OR=1.042216 [1.019674, 1.065256]、p=2.10×10⁻⁴
与流水线五个量完全吻合；且流水线记录的 effect_allele.outcome = C，已正确统一到暴露侧的效应等位 ✅。
这次验证覆盖了什么
复合变异 ID 解析（A0/A1 → other/effect）· UKB-PPP 的 BETA 按 A1 定向 · FinnGen 的 alt=effect 约定 · 谐化的两个分支（等位一致 + 等位翻转） · Wald 比 · OR/CI · p 值 · R²/F。
⇒ MR 主体计算链是今天唯一一件被"证明对了"、而不只是"没发现错"的事。
但这次验证没有覆盖什么
共定位（hyprcoloc/coloc.abf）、LDSC、PheWAS、单细胞富集，只做了代码审阅， 没有独立复算。 手工复算 hyprcoloc 的簇后验概率需要重实现整个算法，成本过高。 这几块的可信度目前只建立在"代码读下来是对的"之上，弱于 MR 主体。
十一、遗留项
config.yaml里的reference/gene_annotation.tsv、reference/ld/EUR、reference/bin/plink在本机均不存在。这几处只被run_all.R单暴露引擎使用，批量筛查主线不碰， 故不影响现有结果；但换课题走run_all.R那条路时会报错。- R9 线的
MICB_MICA因 UKB-PPP 复合 assay 命名 + tar 未下载而跳过（只影响 R9）。
⚠️ 断点续跑用「文件存在」代替「文件完整」——本轮实测踩到
04_ldsc.R 靠 !file.exists(*.sumstats.gz) 判断是否需要重新 munge。本轮重跑被中断两次， 分别留下 2 个和 3 个截断的 .sumstats.gz——文件大小都在正常的 14.8 MB 量级， 肉眼完全看不出问题。若不清理直接续跑，LDSC 会拿残缺数据算出一个看着合理的 rg 矩阵， 不报任何错。
这与本次审计修掉的三个 bug 是同一个模式：用"存在"代替"完整"，用"没报错"代替"跑对了"。 同样写法的还有 01_screen.R 的 *_all.csv 缓存与 05_hyprcoloc.R 的 _parts/。
建议（本轮未改，避免引入新变量）：统一改成「先写临时文件、成功后再 rename」， 这样中断只会留下临时文件，不会污染断点判断。本轮的处置是把 R9/R13 的 sumstats 全部清空重算， 并在跑完后逐个核对行数（各约 1,195,667 行 ≈ HapMap3 全量），确认完整后才采信结果。
（Synapse 下载脚本 tools/ 已按这个模式写：先落 .part，下完才 Move-Item 改名， 中断只会留下 .part，不会留下"看着完整"的半截 tar。）
十二、明确尚未修复的问题
审计不是为了宣布"没问题了"。以下问题已知、已定位、但本轮有意未改，记录在此以免被遗忘：
| # | 问题 | 为什么没改 | 风险 |
|---|---|---|---|
| 1 | 已修（2026-07-21），见下方 | — | |
| 2 | 6 个 smoke test 里 5 个零断言，PASSED 是脚本跑到最后无条件打印的 | 需要为每个测试补预期值 | 它只证明"能跑完不报错"，不证明"算得对"——本文档第九节引用的"4 个 smoke test 全过"，证据力仅限于此 |
| 3 | R/00_setup.R 的全局选项（随机种子、datatable.na.strings）只在 01_screen.R 生效 | 给 05/06/08/10 补上会改变 fread 的缺失值解析，可能连带改动刚验证完的结果 | fread 在不同脚本里对 "."/""/"#NA" 的处理不一致；05/08/10 无显式种子（实测 hyprcoloc 确定性，重跑逐位一致） |
| 4 | Rmd/report_template.Rmd（191 行）未审阅 | 只在 run_all.R 出报告时用，不在结果路径 | 低 |
实测：断点缓存确实被静默复用了——但这次结果没被污染
问题 1 不是理论风险，本轮真的发生了：清 _parts 的命令写在一个被中断的编排脚本里， 从未执行到；后续的可续跑脚本又按「文件存在＝已完成」跳过了它们。结果是 R13 有 31 个、R9 有 23 个 _parts 仍是 07-18 生成的，早于 07-20 的筛查修复。
发现方式：核对产物时间戳——这也是本轮另一个教训（rerun_pool.sh 漏跑阶段 08， 同样是靠时间戳发现的）。判断"跑完了"必须核时间戳，不能凭印象。
三层核实，逐层加硬：
| 核实 | 结果 |
|---|---|
| 缓存里记录的性状数 vs 按当前显著性应有的性状数 | 0 个不符 |
| 区域缓存的 beta/se 是否被解析成字符型 | 161 个缓存全部数值型 ✅ |
| 删掉两个旧缓存（CSF2、APOE）重算，与备份逐字节比对 | 完全一致 ✅ |
⇒ 当前共定位结果可信——这是重算证实的，不是推断。
但这次没出事是运气，不是设计
缓存之所以仍然有效，是因为本轮筛查修复的净效果是**「新增 46 条、消失 0 条」**—— 没有任何蛋白丢掉显著疾病，所以旧的性状组合恰好还成立。
只要哪次重跑让某个蛋白失去一个显著疾病，这批缓存就会静默给出错误的共定位结论， 而且不会有任何报错。
✅ 已修（2026-07-21）：原子写入 + 输入指纹
新增共用缓存层 analysis/_cache_io.R：
- 原子写入：先写同目录
.tmp再rename。中断只会留下.tmp， 绝不会留下"看着完整"的半成品。 - 输入指纹：每份缓存旁写一个
.fp文件记录其输入，输入变了自动失效。 存的是明文 key 而非哈希，出问题时能直接看出是哪一项变了，例如：v1|TNFRSF10A|AMD,DryAMD,WetAMD|5e+05|0.5|0.5
指纹涵盖范围：
| 缓存 | 指纹包含 |
|---|---|
01_screen.R 的 *_all.csv | 结局文件身份（名/大小/mtime）· 病例对照数 · 匹配模式 · 效应量与 P 值编码 · se_from_ci · 暴露文件身份与样本量 |
05_hyprcoloc.R 的 _parts/ | 该蛋白的性状组合（蛋白 + 排序后的疾病）· 窗口 · hyprcoloc 两个阈值 |
验证方式（故意制造失效）：把 TNFRSF10A 的指纹里的疾病集合从 AMD,DryAMD,WetAMD 改成 AMD,DryAMD，再跑 05——
| 检查 | 结果 |
|---|---|
| 是否触发重算 | ✅ 打印「输入已变，重算: TNFRSF10A.csv」，其余 54 个正确跳过 |
| 重算结果与原文件比对 | ✅ 逐字节完全一致（CSF2 同） |
| 指纹是否刷新回正确值 | ✅ |
是否残留 .tmp | ✅ 无 |
⇒ 新机制不改变任何现有数字，只是让下次重跑不再有静默复用的风险。 升级前的既有缓存由 tools/stamp_cache_fingerprints.R 一次性补写指纹， 其合法性依据正是上面那三层核实。
本轮已修但值得单独记一笔的引擎路径问题（run_all.R 侧，不影响现有结果）： R/09_hyprcoloc.R 原用 grepl 做性状名匹配——子串匹配会让 AMD 命中只含 WetAMD/DryAMD 的簇而误判共定位，已改精确匹配；同文件的区域谐化原先 完全没有回文 SNP 防护（A/T、C/G 无法凭等位编码区分正负链，链翻转会被误判成 "等位一致"导致 β 符号错误），已改为默认剔除并告警。
十三、共定位与 LDSC 的独立复算（2026-07-25 补做）
第十节的手工复算只覆盖了 MR 主体。本轮把此前"只审代码未独立复算"的共定位与 LDSC 也补上，脚本在 tools/verify_coloc_pqtl.R、tools/verify_ldsc.R。
13.1 共定位：跨方法交叉验证 ✅
流水线用 hyprcoloc（多性状）。独立用 coloc.abf（不同方法、成对）对头号靶点 TNFRSF10A 的 pQTL × FinnGen AMD 重跑：
| 量 | 独立 coloc.abf | 流水线 hyprcoloc |
|---|---|---|
| 共享 SNP 数 | 5612 | 5609 |
| 共定位后验 | PP.H4 = 0.9994 | best_pp = 0.9219 |
两种方法都强力支持共定位（SNP 数几乎一致），跨方法交叉验证通过——hyprcoloc 的"支持"判定站得住。
13.2 LDSC：发现 rg 的标准误被低估 ⚠️（点估计无误）
独立从 ldsc_result.rds 的 S/V 矩阵重算 AMD 亚型间遗传相关：
| 项 | 结论 |
|---|---|
| rg 点估计 | ✅ 与 GenomicSEM .log 逐位吻合（AMD-WetAMD 0.931、AMD-DryAMD 0.972、WetAMD-DryAMD 0.865） |
| h2 及其 se | ✅ 与 .log 逐位吻合（AMD h2=0.1267, se=0.0225） |
| rg 的 se / z / p | ⚠️ 被低估 |
analysis/04_ldsc.R 用 delta 法从 V 矩阵算 rg 的 se（第 110–124 行）。我独立重算精确复现了它（AMD-WetAMD：独立 0.0237 vs CSV 0.024）——代码实现无误。但它与 GenomicSEM .log 自报的 rg se（0.1594）差 6.6 倍：
- 原因：rg 逼近 1 时 delta 线性化会低估 se；GenomicSEM 的 log se（直接分块 jackknife）更可靠。
- 后果：
genetic_correlation_rg_long.csv里 AMD-WetAMD 的 z=39.22、p=0 是虚高（用 0.1594 则 z≈5.84、p≈5×10⁻⁹）。 - 对结论无影响（rg 点估计正确，AMD 亚型遗传上近乎相同这一结论稳固；即便 z≈5.8 也极显著），但报告的精度被夸大。
✅ 已修复（2026-07-25）
analysis/04_ldsc.R 已改为优先解析 GenomicSEM .log 的原生 jackknife se（缺失才回退 delta），并保留 se_delta 列作透明对照；两个家族的 genetic_correlation_rg_long.csv 已用原生 se 重出（tools/fix_ldsc_rg_se.R）。修正后：AMD↔WetAMD 由 z=39.22/p=0 → z=5.84/p=5.2×10⁻⁹；RetinopathyStrict↔T2D → z=1.53/p=0.13（正确地不显著）。此问题只显著影响**高 rg（近 1）**的配对；近零 rg（如 T2D↔糖网 rg≈0.05）delta 与 jackknife 本就一致（0.044 vs 0.043）。rg 点估计、h2、h2-se 及 REPORT 的 rg 矩阵表均未受影响（它们只含点估计）。