主题
02 · 第 0 步 · 数据验收(输入断言 + LDSC)
执行日期:2026-08-07 状态:✅ 完成(A1/A3 通过;LDSC 与旧结果逐字节相同;A2/A4 因依赖关系移至第 1/2 步之后)
本页上半部分是在运行之前写的
按本轮纪律,「数据 / 版本 / 方法 / 预期结果」必须在跑之前落文,跑完只填「实际结果 / 解读 / 去向」。 ★ 跑前写下的预期与实际不符,一律先按缺陷处理,不许现场解释掉。 理由:事后补填的"预期"没有约束力——任何差异都能事后编出合理解释, 而这类错误在输出上和"没出错"长得一模一样。
结局命名
本页表格用内部 ID(T1D_gcst / T2D_mahajan 等)——它们是文件名与代码里的字符串,实录必须能对得上。论文显示名分别是 Type 1 diabetes / Type 2 diabetes,对照表见总览页。★ 内部 ID 绝不允许出现在图上或论文表格里。
一、这一步做什么、为什么做
v4 把数据验收定为可选的非标准步骤,触发条件只有两个: ① 换了结局数据源 ② 需要估计样本重叠。
本轮两个条件都不满足(数据源与上一轮相同)。
★ 本轮仍跑 LDSC,是用户指定的可复现性验证,不是流程触发。 记录与写作时必须如实这样描述,不能说成"流程要求"。
它验收的是数据源,不是候选蛋白——这一步不淘汰任何蛋白。
二、用什么数据、什么版本
2.1 结局(LDSC 纳入 6 个,读 outcome_manifest.csv 的 release == "R9_dm")
| id | 来源 | 版本 / 编号 | build |
|---|---|---|---|
| Retinopathy | FinnGen | R9 | GRCh38 |
| Maculopathy | FinnGen | R9 | GRCh38 |
| Nephropathy | FinnGen | R9 | GRCh38 |
| Neuropathy | FinnGen | R9 | GRCh38 |
| T1D_gcst | GWAS Catalog(harmonised GWAS-SSF v1.0) | GCST90824163 | GRCh38 |
| T2D_mahajan | DIAMANTE noUKBB | Mahajan 2018b European | GRCh37(LDSC 侧按 rsID 处理,无需 liftover) |
NeovascGlaucoma 不入 LDSC(病例 1,100 太小),与 R13_dm 口径一致。
2.2 参考数据
| 用途 | 路径 | 说明 |
|---|---|---|
| LD scores + HapMap3 位点表 | /mnt/d/mrdata/ldsc/eur_w_ld_chr/ | 实测就位:48 个文件;w_hm3.snplist 17 M |
| 等位频率第三方裁决 | /mnt/d/mrdata/ldpanel/EUR(1000G EUR,GRCh37,按 rsID 匹配) | 用于 EAF 方向断言 |
2.3 代码与环境
- 运行目录
/home/research/mr-pipeline-r9v3(属主research),以research身份运行 - git 基线:
F:\project\mr-pipeline-r9v3提交9a0edc3 config.yaml:resume: false(从零重算)、缓存指纹v3
三、什么检验方法
3.1 输入断言(A 组,无需参考面板以外的数据)
| 脚本 | 查什么 | 为什么这层不可省 |
|---|---|---|
tools/verify_input_sanity.R | 列映射审计(哪一列被当成 effect_allele / other_allele / eaf);Beta 是否等于 log(OR)(防重复取对数);全文件按步长抽样,不只取开头 | ★ 输入整体错位时内部一致性是完美的:方向错误不改变 p 值、不改变 CI 宽度、不改变任何一张图。T1D 源 EAF 整列反转就是这样躲过全部审计的 |
tools/verify_t1d_eaf.R | 用全部工具位点对 1000G EUR 做第三方裁决,判 GCST90824163 的 EAF 列是否整列反转 | 只用非回文位点(排除 A/T、C/G)——等位身份靠字母即可确定,无链歧义 |
tools/verify_eaf_all_outcomes.R | 把上面的裁决扩展到每一个在用结局 | 此前其余结局只有"暴露 eaf ↔ 结局 eaf 相关"这种内部一致性证据,不是外部参照 |
tools/verify_harmonise.R | ① 谐化后两侧效应等位逐条相等;② Wald 比 b == beta.outcome/beta.exposure;③ 回文位点不得残留模糊回文(eaf 0.42–0.58) | 谐化做错 = 因果方向整体反过来,而 p 值/CI/图形全都正常,统计量上完全看不出 |
⚠️ verify_harmonise.R 需要 01_screen.R 的产物,因此排在第 1/2 步之后补跑,本步只跑前三个。
3.2 LDSC(B 组)
- 工具:GenomicSEM 的
munge()+ldsc() - 参考面板:
eur_w_ld_chr(欧裔 LD scores),位点限定 HapMap3(w_hm3.snplist) - 输入按格式分派:
finngen/gcst_ssf(T1D,beta+SE 直给) /mahajan(chr:pos 无 rsID, 用ldpanel/EUR.bimGRCh37 注回 rsID —— 实测 1,178,756 / 1,217,312 = 96.8% 落在 HapMap3) - 输出:观测尺度 h²、遗传相关矩阵 rg
四、★ 预期结果(跑之前写下)
4.1 断言
| 项 | 预期 |
|---|---|
verify_input_sanity.R | 全部通过,非零码退出即缺陷 |
verify_t1d_eaf.R | GCST90824163 判为方向正确(08-05 已修;修复前是"整列反转") |
verify_eaf_all_outcomes.R | 每个结局均判方向正确 |
4.2 LDSC —— 应与旧结果数值一致(同数据、同面板、同代码)
遗传相关 rg(旧值,15 对):
| trait1 | trait2 | rg | se | p |
|---|---|---|---|---|
| Retinopathy | Nephropathy | 0.966 | 0.085 | 6.2e-30 |
| Retinopathy | Maculopathy | 0.946 | 0.102 | 1.3e-20 |
| Maculopathy | Nephropathy | 0.942 | 0.104 | 1.5e-19 |
| Neuropathy | Nephropathy | 0.751 | 0.111 | 1.5e-11 |
| Retinopathy | T2D_mahajan | 0.744 | 0.048 | 5.6e-54 |
| Retinopathy | Neuropathy | 0.734 | 0.103 | 8.6e-13 |
| Nephropathy | T2D_mahajan | 0.693 | 0.055 | 3.1e-36 |
| Maculopathy | T2D_mahajan | 0.680 | 0.056 | 1.6e-34 |
| Maculopathy | Neuropathy | 0.601 | 0.115 | 1.8e-07 |
| Neuropathy | T2D_mahajan | 0.591 | 0.061 | 4.1e-22 |
| Maculopathy | T1D_gcst | 0.454 | 0.061 | 9.9e-14 |
| Retinopathy | T1D_gcst | 0.453 | 0.051 | 2.7e-19 |
| Nephropathy | T1D_gcst | 0.433 | 0.061 | 1.4e-12 |
| Neuropathy | T1D_gcst | 0.359 | 0.074 | 1.1e-06 |
| T1D_gcst | T2D_mahajan | ★ 0.087 | 0.034 | 0.0108 |
观测尺度 h²(旧值):
| trait | h²_obs | se | z |
|---|---|---|---|
| Retinopathy | 0.1875 | 0.0187 | 10.05 |
| Maculopathy | 0.3052 | 0.0426 | 7.16 |
| Neuropathy | 0.2794 | 0.0478 | 5.85 |
| Nephropathy | 0.2717 | 0.0338 | 8.03 |
| T1D_gcst | 0.1682 | 0.0150 | 11.19 |
| T2D_mahajan | 0.1222 | 0.0069 | 17.72 |
★ 最关键的一个数:rg(T1D_gcst, T2D_mahajan) = 0.087。 它是当初退役整个旧 T1D 数据源的依据——旧源 GCST90475661 给出 rg = 0.886, 说明它实为 T2D 主导。这个数必须复现;若明显偏离,说明数据源身份出了问题, 是红旗级缺陷,不是"数值抖动"。
4.3 允许的容差
| 项 | 容差 | 超出即缺陷 |
|---|---|---|
| rg / h² 点估计 | 小数点后三位完全一致 | 是 —— LDSC 无随机性,同输入同面板应完全可复现 |
| munge 保留位点数 | 完全一致 | 是 |
| sumstats 文件数 | 6 | 是 |
⚠️ 若出现不一致,第一嫌疑不是"算法随机",而是:① 输入文件被动过;② 包版本变了;③ 参考面板被动过。 按此顺序排查,不许先假设是正常波动。
4.4 预期产物路径
/home/research/mr-pipeline-r9v3/results/ldsc_R9_dm/
├ genetic_correlation_rg.csv 6×6 矩阵
├ genetic_correlation_rg_long.csv 15 行长表(含 se / z / p)
├ h2_observed.csv 6 行
├ ldsc_result.rds 原始返回对象
├ *_input.tsv 6 个 munge 输入
├ *.sumstats.gz 6 个(每个约 15 M)
└ *_munge.log五、实际结果
5.1 A1 verify_input_sanity.R —— ✅ 通过(0 项失败,退出码 0)
| 检查 | 结果 |
|---|---|
B1 Salem2019 Beta vs log(OR) | n=83,115,最大差 4.86e-06 OK |
B2 各结局 p 与 -log10p 自洽(全文件步长抽样) | 7 个结局全 OK,中位 |Δlog10p| ≤ 3.8e-03 |
B3 暴露 log10(p) 编码 | 取值 [10.80, 5369.50] 全正 = -log10p,与 pval_encoding: neglog10 一致;对照重算相关 0.9952 |
| D 哨兵位点硬断言 | 7/7 OK |
D 组值得单列 —— 它验的是「文件里装的是不是它声称的表型」:
| 数据源 | 位点 | 判定 |
|---|---|---|
| T1D_gcst | rs2476601 (PTPN22) | 风险等位 A → OR 1.907,p=5.8e-184 OK(T1D 头号非 HLA 位点) |
| T1D_gcst | rs689 (INS) | 风险等位 T → OR 2.081,p=2.8e-305 OK |
| T1D_gcst | rs7903146 (TCF7L2) | 阴性对照 p=0.67 OK —— 真 T1D 不应有 T2D 头号信号 |
| T1D_crouch | rs2476601 / rs689 / rs7903146 | OR 1.722 / 1.790 / 阴性 p=0.075 全 OK |
| T2D_mahajan | rs7903146 (TCF7L2) | 风险等位 T → OR 1.350,p=4.5e-238 OK |
★ 那个阴性对照是这层里最有价值的一项:它排除的正是 2026-07-29 踩过的坑 (旧 T1D 源 GCST90475661 实为 T2D 主导)。
5.2 A3 verify_eaf_all_outcomes.R —— ⚠️ 首跑与预期不符 → 调查 → 判为非缺陷 + 修了检查器
这是「跑前预期与实际不符」纪律的一次实际执行,过程留痕
首跑结果:T1D_gcst 判为 ★整列反转,93.7%,cor = −0.9922,退出码 1。 这正是 2026-08-05 修复之前的原始症状。 按纪律先当缺陷处理,不许现场解释,遂调查。
调查结论:不是缺陷,是我写的跑前预期错了。 证据链(全部读代码实证,非推断):
| 位置 | 事实 |
|---|---|
outcome_manifest.csv | T1D_gcst 的 eaf_orientation = **other**(其余 6 个结局为 effect) |
R/adapters/read_tsv.R:59-69 | 声明为 other 时,在读取时做 1-eaf 校正 |
R/adapters/assert_eaf.R:20 | 注释明写「断言在 eaf_orientation 校正之后执行」 |
verify_eaf_all_outcomes.R | 读的是原始文件 |
→ 原始文件确实反转,这是已知且已声明的事实;修复不是改文件,而是在读取时翻转。 我把两个脚本的检查对象搞混了:一个查原始数据源,一个查流水线产物。 正确的预期应为「T1D 原始文件应判反转,且与 manifest 声明一致」。
但调查中发现了一个真问题(已修):
原实现不看 manifest,因此对 T1D_gcst 永远返回 ★未通过 + 退出码 1。 一个永远失败的断言会训练人忽略红色告警——而 2026-07-29 → 08-05 那个 EAF 反转 bug 能活整整一周,靠的正是没人把告警当真。
改为 manifest 感知,判定更严而非更松:
| 原始文件实测 | manifest 声明 | 判定 |
|---|---|---|
| 正常 | effect | OK |
| 整列反转 | other | OK(反转 · 已在 manifest 声明) |
| 整列反转 | effect | ★ 未声明的整列反转 ← 原始 bug 的形状 |
| 正常 | other | ★★ 过度校正:文件正常却声明 other ← 原实现完全看不见的第二种失效 |
★ 第四种是新增的:给一个本来正常的文件误设 other,流水线会把好数据翻错, 原实现对此毫无察觉。
修改后重跑 —— ✅ 通过,退出码 0:
| id | n_cmp | 反转率 | cor_1000G | 原始文件 | manifest 声明 | 判定 |
|---|---|---|---|---|---|---|
| Retinopathy | 1210 | 2.3% | +0.9832 | 正常 | effect | OK |
| Maculopathy | 1210 | 2.3% | +0.9832 | 正常 | effect | OK |
| NeovascGlaucoma | 1210 | 2.3% | +0.9833 | 正常 | effect | OK |
| Neuropathy | 1210 | 2.3% | +0.9833 | 正常 | effect | OK |
| Nephropathy | 1210 | 2.3% | +0.9832 | 正常 | effect | OK |
| T1D_gcst | 1277 | 93.7% | −0.9922 | 整列反转 | other | OK(已声明) |
⚠️ T2D_mahajan 可比对位点为 0(chr:pos 无 rsID),不足以判定 —— 不是通过也不是失败。
产物:results/platform_check/eaf_all_outcomes.csv
5.3 ⚠️ A2 / A4 无法在本步运行(依赖关系,跑前记录里漏了)
| 脚本 | 读什么 | 结论 |
|---|---|---|
verify_t1d_eaf.R | results/screen_R9_dm/T1D_gcst_all.csv | 需 01_screen.R 产物 |
verify_harmonise.R | results/screen_R9_dm/ | 需 01_screen.R 产物 |
★ 这意味着**「数据验收」在本流水线里并不是一个纯前置步骤**:四个断言里只有 verify_input_sanity.R 真正跑在第 1 步之前;另外三个查的是谐化后的流水线产物, 本质上是第 1/2 步的产物验收。跑前记录里我只标了 verify_harmonise.R 有此依赖, 漏了 verify_t1d_eaf.R。
→ A2 / A4 移到第 1/2 步之后执行,届时在 03 页记录。 ★ 校正是否真的生效,要等 A2 才有答案 —— 本步只证明了「原始文件状态与声明一致」。
5.4 B 组 LDSC —— ✅ 与旧结果逐字节相同(退出码 0)
预处理行数(6 个性状全部纳入,与预期一致):
| trait | 预处理行数 | munge 后保留位点 |
|---|---|---|
| Retinopathy | 18,802,626 | 1,191,573 |
| Maculopathy | 18,802,442 | 1,191,573 |
| Neuropathy | 18,801,371 | 1,191,573 |
| Nephropathy | 18,802,440 | 1,191,573 |
| T1D_gcst | 56,442,328 | 1,211,422 |
| T2D_mahajan | 7,666,479 | 1,176,437 |
munge 耗时 14 分 25 秒。
★ 与旧结果的逐项比对(这是本轮"顺带验准确性"的第一份证据)
| 比对项 | 方法 | 结果 |
|---|---|---|
genetic_correlation_rg.csv | cmp 逐字节 | 完全相同 |
genetic_correlation_rg_long.csv | cmp 逐字节 | 完全相同 |
h2_observed.csv | cmp 逐字节 | 完全相同 |
| munge 日志保留位点数 | cmp | 完全相同 |
6 个 .sumstats.gz | 字节大小 | 6/6 相同 |
→ §4.2 预期表里的 15 个 rg、6 个 h² 全部命中,无一偏离。
★ 最关键那一个数复现了:
rg(T1D_gcst, T2D_mahajan) = 0.087 se=0.034 z=2.55 p=0.0108它是当初退役旧 T1D 源(GCST90475661,给出 rg = 0.886 → 实为 T2D 主导)的直接依据。
产物路径:/home/research/mr-pipeline-r9v3/results/ldsc_R9_dm/ 运行日志:results/logs/ldsc_R9_dm_20260807.log
六、如何解读
6.1 LDSC 结果本身说明什么
- 四个并发症之间遗传相关极高(Retinopathy–Nephropathy 0.966、 Retinopathy–Maculopathy 0.946、Maculopathy–Nephropathy 0.942)—— 它们在遗传层面高度共享,这是后面「一个蛋白横跨多个并发症」现象的背景, ⚠️ 但不能据此说蛋白的跨病效应是真的:高 rg 是性状层的事实, 与单个蛋白的因果无关。v4 明确 LDSC 不在候选链上,它只作背景描述。
- 并发症与 T2D 的 rg(0.59–0.74)明显高于与 T1D(0.36–0.45), 与本课题「多数关联反映糖尿病易感性」的主句方向一致,但同样只是背景,不是证据。
rg(T1D, T2D) = 0.087说明当前 T1D 源确实是 T1D 而非 T2D 主导, 数据源身份核验通过。
6.2 这次比对说明什么(方法学价值)
resume: false 从零重算,与三个月前的旧结果逐字节相同。这一条同时证明了三件事:
- 流水线在这一层是确定性的——没有未受控的随机性
- 环境未漂移——R 包版本、参考面板、输入文件都没被动过
- 旧结果在这一层是可信的——这是本轮"顺带验准确性"要拿的东西
⚠️ 但只能推到这一层:LDSC 用的是全量 GWAS,不经过工具选择与谐化。 主链(第 1/2 步)是否同样可复现,要等 01_screen.R 跑完才知道。
6.3 一条不应被忽略的空白
T2D_mahajan 在 EAF 方向核验里可比对位点为 0(chr:pos 无 rsID), 既非通过也非失败。也就是说:七个结局里有一个的等位频率方向从未被外部参照验证过。 它在 LDSC 里能算(按 rsID 注回后 96.8% 落在 HapMap3),但那是另一条路径。 → 记入待办,并在 Limitations 如实写明。
七、进正文还是补充(暂定)
| 内容 | 暂定去向 | 依据 |
|---|---|---|
| 输入断言 | 都不进 | 属可复现性材料;但若断言曾抓出真错(如 T1D EAF 反转),该错本身要进 Limitations 或 Methods 的数据质控段 |
| rg 矩阵 | 补充材料 | 它是背景描述,不是本研究的因果结论;v4 明确 LDSC 不在候选链上 |
rg(T1D,T2D)=0.087 | 正文方法段一句 | 它是"为何选用这一 T1D 源"的直接依据,属数据源选择的正当性说明 |
| h² 表 | 补充材料 | 同 rg |
⚠️ 暂定,第 7 步统一复核。