主题
方向感知药靶 PheWAS：获益 / 风险判读
这是 baseline（导师
Whole story line_v2.docx的 Fig5 与 Discussion 落脚点）有、我们原来唯一没有的一块。 做的过程中实测结果与原计划相反，产出定位因此改了，下面第一节先说清楚这件事。
产物：
| 文件 | 内容 |
|---|---|
analysis/37_drug_target_phewas.R | 分析脚本 |
results/drug_target/benefit_risk_matrix.csv | 长表，836 行逐关联判读 |
results/drug_target/benefit_risk_summary.csv | 逐蛋白汇总 |
results/drug_target/pleiotropy_spectrum.csv | 95 个蛋白的多效性谱系 |
results/drug_target/sparsity_check_summary.csv | 稀疏度确认（放宽阈值重查） |
figures/{pleiotropy_spectrum, benefit_risk_heatmap, benefit_risk_heatmap_MHC} | 三张图，PNG + 矢量 PDF |
一、★ 原计划落空：三个 A 级候选没有脱靶疾病性状
原计划是"三个 A 级候选各出一列，能读出该抑制还是该激动、代价是什么"。实测做不到，因为它们的哨兵 SNP 本身就不多效：
| 候选 | 哨兵 SNP | Bonferroni 显著关联 | 其中非分子性状 | 其中疾病性状 |
|---|---|---|---|---|
| ERMAP | rs11210710 | 6 | 0 | 0 |
| IFNAR1 | rs914142 | 4 | 1（Red cell distribution width） | 0 |
| APOL1 | rs136168 | 2 | 0 | 0 |
| （对照）APOE | rs429358 | 867 | 595 | 63 |
| （对照）AGER | rs204993 | 431 | — | 91 |
按原框架画出来会是三列空白。
不是查询阈值造成的假象
主分析的 PheWAS 查询阈值是 p < 1e-5（比 Bonferroni 阈值 1.06e-08 宽松两个多数量级）。为排除"漏查"，对这三个哨兵 SNP 放宽到 p < 1e-4 单独重查：
| 候选 | 重查返回 | Bonferroni 显著 | 非分子且显著 | 非分子最小 p |
|---|---|---|---|---|
| ERMAP | 10 | 6 | 0 | 2.1e-05 |
| IFNAR1 | 23 | 4 | 1 | 9.9e-09 |
| APOL1 | 7 | 2 | 0 | 9.96e-06 |
| APOE | 1,557 | 896 | 621 | 0 |
放宽后新增的关联全部落在 Bonferroni 阈值之外。稀疏是真实的多效性差异，与 APOE 差两个数量级。
因此改了定位
产出从「A 级候选的获益 / 风险矩阵」改为：
- 多效性谱系（
pleiotropy_spectrum）—— A 级三候选作为低多效性 = 干净靶点的正面安全信号呈现，与 APOE / MHC / ABO 的高多效对照； - 获益 / 风险矩阵（
benefit_risk_heatmap）—— 在确实有疾病性状的那批蛋白上做，方法完整，A 级三个作为空行保留在图上，让"干净"这件事看得见。
这是一个正面结论，不是回避：药靶 PheWAS 的本来用途就是查脱靶负担，查出来"没有"是对候选有利的证据，只是不能写成矩阵。
二、方法
干预方向
对每个蛋白取索引结局（优先并发症，同蛋白多结局取 FDR 最小者），用该结局的 MR 效应定干预方向：
b > 0（蛋白升高→风险升高）→ 策略 = 抑制b < 0（蛋白升高→风险降低）→ 策略 = 激动
95 个蛋白里 31 个索引到并发症，64 个索引到糖尿病本身。与 results/druggability.csv 已有的 mr_direction 交叉核对 7/7 一致（该文件只覆盖 7 个共定位支持的蛋白）。
判读规则
把 PheWAS 的 beta 对齐到 MR 工具的效应等位后，再除以 beta.exposure 换算成"每 +1 SD 蛋白"的尺度（与 b 同一尺度的 Wald 比）：
wald_trait = beta_phewas_aligned / beta.exposure
sign(wald_trait) == sign(b_index) → 干预该靶点也会降低该病风险 → 潜在额外获益
sign(wald_trait) != sign(b_index) → 干预会升高该病风险 → 潜在不良反应两者除以同一个 beta.exposure，所以符号比较不受工具方向影响。
三、三个坑的处理与实测
坑 1 等位方向对齐 —— 实测是空操作，但验证过了才敢这么说
PheWAS 返回的 beta 相对 OpenGWAS 自己的效应等位，必须先与 effect_allele.exposure 对齐。脚本实现了完整的翻转矩阵（同向 / 交换 / 换链 / 换链且交换 / 不匹配），回文位点（A/T、C/G）用 EAF 消歧。
实测结果：
- 93 / 93 个 SNP 两侧等位完全一致（独立于主脚本、直接从两个原始文件各取一次核对）
- 可对齐 15,572 条，无法对齐 217 条（回文且 EAF 缺失或过近 0.5，已剔除）
- 回文位点经 EAF 消歧后
flip = -1的仅 3 条 - 把 flip 全部强制为 +1（即完全忽略对齐），获益 / 风险判读改变 0 行
也就是说 OpenGWAS 与 UKB-PPP 都已对齐到同一参考等位，这一步在本数据上没有实际作用。但结论是"验证后为空操作"，不是"假设它对齐了" —— 换数据源就未必成立，脚本里的断言要保留。
对齐正确性的正面旁证：PheWAS 里命中我们自己结局（DM_MACULOPATHY 等）的 16 条关联，对齐后方向与 MR 估计 16 / 16 同号。
坑 2 性状方向的临床含义
"血细胞计数升高"没有好坏之分，只对疾病类性状做判读。分类办法：
| 类别 | 判据 | Bonferroni 显著条数（非糖尿病相关） |
|---|---|---|
| 分子性状 | ID 前缀 eqtl-/prot-/met- 或 性状名以 "levels" 结尾 或 是 ENSG 号 | 2,087 |
| 疾病（EFO） | dataset_domains 的 EFO 父类含 disease / disorder / cancer | 616 |
| 疾病（FinnGen 推定） | 域为 Unclassified 且 ID 为 finn-* | 220 |
| 定量或未分类 | 其余 | 5,777 |
★ 光按 ID 前缀滤分子性状会漏：
ebi-a-GCST900102xx这批「Galectin-4 levels」「E-selectin levels」是蛋白量性状但不带prot-前缀，所以加了性状名判据。
另外把靶标自身相关（性状名含 diabet，或 ID 含 DIAB/DM_/T1D/T2D）的 275 条单独标出——那不是脱靶，是内部一致性。
坑 3 MHC 区蛋白高多效
MHC 区（chr6:25–34 Mb）蛋白的关联密集但由长程 LD 驱动，没有靶点特异的解释力。正文热图排除，补充图单列。
★ 这里踩过一个静默丢数据的坑：第一版直接 merge
target_decision_table.csv的in_MHC列，但那张表只覆盖 31 个蛋白，其余 64 个的in_MHC是 NA，于是in_MHC == FALSE把它们整批悄悄滤掉了——正文热图只剩 8 个蛋白，且不报错。 改成按哨兵 SNP 坐标独立判定，并与分层表交叉核对（重叠的 29 个蛋白 29/29 一致）。 实际 MHC 蛋白是 20 个，不是分层表里的 11 个。
四、结果
多效性谱系
仅统计 Bonferroni 显著、非分子、非糖尿病相关的关联：
| 分组 | 蛋白数 | 脱靶关联中位数 | 脱靶疾病性状中位数 |
|---|---|---|---|
| A 级（优先） | 3 | 0（范围 0–1） | 0 |
| B 级（共定位支持） | 4 | 11.5（范围 0–535） | 1.5 |
| C 级（MHC / LD 混杂） | 11 | 156（范围 55–327） | 33 |
| D 级（仅 MR，非 MHC） | 11 | 14（范围 0–66） | 0 |
| 未分层（糖尿病索引） | 63 | 15.5–70 | 0–9 |
| 合计：MHC 区 | 20 | 106.5 | 24 |
| 合计：非 MHC | 72 | 12.5 | 0 |
获益 / 风险矩阵
- 836 行 / 49 个蛋白 / 255 个疾病性状
- 潜在额外获益 448 条，潜在不良反应 388 条
- 其中 593 行（20 个蛋白）来自 MHC 区，只进补充图；正文图为 29 个非 MHC 蛋白 + 补进来的 A 级 3 个空行
非 MHC 蛋白里净值最突出的：
| 净获益 | 净风险 | ||
|---|---|---|---|
| ABO | +29（33 获益 / 4 风险） | APOE | −57（3 / 60） |
| IL7R | +21（21 / 0） | IL10 | −21（0 / 21） |
| KLK1 | +16（16 / 0） | APOBR | −9（0 / 9） |
这些是假设性判读，不是药理结论：单哨兵 SNP 工具、未做多效性检验、FinnGen 推定疾病类未经人工核对。ABO / APOE 这类本身就是全基因组著名多效位点，其"净获益 / 净风险"更多反映位点多效而非靶点药理。
五、验收与自检
| 验收项 | 结果 |
|---|---|
| 单元自检：随机抽 3 条手工核对等位方向 | 脚本每次运行都打印 3 条完整明细（含两侧等位、align 判定、flip、原始与对齐后 beta），供人工复核 |
| 阳性对照：IFNAR1 已知药是干扰素激动剂，脱靶应见免疫 / 血液性状 | 部分成立——唯一的非分子显著关联是 Red cell distribution width（血液类），方向为负；免疫类性状未达 Bonferroni |
| 三个 A 级候选各出一列 | 未达成，见第一节：疾病类脱靶为 0，已改为以空行呈现 |
干预方向与 druggability.csv 一致 | 7 / 7 |
| MHC 标记与分层表一致 | 29 / 29 |
| 自结局方向旁证 | 16 / 16 同号 |
脚本内含 6 处 stopifnot 断言：单 SNP 工具、每 rsid 单组 ea/nea、干预方向与既有文件一致、MHC 标记无缺失、MHC 标记与分层表一致、汇总表行数。
六、局限
- 3 个蛋白的工具 SNP 根本不在 PheWAS 结果里：
ERN1、LACTB2、VWC2L。其中 LACTB2 与 VWC2L 属于 13 个并发症特异蛋白。它们的n_offtarget记为 NA 而非 0 —— "没测"与"测了没有"是两回事，谱系图里也不画它们。 - "无脱靶"受限于 OpenGWAS 的覆盖与功率。50,164 个数据集里罕见病、非欧人群性状覆盖有限；只能说"未发现证据"，不能说"已排除脱靶风险"。
- FinnGen 推定疾病类（220 条）未逐条人工核对，其中可能混入用药 / 手术编码类端点。已在
disease_class列单独标记，可随时剔除重算。 - 回文位点靠 EAF 消歧：矩阵中 38 行来自 4 个回文位点蛋白（CDSN 13、TRIM40 20、DDR1 3、CD164 2），已在
palindromic列标出。 - 单哨兵 SNP 工具，做不了 Cochran's Q / MR-Egger 截距，无法区分真脱靶与水平多效。这与全课题的同一条已知局限一致。
相关页面
- A 级候选证据链与两处硬伤
- baseline 文章逐条对比 · Yuan 2023
- 共定位证据分档与 SuSiE 降级
- R9 第二三轮核查（PheWAS 缓存指纹修复的来龙去脉）