主题
14 · 第 6 步 · 描述性注释
日期:2026-08-08 状态:6-1 跑前记录(尚未运行)
★★ 本步的铁律:只加分,不否决
v4 明文:6-1~6-4 只加分不否决。全流程的否决点只有两个 (第 3 步反向 MR、4X 表位/PAV),第 5 步只升降措辞。
第 6 步无权改变任何候选的去留。 → 断言 N2 强制校验 candidate_status.csv 的 retained 集合跑前跑后完全相同。
〇、输入集(用户 2026-08-08 拍板)
只做共定位支持的 11 个蛋白(PP.H4 ≥ 0.8,全部非 MHC)。 MHC 且 PP.H4 < 0.8 的不进入本步,改走「MHC 专用核验」。
四个子步骤 6-1 / 6-2 / 6-3 / 6-4 全部只做这 11 个。
已核对的输入(实测,非推断)
| 蛋白 | 哨兵 SNP | MR 结局数 | 最高 PP.H4 |
|---|---|---|---|
APOE | rs429358 | 3(★ 其中 Reti 在第 3 步被否决) | 0.998 |
NUDT5 | rs10508438 | 1 | 0.977 |
ERMAP | rs11210710 | 1 | 0.967 |
PAM | rs149802978 | 1 | 0.959 |
WARS | rs2273804 | 1 | 0.956 |
GALNT3 | rs2116546 | 1 | 0.948 |
IFNAR1 | rs914142 | 1 | 0.940 |
NOTCH2 | rs2641348 | 2 | 0.897 |
SIGLEC5 | rs1106476 | 2 | 0.887 |
ACRBP | rs7959658 | 1 | 0.867 |
LACTB2 | rs191588099 | 1 | 0.800 |
★ 每个蛋白恰好 1 个工具 SNP,共 11 个唯一 SNP(已断言)。
报告纪律:分母不能消失
采用「只做 11 个」的口径,正文与补充表必须同时写明:
「52 对 / 28 蛋白中,13 对 / 11 蛋白 PP.H4 ≥ 0.8」
未做深挖的 17 个蛋白不等于被排除,须在补充表逐个注明类别: MHC 假设不成立(11)/ H4 主导但未达线(4)/ 提示不同变异(2)。
一、6-1 on-target PheWAS
1.1 文献底本
| # | 来源 | 做法 |
|---|---|---|
| C1R | medRxiv 2026.07.14.26358022 · Methods 第 407–414 行 | GWAS ATLAS,4,756 个性状 / 28 个表型大类,阈值 P < 5×10⁻⁸ |
| C1R | Figure 4 图注 | ★ 方向感知:"associations of risk alleles … triangles up representing positive associations and down representing negative" |
| C1R | Discussion 第 296–298 行 | ★ 自陈局限:"PheWAS based on a cis-pQTL cannot fully recapitulate pharmacological inhibition" |
| 本课题 | analysis/03_phewas.R + 37_drug_target_phewas.R | OpenGWAS,Bonferroni 校正,方向感知获益/风险判读 |
未获取的文献
Schmidt AF, et al. Genetic drug target validation using Mendelian randomisation.Nat Commun 2020;11:3255 · doi 10.1038/s41467-020-16969-0 —— 全文未取到 (PMC 反复返回不相关内容)。故本页不引用其具体表述。 若正文需要方法学权威出处,须先取到原文再写。
1.2 与 C1R 的两处口径差异(须在 Methods 声明)
| C1R | 本课题 | 说明 | |
|---|---|---|---|
| 数据库 | GWAS ATLAS(4,756 性状) | OpenGWAS(数据集数以运行时 provenance 为准) | 二者都是汇总库,覆盖不同 |
| 阈值 | 固定 5×10⁻⁸ | Bonferroni(= 0.05 / 实际检验数) | 见下 |
★ 阈值这次会翻转方向,必须先想清楚
旧仓那次查了 94 个 SNP,Bonferroni 阈值 = 1.06×10⁻⁸(比 5×10⁻⁸ 更严)。 本轮只查 11 个 SNP,检验数降到约 1/9,Bonferroni 阈值会放宽到约 9×10⁻⁸ —— 反而比全基因组显著性 5×10⁻⁸ 还松。
定:p_threshold = min(Bonferroni_11, 5e-8) —— 取两者中更严的, 并在产物中同时保留两列(bonferroni_sig、gws_sig)供读者自行取用。
理由:Bonferroni 分母应与实际做的检验数一致(这是它的定义), 但一条 SNP–性状关联若连全基因组显著都达不到,本身就不可信。两条都要过。
1.3 方法与代码
| 步 | 脚本 | 本轮须改 |
|---|---|---|
| 查询 | analysis/03_phewas.R | 输入集由「所有显著 SNP」改为这 11 个;缓存指纹随之失效重查 |
| 判读 | analysis/37_drug_target_phewas.R | ★ 现读 results/targets/target_decision_table.csv(tier 表,用户已定删除)→ 改读 results/candidate_status.csv 与 coloc/coloc_verdict.csv |
方向感知的含义:PheWAS 给的是「某等位与某性状的关联」。 要读成「干预该靶点会怎样」,必须把等位翻到与 MR 工具同一效应等位, 再结合 implied_drug_action(抑制 / 激动)判断某个脱靶性状是获益还是风险。
coloc_verdict.csv 已有 implied_drug_action 列:这 11 个里 抑制 5 个(APOE NUDT5 WARS SIGLEC5 ACRBP)、 激动/补充 6 个(ERMAP PAM GALNT3 IFNAR1 NOTCH2 LACTB2)。
二、★ 预期(写于运行之前)
2.1 机械可验的强预期
| # | 预期 |
|---|---|
E1 | 查询覆盖 11 个 SNP 全部,含 rs191588099 |
E2 | 产物中每个 SNP 都有「已查询」记录,哪怕 0 条关联 |
E3 | Bonferroni 分母 == 实际 SNP 数 × 运行时数据集数,且写进 phewas_provenance.csv |
E4 | 每条关联恰属疾病 / 定量 / 分子三类之一(分类穷尽且互斥) |
2.2 数量预期(来自旧仓同口径实测)
旧仓在 1.06×10⁻⁸ 下,这 10 个 SNP 的显著关联数:
APOE rs429358 896 ← ★ 极端多效(APOE ε4)
NUDT5 rs10508438 34
PAM rs149802978 34
SIGLEC5 rs1106476 29
WARS rs2273804 27
NOTCH2 rs2641348 25
ACRBP rs7959658 19
ERMAP rs11210710 6
IFNAR1 rs914142 4
GALNT3 rs2116546 2
LACTB2 rs191588099 ★ 旧仓从未查过| # | 预期 | 理由 |
|---|---|---|
E5 | 本轮阈值 5×10⁻⁸ 比旧的 1.06×10⁻⁸ 更松 → 各 SNP 关联数应 ≥ 上表 | 单调性;若变少说明 OpenGWAS 库有变,须记录 |
E6 | APOE 仍远超其余 10 个,量级差 1–2 个数量级 | rs429358 是已知极端多效位点 |
E7 | ★ 疾病类脱靶:ERMAP 与 IFNAR1 接近 0 | 2026-07-31 实测 ERMAP 0 条、IFNAR1 1 条(红细胞分布宽度,属定量性状不参与判读) |
E8 | LACTB2 结果未知 | 首次查询;★ 无论结果如何都必须与「没查」区分开 |
预期与实际不符怎么办
一律先按缺陷处理,不许现场解释。查清是数据变了、代码错了、还是预期本身写错了, 判完再写进本页。
三、验收断言(抓不到即 stop())
3.1 全步通用(N 系列)
| # | 断言 |
|---|---|
N1 | 输入集 == coloc_verdict.csv 中 coloc_tier=="strong" & in_mhc==0 的蛋白集合,且恰为 11 个、11 个唯一哨兵 SNP |
N2 | ★★ 跑前跑后 candidate_status.csv 的 retained 集合完全相同(第 6 步无否决权) |
N3 | 产物取值全 ASCII;列名不含内部流程编号(4a step3 等);不含内部结局 ID(T1D_gcst R9_dm 等) |
N4 | ★★ 外部 API 失败必须 stop(),不得当成阴性 —— 本轮已踩三次(VEP ×1、g:Profiler ×2) |
3.2 6-1 专属(P 系列)
| # | 断言 |
|---|---|
P1 | ★ 每个 SNP 都有「已查询」标记;「查了 0 条」与「没查」必须是不同取值 |
P2 | 等位对齐:PheWAS 的 beta 相对 OpenGWAS 自身 effect allele,须翻到 MR 工具的 effect_allele.exposure 后再取符号;回文位点用 EAF 消歧,消不掉的一律剔除并记数 |
P3 | Bonferroni 分母 == 实际查询 SNP 数 × 运行时数据集数(写入 provenance) |
P4 | 跨仓对账:10 个旧仓覆盖的 SNP,在同一阈值 1.06×10⁻⁸ 下关联数与旧仓比对;不一致必须逐条列出并说明(可能是 OpenGWAS 库变动,不是自动通过) |
P5 | 性状三分类穷尽且互斥:每条关联恰属疾病 / 定量 / 分子之一,未分类数 == 0 |
P6 | implied_drug_action 逐蛋白与 coloc_verdict.csv 完全一致(抑制 5 / 激动 6) |
四、★ 已知坑(本轮必须防住)
来自 2026-07-31 的实测,见记忆 project_drug_target_phewas_20260731:
| # | 坑 | 本轮对策 |
|---|---|---|
| 1 | in_MHC 静默丢数据 —— 直接 merge target_decision_table.csv 的 in_MHC,该表只覆盖 31 个蛋白,其余是 NA,in_MHC == FALSE 把它们整批悄悄滤掉且不报错 | 改用 candidate_status.csv 的 in_mhc(覆盖全部 57 行);断言无 NA |
| 2 | 「没测」≠「测了没有」 —— LACTB2 等三个蛋白的 SNP 根本不在 phewas_all.csv 里,被记成「0 条脱靶」 | 断言 P1;LACTB2 正是这次的当事者 |
| 3 | 只按 ID 前缀滤分子性状会漏 —— ebi-a-GCST900102xx 那批「Galectin-4 levels」不带 prot- 前缀 | 补性状名判据(以 " levels" 结尾);断言 P5 |
| 4 | data.table 的 j 里变量名不能叫 p —— 会被解析成 p 值列,汇总表整个错乱且不报错 | 代码审查时专查 |
| 5 | 缓存指纹必须含 SNP 集合 + 查询阈值 —— 否则换了 SNP 集合还复用旧缓存 | 03_phewas.R 已实现;本轮 SNP 集合由 94→11,指纹必然失效 |
| 6 | 等位对齐上次实测是空操作(93/93 本来就一致)—— 但结论是「验证后为空操作」而非「假设它对齐」 | 断言 P2 保留 |
五、产物(预期)
results/phewas/
phewas_raw.csv 原始查询结果(只在真正查询后写一次)
phewas_all.csv 派生列重算后的交付物
phewas_provenance.csv ★ 查询日期 / SNP 数 / 数据集数 / 阈值 / 包版本
query_coverage.csv ★ 新增:逐 SNP 的「已查询 / 关联数」,区分 0 与未查
results/drug_target/
ontarget_phewas.csv ★ 方向感知判读(获益 / 风险 / 不判读)
ontarget_summary.csv 逐蛋白汇总六、6-2 成药性(★ 改为先跑)
6.1 为什么把 6-2 提到 6-1 前面
用户 2026-08-08 决定。结果不受影响(6-1~6-4 互不为输入、都不否决),但有两个实在好处:
- 成药性给 PheWAS 提供对照基准 —— 若靶点已有上市药,PheWAS 的脱靶谱可与该药已知不良反应谱互相印证。C1R 就是这个写法:先查出
C1R有 conestat alfa / C1-酯酶抑制剂(遗传性血管性水肿),再说 PheWAS 无显著关联。 - 不依赖会过期的 token —— Open Targets 公共 API 无需鉴权;OpenGWAS 的 JWT 有效期只有 14 天(本轮就是因过期而中断)。
★ 澄清:PheWAS 的方向感知不依赖 6-2。用药方向已在 coloc_verdict.csv 的 implied_drug_action(inhibit 5 / agonize_or_supplement 6),由 MR 效应方向推得。
6.2 数据源:Open Targets Platform(GraphQL v4)
查询字段:tractability{label modality value}(小分子 SM / 抗体 AB 可及性分桶)、 drugAndClinicalCandidates{count rows{maxClinicalStage drug{name mechanismsOfAction{rows{actionType}}}}} (药名、最高临床阶段、动作类型)。
方向匹配规则(脚本内 NEED_ACTION):
| MR 推得的用药方向 | 需要的 actionType |
|---|---|
inhibit | INHIBITOR / ANTAGONIST / NEGATIVE MODULATOR / BLOCKER / DEGRADER |
activate(= agonize_or_supplement) | AGONIST / ACTIVATOR / POSITIVE MODULATOR / OPENER |
6.3 ★ 要不要加 DrugBank —— 调研结论:先不加,改用 ChEMBL 兜底
| 项 | 事实 |
|---|---|
| Open Targets 的已知缺口 | 官方文档明载:"the Open Targets Platform only includes data from ChEMBL if a drug has both a disease indication and a human protein target (linked via its MOA). Therefore, there will be drugs in ChEMBL that are not present in the Open Targets Platform." → 确实会漏 |
| DrugBank 的获取成本 | 学术许可免费,但 API 访问需付费许可,无自助免费 API;要用得申请许可 + 下载全量 XML |
| ★ 更省的等价补丁 | 直接查 ChEMBL 开放 REST API(无需鉴权)—— 它正是 OT 的上游,能补上 OT 因「必须有疾病适应症 + MOA」而滤掉的那部分 |
定:本轮先跑 Open Targets;凡结果为 no_known_drug 的蛋白,再用 ChEMBL 开放 API 复核。 两者都空、且确实需要适应症级细节时,再考虑申请 DrugBank 学术许可。
理由:DrugBank 相对 OT 的增量主要在「适应症描述」与「非主靶点关系」, 而我们这一步要回答的是「有没有药 / 什么模态 / 到哪期 / 方向对不对」—— 这四项 OT + ChEMBL 已经覆盖。为这点增量付出许可申请的时间成本不划算。
6.4 ★ 预期(写于运行之前)
机械可验
| # | 预期 |
|---|---|
D1 | 11 个蛋白全部解析到 Ensembl ID(NOT_FOUND 与「查到但无药」必须分开) |
D2 | ★ WARS 须经别名表命中 WARS1(HGNC 现行符号)—— 与富集那次栽的是同一个改名问题 |
D3 | 用药方向与 coloc_verdict.csv 一致:inhibit 5 / activate 6 |
内容预期(★ 最要紧的一条)
| # | 预期 | 依据 |
|---|---|---|
D4 | ★★ IFNAR1 会查到已上市药(anifrolumab,抗 IFNAR1 单抗,SLE,Phase 4),但 dir_match 应为 direction_opposite_or_unknown | 我们推得 IFNAR1 需激动/补充(风险降低方向),而现有药是拮抗——方向正好相反 |
D5 | NOTCH2 可能查到抑制类药物(γ-分泌酶抑制剂 / 抗 NOTCH 抗体),方向同样相反 | 我们推得需激动 |
D6 | 其余多数蛋白(ACRBP ERMAP GALNT3 LACTB2 NUDT5 PAM SIGLEC5 WARS)预期 no_known_drug | 均非经典药靶 |
D7 | tractability 多数有 AB 或 SM 分桶(预测性可及性,不等于有药) | OT 的 tractability 是预测分桶 |
D4 若成立,是本步最有价值的产出
它意味着已有的 IFNAR1 药物走的是相反方向 —— 直接用于糖尿病并发症会加重而非改善。 这既是安全性提示,也说明该靶点若要开发需要全新的激动/补充策略,不能重定位现有药。
⚠️ 但必须同时写明:cis-pQTL 推得的方向是「血浆蛋白水平↑/↓ 与疾病风险的关系」, 不等于药理干预的净效应(C1R 自陈的同一条局限:"PheWAS based on a cis-pQTL cannot fully recapitulate pharmacological inhibition")。
6.5 验收断言
N1 产物 11 行、蛋白不重复 · N3 全 ASCII + 无内部编号(本轮已修三处: 可成药/方向一致 等中文取值、release 列的 R9_dm)· N4 网络失败即 exit(1)
D1 全部解析到 Ensembl · D2 dir_match 取值域封闭 · D3 方向计数 5/6
6.6 运行结果
运行:python3 analysis/07_druggability.py 产物:results/druggability.csv(11 行 × 17 列) 断言:12 条全 PASS(含修复后新增的 D4/D5/D6)
| 蛋白 | 用药方向 | 小分子 | 抗体 | 药物数 | 已上市 | 最高期 | 方向判定 |
|---|---|---|---|---|---|---|---|
IFNAR1 | activate | no | yes | 12 | ★ 12 | 4 | direction_consistent |
NOTCH2 | activate | yes | yes | 1 | 0 | 2 | direction_opposite_or_unknown |
APOE | inhibit | yes | yes | 0 | 0 | 0 | no_known_drug |
PAM | activate | yes | yes | 0 | 0 | 0 | no_known_drug |
NUDT5 | inhibit | yes | no | 0 | 0 | 0 | no_known_drug |
WARS | inhibit | yes | no | 0 | 0 | 0 | no_known_drug |
ACRBP | inhibit | no | yes | 0 | 0 | 0 | no_known_drug |
ERMAP | activate | no | yes | 0 | 0 | 0 | no_known_drug |
GALNT3 | activate | no | yes | 0 | 0 | 0 | no_known_drug |
SIGLEC5 | inhibit | no | yes | 0 | 0 | 0 | no_known_drug |
LACTB2 | activate | no | no | 0 | 0 | 0 | no_known_drug |
可及性(tractability,预测分桶,不等于有药):抗体可及 8/11、小分子可及 5/11、 两者皆无 1/11(LACTB2)。
★★ 预期对账:D4 预测错了,但错得有价值
我跑前写的是「IFNAR1 会查到 anifrolumab(拮抗),故方向相反」。实际相反:
INTERFERON BETA-1B / ALFA-2B / ALFA-N3 / ALFACON-1 /
PEGINTERFERON BETA-1A / ALFA-2A / ALFA-2B /
ALBINTERFERON ALFA-2B / ROPEGINTERFERON ALFA-2B /
INTERFERON ALFA-2A / INTERFERON BETA-1A AGONIST / POSITIVE MODULATOR 已上市 x11
ANIFROLUMAB ANTAGONIST 已上市 x1我只想到了 anifrolumab,忘了干扰素本身就是 IFNAR1 的激动剂。IFNAR1 两个方向都有已上市药,direction_consistent 是对的。
★ 这个更正把结论从「现有药方向不对」翻转为「方向一致的已上市药有 11 个」—— 对 IFNAR1 是利好,且指向药物重定位的可能性。
但不能就此说「可以重定位干扰素」
三条必须同时写:
- cis-pQTL 推得的方向 ≠ 药理干预净效应(C1R 自陈的同一条局限)
- 干扰素类有明确且严重的不良反应谱,
IFNAR1通路双向都有临床用途(anifrolumab 治 SLE) - 我们的
IFNAR1证据来自血浆蛋白水平,与受体激动的组织效应不是一回事
★ 抓到的缺陷:APPROVAL 被静默算成 0 期
第一次运行时 IFNAR1 显示 drugs=12 但 maxphase=0。查原始 API 返回:
maxClinicalStage = 'APPROVAL' ← 而 PHASE_NUM 里没有这个键
PHASE_NUM.get('APPROVAL', 0) → 0 ← 12 个已上市药全被记成 0 期又是「.get(x, 默认值) 吞掉未知枚举」这一类。修复两处: ① 补齐 APPROVAL/PRECLINICAL;② 未知取值不再默认 0,改为记账后由断言 D4 报错。 另加 D5(有已上市药者 max_phase 必须 == 4)防止再次退化, 并新增四列:n_approved_drugs / n_approved_direction_match / n_approved_direction_opposite / approved_action_types。
顺带修掉的 N3 违规
产物原先含中文取值(可成药 / 方向一致 / 无药)与内部版本 ID(release=R9_dm), 已全部改为 ASCII 并删列。
6.7 ChEMBL 复核:★ 未完成,EBI 服务当前降级
按 §6.3 的承诺,对 9 个 no_known_drug 靶点做 ChEMBL 复核 (analysis/08_chembl_crosscheck.py)。当前跑不通:
| 路径 | 结果 |
|---|---|
本机 WSL curl | rc=28 超时 / 返回非 JSON |
myserver curl | http=000 / 返回 EBI 的 HTML 页面 |
| WebFetch | HTTP 500 |
三条路径都不通 → 是 EBI 端的问题,不是我们的网络。 (间歇性:GALNT3→CHEMBL4523291、NUDT5→CHEMBL4105713 曾成功解析。)
最终结果(9 个靶点):
ACRBP / APOE / ERMAP / PAM / SIGLEC5 / WARS -> 查询失败 (5/9)
GALNT3 -> CHEMBL4523291 机制记录 0 条
NUDT5 -> CHEMBL4105713 机制记录 0 条
LACTB2 -> ChEMBL 无对应 target
[FAIL] N4 无网络层失败未完成就是未完成
脚本按 N4 硬失败,不会产出「ChEMBL 也没有」的假阴性。 在复核完成前,那 9 个只能写成「Open Targets 未见已知药物」, 不能写「无药可用」 —— 因为 OT 官方明载它会漏掉「无疾病适应症或未连 MOA」的药。
→ 待 EBI 服务恢复后重跑本脚本;DrugBank 是否需要,等复核结果出来再判。
★ 又抓到一个产物设计缺陷:先写盘、后断言
原脚本是「先写 CSV → 再跑断言」。本次断言正确地 FAIL 了, 但磁盘上已经留下一份看起来完整、实则无效的产物 —— 以后谁捡起它都看不出问题。
改法:断言全过才写正式产物;未过则写成 .partial 并以非零码退出, 让无效产物在文件名上就自曝。本次那份已改名 results/druggability_chembl_crosscheck.csv.partial(未删除,留作排查)。
★ 这条应升级为通用约定,与既有三条并列:
④ 断言未全过,不得写出正式产物(写
.partial并非零退出)
七、6-1 运行结果
7.1 查询(已完成)
运行:Rscript analysis/03_phewas.R only=coloc_strong溯源:11 SNP × 50,164 数据集 = 551,804 次检验;查询阈值 1e-5; 主判据 5e-8(固定);Bonferroni 辅助 9.06e-08
★ 首次尝试因旧 token 过期返回 401,脚本按 N4 硬失败退出(exit=1), 未产出任何「0 条脱靶」的假阴性 —— 断言按设计生效。Token 已续期(至 2026-08-22)。
返回 1,659 条关联(5e-8 显著 1,157 条 / Bonferroni 显著 1,200 条)。
| 蛋白 | 哨兵 | MAF | 返回 | 5e-8 显著 | 旧仓(1.06e-08) |
|---|---|---|---|---|---|
APOE | rs429358 | 0.156 | 1262 | ★ 951 | 896 |
PAM | rs149802978 | 0.054 | 88 | 42 | 34 |
NUDT5 | rs10508438 | 0.336 | 68 | 38 | 34 |
WARS | rs2273804 | 0.260 | 53 | 32 | 27 |
SIGLEC5 | rs1106476 | 0.113 | 59 | 32 | 29 |
NOTCH2 | rs2641348 | 0.111 | 64 | 30 | 25 |
ACRBP | rs7959658 | 0.361 | 41 | 20 | 19 |
ERMAP | rs11210710 | 0.476 | 6 | 6 | 6 |
IFNAR1 | rs914142 | 0.271 | 9 | 4 | 4 |
GALNT3 | rs2116546 | 0.298 | 9 | 2 | 2 |
LACTB2 | rs191588099 | ★ 0.0071 | 0 | 0 | ★ 旧仓从未查过 |
7.2 ★ 预期对账
| # | 预期 | 实际 |
|---|---|---|
E1 | 11 个 SNP 全查,含 rs191588099 | ✅ |
E2 | 每 SNP 都有「已查询」记录,哪怕 0 关联 | ✅ 覆盖表 11 行,LACTB2 = queried + 0 条 |
E3 | Bonferroni 分母 = SNP 数 × 数据集数 | ✅ 11 × 50,164 = 551,804 |
E5 | 各 SNP 关联数 ≥ 旧仓(本轮阈值更松) | ✅ 10/10 全部 ≥(见上表右两列) |
E6 | APOE 远超其余 1–2 个数量级 | ✅ 951 vs 次高 42,约 23 倍 |
E7 | 疾病类脱靶 ERMAP/IFNAR1 接近 0 | ⏸ 待方向感知判读(需性状分类) |
E8 | LACTB2 结果未知 | ✅ 返回 0 条 —— 但见下 |
其他:别名 rsID 按坐标归一 3 → 0 条未能归一;蛋白归属缺失 0 条。
7.3 ★★ LACTB2 的 0 条不能读成「干净靶点」
rs191588099 的 MAF = 0.0071(0.71%),比次低的 PAM(0.054)还稀有一个数量级, 是 11 个里唯一的罕见变异,也是唯一返回 0 条的。
OpenGWAS 的 phewas 只返回该变异在各数据集中达阈的关联; 变异若根本不在某数据集里,同样不返回。而大量 GWAS 汇总统计做过 MAF>1% 过滤。
这是「没测 ≠ 测了没有」的第三层
| 层 | 问题 | 状态 |
|---|---|---|
| ① | SNP 根本没被提交查询 | 2026-07-31 踩过,已由 query_coverage.csv 修掉 |
| ② | 查询失败被当成阴性 | 已由 N4 修掉(本轮 401 即验证) |
| ③ | ★ 查了,但多数数据集里可能没有这个变异 | 覆盖表记录不了这一层 |
→ LACTB2 只能写成「该罕见变异在 OpenGWAS 中未见达阈关联, 多效性无法评估」,不得写成「无脱靶」。
跨 11 个蛋白的 MAF 与关联数 Spearman 相关仅 −0.241(APOE MAF 0.156 却有 1262 条), 所以不是 MAF 普遍驱动关联数;LACTB2 是单点的稀有性问题,不是全局趋势。
7.4 ⏸ 方向感知判读:待决
余下一半(获益/风险判读、疾病 vs 定量 vs 分子的三分类,即 E4/E7) 由 37_drug_target_phewas.R 完成,但它有两处依赖不可用:
| 依赖 | 问题 |
|---|---|
results/targets/target_decision_table.csv | ★ tier 表,用户已定删除,v3 不存在 → 改读 candidate_status.csv |
results/phewas/dataset_domains.csv | ★★ v3 没有,且全流水线没有任何脚本生成它 —— 旧仓那份是临时手工产物,domain_source 列还含中文(违反 N3) |
性状三分类必须先解决。 用户 2026-08-08 定:重建为可重跑脚本。
7.5 性状三分类:新建 48_trait_domains.R
输入 OpenGWAS 官方元数据(gwasinfo,1,500 个数据集,带指纹缓存) 输出 results/phewas/dataset_meta.csv · dataset_domains.csv · dataset_domains_residual_review.csv
| 类 | 定义 | 数量 |
|---|---|---|
disease | 病例-对照 / 疾病性状 —— 方向可读成风险升高或降低 | 263 |
quantitative | 临床或人体测量的连续性状(HbA1c、血压、BMI…)—— 方向有临床含义但需逐个解释 | 698 |
molecular | 组学平台读数(蛋白/代谢物/表达/脑影像)—— 不参与获益/风险判读 | 539 |
判定级联(来源标签与判定路径同步写,见下方教训):
1. id 前缀 ∈ {prot, met, eqtl, ubm} -> molecular (539)
2. subcategory == "Protein" -> molecular
3. ★ ncase > 0(病例-对照,最硬信号) -> disease (248)
4. category ∈ {Disease, Binary} -> disease (122 含定量)
5. subcategory 明示疾病系统 -> disease (5)
6. category ∈ {Continuous, Risk factor, ...} -> quantitative
7. ★ 末位兜底:性状名明示疾病 -> disease (9)
8. 其余 -> quantitative (577)断言 7 条全 PASS(T1–T7 + N3)。
★★ 推翻一条旧规则
旧记忆写「性状名以 levels 结尾 ⇒ 分子性状」。本轮实测该规则过度捕获:
Hemoglobin A1c levels · Glucose levels · Calcium levels
LDL cholesterol levels · Total testosterone levels在一个糖尿病研究里把 HbA1c 与血糖剔出判读是严重错误。 已废弃该规则, 改以 ncase + category + ID 前缀为准,并立 T4 永久防守。
★ 首版的第二个错:category 不能用来判平台来源
首版用 category 匹配 metabolit|protein 判分子,误伤 4 个: ukb-d-30740(血糖)、ukb-d-30750(HbA1c)是 UK Biobank 血生化字段, OpenGWAS 却把 category 标成 Metabolites。由 T4 抓出。 改为平台来源只认数据集 ID 前缀。
★★ 我误报了一次,根因是溯源列本身写错了
我看到 domain_source == "fallback_quantitative" 有 628 个,报了 「628 个疾病可能被当成定量」。这是假警报。
实际分类是对的(Alzheimer's、T2D、Early AMD、Cataracts 全部正确判为 disease, 走的是 ncase 规则)。错的是 domain_source 这一列 —— 它在级联之外单独计算, 把 88 个「靠 ncase 判为疾病」的行标成了 fallback_quantitative。
一个不反映真实判定路径的溯源列,比没有更糟 —— 它会主动误导。
已改为判定与来源在同一次赋值里写,并立 T6 永久防守。
T7:残余风险已量化并清零
补了 T7 度量「兜底判为定量的里面还混着多少疑似疾病」(这是低估疾病脱靶的方向, 对安全性筛查不利)。首次测得 9 个(CAD (Firth correction)、T2D (SPA correction) 这类同一疾病的不同统计校正版本),加末位兜底规则后 降至 0。
⚠️ 但关键词法本身不可靠 —— 本轮它就漏过 Early age-related macular degeneration (那个是靠 ncase 救回来的)。故该规则只放在最后一位,且单列来源标签供复核。
7.6 方向感知判读:新建 49_ontarget_phewas.R
未沿用 37_drug_target_phewas.R,改为新脚本,原因两条:
- 37 读
results/targets/target_decision_table.csv(tier表,已定删除,且只覆盖 31 个蛋白 —— 直接 merge 会让其余蛋白的in_MHC变 NA 后被整批静默滤掉,第一版就踩过) - ★★ 37 内含一条已被本轮推翻的性状再分类规则:
grepl("\\blevels$", trait) => molecular
移植了 37 最有价值的部分:等位翻转矩阵、回文位点用 EAF 消歧、随机抽样自检。
等位对齐:1,641 条可对齐 / 18 条无法对齐(回文不可消歧或碱基不匹配,已剔除)。 对齐旁证:命中自身结局 1 条,方向与 MR 同号 1/1。
断言 12 条全 PASS(N1–N3 + O1–O9),其中 N2 校验 跑前跑后 retained 集合完全相同 —— 第 6 步无否决权,已在代码里落实。
★ 判读前抓到的两处「疾病类」渗漏
| # | 现象 | 根因 | 修法 |
|---|---|---|---|
| 1 | SIGLEC5/NOTCH2/APOE 的「疾病脱靶」里出现 basophil / monocyte / lymphocyte cell count | 这三个的 category 是 Continuous,但 subcategory 是 Immune system;而我把 subcategory 关键词规则排在 Continuous 之前,immune 命中了。「Immune system」是身体系统标签,不是疾病标志 | 连续性状规则前移;关键词由裸 immune 改为 autoimmune;立 T8 |
| 2 | PAM 出现 3 条「脱靶获益」:metformin 用药码、ICD10 E11.9 | on_target_axis 只查了 diabet 字样,漏掉 UKB 的 ICD 编码型与用药编码型糖尿病代理 | 补 E10/E11 ICD 正则与 16 种降糖药名;立 O8/O9 |
⚠️ 反向也要守住:APOE 的他汀/依折麦布用药码是真脱靶(反映脂质轴),不能一起剔除。
修复前后:PAM 疾病脱靶 3 → 0;SIGLEC5 2 → 0;NUDT5 2 → 1;NOTCH2 4 → 3。
八、★ 第 6 步当前结论:11 个蛋白
| 蛋白 | 用药方向 | 索引结局 | OR | 显著脱靶 | 疾病类 | 定量 | 分子 | 获益 | 风险 | 已上市药 |
|---|---|---|---|---|---|---|---|---|---|---|
APOE | inhibit | 肾病 | 1.151 | 926 | 131 | 469 | 326 | 25 | ★ 106 | 0 |
NUDT5 | inhibit | 视网膜 | 3.080 | 29 | 1 | 24 | 4 | 1 | 0 | 0 |
SIGLEC5 | inhibit | 黄斑 | 1.105 | 29 | 0 | 23 | 6 | 0 | 0 | 0 |
WARS | inhibit | 视网膜 | 1.385 | 28 | 0 | 20 | 8 | 0 | 0 | 0 |
NOTCH2 | activate | 视网膜 | 0.598 | 24 | 3 | 15 | 6 | 0 | 3 | 0 |
ACRBP | inhibit | 视网膜 | 1.450 | 20 | 0 | 8 | 12 | 0 | 0 | 0 |
PAM | activate | 视网膜 | 0.888 | 18 | 0 | 13 | 5 | 0 | 0 | 0 |
ERMAP | activate | 黄斑 | 0.329 | 6 | 0 | 0 | 6 | 0 | 0 | 0 |
IFNAR1 | activate | 黄斑 | 0.782 | 4 | 0 | 1 | 3 | 0 | 0 | ★ 12 |
GALNT3 | activate | 视网膜 | 0.846 | 2 | 0 | 0 | 2 | 0 | 0 | 0 |
LACTB2 | activate | 视网膜 | 0.674 | NA | NA | NA | NA | NA | NA | 0 |
8.1 三条主要发现
① 8/10 可评估的蛋白疾病类脱靶 = 0。 SIGLEC5 WARS ACRBP PAM ERMAP IFNAR1 GALNT3 全为 0,NUDT5 仅 1 条(FEV1/FVC < 0.7,判为获益)。 ★ 其中 ERMAP 与 IFNAR1 的 0 复现了 2026-07-31 的独立结果(预期 E7 命中)。
② APOE 是唯一的高多效离群点:926 条显著脱靶、131 条疾病类、106 条判为风险 (阿尔茨海默病、冠心病、脂质轴)。这是 rs429358(APOE ε4)的已知特性,不是新发现。
③ NOTCH2 有 3 条一致的风险信号,全部是克罗恩病(wald 0.705–0.869,p 最小 2.8e-11)。 我们推得 NOTCH2 需激动,而激动会升高克罗恩病风险 —— 这是本步唯一一条 「靶点特异、方向明确」的安全性提示,值得单独写。
8.2 ⚠️ 两条必须随结果一起写的局限
① LACTB2 不可评估,不是「干净」。rs191588099 的 MAF=0.0071, OpenGWAS 返回 0 条 —— 见 §7.3 的「第三层」。
② ncase 判疾病无法区分「真疾病」与「二元问卷项」。APOE 的 131 条里约 29 条是问卷/自报项(Spread type: Butter、 Illnesses of mother:、Major dietary changes 等)。
⚠️ 而我为量化它写的关键词筛查本身也过度捕获 —— 它把 Nonalcoholic fatty liver disease 和 Delirium, not induced by alcohol... 按 alcohol 误命中了。关键词法在两个方向都不可靠。
→ 故不再加关键词补丁,而是如实记录:该局限只影响 APOE (其余 10 个蛋白的疾病类脱靶分别是 3、1、0×8,逐条看过,无问卷项)。
九、6-4 组织 / 单细胞(跑前记录)
状态:尚未运行。范围:11 个蛋白(用户 2026-08-08 定)
9.1 为什么 6-4 做 11 个而 6-3 只做 IFNAR1
用户问「C1R 最后聚到一个蛋白,我们是不是也只做 IFNAR1」。核对原文,它是两层结构:
比较层(4 个都做):"Tissue expression profiling of the four prioritized candidates…" "We next evaluated the translational potential of the four candidates…" "C1R emerged as the most compelling target among the four candidates"
深挖层(只有 C1R):Figure 5 衰老相关疾病表达、UKB 观察性验证、实验验证规划
→ 聚焦到一个是比较之后的结论,不是分析的起点。 它的组织表达(FUMA)做的是 4 个。
我们的比较层其实比 C1R 还厚一层(共定位 11 → 成药性 11 → PheWAS 11),结论已出: IFNAR1 是唯一「共定位 strong + 已上市药方向一致 + 脱靶谱干净」三样齐全的。
但 6-4 仍须做 11 个 —— 一条具体理由
阴性结果只有在有对照时才可解释。
上一轮实测 IFNAR1 在视网膜 eQTL 层完全阴性(见 9.4)。若只跑它, 拿到阴性后分不清是「视网膜不是相关组织」还是「这套数据在该位点没功效」。 11 个一起跑、别的蛋白(如 ERMAP)在同一套数据里出信号, IFNAR1 的阴性才是有意义的阴性。这个内部对照白送。
→ 写作时:11 个的表进补充材料,正文写「在 11 个候选中比较后 IFNAR1 因⋯⋯被选为主推」, 让「聚焦到一个」是读者跟着推出来的。
9.2 文献底本
| # | 来源 | 内容 |
|---|---|---|
| R1 | Advani J, Mehta PS, et al. QTL mapping of human retina DNA methylation identifies 87 gene-epigenome interactions in age-related macular degeneration. Nat Commun 2024;15:1972 · PMID 38438351 · doi 10.1038/s41467-024-46063-8 | 视网膜 eQTL 403 眼 / mQTL 152 眼 / eQTM;Zenodo 开放(record 10569726)。★ 本地有 PDF:F:\project\s41467-024-46063-8.pdf |
| R2 | Single-cell atlas of the transcriptome and chromatin accessibility in the human retina. Nat Genet 2025 · PMID 41578023 · doi 10.1038/s41588-025-02454-1(通讯 Rui Chen) | HRCA:约 390 万细胞 / 125 供者 / >130 细胞类型,CELLxGENE 开放 |
| R3 | C1R(medRxiv 2026.07.14.26358022)Methods "Gene expression analyses" | 用 FUMA GENE2FUNC(基于 GTEx 的通用组织) |
★ 我们比 C1R 强的地方要写进 Discussion:他们用 GTEx 通用组织做描述性表达定位; 我们用视网膜专属 eQTL 做 MR + 共定位(因果层),再叠单细胞定位(描述层)。 ⚠️ 两层性质不同,不能混着说 —— 见记忆 project_mr_beginner_skill_20260805 「单细胞分 sc-eQTL MR(因果)vs 表达定位(描述)」。
9.3 数据现状(★ 一处违规必须先处理)
| 数据 | 位置 | 状态 |
|---|---|---|
| HRCA 单细胞 | D:\mrdata\singlecell\hrca_allcells.h5ad(2.5 GB) | ✅ 合规 |
| Advani 视网膜 eQTL | /home/research/multiomics-mr/05_retina/data/eQTL_retina.txt(384 MB / 202 万行) | ★ 不在 D 盘,违反「原始数据一律存 D 盘」 |
★ D 盘现有的 D:\mrdata\eqtl_retina_strunz\ 是另一个来源(Strunz,FastQTL nominal 分染色体), 不是 Advani,不能混为一谈。
本地 HRCA 是子集,不是全图谱
本地 h5ad 为 265,767 细胞 / 18 个细胞类型(含 RPE), 是 CELLxGENE 该 collection 下的一个数据集,不是 R2 描述的 390 万细胞全图谱。 Methods 里必须写清用的是哪个子集,不能引用全图谱的规模。
9.4 ★ 上一轮结果(实测复核,非记忆引用)
multiomics-mr/05_retina/results_eqtl_r9dm/ 只覆盖了 3 个基因:
| 基因 | 结局 | 方法 | b | p | PP.H4 |
|---|---|---|---|---|---|
ERMAP | 黄斑 | Wald ratio | ★ +0.0929 | 2.52e-05 | ★ 0.9488 |
ERMAP | 视网膜 | Wald ratio | +0.0058 | 0.661 | 0.0081 |
IFNAR1 | 黄斑 | IVW | +0.0031 | 0.790 | 0.0189 |
IFNAR1 | 视网膜 | IVW | −0.0001 | 0.975 | 0.0098 |
APOL1 | — | — | — | — | 0 条记录(非 eGene) |
★★ ERMAP 的方向冲突已由数据坐实(不是凭记忆): 视网膜 eQTL 层 b = +0.0929(视网膜表达↑ → 黄斑病变风险↑), 而血浆层 OR = 0.329(血浆 ERMAP↑ → 风险↓)。两层方向相反。
⚠️ 但记忆 project_multiomics_mr 载明:跨层往往不是同一个变异 (实测 r² 仅 0.095 / 0.031),故「方向相反」不构成对血浆层发现的反对。 这一条必须随结果一起写,否则会被读成自相矛盾。
9.5 要改的两处接线(★ 与前面同一类问题,已是第三次)
| 脚本 | 现状 | 改法 |
|---|---|---|
multiomics-mr/05_retina/05d_retina_eqtl_r9dm.R | R9 <- "/home/research/mr-pipeline-r9"(旧仓);靶点读 common/targets_r9dm.csv(只有 3 个) | 指向 mr-pipeline-r9v3;靶点改由 coloc_verdict.csv 生成 11 个 |
mr-pipeline-r9v3/analysis/_singlecell_io.py | load_targets() 读 results/network_R9_dm/target_list.csv(tier 表,v3 不存在) | 改读 coloc_verdict.csv |
9.6 ★ 预期(写于运行之前)
| # | 预期 | 依据 |
|---|---|---|
S1 | Advani 只覆盖 9,408 个基因(非全转录组)→ 11 个里会有若干个不是 eGene(0 条记录) | 记忆 reference_retina_eqtl 明载 |
S2 | ★ ERMAP 黄斑 PP.H4 复现 ≈ 0.949、b 为正;IFNAR1 两个结局均阴性 | 9.4 实测 |
S3 | APOL1 不在本轮 11 个里(4X 已否决),不会出现 | 第 4X 步 |
S4 | 单细胞:IFNAR1 曾测得伞状神经节 74.4% / AUC 0.757;ERMAP AUC 0.536 = 无富集 | targets_r9dm.csv 的 note 列 |
S5 | 至少 1 个非 ERMAP 的蛋白在视网膜层有信号 —— 否则 IFNAR1 的阴性无法解释(见 9.1) | 本步设计目的 |
9.7 验收断言
| # | 断言 |
|---|---|
N1 | 输入集 = coloc_verdict.csv 的 strong 非 MHC,恰 11 个 |
N2 | ★★ 跑前跑后 candidate_status.csv 的 retained 集合完全相同 |
N3 | 产物全 ASCII、无内部流程编号与内部结局 ID |
N4 | ★ 输入文件缺失 / 读取失败必须 stop(),不得当成「该基因无记录」 |
R1 | ★ ENSG 号一律取 druggability.csv 的 ensembl 列,并用坐标与哨兵 SNP 交叉核实同一 cis 窗口 —— 2026-07-31 正是 ENSG 写错导致 ERMAP 被误判为「做不了」 |
R2 | ★★ 「非 eGene(查了,0 条)」与「没查」必须是不同取值 —— 与 6-1 的 query_coverage 同一约定 |
R3 | ★ 复现对照:ERMAP/IFNAR1 的 4 行结果与 9.4 的数值一致(容差按浮点) |
R4 | 每个有信号的基因都必须同时记录血浆层方向与视网膜层方向,并标注跨层是否同一变异(r²) |
R5 | 单细胞:基因符号解析失败(missing)必须显式列出,不得静默当成「不表达」 |
9.8 数据已入 D 盘(2026-08-08)
D:\mrdata\eqtl_retina_advani2024\
eQTL_retina.txt 384 MB md5 5a696f5ba26dfadd83dd09dfa9cf9647
mQTL_retina.txt 522 MB md5 d49613b66fbcc4e53b2a4e957803840f
eQTM_retina.txt 1.8 MB md5 a7570df143406a26ddb73f944038d53d三个文件 md5 与 WSL 源文件全部一致。50_retina_eqtl.R 已改为从 D 盘读取。
⚠️ EUR.frq(1000G 派生的频率文件,436 MB)仍在 multiomics-mr/05_retina/data/, 它是派生数据不是原始数据,暂不纳入 D 盘规则。
十、6-4 视网膜 eQTL 运行结果
脚本 analysis/50_retina_eqtl.R(新建) 产物 results/retina_eqtl/范围 11 个蛋白
归类更正:这一层是因果层,不是描述性注释
v4 把第 6 步定义为「描述性注释」,但本脚本做的是 「视网膜 eQTL → 眼病」的 MR + 共定位,属因果方法。
→ 写 Methods 时必须单列为「组织层因果验证」,不得并入描述性注释。 N2 断言仍锁住「不改变任何候选去留」,所以它依然无否决权 —— 但方法归类要正确。
★ 对比:单细胞那一层是真正的描述(只查表达定位,不做 MR、不做共定位)。 两层性质不同,见记忆 project_mr_beginner_skill_20260805。
10.1 eGene 覆盖(★ 已回答 S1)
| 基因 | ENSG | 记录数 | 显著 | 视网膜 eGene |
|---|---|---|---|---|
GALNT3 | ENSG00000115339 | 924 | ★ 782 | ✔ |
NUDT5 | ENSG00000165609 | 278 | 241 | ✔ |
IFNAR1 | ENSG00000142166 | 172 | 167 | ✔ |
LACTB2 | ENSG00000147592 | 147 | 147 | ✔ |
ERMAP | ENSG00000164010 | 136 | 136 | ✔ |
WARS | ENSG00000140105 | 96 | 96 | ✔ |
ACRBP | ENSG00000111644 | 25 | 25 | ✔ |
APOE | ENSG00000130203 | 4 | 4 | ✔ |
NOTCH2 PAM SIGLEC5 | — | 0 | 0 | ✘ 非 eGene |
8/11 是视网膜 eGene(S1 预期命中)。R1b 坐标交叉核实通过 (哨兵与 eQTL 基因最大距离 0.952 Mb,均在 cis 窗口内)。 工具:F≥10 过滤后 1,279 条(最小 F=15.4,中位 27.2)。
★ 注意 LACTB2:它在 6-1 PheWAS 里因罕见变异(MAF 0.0071)不可评估, 但在这一层有 147 条记录、可评估 —— 两层的可评估性不同,须分别叙述。
10.2 MR + 共定位结果(全 12 行)
| 基因 | 结局 | nsnp | b | p | FDR | PP.H4 | vs 血浆 |
|---|---|---|---|---|---|---|---|
ERMAP | 黄斑 | 1 | +0.0929 | 2.5e-05 | 1.5e-04 | ★ 0.9488 | 相反 |
APOE | 黄斑 | 1 | −0.0079 | 5.3e-03 | 0.016 | NA | 相反 |
WARS | 黄斑 | 1 | −0.0068 | 0.023 | 0.045 | 0.1735 | 血浆不显著 |
NUDT5 | 黄斑 | 7 | −0.0081 | 0.20 | 0.30 | 0.0131 | — |
GALNT3 | 黄斑 | 6 | +0.0077 | 0.57 | 0.69 | 0.0161 | — |
IFNAR1 | 黄斑 | 3 | +0.0031 | 0.79 | 0.79 | 0.0189 | 相反 |
WARS | 视网膜 | 1 | −0.0077 | 2.0e-05 | 1.2e-04 | ★★ 0.8086 | 相反 |
NUDT5 | 视网膜 | 7 | −0.0201 | 8.7e-04 | 2.6e-03 | 0.0346 | 相反 |
APOE | 视网膜 | 1 | −0.0030 | 0.083 | 0.17 | NA | 相反 |
GALNT3 | 视网膜 | 6 | +0.0113 | 0.33 | 0.49 | 0.0089 | — |
ERMAP | 视网膜 | 1 | +0.0058 | 0.66 | 0.79 | 0.0081 | — |
IFNAR1 | 视网膜 | 3 | −0.0001 | 0.97 | 0.97 | 0.0098 | — |
steiger_ok 全部 TRUE。
10.3 ★★ 三条发现
① 新增一个强共定位:WARS × 视网膜病变 PP.H4 = 0.809、p = 2.0e-05。 上一轮只跑 3 个基因,从未见过它。若只跑 IFNAR1 会完全错过。
② ERMAP 黄斑 PP.H4 = 0.9488 与旧仓数值完全复现(R3b 通过), 且 b = +0.0929 与血浆 OR 0.329 方向相反,由数据坐实(非记忆引用)。
③ IFNAR1 在这一层彻底阴性:黄斑 p=0.79 / PP.H4=0.019, 视网膜 p=0.97 / PP.H4=0.010。
这个阴性只有靠对照才立得住 —— S5 的意义
同一套数据、同一套流程里,ERMAP 出 0.949、WARS 出 0.809, GALNT3 有 782 条显著 eQTL 记录。 → IFNAR1 的阴性不是数据没功效,是真阴性。
若只跑 IFNAR1,交付的将是「两个结局全阴性、无任何同数据集对照」, 分不清「视网膜不是相关组织」与「这套数据在该位点没功效」。
10.4 ⚠️ 跨层方向必须带的限定
ERMAP(黄斑)、WARS(视网膜)、NUDT5(视网膜)、APOE 的视网膜层方向 都与血浆层相反。
这不构成对血浆层发现的反对 —— 本课题实测跨层往往不是同一个变异 (r² 仅 0.095 / 0.031,见记忆 project_phase2_multiomics_r9dm_20260801)。 产物已加 direction_vs_plasma 列,解读时必须同时给出这一条,否则读起来像自相矛盾。
十一、★ 数据版本核对(2026-08-08)
用户要求核实「现在用的数据是不是最新最全的」。逐层查证:
| 层 | 我们用的 | 当前最新 | 判断 |
|---|---|---|---|
| 血浆 pQTL(主线暴露) | UKB-PPP discovery n=34,557 | full n=54,219;2025-10 Olink+SomaScan 整合图谱(>90,000) | ⚠️ 故意不用最大的(discovery 有留出复制;full 原文自称 putative),已在定稿决定里 |
| 视网膜 eQTL | Advani 2024,403 眼 | 同一个 | ✅ 已是最大的开放数据。EyeGEx 406 眼受控访问且 312/406 是病例;mega-analysis 2020 仅 311 眼 |
| 单细胞 | scRNA 265,767 / 18 类 | snRNA 3,177,310 / 31 类 | ★ 差 12 倍 —— 已换,见 11.1 |
| 血液 eQTL | eQTLGen Phase 1 · 31,684 | ★ Phase 2 · 43,301(medRxiv 2026-02 · PMID 42051578) | ★ 落后 37%,属 6-3 范围 |
| sQTL / 多组织 eQTL | GTEx v8 | ★ GTEx v10(2024-11,eQTL 样本 +23%) | ★ 落后一个大版本,属 6-3 范围 |
| 代谢物 · 血液 mQTL | Nightingale meta_EUR · GoDMC | 未核 | 本轮不在范围内 |
★ 两处落后都在 6-3 跨组学层,后果是功效略低而非方向改变。 ⚠️ eQTLGen Phase 2 的 cis-eQTL 汇总统计能否公开下载我未查实(官网 Resources 页 404)。
11.1 单细胞换全集:查清了「更全」到底全在哪
同一 CELLxGENE collection(doi 10.1038/s41588-025-02454-1)下 8 个数据集,合计 551 万细胞。 两个 all-cells 版本对比:
| 细胞数 | 类型数 | 大小 | |
|---|---|---|---|
scRNA-seq all cells(原用) | 265,767 | 18 | 2.54 GB |
snRNA-seq all cells(已下载) | 3,177,310 | 31 | 37.67 GB |
逐类型差集核对结论:
- ✅ RPE、小胶质、Müller、星形胶质,两个版本都有 —— 关键大类没丢
- 多出的 15 类几乎全是双极细胞与节细胞的精细亚型(
diffuse bipolar 1/2/3a/3b/4/6、ON/OFF midget ganglion…),对本课题不关键 - ★ 真正收益在细胞数(12 倍) —— 小胶质这类稀有细胞在 265k 里可能只有几百个, 在 3.18M 里是几千,对检出功效与 AUC 富集计算差别是实的
换数据集也解决不了的硬伤
两个版本都没有内皮细胞与周细胞。 糖尿病视网膜病变是血管病, 这一层缺失必须写进 Limitations。
下载:wget -c 直连 11.1 MB/s(并发 4 路仅 8 MB/s,瓶颈在上行带宽,不值得分块)。
十二、当前范围审计(★ 逐产物点过,非印象)
| 步骤 | 蛋白数 |
|---|---|
| 第 2 步 MR 显著 | 31 |
| 第 3 步 反向 MR | 95(含未显著者) |
| 4b 换队列 | 66 |
| 4c 两侧都换 | 12 |
| 第 5 步 共定位 | 28 |
| 6-1 PheWAS 查询 / 判读 | 11 |
| 6-2 成药性 | 11 |
| 6-4 视网膜 eQTL | 11 |
| 6-3 跨组学 | 11(2026-08-08 由 IFNAR1 改定) |
★ 没有任何一步只做了 1 个蛋白。
十三、★ Voigt 2026 单细胞 eQTL 原文核读(2026-08-08)
出处:bioRxiv 10.64898/2026.03.30.714946,2026-04-01 posted,未经同行评议。 读的是本地 PDF F:\project\nihpp-2026.03.30.714946v1.pdf, pdftotext -layout 抽出 842 行,逐行读完全部 28 页(含 SI Table 1/2 与全部 46 条参考文献)。
13.1 它到底是什么
| 项 | 值 |
|---|---|
| 设计 | 单细胞(snRNA-seq)cis-eQTL,非 bulk |
| 组织 | 黄斑区的视网膜 + RPE + 脉络膜 |
| 供体 | 新测 88 供体 / 89 眼 + 既往 37 例,合计 122 眼 |
| 进入 eQTL 的供体 | 视网膜 78 人;RPE/脉络膜 108 人 |
| 细胞核 | 432,415(视网膜 106,000 / RPE-脉络膜 326,415) |
| 基因分型 | Illumina Infinium Genome Diversity Array-8 v1.0,1.8M 位点 |
| 填补 | Michigan Imputation Server · MiniMac4 · Eagle · 1000G Phase3 · Rsq > 0.3 |
| eQTL 方法 | MatrixEQTL 线性回归;协变量 = 10 个 PEER 因子 + 性别 + 实验批次 + 表达 PC |
| cis 窗口 | 基因 ±1 Mb |
| 显著性 | FDR < 0.05 |
| 产出 | 8,393 typed + 91,987 imputed eSNP,10,298 个 eGene |
13.2 ★★ 它确实有内皮和周细胞 —— 这正是 HRCA 的窟窿
SI Table 2 列出进入 eQTL 分析的 22 个细胞类型及供体数:
| 类别 | 细胞类型(供体数) |
|---|---|
| ★ 血管 | Choriocapillaris 79 · Artery 75 · Vein 81 · Pericyte 79 |
| RPE | RPE 80 |
| 脉络膜其他 | Fibroblast 103 · Melanocyte 92 · Schwann 76 |
| 免疫 | T-Cell 98 · Macrophage-res 93 · Macrophage-inf 83 · Dendritic 76 · B-Cell 59 |
| 视网膜 | Amacrine 74 · Cone 72 · Rod 72 · OFF-BC 71 · ON-BC 71 · Rod-BC 71 · Horizontal 69 · Müller 69 · RGC 69 |
HRCA 两个版本都没有内皮细胞与周细胞,而糖网是血管病 —— 这篇是目前唯一能补上这一格的人眼数据。
13.3 ★ 但它现在不能直接拿来用,四条硬障碍
| # | 障碍 | 依据(原文) |
|---|---|---|
| 1 | ★★ 没有汇总统计 | 全文没有 Data Availability 章节,没有任何 accession 号。讨论结尾只有一句「all raw and processed gene expression data as well as all genotype data can be downloaded at publicly available data repositories」,未说 eQTL 汇总统计。SI File 2–5 是差异表达结果,不是 eQTL |
| 2 | ★ MAF ≥ 5% 才测 | Methods:「a minor allele frequency of at least 5%」→ 低频哨兵直接出局(LACTB2 MAF=0.0071 必然测不到) |
| 3 | 只报 FDR<0.05 的显著对 | 共定位需要全量名义 p —— 没有全量就做不了 coloc,只能做「查表」 |
| 4 | 样本量远小于 Advani | 视网膜 78 人 vs Advani 403 眼 |
13.4 还有三条会影响解释
- 细胞类型供体数 < 50 的被整类排除;且供体在该类型需 ≥ 10 个细胞
- 疾病谱是 AMD(89 眼里 40 例 AMD),不是糖尿病;供体高龄(对照组 38 人 > 70 岁)
- 组织限定黄斑,与我们「黄斑病变」结局对得上,与「视网膜病变」结局对不齐
13.5 ★ 补充材料实查(media-1.xlsx / media-2.xlsx,2026-08-08)
是什么:各 31 个 sheet = 31 个细胞类型(比 SI Table 2 的 22 类更细,含双极细胞亚型)。 列结构完全相同:
gene · pct_cluster_exp · donor_cell_cutoff · group_1_donors · group_2_donors
logFC_pseudobulk · logCPM_pseudobulk · F_pseudobulk · PValue_pseudobulk · FDR_pseudobulk→ 是 pseudobulk 差异表达(SI File 2–4 中的两个),不是 eQTL 汇总统计。media-1 为「10 供体 vs 18 供体」、media-2 为「4 供体 vs 18 供体」; ⚠️ 文件内没有比较标签,无法从文件本身断定各自对应干性 AMD / GA / MNV 中的哪一个。
★ 再次印证 13.3 第 1 条:eQTL 结果确实没有公开。
11 个候选基因的实测覆盖
| 基因 | 出现于几个细胞类型 | 四类血管细胞 |
|---|---|---|
APOE ERMAP IFNAR1 LACTB2 NOTCH2 NUDT5 PAM WARS | 31/31 | ✅ 全覆盖 |
GALNT3 | 29/31 | ✅ 全覆盖 |
ACRBP | 15/31 | ❌ 无 |
SIGLEC5 | 1/31 | ❌ 无 |
★ ACRBP / SIGLEC5 的缺席是**「未进入检验」不是「不表达」**(未过表达过滤), 用的时候必须三态记账,不得写成阴性。
★★ WARS 在这份数据里就叫 WARS,不是 WARS1 —— 与 HRCA 相反。 别把别名规则写死,要按数据集分别处理。
血管细胞里的表达比例(pct_cluster_exp)
| 基因 | 脉络膜毛细血管 | 动脉 | 静脉 | 周细胞 |
|---|---|---|---|---|
PAM | 0.369 | 0.405 | 0.431 | 0.278 |
WARS | 0.170 | 0.215 | 0.207 | 0.074 |
NOTCH2 | 0.065 | 0.081 | 0.086 | 0.173 |
IFNAR1 | 0.075 | 0.098 | 0.079 | 0.077 |
APOE | 0.069 | 0.078 | 0.101 | 0.090 |
NUDT5 | 0.036 | 0.050 | 0.044 | 0.059 |
GALNT3 | 0.018 | 0.016 | 0.042 | 0.013 |
ERMAP | 0.017 | 0.019 | 0.018 | 0.020 |
LACTB2 | 0.016 | 0.016 | 0.022 | 0.018 |
⚠️ 差异表达那几列对我们没用:11 个基因在四类血管细胞里 FDR 几乎全是 0.9999 (唯一例外 LACTB2 × 周细胞 = 0.978),即 AMD vs 对照无差异 —— 本来也不该有,我们的病因是糖尿病不是 AMD。
→ 有用的只有 pct_cluster_exp 一列,可作 6-4b 血管细胞那一格的补充。
13.6 结论(不要越过)
① 它替代不了 Advani 做 6-4a。 6-4a 是 MR + 共定位,需要全量汇总统计和更大样本。
② 它可能补 6-4b 的一格,但仅限描述性表达定位,且必须先拿到数据。
③ 现在的正确动作是联系作者 / 等正式发表拿 accession,不是把它当成已有数据层。 在拿到之前,Limitations 里「无内皮/周细胞」这句话照写不误。
十四、6-3 跨组学定向验证(跑前记录 · 尚未运行)
★ 范围:11 个蛋白(2026-08-08 用户改定,原定只做
IFNAR1)。 理由与 6-4a 相同 —— 这是因果层(MR + 共定位),单个蛋白的阴性无法解释。
14.1 这一层要回答什么
第 2 步的暴露是血浆蛋白水平。6-3 换一个分子层级问同一个因果问题:
同一个基因,转录本水平(eQTL)、剪接水平(sQTL)的遗传预测值, 是否也与糖尿病并发症相关,方向是否一致?
★ 它是加分项,不是否决点(v4 §六)。跨层方向相反不构成反对 —— 实测跨层往往不是同一个变异(旧仓 cross_layer_ld.csv 实测 r² 仅 0.095 / 0.031)。
14.2 数据(★ 2026-08-08 已实测核实版本)
| 层 | 用什么 | 版本判定 |
|---|---|---|
| 血液 eQTL | eQTLGen Phase 1 cis-eQTL-FDR0.05.txt.gz(31,684 人) | ❌ Phase 2(43,301)汇总统计不可下载。官网 cis-eqtls.html 只挂 2019-12-11 的 Phase 1 文件;Phase II 站点 Resources 页明写只提供 Phase I → 只能用 Phase 1,写进 Limitations |
| 多组织 sQTL | GTEx v10 GTEx_Analysis_v10_sQTL.tar(1.96 GB) | ✅ 可公开下载,非 requester-pays:gs://adult-gtex/bulk-qtl/v10/single-tissue-cis-qtl/。★ 应从现用的 v8 升级 |
| (备选)多组织 eQTL | GTEx_Analysis_v10_eQTL.tar(2.56 GB) | ✅ 同上;另有 susie-qtl/ 精细定位可用于共定位 |
requester-pays 只限 all_associations(全量名义 p:eQTL 260 GB / sQTL 860 GB)。 ⚠️ 这会限制共定位 —— 开放桶里是显著对,做 coloc 需要区域全量。 处理办法见 14.5「已知障碍 B」。
⚠️ 开跑前必做的数据搬迁:eQTLGen Phase 1 现在在 WSL /home/research/multiomics-mr/data/eqtlgen/(cis-eQTL-FDR0.05.txt.gz 322 MB、 SNP_AF.txt.gz 240 MB),不在 D 盘,违反「原始数据一律存 D 盘」→ 先移到 D:\mrdata\eqtlgen_phase1\ 并校验 md5。
14.3 ★ 预期(写于运行之前,与实际不符一律先按缺陷处理)
| # | 预期 | 依据 |
|---|---|---|
| E1 | 11 个基因里 8–11 个在 eQTLGen 有 cis-eQTL | 旧仓实测 IFNAR1 799 条、ERMAP 1354 条;eQTLGen 覆盖 88.6% 基因 |
| E2 | ★ LACTB2 很可能不可评估 | 其哨兵 MAF=0.0071,eQTLGen 用的是常见变异 |
| E3 | ★ WARS 必须用 WARS1 查 | HGNC 已改名;这是本轮踩过的坑 |
| E4 | 跨层方向至少有 2–3 个相反 | 旧仓 ERMAP 血浆 vs 视网膜已相反;跨层非同一变异 |
| E5 | IFNAR1 血液 eQTL 层有信号(与视网膜阴性形成对照) | 旧仓 07_r9dm 已跑过 IFNAR1,可直接对账 |
| E6 | GTEx v10 里没有视网膜组织 | GTEx 不含眼组织 —— sQTL 层只能用血/肾等替代组织,这是结构性限制不是缺陷 |
★ E6 要特别小心:别把「GTEx 无视网膜」写成阴性结果,那是数据集本身不含该组织。
14.4 验收断言(抓不到即 stop(),且 ★ 断言未全过不得写正式产物)
| ID | 断言 |
|---|---|
N1 | 输入集恰 11 个,且与 coloc_verdict.csv 的 strong 非 MHC 集合完全相同 |
N2 | ★ 跑前跑后 retained 集合完全相同(6-3 不得剔除任何候选) |
N3 | 产物全 ASCII,无中文取值、无内部流程编号(★ 旧仓 07_r9dm 的 可分析/糖尿病性黄斑病变 是反面教材) |
N4 | 外部依赖失败必须 stop() —— 文件读不到 / 解压失败 / 0 行,一律不许当成「没有 eQTL」 |
X1 | ★ 符号解析必须显式记账:每个基因输出 lookup_status ∈ {queried, symbol_unresolved},WARS→WARS1 的改名要留痕 |
X2 | ★ 覆盖表必须区分三态:no_cis_eqtl(查了没有)/ not_queried(没查)/ below_maf_threshold(低频不可评估) |
X3 | 参照等位必须与血浆层对齐后才比较方向 —— 未对齐的方向比较一律作废 |
X4 | GTEx 版本号写进产物元数据(v10),且断言不是 v8 |
X5 | Steiger filtering 的 units 与 prevalence 必须齐全,缺失时显式跳过并记录,不得静默通过 |
X6 | ★ 共定位若因「只有显著对、无区域全量」做不了,产物必须写 coloc_status = not_assessable_no_full_sumstats,不得留空或写 0 |
14.5 ★ 已知障碍(现在就写下来,别到跑的时候才发现)
A. eQTLGen Phase 1 手头这份是 FDR<0.05 的显著对,做不了共定位 ★ full 版已实测可下(2026-08-08 curl -I 验证):
https://download.gcc.rug.nl/downloads/eqtlgen/cis-eqtl/
2019-12-11-cis-eQTLsFDR-ProbeLevel-CohortInfoRemoved-BonferroniAdded.txt.gz
HTTP 200 · content-length 3,880,187,026(3.88 GB)· last-modified 2020-09-11→ 要做血液层共定位就必须补下这 3.88 GB。
B. GTEx 开放桶只有显著对,但 SuSiE 精细定位也开放 全量 all_associations 在 requester-pays 桶(eQTL 260 GB / sQTL 860 GB),不下。 ★ 实测 susie-qtl/ 在开放桶且很小:
gs://adult-gtex/bulk-qtl/v10/susie-qtl/
GTEx_v10_SuSiE_eQTL.tar 174 MB
GTEx_v10_SuSiE_sQTL.tar 169 MB→ 走 SuSiE 路线:用可信集(credible set)做共定位,不需要区域全量。 开跑前只需确认其列结构与组织清单。
C. GTEx 无眼组织 —— 见 E6,写 Limitations,不当阴性。
D. 旧仓 06_sqtl_smr 是孤儿层 —— 只有 6 条染色体、走的是 R13 线、无视网膜。 ★ 不得复用其结论,要跑就重跑。
14.6 ★ 下载完成后的实查(2026-08-08,开跑前必读)
已落盘并三重校验:
D:\mrdata\eqtlgen_phase1\
2019-12-11-cis-eQTLsFDR-ProbeLevel-...txt.gz 3,880,187,026 字节(= 源站 content-length)
md5 14a8bd98bcab93ed83f03aeb1add0bc4(与服务器端一致)
gzip -t 通过
cis-eQTL-FDR0.05.txt.gz 322 MB md5 3073e2f39d0847692c053949e85723d9
SNP_AF.txt.gz 240 MB md5 8ee9fb9476b4f8c6c97ceadc8b16733c下载走 myserver 中转:本地直连仅 113 KB/s(需 9.5 小时), myserver→荷兰源站 6.2 MB/s(4 分钟),回传 3.4 MB/s(19 分钟)。
两列结构上的硬约束(会直接改写脚本设计):
| # | 事实 | 后果 |
|---|---|---|
| ★ 1 | 文件只有 Zscore,没有 beta / se,也没有等位频率 | 必须用 SNP_AF.txt.gz 的 MAF 换算:beta = Z / sqrt(2p(1−p)(N+Z²))、se = 1 / sqrt(2p(1−p)(N+Z²))(Zhu 2016 SMR 近似)。换算失败的 SNP 必须显式剔除并计数,不得静默当 0 |
| ★★ 2 | SNPPos 是 GRCh37/hg19,我们主线是 GRCh38 | 坐标不能直接比。按 v4 §1「统一 build 必须最先」,本层要么按 rsID 匹配、要么先 liftover;两种都必须留下匹配率记账 |
→ 据此补两条断言:
| ID | 断言 |
|---|---|
X7 | ★ Z→beta/se 换算的输入(MAF、N)缺失时必须剔除并计数;剔除率 > 10% 即 stop() |
X8 | ★★ 产物元数据必须写明本层 build = GRCh37、对接方式 = rsID;各结局 rsID 命中率入表 |
★ X8 阈值已从计划里的「<80% 即 stop」改为「<50% stop,<80% 告警」 —— 写计划时我把它当成「同源数据的匹配率」,实际上这里是 eQTLGen(1000G 填补、HapMap 时代 rsID)与 FinnGen R9 的 rsID 交集, 天然达不到 80%,用 80% 会把正常情况判成失败。 改成双阈值:低于 50% 认定断链 stop(),50–80% 打印告警并要求在结果里说明。
14.4b ★ 写脚本时新增/偏离计划的地方(如实记账)
| # | 改动 | 原因 |
|---|---|---|
| 1 | 新增断言 O1:结局清单与病例数从 outcome_manifest.csv 读,不硬编码 | 主线 01_screen.R 早有此约定;50_retina_eqtl.R 因只跑两个眼病而硬编码,本层跑四个结局不宜照抄 |
| 2 | 结局范围 = 候选实际涉及的全部 4 个并发症(视网膜/黄斑/肾病/神经病变),不限于索引结局 | 全血是系统性组织,没有理由只对眼病测;且四个一起跑才有跨结局对照 |
| 3 | cis 窗口 ±500 kb(不是 eQTLGen 自带的 ±1 Mb) | 与主线第 1 步口径一致 |
| 4 | ★ 明确不使用 mr_funcs.R 里的 OUTCOMES 常量 | 那是 R13 线的结局表,R13 暂停中 |
14.7 执行顺序(每步停下等点头)
① 数据搬迁 + md5:eQTLGen Phase 1 → D:\mrdata\eqtlgen_phase1\
② 决定是否补下 eQTLGen full 版(取决于要不要做血液层共定位)
③ 下 GTEx v10 sQTL(1.96GB) + eQTL(2.56GB) + 探 susie-qtl 格式 → D:\mrdata\gtex_v10\
④ 写 51_crossomics_eqtl.R:11 个基因 × 血液 eQTL → MR + 共定位(若拿到 full)
⑤ 写 52_crossomics_sqtl.R:11 个基因 × GTEx v10 sQTL
⑥ 汇总跨层方向表(★ 必须先对齐参照等位)十五、6-3 运行结果(2026-08-08)
脚本 analysis/51_crossomics_eqtl.R,10.2 分钟,19 条断言全部 PASS。 产物 results/crossomics_eqtl/:crossomics_eqtl_mr.csv(40 行)· crossomics_eqtl_mr_all_methods.csv · blood_egene_coverage.csv · crossomics_eqtl_metadata.csv。
15.1 ★★ 跑出来的第一个真缺陷:SIGLEC5 被误判成「无全血 eQTL」
首跑按 druggability.csv 的 ENSG 查,SIGLEC5 = ENSG00000268500,抽到 0 行, 被 X2 记成 no_cis_eqtl。按基因符号复查:
eQTLGen 里 SIGLEC5 = ENSG00000105501 9,235 条 cis 记录 / 1,058 条显著**是注释版本不一致,不是生物学上不表达。**11 个基因逐个按符号交叉核过, 只有 SIGLEC5 一个受影响,其余 10 个双键一致。
★ 教训:
X2三态记账能分「没查 vs 查了 0 条」,但分不出「用错键去查」。 「查询覆盖」这件事有两层,键本身对不对是更外面的一层。
修法:ENSG + 基因符号(含 WARS1→WARS 别名)双键并集抽取,并新增三条断言:
| ID | 断言 |
|---|---|
X9a | 抽取必须用双键(键数 ≥ 2× 靶点数) |
X9b | 抽出的每一行都必须能回填到靶点(未回填 > 0 即 stop(),防 awk 键与 R 键不一致) |
X9c | 源库 ENSG 与 druggability.csv 不一致的靶点必须显式列出 |
同时把「本层出不来」拆成两种,不再混为一谈:
有 cis 记录但无显著 eQTL→ 真阴性,可解释(本轮只有APOE)双键均查不到→ 不可评估,不是阴性(本轮为空)
15.2 我在同一天因为同一个原因崩了两次
data.table(key = ...) —— key 是 data.table() 的形参, 写成列名会被当成「设主键」,报 some columns are not in the data.table。
第一次在元数据表(改名 field),第二次在 keymap(改名 lookup_key)。 两次都崩在 fwrite 之前,没留下半成品 —— 这是「断言/产物顺序」那条纪律的直接收益。
★ 另外:parse() 语法检查查不出这类错(语法合法,是运行期参数名冲突)。
15.3 覆盖(10/11 可分析)
| 基因 | 源库 ENSG | 匹配键 | cis 被测 | FDR<0.05 | 工具(clump 后) |
|---|---|---|---|---|---|
PAM | ENSG00000145730 | ensembl | 7,098 | 5,032 | 50 |
SIGLEC5 | ENSG00000105501 | ★ symbol | 9,235 | 1,058 | 64 |
LACTB2 | ENSG00000147592 | ensembl | 6,043 | 1,444 | 21 |
WARS | ENSG00000140105 | ensembl | 7,046 | 1,399 | 55 |
ERMAP | ENSG00000164010 | ensembl | 5,698 | 1,354 | 38 |
IFNAR1 | ENSG00000142166 | ensembl | 6,043 | 799 | 39 |
NUDT5 | ENSG00000165609 | ensembl | 7,619 | 372 | 18 |
NOTCH2 | ENSG00000134250 | ensembl | 3,671 | 212 | 9 |
GALNT3 | ENSG00000115339 | ensembl | 5,898 | 210 | 2 |
ACRBP | ENSG00000111644 | ensembl | 6,667 | 95 | 4 |
APOE | ENSG00000130203 | ensembl | 7,076 | 0 | — |
流程量:抽 72,094 行 → cis ±500kb 35,070 → FDR<0.05 9,263 → 频率对不齐 0 条 → F≥10 全留(最小 F=18.2,中位 108.5)。 rsID 命中率最低 90.8%。
⚠️ GALNT3(2 条)、ACRBP(4 条)工具数过少,其结果只能当探索性看。
15.4 ★★ 三条发现
① NUDT5 × 视网膜病变 —— 唯一的双层共定位,且同向
| 层 | PP.H4 | 效应 | 方向 |
|---|---|---|---|
| 血浆蛋白(第 5 步) | 0.9773 | OR 3.080 | risk_increasing |
| 全血转录本(本层) | 0.9821 | b = +0.111(p=1.8e-4,FDR=8.9e-4) | risk_increasing |
11 个候选里只有它出现「两个独立分子层、同一结局、共定位都 >0.97、方向一致」。 生物学上也自洽:转录本↑ → 蛋白↑ → 视网膜病变风险↑。
② IFNAR1 × 神经病变 PP.H4 = 0.946(p=7.4e-5)—— 一个新结局上的强信号
★ 神经病变不是 IFNAR1 的索引结局(血浆层在该结局上不显著)。 这是本层独立冒出来的,需要单独判断是发现还是巧合。
③ IFNAR1 在自己的索引结局(黄斑)上:PP.H4 仅 0.235,且方向与血浆相反
按既定规则,PP.H4 低 = 不是同一个变异 → 「方向相反」不构成反对。 WARS 亦然(视网膜 0.716 / 神经 0.722,均未达 0.8)。
④ 被救回来的 SIGLEC5 没有改变主结论:黄斑 p=1.1e-3、视网膜 p=1.7e-3, 方向都与血浆一致,但 PP.H4 只有 0.026 / 0.015 —— MR 显著而不共定位, 提示是 LD 混杂而非共享因果变异。它的加入只让 FDR 分母从 9 变 10,改动很小。
15.5 必写进 Limitations
- eQTLGen Phase 2(43,301 人)已发表但汇总统计未公开,本层用 Phase 1(31,684 人)
- 暴露侧 build 为 GRCh37;cis 窗口在其自身坐标内计算,与结局一律按 rsID 对接,全程不跨 build
beta/se由Zscore+ 等位频率换算(SMR 近似),不是原始回归系数- 跨层往往不是同一变异(本课题实测 r² 仅 0.095 / 0.031),方向相反不构成反对
十六、6-4b 单细胞表达定位运行结果(2026-08-08)
脚本 analysis/52_singlecell_expression.py,11.0 分钟,14 条断言全部 PASS。 产物 results/singlecell/:singlecell_expression.csv(341 行 = 11 基因 × 31 类)· singlecell_gene_resolution.csv · singlecell_metadata.csv。
图谱:hrca_snrna_allcells.h5ad,3,177,310 细胞 × 35,475 基因 / 31 类, 非零 8,248,015,746(每行均 2,596)。用的是 X(归一化值),不是原始计数。
16.1 ★ SIGLEC5 第二次踩同一个坑,被双键拦下
图谱 var 里没有 ENSG00000268500(druggability.csv 给的), 有的是 ENSG00000105501。若只用主 ENSG,它会第二次被静默报成「不表达」。
resolved_by = ensembl_alternate 已逐行记进 singlecell_gene_resolution.csv。 → 跨数据集的基因标识必须多路回退 + 逐个留痕,这已是同一天第二次。
16.2 ★★ AUC 首跑全错:平局按行序破了
首跑用 argsort 给 1..n 的唯一秩,等于对并列的 0 值按行号排先后。 而本图谱的细胞按类型排序存放:
| 细胞类型 | n | 行号相对位置 |
|---|---|---|
| retinal pigment epithelial cell | 863 | [1.000, 1.000](全在文件最末) |
| microglial cell | 4,894 | [0.998, 1.000] |
| Mueller cell | 221,612 | [0.924, 0.994] |
| retinal rod cell | 1,066,056 | [0.523, 0.859] |
于是排在末尾的细胞类型凭空拿到高秩。露馅的数字: SIGLEC5 × RPE 表达比例 0.0% 却得到 AUC 0.996。
- 受影响:只有
auc_vs_rest。pct_expressing/mean_expression不用秩,不受影响 - 修法:
scipy.stats.rankdata(method="average") - 旧产物改名
.badauc-171100留档,未删
新增断言 SC5:
一个细胞都不表达的类型,AUC 不可能 > 0.5 (平均秩下全零组 AUC = 0.5 × 对照组中同为 0 的比例 ≤ 0.5)
★ 这条若首跑就写,当场就会被拦下。「不可能的数字」值得写成断言,而不是靠眼看。
16.3 结果(修正后)
| 基因 | AUC 最高的细胞类型 | AUC | 表达比例 | 次高 |
|---|---|---|---|---|
PAM | GABA 能无长突细胞 | 0.866 | 89.4% | OFF 伞状神经节 0.837(97.8%) |
APOE | Müller 细胞 | 0.794 | 66.9% | 小胶质 0.697(51.5%) |
NOTCH2 | Müller 细胞 | 0.767 | 61.2% | 星形胶质 0.746(60.7%) |
NUDT5 | OFF 伞状神经节 | 0.672 | 59.7% | ON 中央凹侏儒神经节 0.623 |
ERMAP | 视锥细胞 | 0.661 | 35.5% | S 视锥 0.575 |
GALNT3 | Müller 细胞 | 0.651 | 31.1% | 星形胶质 0.539 |
WARS | OFF 伞状神经节 | 0.591 | 33.5% | 视锥 0.558 |
IFNAR1 | OFF 伞状神经节 | 0.560 | 36.2% | H1 水平细胞 0.553 |
LACTB2 | OFF 伞状神经节 | 0.532 | 13.8% | S 视锥 0.523 |
SIGLEC5 | 小胶质 | 0.516 | 3.7% | 侵入型侏儒双极 0.504 |
ACRBP | 弥散双极 3b | 0.505 | 1.6% | — 无富集 |
16.4 与旧版(scRNA 子集 265,767 细胞 / 18 类)的差异
| 基因 | 旧值 | 新值 | 说明 |
|---|---|---|---|
IFNAR1 | 伞状神经节 74.4% / AUC 0.757 | 36.2% / 0.560 | 细胞数 12 倍、snRNA vs scRNA,数值变化在预期内 |
ERMAP | AUC 0.536 = 无富集 | 视锥 0.661 | ★ 全集下有轻度视锥富集,旧结论需更新 |
★ ERMAP 的视锥定位与它 黄斑视网膜 eQTL 共定位 0.949 在解剖上自洽(视锥密集于黄斑)。
16.5 必写进 Limitations
- ⚠️ 本图谱无内皮细胞、无周细胞(
majorclass仅 10 类),而糖网是血管病 —— 换数据集也补不上(Voigt 2026 有,但数据未公开,见 §13) - 本层是描述性表达定位,不是因果证据;不得据此排除或提升任何候选
- AUC 是「该类 vs 其余全部」的判别力,受该类细胞数影响;
ACRBP/SIGLEC5整体表达极低, 其 AUC 接近 0.5 应理解为信息量不足,不是「确证无富集」
十七、6-5 组织特异 eQTL + 6-6 sQTL(跑前记录 · 尚未运行)
2026-08-08 用户拍板方案 A:做 6-5 与 6-6,mQTL 层砍掉(理由见 17.2)。 代谢组 / 脂质组不做(范畴不同,见 17.1)。
17.1 为什么还要加这一层:组织不对称,不是缺组学门类
先厘清一个被混在一起的范畴问题:
| 类 | 是什么 | 归属 |
|---|---|---|
| A. 同一基因的其他分子层 | eQTL / sQTL / mQTL —— 问「这个基因在别的层面也有因果信号吗」 | ✅ 就是 6-3 / 6-4 |
| B. 代谢物 / 脂质 | 不是这 11 个基因,是下游表型或独立暴露 | ❌ 属于 MVMR / 两步中介,或另一篇 |
★ 已核:药靶多组学 MR 文献用的一律是 pQTL + eQTL + sQTL + mQTL,全部是同一基因 (GSTM4 偏头痛、 BTN3A2 肾结石、 CPXM1 骨质疏松、 胶质瘤脑多组学)。 没有一篇把代谢组当作多组学验证层 —— 代谢物回答不了「是不是这个基因」。
且糖网方向的代谢/脂质 MR 已饱和(179 脂质种 × DR、脂质组+炎症因子中介、 血液代谢物 × DR、肠菌-代谢物轴 × DN 均已发表),做了也是重复。
真正的缺口是组织:
| 结局 | 对应组织 eQTL | 状态 |
|---|---|---|
| 视网膜病变 / 黄斑病变 | 视网膜(Advani 2024,403 眼) | ✅ 6-4a 已做 |
| 肾病 | 肾组织 | ❌ 缺 |
| 神经病变 | 神经组织 | ❌ 缺 |
这个缺口卡着两条已有结果:
APOE索引结局是肾病(血浆 coloc 0.998),从未在肾组织验证IFNAR1全血 eQTL 的 PP.H4=0.946 出在神经病变(非索引结局,§15 标注「需单独判断是发现还是巧合」)—— 胫神经组织正好能判
17.2 数据源选型(★ 全部实测,非看介绍页)
选定:eQTL Catalogue,最新 r8 pre-release(2026-01), 收录 GTEx v8,统一重处理,GRCh38(r5 起用 1000G 30x on GRCh38,已无需 liftover)。
胫神经 nerve_tibial n=532 QTD000286.all.tsv.gz 3.7 GB + .tbi
肾皮质 kidney_cortex n= 73 QTD000261.all.tsv.gz 2.4 GB + .tbi
全血 blood n=670 QTD000356.all.tsv.gz 2.6 GB + .tbi
胫神经 sQTL leafcutter QTD000290.cc.tsv.gz 916 MB + .tbi
肾皮质 sQTL leafcutter QTD000265.cc.tsv.gz 44 MB + .tbi
全血 sQTL leafcutter QTD000360.cc.tsv.gz 575 MB + .tbi
+ 每个数据集的 .permuted.tsv.gz(约 2 MB,列出全部被测性状)
合计约 10.2 GBFTP 根:https://ftp.ebi.ac.uk/pub/databases/spot/eQTL/sumstats/QTS000015/<QTD*>/ 索引表:github.com/eQTL-Catalogue/eQTL-Catalogue-resources → tabix/tabix_ftp_paths.tsv
★★ 下面这段结论已被 §19.1 部分推翻,保留原文以便对照
当时只看了变异层的 p 分布就下了「不是显著性过滤」的判断。正确说法是: 变异层没筛(零分布完整,故共定位有效),但性状层就是按显著性筛的(.cc 只含 FDR 显著的分子性状,保留比例 0.58%–9.9%)。以 §19.1 的表为准。
★ 实测澄清:.cc 不是显著性过滤。 取头部 3 MB 解压看 pvalue 分布:
| 文件 | n | min | median | max |
|---|---|---|---|---|
QTD000261.all(肾 ge) | 19,998 | 3.95e-06 | 0.489 | 0.99999 |
QTD000261.cc(肾 ge) | 19,998 | 7.15e-12 | 0.438 | 0.99998 |
QTD000290.cc(神经 sQTL) | 19,998 | 3.05e-18 | 0.466 | 0.999999 |
→ .cc 是按分子性状筛(只保留有信号的 trait,但该 trait cis 窗内的变异是全的), 不是按变异筛。因此 .cc 可以做共定位,但缺席某基因时含义是 「该性状无可信信号」而非「未测」,必须靠 .permuted 区分。
★ 被排除的三个选项(记账,免得以后重查)
| 数据源 | 排除理由(实测) |
|---|---|
| Susztak 肾 eQTL n=686 | 只有显著对(Kidney_eQTL_Meta_S686_Significant.q0.01.txt.gz,FDR<0.01),无全量 → 做不了 coloc。与 Voigt 2026 同一个坑 |
| GTEx 原始 tar(v10 开放桶) | 开放桶 2.56 GB 只有显著对;全量 all_associations 在 requester-pays(eQTL 260 GB / sQTL 860 GB) |
| GoDMC mQTL | ★★ README 逐字读:"every study performed a full analysis... returning only associations at a threshold of p<1e-5" → assoc_meta_all.csv.gz(5.9 GB)是 p<1e-5 候选表,不是全量,零分布被截断 → coloc.abf 不可用。且 build37、SNP 标识为 CHR:POS:SNP 非 rsID、n=27,750、450k+EPIC。只能做 MR,而单 cis 工具的 Wald p 恒等于结局 p → 加分近乎为零,砍掉 |
★ Susztak 仍留作后备:n=686 vs GTEx 肾 n=73,可做「该基因在肾里究竟有无 cis-eQTL」 的存在性检验(不需要全量)。是否启用见 17.4 的 E3。
17.3 ★ 预期(写于运行之前;与实际不符一律先按缺陷处理,不许现场解释)
| # | 预期 | 依据 |
|---|---|---|
E1 | 胫神经(n=532,与视网膜 Advani 403 同量级)11 个基因中 6–9 个达 FDR<0.05 显著 | 6-4a 视网膜实测「10 可分析 / 部分显著」的量级 |
E2 | ★★ 肾皮质(n=73)显著者 ≤3 个 | .cc 仅 50 MB vs 神经 1.2 GB(≈1:24),已提示有信号的性状极少 |
E3 | ★★ APOE 是肾层的功效对照 —— APOE 在肾脏高表达,若连它都测不出 cis-eQTL,即判定 n=73 完全无功效,该层所有阴性一律记 underpowered,不得写 no_effect | 「阴性只有靠对照才可解释」(第 9 节教训①) |
E4 | ★★ 本层最有价值的单条:IFNAR1 × 神经病变在胫神经是否复现 §15 的 PP.H4=0.946 | 非索引结局信号需独立层判真伪 |
E5 | GTEx 全血(n=670)是流水线正对照:应能复现 eQTLGen(n=31,684)已确认信号中的至少 1 个,但功效低得多 | 同组织不同队列;全 0 即接线错误而非生物学 |
E6 | ★ SIGLEC5 的 ENSG 又会不一致(druggability.csv 给 ENSG00000268500) | 同一坑第三次(eQTLGen …105501、HRCA …105501) |
E7 | rsID 命中率 >95% | 本层 GRCh38 且 1000G 30x,与 FinnGen R9 同 build(eQTLGen 因跨 build 只有 90.8%) |
E8 | sQTL 层部分基因缺席 .cc | .cc 按性状筛;须靠 .permuted 定三态 |
E9 | ★ GTEx 无眼组织 —— 眼病结局在本层只能用血/神经替代,这是结构性限制不是阴性结果 | 与 6-3 的 E6 同 |
E10 | LACTB2 大概率仍不可评估 | 哨兵 MAF=0.0071 |
17.4 验收断言(抓不到即 stop();★ 未全过只准写 .partial 且非零退出)
| ID | 断言 |
|---|---|
N1 | 输入集恰 11,且与 coloc_verdict.csv 的 strong 非 MHC 集合完全相同 |
N2 | ★ 跑前跑后 retained 集合完全相同(本层不得剔除任何候选) |
N3 | 产物全 ASCII,无中文取值、无内部流程编号 |
N4 | 外部依赖失败必须 stop() —— 下载失败 / tabix 报错 / 解压失败,一律不许当成「没有 eQTL」 |
K1 | ★ 基因标识 ENSG → 备用 ENSG → 符号 → 别名 四路回退,逐行写 resolved_by;解析失败必须显式列出(E6) |
K2 | ★★ 覆盖表四态:queried / no_signal_in_cc / not_tested(查 .permuted 确认)/ not_queried。四者不得合并 |
K3 | ★★ 肾层每条阴性必须带 negative_reason ∈ {no_cis_eqtl, underpowered_n73},由 E3 的 APOE 对照决定取值;不得单写 no_effect |
K4 | 元数据须写 exposure_build=GRCh38 且 == 结局 build;rsID 命中率入表,<80% 即 stop()(本层同 build,无 eQTLGen 那种天然折损,故不设双阈值) |
K5 | ★ 正对照断言:GTEx 全血层复现 eQTLGen 已显著基因数 ≥1;全 0 即 stop()(E5) |
K6 | 参照等位与血浆层对齐后才比较方向;未对齐的方向比较一律作废 |
K7 | 共定位输入来源入表:ge 层 coloc_input=all,sQTL 层 coloc_input=cc_subset |
K8 | tabix 切片行数 >0 才算 queried;0 行必须回查 .permuted 再定性(K2) |
K9 | ★ sQTL 一基因多内含子簇:簇层做 FDR,每基因取 PP.H4 最大簇并记录 n_clusters;不得只报最大值而不报簇数 |
K10 | 每个下载文件 md5 与服务器端一致、gzip -t 通过、字节数 == 源站 content-length |
K11 | ★ 三个组织同一套代码路径跑(血/肾/神经只换数据集参数),避免「组织间差异其实是脚本差异」 |
17.5 已知障碍(现在写下来,别到跑的时候才发现)
| # | 障碍 | 处置 |
|---|---|---|
| A | ★★ 肾皮质仅 n=73,GTEx 最小组织之一 | E3 的 APOE 功效对照 + K3 的 negative_reason;Limitations 必写 |
| B | ★ sQTL 只有 .cc 没有 .all | .permuted 定三态(K2/K8) |
| C | GTEx 无眼组织 | E9,写 Limitations,不当阴性 |
| D | ★ 两台服务器可用空间仅 14 GB / 16 GB,D 盘 6.1 TB 富余 | 分批下载→回传→删,不可一次拉完 10.2 GB |
| E | ★★ leafcutter 的 molecular_trait_id 是内含子簇(如 1:829104:841200:clu_46823_+),一基因多簇 | 按 gene_id 聚合;K9 定规则 |
| F | 旧仓 06_sqtl_smr 是孤儿层(6 条染色体 / R13 线 / 无视网膜) | ★ 不得复用其结论,重跑 |
17.6 执行顺序(★ 每步停下等点头)
① 分批下载 → D:\mrdata\eqtl_catalogue\ + md5 + gzip -t (走 myserver 中转)
② 格式实查:列结构 / SIGLEC5 解析 / 11 基因覆盖 → 回填本节
③ 写 53_tissue_eqtl.R(6-5:血·肾·神经 × 11 基因 → MR + coloc)
④ 跑 → 写结果记录
⑤ 写 54_sqtl.R(6-6:同三组织 leafcutter)
⑥ 跑 → 写结果记录
⑦ 汇总跨组织/跨层方向表(★ 必须先对齐参照等位)七要素归属:本层与 6-3 / 6-4a 同属因果层(MR + 共定位), Methods 里与 6-4b(描述性查表)分开归类;结果进补充材料, 仅 APOE/肾病与 IFNAR1/神经病变两条若成立则进正文。
17.7 ★ 数据落盘与格式实查(2026-08-08 19:25,写脚本之前)
落盘校验:18 个文件、10.94 GB、22 分钟,fail=0,无 .BAD、无残留段。
D:\mrdata\eqtl_catalogue\ 每文件:字节数 == 源站 content-length ✓ gzip -t ✓
MD5SUMS.txt 已生成(★ EBI 不发布官方 md5,
故校验依据是「字节数 + gzip -t」,不是与源站校验和比对)下载方式:本地代理 10808 直连,400 MB 分段 × 10 路并发,聚合 10.2 MB/s。 ★ 实测三方对比:本地代理 1.06 MB/s、本地直连 0.57 MB/s、myserver→EBI 1.01 MB/s —— 瓶颈在 EBI 侧,中转无收益。与 eQTLGen 那次(荷兰源站直连仅 113 KB/s, 必须走中转)结论相反,不可照搬上次经验。 EBI 限单连接不限单客户端,故并行有效。
WSL 无 tabix/bgzip → 不安装,改用与 6-3 相同的 zcat | awk 流式过滤。
11 个基因在三组织的被测与显著情况(.permuted 的 p_beta)
| 基因 | 胫神经 n=532 | 肾皮质 n=73 | 全血 n=670 |
|---|---|---|---|
ACRBP | ★ 7.15e-04 | 0.252 | 0.0527 |
APOE | ★ 5.46e-17 | 0.397 | 0.567 |
ERMAP | 0.0568 | 0.273 | 0.705 |
GALNT3 | ★ 5.32e-41 | 0.108 | 0.410 |
IFNAR1 | ★ 1.34e-11 | 0.200 | ★ 0.0392 |
LACTB2 | ★ 1.32e-04 | 0.622 | ★ 1.54e-04 |
NOTCH2 | ★ 1.96e-02 | 0.796 | 0.967 |
NUDT5 | ★ 3.26e-09 | 0.931 | ★ 0.0492 |
PAM | 0.267 | 0.984 | ★ 9.45e-63 |
SIGLEC5 | ★ 3.05e-05 | 0.657 | ★ 2.85e-27 |
WARS | 0.264 | 0.543 | ★ 2.81e-26 |
| p_beta<0.05 | 8/11 | ★★ 0/11 | 6/11 |
| 被测性状总数 | 23,788 | 24,310 | 16,701 |
跑前预期对账(★ 逐条,含我预测错的)
| # | 预期 | 实际 | 判定 |
|---|---|---|---|
E1 | 胫神经 6–9 个显著 | 8/11 | ✅ 命中 |
E2 | 肾皮质显著者 ≤3 | 0/11 | ✅ 命中(比预期更极端) |
E3 | APOE 作肾层功效对照 | 肾 p_beta=0.397 不显著;同基因胫神经 5.46e-17 | ★★ 触发:判定肾皮质 n=73 完全无功效 |
E5 | 全血能复现 eQTLGen 已确认信号 ≥1 | 6/11(含 IFNAR1/NUDT5/WARS/PAM/SIGLEC5/LACTB2) | ✅ 正对照成立 |
E6 | ★ SIGLEC5 的 ENSG 会第三次不一致 | 错了 —— eQTL Catalogue 用的就是 ENSG00000268500,与 druggability.csv 一致,备用键未启用 | ❌ 我预测错 |
★ E6 记为预测失败。多路回退仍保留(成本近零的保险),但 「跨数据集 ENSG 必然不一致」这个推广是过头的 —— eQTLGen 与 HRCA 不一致, 不代表所有数据源都不一致。
★★ 由实查得出、会改变 6-5 交付物的结论
肾皮质层 0/11 显著,且 APOE 功效对照未通过。 注意肾皮质被测性状 24,310 个,比胫神经的 23,788 还多 —— 不是「测得少」, 是纯粹没功效(n=73)。
→ 后果:在无边际 eQTL 信号的基因上做共定位近乎无信息(PP.H0/H3 主导)。 APOE × 肾病这一问,GTEx 肾皮质 n=73 回答不了。 按 E3 跑前写死的规则,该层全部结果只能记 underpowered,不得写 no_effect。
→ 胫神经层健康:8/11 显著,含 IFNAR1 p_beta=1.34e-11 → ★ E4(IFNAR1 × 神经病变那个 0.946)是可以回答的,这是本层的主要价值所在。
17.8 ★★ Susztak 肾存在性检验 —— 它推翻了我自己写的 E3
用户 2026-08-08 选方案 B:GTEx 肾照跑 + 补 Susztak n=686 做存在性检验。
数据:Kidney_eQTL_Meta_S686_Significant.q0.01.txt.gz 28,613,297 字节 · gzip -t 通过 · 1,179,179 条显著 cis 对(FDR<0.01) → D:\mrdata\kidney_eqtl_susztak\ 出处:Liu H, et al.(Susztak lab),4 项研究 meta,n=686。只有显著对,无全量 → 不能做 coloc。
★ 下载坑:页面给的
https://figshare.com/ndownloader/files/33957947恒返 HTTP 202(0 字节),加 UA+cookie 后变 403;两个不同出口 IP 均如此。 走 figshare API 的 302 才暴露真实主机ndownloader.figshare.com—— 换这个立刻 200,4.15 MB/s。 ⚠️ 首次那个 0 字节文件被gzip -t当场抓住。若无校验,它会被读成「肾里没有 eQTL 数据」。
11 个基因在肾组织(n=686)的显著 cis 对数
| 基因 | 显著对数 | 最小 P | 最强 SNP |
|---|---|---|---|
PAM | 165 | 4.471e-16 | rs385827 |
WARS | 131 | 1.669e-16 | rs941926 |
LACTB2 | 96 | 2.501e-05 | rs201659904 |
IFNAR1 | 93 | 6.881e-12 | rs2040109 |
GALNT3 | 3 | 8.001e-06 | rs1432275 |
NUDT5 | 2 | 1.842e-05 | rs10906083 |
SIGLEC5 | 1 | 2.188e-05 | rs8104955 |
ACRBP | 0 | — | 无显著 cis-eQTL |
APOE | ★★ 0 | — | 无显著 cis-eQTL |
ERMAP | 0 | — | 无显著 cis-eQTL |
NOTCH2 | 0 | — | 无显著 cis-eQTL |
★ APOE 的 0 已核实不是区域缺失:同一 chr19:44.5–45.5 Mb 窗口内 APOC2 155 条、TOMM40 2 条、NECTIN2 1 条,区域覆盖正常。
★★ 因此 E3 是选错了功效对照,必须更正
跑前我把 APOE 定为肾层的功效对照(「若连 APOE 都测不出,即判 n=73 无功效」)。 这个设计是错的 —— 现在知道 APOE 在肾里即便 n=686 也没有显著 cis-eQTL, 它在 GTEx n=73 的阴性与真实缺失完全一致,因此不能用来诊断功效。
正确的功效对照应是在 n=686 确有强信号、而在 n=73 落空的基因:
| 对照基因 | 肾 n=686 | GTEx 肾 n=73 (p_beta) | 读数 |
|---|---|---|---|
PAM | 165 条,p=4.47e-16 | 0.984 | 落空 |
WARS | 131 条,p=1.67e-16 | 0.543 | 落空 |
IFNAR1 | 93 条,p=6.88e-12 | 0.200 | 落空 |
→ 三个在 n=686 极显著的基因在 n=73 全部落空 —— 这才是 GTEx 肾皮质 n=73 功效不足的证据,结论与原先一致,但依据换了。
对 APOE × 肾病的实质结论
血浆层 APOE × 肾病共定位 0.998,但肾组织层没有可检出的 cis-eQTL(n=686, FDR<0.01) → 该血浆信号没有肾脏转录本层面的对应物。
⚠️ 措辞边界:Susztak 只发布显著对,故「FDR<0.01 无显著对」不等于「不存在 eQTL」。 但 n=686 且邻近基因均被检出,该阴性可写,措辞用 「未检出显著 cis-eQTL」而非「无 eQTL」。
17.9 ★ 另两个写脚本前查实的列语义坑
| # | 事实(实测) | 后果 |
|---|---|---|
| ★★ 1 | maf 不是效应等位频率。ac/an 才是 alt(=效应)等位频率,可 >0.5;maf = min(ac/an, 1−ac/an)。实测 199,999 行中 41,554 行(21%) 两者不等,恰为 ac/an>0.5 者(例:ac/an=0.80137 vs maf=0.19863) | 必须用 ac/an 算 eaf。用 maf 会让 21% 的变异频率翻转,MR 与共定位静默偏掉 |
| 2 | .all 自带 beta / se / rsid,rsid 缺失 0 条 | 不需要 Z→beta 换算(与 6-3 的 eQTLGen 不同);按 rsid 对接结局 |
基因坐标(.all 无基因位置列)取自 Ensembl REST,GRCh38,写死进脚本并加断言交叉核对:
ACRBP chr12 6,638,075-6,647,441 (-) APOE chr19 44,903,787-44,909,396 (+)
ERMAP chr1 42,817,078-42,845,920 (+) GALNT3 chr2 165,747,339-165,846,201(-)
IFNAR1 chr21 33,324,387-33,359,864 (+) LACTB2 chr8 70,619,104-70,669,306 (-)
NOTCH2 chr1 119,909,256-120,100,779(-) NUDT5 chr10 12,164,844-12,196,155 (-)
PAM chr5 102,753,981-103,031,149(+) SIGLEC5 chr19 51,610,948-51,630,401 (-)
WARS1 chr14 100,333,782-100,376,805(-)★ ENSG00000140105 的 Ensembl display_name 是 WARS1,再次印证改名。
17.10 ★ 写 53_tissue_eqtl.R 时新增/偏离计划的地方(如实记账)
| # | 改动 | 原因 |
|---|---|---|
| 1 | ★★ E3 的功效对照由 APOE 换成「Susztak n=686 显著对 ≥50 条」的基因(实为 PAM 165 / WARS 131 / LACTB2 96 / IFNAR1 93) | 见 §17.8:APOE 在 n=686 也无显著 cis-eQTL,用它诊断功效是范畴错误 |
| 2 | ★ negative_reason 由两值扩为三值:no_cis_eqtl_even_at_n686 / underpowered_n73 / no_effect | 「肾里本来就没有」与「n=73 测不出」是不同结论,不能合并 |
| 3 | ★★ 共定位输入用全 cis 区(含不显著变异),而非 F≥10 工具集 | 50_retina_eqtl.R 当时只有显著对可用;本层有 .all 全量,给 coloc.abf 完整零分布才正确。产物列 coloc_input=all_nominal_cis 留痕 |
| 4 | ★★ eaf 用 ac/an 而非 maf | §17.9:21% 的行两者不等,用 maf 会静默翻转频率 |
| 5 | 基因坐标写死(Ensembl GRCh38),加断言 K12b:TSS 必须落在该基因变异跨度内且染色体一致 | 目录窗为 TSS±1 Mb,故 TSS 必在跨度内;写死使脚本不依赖网络可复现 |
| 6 | K4 只设单阈值,不像 6-3 那样用 50/80 双阈值 | 本层与结局同为 GRCh38 且同 1000G 30x,无跨 build 的天然折损。★ 阈值最初写 80%,实测命中率 78.7% 反而把断言打挂了 —— 根因是我的预期 E7 错了(见 §18.8),查清后阈值定为 70%,并在断言消息里写明实测基线 |
| 7 | 不做 Z→beta 换算 | .all 自带 beta/se/rsid(与 6-3 的 eQTLGen 不同) |
| 8 | 断言分 hard() / chk() 两级 | 输入集类断言(N1)立即 stop();结果类断言累积后只写 .partial 且非零退出 |
预期运行开销:流式过滤 9.3 GB(三个 .all)+ 3 组织 × 11 基因 × 4 结局 = 132 次共定位。 提取结果会缓存为 _raw_<dataset>_11genes.tsv,重跑可复用。
17.11 ★★ 首跑当场抓到的两个缺陷(2026-08-08 20:16,运行中发现)
缺陷 1:结局集漂移 —— 丢掉了 Neuropathy,而它正是本层最有价值的一问
首跑日志第 4 行:
[PASS] O1 结局清单来自 manifest,共 3 个:Retinopathy, Maculopathy, Nephropathy**少了 Neuropathy。**根因是我写成
r
DIS <- sort(unique(V[coloc_tier == "strong" & in_mhc == 0, disease])) # ← 错而 6-3 用的是 sort(unique(V$disease))。实测:
| 取法 | 得到的结局 |
|---|---|
coloc_verdict 全部 disease | Maculopathy, Nephropathy, Neuropathy, Retinopathy ← 6-3 用这个 |
按 coloc_tier=="strong" 过滤 | Maculopathy, Nephropathy, Retinopathy ← 我错用了这个 |
Neuropathy 只出现在非 strong 的候选对里,按 strong 过滤会把它整个丢掉 —— 于是 E4(IFNAR1 × 神经病变的 PP.H4=0.946 能否在胫神经复现)根本跑不到, 而胫神经恰恰是为神经病变才引入的组织。
修法:改回 sort(unique(V$disease)),并新增硬断言 O2: 结局集必须与 6-3 的产物完全一致,否则 stop()。防同类漂移再犯。
★ 教训:「输入集/结局集从哪来」本身就该有断言。
O1只断言了「来自 manifest 且数量对得上」,却没断言这个集合是对的—— 一个自洽但错误的集合照样能 PASS。
缺陷 2:缓存复用判据是「文件存在且 >1000 字节」
r
if (file.exists(raw) && file.size(raw) > 1000) return(raw) # ← 中途被杀会静默复用截断文件这正是本项目记录在案的老 bug「断点续跑:文件存在=已完成」。 修法:先写 .building 再 file.rename() 原子落地,只认落地后的名字。
★ 处置决策:不杀当前进程
当前进程正在做最贵的一步(流式过滤 9.3 GB 建三个缓存)。 此刻杀掉会留下截断的 _raw_*.tsv,而修复前的复用判据会静默采用它。 故:让它跑完 → 缓存完整落地 → 用修好的脚本热缓存重跑(提取阶段秒过)。 首跑那份 3 结局的产物会被重跑覆盖。
十八、6-5 运行结果(定版,2026-08-08 22:03–22:19,22 条断言全 PASS)
耗时 15.9 分(缓存热)。产物 results/tissue_eqtl/,无 .partial,[FAIL] 计数 0。
exposure eQTL Catalogue r8 pre-release / GTEx v8 / GRCh38
cis ±500 kb · eaf = ac/an · coloc 输入 = 全 cis 区名义统计(coloc_input=all_nominal_cis)
结局对接 rsid 精确匹配 · 命中 31,972 / 40,637 = 78.7%(四结局完全一致)
工具 cis 变异 108,180 → F>=10 工具 5,586(最小 F=10.0,中位 F=23.2)
产物规模 132 行 = 11 基因 × 4 结局 × 3 组织★ 本节数字是第三版(定版),前两版已改名留档、请勿引用
| 版本 | 时间 | 差别 | 为什么废弃 |
|---|---|---|---|
.partial-3outcomes-212431 | 20:42 | 只有 3 个结局 | 结局集漂移,丢了 Neuropathy(见 §17.11) |
.commasplit-220346 | 21:44 | 4 结局,但 rsID 拆逗号对接 | ★ 自创无先例的做法,已撤回(见 §18.6) |
| 无后缀(定版) | 22:19 | 4 结局 + 精确 rsID 匹配 | ✅ 以此为准 |
18.1 ★★ PP.H4 ≥ 0.5 的全部结果(只有 3 条,全在全血)
| 组织 | 基因 | 结局 | PP.H4 | nsnp | b | p | vs 血浆 |
|---|---|---|---|---|---|---|---|
| blood | WARS | Retinopathy | 0.9252 | 8 | −0.1222 | 2.59e-06 | opposite |
| blood | IFNAR1 | Neuropathy | 0.8421 | 2 | +0.4187 | 0.262 | 不可比 |
| blood | WARS | Neuropathy | 0.6601 | 8 | −0.2507 | 4.25e-10 | 不可比 |
**胫神经与肾皮质无一条 ≥0.5。**各自最高:
| 组织 | Top 3 |
|---|---|
| nerve_tibial | ACRBP×Reti 0.318 · GALNT3×Reti 0.318 · SIGLEC5×Neuro 0.272 |
| kidney_cortex | NOTCH2×Reti 0.084 · NOTCH2×Macu 0.075 · WARS×Reti 0.064 |
18.2 ★ E4 有了答案,而且是两半
跑前 E4 问:6-3 全血 eQTLGen 给出的 IFNAR1×神经病变 PP.H4=0.946 是发现还是巧合?
| 层 | 数据 | PP.H4 | 读数 |
|---|---|---|---|
| 6-3 全血 | eQTLGen n=31,684 | 0.946 | 原始发现 |
| 6-5 全血 | GTEx n=670(独立队列) | ★ 0.8421 | 复现了 |
| 6-5 胫神经 | GTEx n=532 | ★ 0.0584 | 不外推到神经组织 |
| 6-5 肾皮质 | GTEx n=73 | 0.0366 | 该层整体功效不足,不作解读 |
★★ 两条都重要:
- 不是 eQTLGen 的偶然 —— 换独立队列仍 0.84
- 但它是「血液层」现象,不是「神经组织」现象 —— 且这是有功效的阴性:
IFNAR1在胫神经是强 eGene(p_beta=1.34e-11),不是测不出
18.3 ★★★ NUDT5 的头号结论被这一层挑战
计划里写着「NUDT5 是唯一两个独立分子层、共定位均 >0.97、方向一致的候选」。
| 层 | 数据 | NUDT5 × Retinopathy PP.H4 |
|---|---|---|
| 血浆 pQTL | UKB-PPP | 0.9773 |
| 6-3 全血 | eQTLGen n=31,684 | 0.9821 |
| ★ 6-5 全血 | GTEx n=670 | ★★ 0.0328 |
| 6-5 胫神经 | GTEx n=532 | 0.2479 |
| 6-5 肾皮质 | GTEx n=73 | 0.0418 |
同为全血,两个独立数据集给出 0.982 与 0.033。
功效解释成立吗?部分成立,但不够干净:
GTEx 全血 p_beta | 全血 cis 变异 / p<0.05 | GTEx 全血 PP.H4 | |
|---|---|---|---|
WARS | 2.81e-26(强 eGene) | 3,609 / 798 | 0.925 |
IFNAR1 | 0.0392(弱 eGene) | 3,099 / 257 | 0.842 |
NUDT5 | 0.0492(弱 eGene) | 4,821 / 358 | 0.033 |
→ 样本量差 48 倍(670 vs 31,684),NUDT5 在 GTEx 全血只是勉强的 eGene, 功效不足是合理解释;但 IFNAR1 的 eGene 强度相当(0.0392 vs 0.0492)却出了 0.840, 所以不能把功效当成完整解释。
★ 写作必须处理这一条,不得回避:NUDT5 的「两层共定位」证据现在有一个 独立血液队列的不复现。要么补更有功效的血液 eQTL 队列,要么在正文明确写出这个不一致。
18.4 ★ APOE × 肾病:三层全阴,且原因已判定
| 组织 | PP.H4 | p_beta | negative_reason | Susztak n=686 显著对 |
|---|---|---|---|---|
| kidney_cortex | 0.0499 | 0.397 | ★ no_cis_eqtl_even_at_n686 | 0 条 |
| blood | 0.0295 | 0.567 | — | 0 |
| nerve_tibial | 0.0050 | 5.46e-17 | — | 0 |
→ 血浆层 APOE×肾病 coloc 0.998,但肾组织层在 n=686 都检不出 cis-eQTL。 措辞用「未检出显著 cis-eQTL」,不得写「无 eQTL」(Susztak 只发布显著对)。
★ 注意胫神经那一行:APOE 在胫神经是极强 eGene(p_beta=5.46e-17)却仍得 0.0050 —— 这是有功效的阴性,比肾层那条更有信息量。
18.5 MR FDR<0.05 的 4 条(★ 与共定位分开读)
有 MR 估计的行 112 / 132(其余为该基因在该组织无合格工具)。
| 组织 | 基因 | 结局 | nsnp | b | p | FDR | PP.H4 |
|---|---|---|---|---|---|---|---|
| blood | WARS | Neuropathy | 8 | −0.2507 | 4.25e-10 | 3.83e-09 | 0.660 |
| blood | WARS | Retinopathy | 8 | −0.1222 | 2.59e-06 | 2.33e-05 | 0.925 |
| nerve | NUDT5 | Retinopathy | — | −0.0983 | 5.77e-05 | 6.35e-04 | 0.248 |
| kidney | ACRBP | Retinopathy | — | −0.0529 | 1.45e-03 | 1.16e-02 | 0.010 |
⚠️ MR 显著而 PP.H4 低(如 ACRBP 0.010)提示 LD 混杂而非共享因果变异。
★ 与拆逗号版相比少了 2 条 —— 两条的原因完全不同,不能一起说
| 掉出的行 | 真实原因 |
|---|---|
ERMAP × Maculopathy(血) | ★ 不是变不显著,是没工具了。ERMAP 在全血过完 F≥10 + clump + Steiger 只剩 1 个工具,而那个工具的 FinnGen rsids 恰好是逗号打包的 → 精确匹配下 nsnp 为空,四个结局全部变成「MR 不可评估」 |
SIGLEC5 × Maculopathy(血) | MR 本身一个字没变(nsnp=15、b=0.103481、pval=0.0062859 逐字节相同),只是 ERMAP 退出后 BH 的秩由第 6 变第 5、检验数 132→131 → FDR 0.0314 → 0.0566 跨过 0.05 |
全表只有 ERMAP 一个基因掉了工具(有 MR 估计的行 116 → 112,正好是它的 4 个结局), 其余 128 行的 nsnp/b/pval 全部逐字节相同。
读法(与事实分开写):这条掉出去在证据上不亏 —— ERMAP 全血共定位 PP.H4 只有 0.060,且单工具 Wald 下 MR 的 p 恒等于结局 GWAS 的 p,那条「显著」本就没有共定位 支撑、也不携带独立信息。但撤回的理由是政策性的(口径一致 + 无先例),不是因为结果更好看 —— 写 Methods 时这两件事不能混为一谈。
18.6 ★★ rsID 对接口径:拆逗号已撤回,全流程统一为精确匹配
问题是怎么冒出来的:本层的 K4 断言(rsID 命中率)实测只有 78.7%,远低于我预期的
95%。查因发现 FinnGen R9 的
rsids列并非单一 rsID:
| 形态 | 占比 |
|---|---|
单个 rsNNN | 94.71% |
| 逗号分隔多 rsID(同一变异的 dbSNP 别名,几乎全为 indel) | 1.92% |
| 空值 | 6.77%(其中 chr23 空值率 100%) |
我一度在 53_tissue_eqtl.R 内覆盖 extract_finngen 做「拆逗号」,产出即 .commasplit-220346。该做法已撤回,四条依据:
| # | 依据 |
|---|---|
| 1 | 无任何参照实现这么做 —— 导师版 1_r9_mr_diabetes.Rmd、主线 R/adapters/read_tsv.R + R/05_harmonise.R、共享的 mr_funcs.R::extract_finngen,三处全是精确匹配 |
| 2 | 拆逗号是半吊子 —— 不治 6.77% 的空 rsID、不治整条 chrX、不治覆盖不到的低频变异;真正的标准做法是 chr:pos:等位 对接 |
| 3 | 被影响者几乎全是 indel,且暴露侧同一变异会挂两个 rsID(实测 1:42324996:C:CT 同挂 rs35963804 / rs397806081,beta 相同)→ 必须靠一个任意的 break 防止同一行被当成两个独立工具 |
| 4 | ★ 共定位侧实测零收益:132 行中 PP.H4 最大只动 0.0071(WARS×Reti 0.9181→0.9252),无一条跨过 0.8 判定线 |
⚠️ 我上次那句「实测零收益」范围说小了,在此更正:零收益只覆盖共定位 PP.H4。 MR 侧不是零 —— 它决定了
ERMAP在全血有没有工具(见 §18.5 的 tip)。 事实与读法要分开写,不能拿「结果差不多」当撤回的理由。
★ 现在全流程的口径(已逐脚本实查,不是凭印象)
| 层 | 取 FinnGen 的实现 | 匹配方式 |
|---|---|---|
| 第 1–5 步主线 | R/adapters/read_tsv.R 的 grep -Fw -f → R/05_harmonise.R 的 out_std[SNP %in% snps] | 精确 |
6-3 全血 eQTL 51_crossomics_eqtl.R | source(mr_funcs.R) → extract_finngen | 精确 |
6-4a 视网膜 50_retina_eqtl.R | 同上 | 精确 |
6-5 组织 eQTL 53_tissue_eqtl.R | 同上(覆盖已删,只留注释说明为何不做) | 精确 |
★ 两条路径都是精确匹配,但由两处不同实现达成:grep -Fw 会把逗号行从磁盘上捞出来, 随后的 %in% 再把它们丢掉 —— 净效果与 awk 的 ($5 in a) 相同。这一点已实测核对, 但它意味着「一致」是巧合式的一致而非设计出来的,第 7 步定稿时应统一到一个函数。
因此 6-3 / 6-4a 不需要重跑(实测:6-3 只差 1 个 SNP、p 值排名第 13,474; 6-4a 差 58 个、含排名第 7 的,但按统一政策属已声明的入选标准而非缺陷)。 ⚠️ 可选:6-4a 另跑一次拆逗号版作敏感性分析(不是修复),不阻塞后续。
★ 写 Methods 时必须把「工具按 rsID 与结局对接」写成预先设定的入选标准并给出流失数字 (STROBE-MR 6d "Explain how missing data were addressed";⚠️ 10a 原文是 "numbers of individuals",不覆盖变异层面,不得拿它当依据)。详见待办计划 §4.4b。
18.7 候选格局(本层之后)
(定版数值;各格取该基因在该组织四个结局中的最大 PP.H4)
| 蛋白 | 血浆 | 视网膜 | 全血(eQTLGen n=31,684) | ★全血(GTEx n=670) | 胫神经(n=532) | 肾(n=73) |
|---|---|---|---|---|---|---|
WARS | 0.956 | 0.809 | 0.716 | ★ 0.925 | 0.055 | 0.064 |
IFNAR1 | 0.940 | 阴性 | 0.946(神经) | ★ 0.842(神经) | 0.058 | 0.057 |
NUDT5 | 0.977 | 0.035 | 0.982 | ★ 0.033 | 0.248 | 0.042 |
ERMAP | 0.967 | ★ 0.949(黄斑,反向) | 0.010 | 0.060 | — | — |
APOE | 0.998(肾) | — | — | 0.032 | 0.014 | 0.063 |
★★ WARS 现在是层数最多的候选(血浆 + 视网膜 + 两个独立全血队列均 >0.7), 反而不是 NUDT5。这与第五版计划里「NUDT5 单独讲」的判断相冲突,需重新评估。
18.8 ★★ 本轮的三次预测失败(如实记账)
跑前记录(§17.3)写下的预期与实际不符时,一律先按缺陷处理,不许现场解释。本轮三次:
| 预期 | 实际 | 根因 |
|---|---|---|
E6 SIGLEC5 的 ENSG 会第三次跨数据集不一致 | 一致 | 把「跨数据集 ENSG 必然不一致」推广过头了 |
E7 rsID 命中率 >95%(同 build 同 1000G 30x) | ★ 78.7% | ★★ 默认了「同 build = 同变异集」。实为 GTEx 深度 WGS 能测到 MAF~1% 的变异,而 FinnGen R9 覆盖不到 —— build 一致只保证坐标可比,不保证变异集相同 |
| 拆逗号后命中率升到 83.9% | 81.1% | 把别名重复(同一变异的两个 rsID)当成了真实回收 |
★ 其中 E7 直接导致 K4 的阈值一开始写成 80% 而把自己的断言打挂(见 §17.10 第 6 行)。 教训:断言阈值应当来自实测基线,而不是来自我的先验估计。
十九、6-6 sQTL(跑前记录 · 尚未运行,2026-08-08 22:5x)
★ 本节全部写于写代码之前。预期与实际不符时一律先按缺陷处理,不许现场解释。
19.1 ★★ 写脚本前的数据实况核查 —— 推翻了两条计划前提
前提一(已推翻):「.cc 不是显著性过滤」
计划 §一之三与本页 §17.2 都写着:.cc 不是显著性过滤(p 中位数 0.44、max 0.9999), 是按分子性状筛。这句话对了一半、错了一半,实测三个组织:
| 组织 | 被检验的内含子簇 | .cc 保留 | 保留者 p_beta 上界 | 未保留者 p_beta 下界 | 例外 |
|---|---|---|---|---|---|
| 胫神经 n=532 | 47,482 | 4,716(9.9%) | 1.02e-03 | 1.07e-31 | 38 |
| 全血 n=670 | 36,038 | 2,878(8.0%) | 8.13e-04 | 5.60e-17 | 20 |
| 肾皮质 n=73 | 53,122 | 308(0.58%) | 5.76e-05 | 8.38e-11 | 1 |
正确的说法:
- ✅ 对:被保留的性状里,所有变异都在(含大量 p≈1 的空值)→ 零分布未被截断 → 共定位在这些性状上有效
- ❌ 错:性状的保留本身就是按显著性的 ——
.cc只含分子性状层 FDR 显著的簇 (sCluster),保留比例随功效变化(肾 n=73 只有 0.58%,神经 n=532 有 9.9%)
→ 后果:.cc 里没有的基因,本层做不了共定位,而且这不是"没测",是"测了不显著"。
前提二(已推翻):「靠 .permuted 做四态覆盖账」
.permuted 没有 gene_id 列(只有 molecular_trait_object_id / molecular_trait_id / n_traits / p_beta),因此无法直接回答「基因 X 的内含子被测过吗」。
19.2 ★★★ 差点造出一条假结果:两种基因映射法在 NUDT5 上打架
.permuted 里的 molecular_trait_id 形如 10:12167893:12169258:clu_23877_+,自带坐标 —— 于是很自然会想「用坐标重叠把簇映射到基因」。实测这么做会出事:
| 方法 | NUDT5 × 胫神经的结论 |
|---|---|
坐标重叠(内含子落在 NUDT5 基因体 10:12,164,844–12,196,155 内) | 5 个簇,最佳 p_beta = 2.92e-07 → 看着像「NUDT5 有很强的胫神经 sQTL」 |
★ gene_id 直连(数据生产方自己的注释) | .cc 里根本没有 NUDT5 |
查那个簇到底属于谁:
簇 clu_23877_+ 内含子 10:12167893:12169258
.cc 里的 gene_id = ENSG00000065665 → SEC61A2(与 NUDT5 在 chr10 上重叠的另一个基因)
NUDT5 自身 = ENSG00000165609★★ 若按坐标归因,就会凭空造出「NUDT5 有胫神经 sQTL」 —— 而 NUDT5 恰恰是 §18.3 里争议中心的那个蛋白,这条假阳性会直接写进正文。
定下的规则:
- 归因只用
gene_id(数据生产方的注释),坐标法一律不得用于归因 - 坐标法只用于覆盖账(判「该区间有没有簇被检验过」),且产出列必须带
_positional后缀,并由断言禁止其进入分析列 - ★ 这条限制要写进 Limitations:重叠基因的剪接归因依赖 leafcutter 的注释选择
19.3 ★★ 本层的可分析面比计划预想的小得多(实查,非估计)
按 gene_id 直连,11 个基因在三个组织 .cc 中的命中:
| 组织 | 可分析的基因 | 行数 / 簇 / 内含子 |
|---|---|---|
| 胫神经 n=532 | WARS · ERMAP · PAM | 23,436/1/3 · 6,223/1/1 · 23,158/2/3 |
| 全血 n=670 | WARS · ERMAP · PAM · ACRBP | 7,787/1/1 · 12,508/1/2 · 7,781/1/1 · 7,685/1/1 |
| 肾皮质 n=73 | WARS | 5,228/1/1 |
可分析的 基因×组织 组合共 8 个(满格是 33 个)。
★★★ 必须提前说明:6-6 回答不了 §18.3 的 NUDT5 争议
NUDT5 与 IFNAR1 在任何组织的 sQTL .cc 里都不存在 —— 这不是接线错,是它们没有显著的剪接 QTL。 所以本层无法为「NUDT5 全血 0.982 vs 0.033」这条不一致提供任何新证据。 跑之前就要把这个预期讲清楚,不能跑完再说。
★ 但本层能给一条有功效的阴性(这正是它的价值): IFNAR1 的内含子在三个组织都被检验过,p_beta 分别为 0.547 / 0.445 / 0.495 → 可以写「IFNAR1 的遗传信号在表达层而非剪接层」, 这是对 6-3/6-5「全血 eQTL 共定位 0.946 / 0.842」的一个机制层面的界定,不是重复。
19.4 ★ 预期(写于运行之前)
⚠️ 区分「已实查确定」与「预测」—— 只有后者才算预期命中/落空:
| # | 预期 | 性质 |
|---|---|---|
S1 | 可分析组合恰 8 个(见 §19.3 表) | 已实查确定,写成硬断言 Q5 |
S2 | NUDT5/IFNAR1 在三组织全部为「测了但不显著」 | 已实查确定,写成断言 Q7 |
S3 | WARS 至少在一个组织出 PP.H4 ≥ 0.5 | ★ 预测。理由:它 eQTL 层全血 0.925,且 sQTL p_beta=5.93e-92(全血)/1.28e-95(神经)极强 |
S4 | ERMAP 在神经或全血出 PP.H4 ≥ 0.5 | ★ 预测。其 sQTL p_beta 8.42e-17 / 6.69e-18,若成立会是本层第二条结论 |
S5 | 肾层(只有 WARS)出不了 ≥0.5 | ★ 预测,理由同 6-5:n=73 功效不足 |
S6 | rsID 命中率落在 70–85% | ★ 预测。6-5 实测 78.7%,本层变异集不同故给区间不给点值 |
S7 | 会出现「MR 显著但 PP.H4 低」的行 | ★ 预测。可分析组合只有 8 个 → FDR 分母极小,须按 LD 混杂读,不得当阳性 |
19.5 验收断言(抓不到即 stop();★ 未全过只准写 .partial 且非零退出)
沿用 6-5 的 N1–N4 / O1–O2 / K12a–K12b(同一批写死的 GRCh38 坐标),新增:
| 断言 | 内容 |
|---|---|
Q1 | ★★ 归因只用 gene_id:任何来自坐标法的列必须带 _positional 后缀,且不得出现在分析列中 |
Q2 | 覆盖表 33 行(11 基因 × 3 组织),四态齐全无 not_queried |
Q3 | ★ coloc_input 全部 = cc_subset(不是 6-5 的 all_nominal_cis),元数据同步记 |
Q4 | ★ 分析单元留痕:每行必须带 n_clusters / n_introns / lead_intron |
Q5 | ★★ 硬断言:可分析的 基因×组织 组合恰为 §19.3 那 8 个,接线错立刻暴露 |
Q6 | rsID 命中率 ≥ 70%(沿用 6-5 的实测基线阈值,不再凭先验拍) |
Q7 | ★★ NUDT5/IFNAR1 全部组织的 coverage_status 必须是 no_significant_scluster,negative_reason = tested_not_significant —— 防止被静默记成 not_tested |
Q8 | 每个 基因×组织×结局 唯一一行 |
19.6 与 6-5 的实现差异(三处)
| # | 差异 | 处理 |
|---|---|---|
| 1 | 只有 .cc 没有 .all | 共定位输入 = cc_subset;缺席基因靠 .permuted + 坐标法记「测了不显著」 |
| 2 | 分析单元是内含子不是基因 | MR/coloc 按内含子跑,按 gene_id 聚合取 PP.H4 最大者,留痕 n_clusters/n_introns/lead_intron;FDR 在内含子层做 |
| 3 | 基因归因 | ★ 只用 gene_id,坐标法列带 _positional 后缀且被 Q1 隔离 |
⚠️ 旧仓 multiomics-mr/06_sqtl_smr 是孤儿层(6 条染色体 / R13 线 / 无视网膜),不得复用其结论。
二十、6-6 sQTL 运行结果(2026-08-08 23:06–23:17,23 条断言全 PASS)
耗时 10.9 分。产物 results/sqtl/,无 .partial,[FAIL] 计数 0。
exposure eQTL Catalogue r8 pre-release / GTEx v8 leafcutter sQTL / GRCh38
分析单元 内含子(molecular_trait_id)→ 按 gene_id 聚合,每基因取 PP.H4 最大的内含子
归因 ★ 只用 gene_id;坐标法列带 _positional 后缀且被断言 Q1 隔离
共定位输入 cc_subset(本层无 .all)· eaf = ac/an · 结局按 rsid 精确匹配
命中率 78.2%(四结局一致;6-5 基线 78.7%)
产物规模 32 行 = 8 个可分析组合 × 4 结局(覆盖表 33 格)20.1 ★★★ 最重要的一条:WARS 的剪接信号三个组织全部共定位
| 结局 | 胫神经 n=532 | 肾皮质 n=73 | 全血 n=670 |
|---|---|---|---|
| Retinopathy | ★★ 0.940 | 0.767 | 0.714 |
| Neuropathy | 0.722 | 0.714 | 0.710 |
| Nephropathy | 0.357 | 0.343 | 0.349 |
| Maculopathy | 0.274 | 0.135 | 0.121 |
★ 三个组织命中的是同一个剪接事件:主导内含子的供体位点完全相同 (14:100369258),受体位点略有差异(…:100374136 神经 / …:100375283 肾 / …:100376260 全血)。leafcutter 的簇编号按数据集各自命名故不同 (clu_59946_- / clu_70469_- / clu_44468_-),不是三个独立事件。
★★ 方向:剪接层与血浆层同向,而表达层反向。WARS×Retinopathy 血浆 OR 1.385(b>0),本层三组织 b 均为 +(direction_vs_plasma = same); 而 6-5 全血 eQTL 层是 b=−0.122(opposite)。 ⚠️ 但跨层未必是同一变异,这条只能作为一致性的加分,不能反过来说表达层错了。
20.2 ★★ 两条预测落空(如实记账)
| 预期 | 实际 | 根因 |
|---|---|---|
S4 ERMAP 在神经或全血出 PP.H4 ≥ 0.5 | ❌ 最高 0.374(神经×黄斑) | 我把「sQTL p_beta 极强(8.4e-17)」当成了「共定位会强」。分子层信号强 ≠ 与疾病共享因果变异 —— 这正是共定位要回答的问题,用 p_beta 去预测它属于推理短路 |
S5 肾层(n=73)出不了 ≥0.5 | ❌ 0.767(WARS×Reti) | ★★ 我把 6-5 得到的「肾 n=73 功效不足」笼统外推到了本层。实为:n 小是平均意义上的功效不足,单个极强的分子性状照样能共定位(WARS 肾 sQTL p_beta=1.3e-11)。「该层功效不足」不能用来预判该层里最强的那个信号 |
命中的三条:S3(WARS 至少一组织 ≥0.5 —— 三个都过)· S6(命中率 70–85%,实得 78.2%)· S7(存在 MR 显著而 PP.H4 低的行 —— ERMAP×Reti fdr 0.023/PP.H4 0.020、 PAM×Neph fdr 0.044/PP.H4 0.034)。 S1/S2 属跑前实查确定,由硬断言 Q5/Q7 守住。
20.3 ★ IFNAR1 与 NUDT5:有功效的阴性,不是没数据
| 基因 | 胫神经 | 全血 | 肾皮质 |
|---|---|---|---|
IFNAR1 | 3 簇被检验,最佳 p_beta 0.547 | 1 簇,0.445 | 3 簇,0.495 |
NUDT5 | 5 簇,2.92e-07(★ 见下) | 4 簇,0.23 | 4 簇,0.0496 |
三个组织全部记 coverage_status = no_significant_scluster、 negative_reason = tested_not_significant(断言 Q7 守住,不得写成「未测」)。
可写的结论:IFNAR1 在 6-3/6-5 有明确的表达层共定位(全血 eQTLGen 0.946 / GTEx 0.842),而剪接层三个组织全阴 → 其遗传信号定位在表达调控而非剪接。 这是对前两层的机制层面界定,不是重复。
⚠️
NUDT5胫神经那个 2.92e-07 是_positional列(坐标重叠),不是它的 —— 该簇的gene_id是SEC61A2。详见 §19.2。产物里这个数只出现在覆盖表, 断言Q1保证它进不了分析表。
★ NUDT5 的 0.982 vs 0.033 之争,本层给不出任何新证据(跑前 §19.3 已声明)。
20.4 MR FDR<0.05 的 14 条(★ 与共定位分开读)
前 6 条(按 FDR):
| 组织 | 基因 | 结局 | nsnp | b | FDR | PP.H4 |
|---|---|---|---|---|---|---|
| kidney | WARS | Retinopathy | 1 | +0.0899 | 1.81e-05 | 0.767 |
| nerve | WARS | Neuropathy | 15 | −0.0998 | 1.03e-04 | 0.722 |
| nerve | PAM | Retinopathy | 8 | −0.0645 | 5.29e-04 | 0.198 |
| kidney | WARS | Neuropathy | 1 | +0.1348 | 7.25e-04 | 0.714 |
| nerve | WARS | Retinopathy | 5 | +0.1242 | 2.44e-03 | 0.940 |
| kidney | WARS | Nephropathy | 1 | +0.0973 | 3.11e-03 | 0.343 |
⚠️ 肾层四格全是 nsnp=1 → 单工具 Wald,MR 的 p 恒等于结局 GWAS 的 p, 不携带独立信息。肾层能立住的是共定位(0.767 / 0.714),不是那个 MR p 值。
⚠️ MR 显著而 PP.H4 低的两条(ERMAP×Reti 0.020、PAM×Neph 0.034)按 LD 混杂读, 不得当阳性 —— S7 预判到了。
20.5 覆盖账(33 格)
| 状态 | 格数 | 含义 |
|---|---|---|
analyzable | 8 | 有 FDR 显著的剪接簇 → 做了 MR + 共定位 |
no_significant_scluster | 20 | ★ 测了但不显著(tested_not_significant) |
not_tested | 5 | 该基因无任何内含子簇被注释/检验:SIGLEC5 三组织全无、LACTB2 神经与肾无 |
★ SIGLEC5 在三个组织都没有内含子簇 —— 与它在 6-4a/6-5 屡次「非 eGene」一致, 写 Limitations 时应合并成一句:该蛋白在多个分子层均缺可用工具。
20.6 ★ 实现上值得记的两点
- 同一基因在不同结局可能选中不同内含子(
lead_intron列留痕)。 如WARS在胫神经:Reti/Macu 选14:100361921:100375283,Neuro/Neph 选14:100369258:100374136—— 同簇不同内含子。这是设计使然(按 PP.H4 取最大), 但报告时必须带lead_intron,否则读者无法复现是哪条。 .cc无.all对应物,故coloc_input记cc_subset(断言Q3)。 与 6-5 的all_nominal_cis不是一回事,跨层比较时不能混。
20.7 候选格局(六层之后)
| 蛋白 | 血浆 | 视网膜 | 全血 eQTLGen | 全血 GTEx | 胫神经 eQTL | 肾 eQTL | ★ sQTL(最高) |
|---|---|---|---|---|---|---|---|
WARS | 0.956 | 0.809 | 0.716 | 0.925 | 0.055 | 0.064 | ★★ 0.940(神经)三组织全 ≥0.71 |
IFNAR1 | 0.940 | 阴性 | 0.946 | 0.842 | 0.058 | 0.057 | 三组织全阴(有功效) |
NUDT5 | 0.977 | 0.035 | 0.982 | 0.033 | 0.248 | 0.042 | 无显著剪接簇 |
ERMAP | 0.967 | 0.949 | 0.010 | 0.060 | — | — | 0.374 |
APOE | 0.998(肾) | — | — | 0.032 | 0.014 | 0.063 | 无显著剪接簇 |
★★ WARS 的领先进一步扩大:血浆 + 视网膜 + 两个全血 eQTL 队列 + 三个组织的剪接层, 是唯一在两种分子机制(表达与剪接)上都有共定位的候选。
相关:共定位方法选型 · 第 5 步实录 · 分析流程定稿 v4 · 第 7 步定稿与作图计划