Skip to content

R9 第一轮定向排查(2026-07-30 晚)

承接 R9 方法学核查与共定位补齐

方法:拿当天查出的 12 类缺陷当检查表,去其余脚本里找同类,不做逐行通读。

为什么不通读

本轮全部缺陷(12 + 新增 5)无一是靠通读代码发现的——全部来自 「换数据实跑」「结果对不上时追根因」「拿已知缺陷模式定向排查」三条路径。

一、排查结果总览

类别结论
03_phewas.R 缓存无指纹确认缺陷,已修并重查
03_phewas.R glob 扫进退役源确认缺陷,已修
区域读取器三份副本确认,只修真 bug,未合并(理由见下)
OpenGWAS 别名 rsID新查出,已修
01_screen.R 缓存指纹✅ 合格,无需改
04_ldsc.R 指纹机制✅ 合格
smoke tests 断言✅ 有实质断言含阳性对照
08/10/18 在 R9 副本✅ 从未产出,属"跑不起来"非"跑错了"

二、① PheWAS 缓存无指纹(唯一影响已发表数字)

r
if (!FORCE && file.exists(CACHE))   # 原:纯文件存在判断

证据链

  • 溯源 phewas_provenance.csvreused_cache=TRUE → 2026-07-29 19:16 那次没真查 OpenGWAS
  • 缓存含 7 个已不再显著的 SNP → 证明它是换 T1D 源之前的 SNP 集合建的
  • 当前 94 个显著 SNP 中 39 个无 PheWAS 数据;而 ieugwasr::phewas 只返回达阈关联, 缺席无法区分「查过无关联」与「根本没查」

修复:指纹取「排序后 SNP 集合的 md5 + 查询阈值」。上线即报:

⚠ 缓存存在但指纹不一致(显著 SNP 集合或阈值已变),重新查询

重查结果:关联数 10,976 → 16,141(+47%),Bonferroni 显著 9,359 条。 此前有近半脱靶关联从未被查询过。

三、② glob 扫进退役源

Sys.glob("results/screen_*/*_significant.csv") 会连 screen_R9_retired (已退役的错误 T1D 源 GCST90475661,实为 T2D 主导)和 screen_R9_sens(敏感性替代源) 一起扫入。这两者是对照用的,其显著 SNP 不该进入主线脱靶分析。

此前未发作,仅因 03 跑在这两个目录创建之前。修复后日志实证:

已排除对照用目录 3 个文件: R9_retired, R9_sens

四、④ 新缺陷:OpenGWAS 用合并/别名 rsID 返回

94 个查询 SNP 中 9 个以别名 rsID 返回(如 rs386699994rs1234415383)。 按 rsid 关联蛋白时匹配不上,protein 变空,52 条关联被静默排除在所有蛋白的 脱靶负担之外,全程无任何提示。

9 个别名按 chr:position 全部唯一命中: APOBR、ITGB7、MOG、CDSN、FGF21、CFB、SIGLEC5、MLN、TREML2。 A 级候选 ERMAP / IFNAR1 / APOL1 不在其中。

修复经过三版才对,教训值得记

  1. 第一版anyNA() 判空 → fread 默认 na.strings="NA",CSV 里空的 protein 字段被读成 "" 而非 NA回填从未触发。 → 判空必须同时认 NA 与空串。
  2. 第二版只补 protein 不归一 rsidphewas_summary(rsid, protein) 聚合, 同一变异以别名和原名各占一行,该蛋白多效负担被重复计数,并把 14_figure_downstream.Rfactor(levels = protein) 撞成 factor level duplicated 报错。
  3. 第三版才正确:按坐标把别名 rsID 改写回规范 rsID 令同一变异合并, 再补归属,再按 (rsid, id) 去重。

★第 2 版的错误若不是 14 崩了根本发现不了——重复计数在数字上看不出异常。 崩溃是好事,静默的错误才可怕。

五、③ 为什么没有合并三份区域读取器

副本chrX 修复空值归 NULLgene 名兜底缓存
_region_io.R
08_replicate_iamdgc.R:39✅(本次)
10_eqtl_coloc.R:54✅(本次)

08/10 在 R9 副本里跑不起来(依赖 R13 的 target_list/verdict,本仓库无 screen_R13_*/network_R13_*),改动无法在此验证,而它们正是 R13 主线的关键路径

把未经测试的重构推上关键路径,风险大于收益——当天已有一次教训 (resolve_tar 只改一处导致 MICB_MICA 被跳过)。 故本次只修真正的数据缺陷(染色体编码),并留 TODO(R13 重跑时) 合并标记。

六、验证干净的部分

  • 01_screen.R 缓存指纹合格:涵盖结局文件名/大小/mtime、病例对照数、匹配模式、 效应与 P 值编码、se_from_ci、暴露文件身份与样本量。换源能正确失效—— 这解释了为何 T1D_gcst_all.csv 单独重算而其他结局未动。
  • 04_ldsc.R 有 fp 指纹机制(含 file.size + mtime)。
  • smoke tests 有实质断言m1 三条命名断言,m3 含阳性对照 (合成同因果变异区域应 PP.H4 >= 0.8)。记忆中"零断言"问题确已修复。

七、操作层面的陷阱(我自己也踩了)

; 串联会掩盖失败

本轮下游链条用 ; 串联,14_figure_downstream.R 崩溃后仍照常打印 CHAIN_DONE判定成败必须看 Execution halted 计数,不能看结束标记—— 与 $? 被外层吞掉给出假 0 是同一类陷阱。

八、三层对照:修复大幅改善数据完整性,但结论未改变

第 1 层:PheWAS 关联

指标修复前修复后
关联条数10,97615,789(+43.9%)
Bonferroni 显著6,4079,107
覆盖唯一 rsid6292 / 94
蛋白归属缺失18 条0 条
reused_cacheTRUE(假复用)FALSE(真查询)

中间过程:完整查询得 16,141 条,别名归一后去重删 352 条(同变异同数据集), 16,141 − 352 = 15,789,数字自洽。

第 2 层:各蛋白多效负担

  • phewas_summary 行数 62 → 92
  • 新增 37 个蛋白(此前从未被查询):APOBR、APOH、CD28、CD40LG、FGF21、IL10、 KIT、KLK1、MLN、OSM、VEGFB 等
  • 退出 4 项:空串(旧版无归属的孤儿行)+ CDC27、DNAJB6、GRP —— 这三个主分析里根本不显著,只在退役源里显著,其退出正是修复 ② 生效的直接证据
  • 55 个共有蛋白中仅 3 个负担有变:CFB +3、CDSN +3、SIGLEC5 +2

第 3 层:靶点分层(决定结论是否改变)

修复前修复后
分层 A / B / C / D3 / 4 / 11 / 13完全相同
A 级候选ERMAP、IFNAR1、APOL1无变化
phewas_flag 变动0 个

为什么没变:新增覆盖的 37 个蛋白绝大多数不在 31 个靶点候选表内;负担有变的 3 个中 CFB 已是 C 级(MHC 区)且本就标注"★高多效"、CDSN 压根不在靶点表、 SIGLEC5 是 D 级且 16→18 未跨过分级阈值。

如何理解"修了但结论没变"

这不意味着修复可有可无。修复之前,我们无从判断结论是否稳健;修完才敢说稳健。 而且旧数据的不完整恰好不在决策路径上,是运气而非设计使然—— 若不修,下次任何一次重跑(显著集合一变)都会把错误放大。

PPT 中引用的脱靶数字经核对全部未变(APOE 642、NOTCH2 20、APOL1 2、ERMAP 6、 GALNT3 2、IFNAR1 4、WARS 24),故汇报材料的多效负担相关内容无需改动。

九、尚未排查

  • 第二轮:数字可追溯性(PPT / docs 里每个数字反查到产生它的文件与代码行)
  • 第三轮:关键结论脱离流水线独立复算
  • 反向 MR / Steiger 实现未逐行看
  • 02 / 19 / 20 / 21 / 23 / 29 / 31 / 33 未逐个查

个人科研与运维文档 · 内容持续修订