主题
12 · 第 4.5 步 · 通路富集
状态：✅ 完成（§0–§4 写于运行之前，§5 为运行后回填） 权限：纯描述，不筛选、不否决 —— 它不改变任何蛋白的去留
〇、输入集：28 个存活蛋白
| 项 | 值 |
|---|---|
| 输入 | results/candidate_status.csv 中 current_status == "retained" 的蛋白 |
| 蛋白数 | 28 |
| 其中 MHC | 11 |
| 其中非 MHC | 17 |
★ 换了输入来源，原因要写清楚
原先 29_enrichment.R 读的是 32_targets_table.R 产的 target_decision_table.csv，取其 step4_survivor 列。
那是旧 v3 口径：step4_survivor = any(repl_class != "opposite") —— 只看 4b 是否方向相反，既不含第 3 步反向 MR 否决的 4 对，也不含 4X 否决的 APOL1。 照它跑，富集输入会比 28 多，且与候选状态总表口径不一致。
→ 现直接读 candidate_status.csv，口径唯一。
依据：C1R 预印本（medRxiv 2026.07.14.26358022）在外部队列验证之后才做通路富集 （"the validated proteins"）。v4 据此规定富集输入 = 4X 之后的存活集。
一、分组：必须按 MHC 拆
| 组名（ASCII） | 蛋白数 |
|---|---|
retained_all | 28 |
retained_no_mhc | 17 |
retained_mhc | 11 |
★★ 为什么 MHC 必须单独跑
MHC 区基因在功能注释上高度富集于免疫通路。与非 MHC 混在一起跑， 会得出「免疫机制驱动糖尿病并发症」这种由 LD 结构造成的假结论。
旧实测（本课题自己的数据）：31 个一起跑，immune response p = 7.2e-07； 只用非 MHC 的 20 个跑，信号消失（最好 p = 0.047）。
→ 「含 MHC vs 不含 MHC」的对比是本课题「MHC 必须单独处理」的核心证据，必须保留。
★ 本轮不做「MR 显著 31」对照组（用户 2026-08-07 决定）。
二、方法
- 工具：g:Profiler REST API（Raudvere et al., NAR 2019）
- 数据源：
GO:BPGO:MFGO:CCKEGGREAC - 多重检验：g:SCS（g:Profiler 自带方法），阈值 0.05
- ★
MICB_MICA会被拆成MICB+MICA，故基因数 ≥ 蛋白数
为什么不用 clusterProfiler
本机 R4.3 下需从源码编译 Bioconductor 依赖，耗时且易半装。 g:Profiler 是同类可引用工具，REST 接口无需编译。
三、★ 预期（写于运行之前）
| # | 预期 | 理由 |
|---|---|---|
| 1 | retained_all（28）会富集到免疫/抗原呈递类条目 | 11/28 是 MHC |
| 2 | ★ retained_no_mhc（17）免疫信号显著减弱或消失 | 旧实测 31→20 时 p 从 7.2e-07 掉到 0.047 |
| 3 | retained_mhc（11）强免疫富集 | MHC 基因本身的注释属性 |
| 4 | 三组都可能条目很少甚至为空 | 输入基因数少（11–29），g:SCS 较保守 |
★ 预期 2 是本步的核心。若非 MHC 组仍有强免疫富集， 说明这批候选的免疫信号不只来自 MHC，那是个需要单独解释的发现， 不许当作「和以前一样」轻轻带过。
3.1 验收断言
| # | 断言 |
|---|---|
| E1 | 输入 = 28 个存活蛋白 |
| E2 | MHC 分组无遗漏无重叠（11 + 17 = 28） |
| E3 | 产物取值全 ASCII（★ 分组名不得用中文——原脚本的中文组名会写进 group 列） |
| E4 | 列名无内部流程编号 |
| E5 | ★ 富集未改动候选状态表 —— 它是描述性步骤，无权改变去留 |
四、产物（预期）
results/enrichment/enrichment_gprofiler.csv（列：source term_id term p_adjterm_size query_size intersection_size group）
五、运行结果
执行日期：2026-08-07 状态：✅ 完成 产物：results/enrichment/enrichment_gprofiler.csv（25 条） 断言：E1–E5 全 PASS
5.1 ★ 结果按数据源分开列（不混排）
为什么按源分列
各源注释到的基因数与背景都不同（KEGG 只覆盖我们 28 个中的 12 个）， 本体结构也不同，分开看才对得上（详见 §5.6）。
★ 不是因为 p 不可比 —— 那是我先前的错误说法，已在 §5.6 更正。
retained_all（28 蛋白）
| 源 | 测到 | g:SCS 阈值 | 显著 | 条目 |
|---|---|---|---|---|
| GO:BP | 1,206 | 3.3e-05 | 4 | regulation of immune response 8.1e-03 (8/970) · T cell mediated immunity 1.8e-02 (4/131) · regulation of T cell mediated cytotoxicity 2.6e-02 (3/46) · T cell mediated cytotoxicity 4.2e-02 (3/54) |
| GO:CC | 170 | 2.6e-04 | 3 | external side of plasma membrane 1.1e-02 · vesicle 1.5e-02 · extracellular region 2.6e-02 |
| KEGG | 61 | 8.2e-04 | 1 | Natural killer cell mediated cytotoxicity 3.5e-02 (3/125) |
| REAC | 116 | 1.5e-04 | 1 | Butyrophilin (BTN) family interactions 4.2e-02 (2/12) |
| GO:MF | 170 | 1.5e-04 | 0 | 最强未达标：heparin binding p=0.314 |
★ retained_no_mhc（17 蛋白）—— 免疫信号彻底不存在
| 源 | 测到 | 阈值 | 显著 | 条目 / 最强未达标项 |
|---|---|---|---|---|
| GO:BP | 894 | 3.8e-05 | ★ 0 | ★ 最强仅 regulation of cellular response to VLDL particle stimulus p=0.878 |
| GO:CC | 132 | 4.0e-04 | 5 | cytoplasmic vesicle 5.9e-03 · intracellular vesicle 6.1e-03 · cytoplasmic vesicle membrane 8.3e-03 · vesicle membrane 9.0e-03 · vesicle 2.2e-02 |
| KEGG | 35 | 9.9e-04 | ★ 0 | 最强 Hepatitis C p=0.23 |
| REAC | 77 | 1.8e-04 | 0 | 最强 Visual phototransduction p=1 |
| GO:MF | 128 | 1.9e-04 | 0 | 最强 metal chelating activity p=0.18 |
★★ 这个阴性有多硬
非 MHC 组的 GO:BP 测了 894 条，显著 0 条，最强的 p 也只有 0.878 —— 离显著差得极远，不是"擦边没过"。KEGG 同样 0 条（最强 p=0.23）。
仅有的 5 条全是 GO:CC 囊泡定位，属细胞组分而非生物学过程， 不构成任何机制主张。
retained_mhc（11 蛋白）
| 源 | 测到 | 阈值 | 显著 | 代表条目 |
|---|---|---|---|---|
| GO:BP | 609 | 4.5e-05 | 9 | regulation of immune response 8.0e-04 · positive regulation of immune system process 1.8e-03 · immune response 3.5e-03 · T cell mediated immunity 2.2e-02 …全部免疫类 |
| REAC | 56 | — | 1 | Butyrophilin (BTN) family interactions |
| GO:MF | 79 | 2.5e-04 | 1 | signaling receptor binding |
| GO:CC / KEGG | 97 / 36 | — | 0 / 0 | — |
5.2 ★★ 预期对账：4/4 命中
| §3 预期 | 实际 | |
|---|---|---|
| 1. 全部组富集到免疫类条目 | GO:BP 4 条全免疫 + KEGG NK 细胞毒性 + REAC Butyrophilin | ✅ |
| 2. 非 MHC 组免疫信号显著减弱或消失 | ★★ GO:BP 0/894，最强 p=0.878；KEGG 0/35 | ✅ 远超预期 |
| 3. MHC 组强免疫富集 | GO:BP 9 条全免疫类 | ✅ |
| 4. 条目可能很少 | 25/3,866 显著 | ✅ |
★★ 本课题「MHC 必须单独处理」最干净的一次证明
| 旧实测（31 → 20 蛋白） | 本轮（28 → 17 蛋白） | |
|---|---|---|
| 含 MHC | immune response p=7.2e-07 | regulation of immune response 8.1e-03 |
| 不含 MHC | 减弱，最好 p=0.047（边缘免疫信号仍在） | ★ GO:BP 0/894，最强 p=0.878 |
→ 这批候选的免疫富集完全由 MHC 区蛋白驱动。 不分组直接拿 28 个跑，会得出「T 细胞 / NK 细胞毒性驱动糖尿病并发症」—— 那是 MHC 区 LD 结构造成的假结论。
富集只是描述，不能当机制证据
- 富集不改变任何蛋白的去留（E5 断言校验）
- 命中数都极小：
Butyrophilin只 2/12、NK cell mediated cytotoxicity只 3/125 - MHC 组的免疫富集本身就是 MHC 基因的注释属性，不是本研究的发现
5.3 关于 KEGG：做了，本轮有 1 条
查询里的数据源始终是 GO:BP GO:MF GO:CC KEGG REAC 五个， 响应元数据也回报了这五个 —— KEGG 从未被排除。
retained_all 命中 1 条 KEGG：Natural killer cell mediated cytotoxicity （3/125，p=0.0351）。非 MHC 组与 MHC 组各 0 条。
★ 首版（用基因符号提交、且丢了 2 个基因）时 KEGG 一条都没有 —— 修正映射后才出现。这说明当时的「KEGG 无富集」是缺陷的产物，不是真结论。
5.4 ★★ 跑中查出的两个缺陷
缺陷一：网络失败被当成阴性结果
| 项 | 内容 |
|---|---|
| 现象 | 三组全部打印「无显著富集」，脚本正常退出（exit 0） |
| 真相 | 本机 WSL 根本连不上 biit.cs.ut.ee —— 直连 90s 超时、走宿主 10808 代理 45s 超时；同一请求在 myserver 上 1 秒返回 |
| 为什么没被断言拦住 | 断言块写在 if (length(all_res)) 里，all_res 为空时整块被跳过 |
| 后果 | 会把「网络故障」写成「这批蛋白无通路富集」这一实质性阴性结论 |
★ 这与第 4X 步的 VEP 静默失败 是同一类事故，本项目第三次遇到。
修法：查询/响应分离——脚本只写 _query_<group>.json、读 _response_<group>.json； 响应缺失／为空／非法 JSON／无 result 字段 → 一律 stop()； 只有 result 存在且长度为 0 才算合法阴性；断言块移出条件块。
缺陷二：★★ 用基因符号提交，静默丢了 2 个基因
| 基因 | 状态 | 原因 |
|---|---|---|
WARS | failed | HGNC 已更名为 WARS1，g:Profiler 不认旧符号 |
SIGLEC5 | ambiguous | 映射到两个 ENSG（ENSG00000268500 有 40 条 GO 注释、ENSG00000105501 有 0 条），g:Profiler 直接弃用 |
后果：非 MHC 组实际只分析了 15 个基因、不是 17。 两个丢掉的基因都与免疫相关（SIGLEC5 是免疫受体）， 直接影响「非 MHC 组有无免疫富集」这个核心结论。
★ E1/E2 查不出来 —— 它们只校验提交数量，不校验实际映射。
修法：改用 Ensembl gene ID 提交（取自 protein_info.csv，28 个存活蛋白全部有）， 从根上消除符号歧义与更名问题；新增断言 E6：failed 与 ambiguous 必须为空 且映射数 == 提交数。修后三组均 28/28 · 17/17 · 11/11 无损。
MICB_MICA 的处理也随之变了
符号提交时它被拆成 MICB + MICA 两个基因； 改 ENSG 后由 protein_info.csv 给出的 ENSG00000204520（MICA） 代表。
这不是将就：第 4X 步实测该探针的 6 个连锁 PAV 全在 MICA、无一在 MICB， 以 MICA 代表有数据支持。故基因数 = 蛋白数 = 28。
5.5 顺带修掉的约定违规
原脚本的分组名是中文（第4步存活-全部 等），会写进产物的 group 列， 违反「CSV 取值一律 ASCII」。已改为 retained_all / retained_no_mhc / retained_mhc， 并加断言 E3 把关。
5.6 ★★ 更正：KEGG 与 GO 的 p 是可比的（我先前说错了）
本节更正 2026-08-07 早先的错误说法
我先前写过「KEGG p=0.035 与 GO:BP p=0.035 证据强度完全不同，因为校正负担不同」。 这是错的。
实测：手算超几何原始 p 与 g:Profiler 返回值逐条比对， 比值在每个「组 × 源」内是常数——即返回的 p_value 已经是校正后的 p：
| 源 / 组 | 原始超几何 p | 返回 p | 倍数 |
|---|---|---|---|
| GO:BP / mhc | 7.18e-07 | 7.98e-04 | ×1111.1 |
| GO:BP / mhc | 1.60e-06 | 1.77e-03 | ×1111.1 |
| GO:BP / all | 5.36e-06 | 8.12e-03 | ×1515.2 |
| GO:CC / no_mhc | 4.79e-05 | 5.92e-03 | ×123.8 |
| REAC / all | 2.26e-05 | 5.26e-03 | ×232.6 |
且 0.05 / 1111.1 = 4.5e-05 正是报告的 retained_mhc / GO:BP 阈值； 0.05 / 1515.2 = 3.3e-05 正是 retained_all / GO:BP 的阈值。
→ g:SCS 的校正因子已经把各源的检验负担吸收掉了， 两个校正后 p 落在同一尺度上，作为校正后 p 可比。
★ 因此 gscs_threshold 列是信息性的（说明该源的负担）， 筛显著直接用 p_adj < 0.05 即可，不需要它。
真正仍然成立的差异（这才是分面的理由）：
| 源 | 词条数 | g:SCS 阈值 | 背景基因数 | 我们 28 个基因中被注释的 |
|---|---|---|---|---|
| GO:BP | 24,546 | 3.3e-05 | 20,972 | 23 |
| GO:MF | 10,123 | 1.48e-04 | 20,208 | — |
| GO:CC | 4,069 | 2.62e-04 | 22,155 | 24 |
| KEGG | ★ 586 | ★ 8.2e-04 | ★ 8,716 | ★ 12 |
| REAC | 2,849 | 1.52e-04 | 11,056 | 16 |
三条真实差异
- ★ 分析的根本不是同一批基因 —— KEGG 只注释了我们 28 个中的 12 个， GO:BP 注释了 23 个。那条 KEGG 结果实际只用了 12 个基因
- 背景不同（8,716 vs 20,972）
- 本体结构不同 —— GO 条目高度嵌套冗余（
immune response⊃regulation of immune response⊃ …），KEGG 通路相对独立； 混排会让 GO 的冗余条目挤占版面
→ 分面作图仍然该做，但理由是 「基因覆盖与本体结构不同、且便于阅读」， 不是「p 不可比」。这个区别必须写清楚，别再说错。
5.7 与常规 KEGG 分析（clusterProfiler / DAVID）的差别
| 项 | 本轮（g:Profiler） | 常规做法 |
|---|---|---|
| 多重检验 | g:SCS（g:Profiler 自有方法） | 多为 BH-FDR |
| 校正范围 | 按源分别校正 | 通常单源单独跑，各自 FDR |
| 背景 | domain_scope: annotated（该源有注释的基因） | 常用全基因组或自定义 universe |
| KEGG 版本 | g:Profiler 快照 e114_eg62_p19_27110d83 | clusterProfiler 走 KEGG REST 取当前版本 |
→ 同一批基因在两种流程下结果不同是预期内的，不是谁算错了。 写 Methods 时必须写明用的是 g:Profiler + g:SCS + annotated 背景 + 版本号。
5.8 产物布局：按源拆开，每个源都有自己的全表
results/enrichment/
├── enrichment_gprofiler_all.csv 全表 3,866 条（含不显著，带 significant 列）
├── enrichment_gprofiler.csv 仅显著 25 条
├── enrichment_source_metadata.csv 各源词条数 / g:SCS 阈值 / 背景 15 行
└── by_source/ ★ 15 个文件 = 3 组 x 5 源
├── retained_all__GO_BP.csv 1,206 条（显著 4，阈值 3.3e-05）
├── retained_all__GO_CC.csv 170 条（显著 3，阈值 2.6e-04）
├── retained_all__GO_MF.csv 170 条（显著 0，阈值 1.5e-04）
├── retained_all__KEGG.csv ★ 61 条（显著 1，阈值 8.2e-04）
├── retained_all__REAC.csv 116 条（显著 1，阈值 1.5e-04）
├── retained_no_mhc__GO_BP.csv 894 条（显著 0，阈值 3.8e-05）
├── retained_no_mhc__GO_CC.csv 132 条（显著 5，阈值 4.0e-04）
├── retained_no_mhc__GO_MF.csv 128 条（显著 0，阈值 1.9e-04）
├── retained_no_mhc__KEGG.csv ★ 35 条（显著 0，阈值 9.9e-04）
├── retained_no_mhc__REAC.csv 77 条（显著 0，阈值 1.8e-04）
├── retained_mhc__GO_BP.csv 609 条（显著 9，阈值 4.5e-05）
├── retained_mhc__GO_CC.csv 97 条（显著 0）
├── retained_mhc__GO_MF.csv 79 条（显著 1）
├── retained_mhc__KEGG.csv ★ 36 条（显著 0）
└── retained_mhc__REAC.csv 56 条（显著 1）★ by_source/ 里的每个文件都是该源的完整结果（含不显著）， 并自带 gscs_threshold 列 —— 作图脚本直接读它，不必再去拼阈值。
为什么必须拆成独立文件，而不是靠总表加 source 列筛
两条都是实际发生过的：
- 总表把 KEGG 与 GO 放在一起，看的人会顺手跨源比 p —— 而它们不可比（KEGG 586 词条阈值 8.2e-04；GO:BP 24,546 词条阈值 3.3e-05）
- 主表只有显著项，「某个源的原始结果在哪」要靠猜 —— 例如 KEGG 在主表里只有 1 行，会被误读成「KEGG 只测了 1 条」， 实际测了 61 条（
retained_all__KEGG.csv）
新增断言 E9 强制每个「组 × 源」组合都有结果。
KEGG 原始结果长这样（retained_all__KEGG.csv，61 条按 p 排序）：
| term | p_adj | significant | hit |
|---|---|---|---|
| Natural killer cell mediated cytotoxicity | 0.0351 | TRUE | 3/125 |
| Kaposi sarcoma-associated herpesvirus infection | 0.127 | FALSE | 3/195 |
| Human papillomavirus infection | 0.564 | FALSE | 3/331 |
| PI3K-Akt signaling pathway | 0.699 | FALSE | 3/358 |
| …（其余 57 条 p 多为 1） |
非 MHC 组的 KEGG（retained_no_mhc__KEGG.csv，35 条）最强也只有 Hepatitis C p=0.23、JAK-STAT signaling pathway p=0.251 —— 无一显著。
六、★★ 作图规则（第 7 步必须遵守，不是建议）
三条硬约束
- 按数据源分面（facet），不做跨源混排的排名图。 理由（见 §5.6 更正）：不是 p 不可比 —— 而是各源注释到的基因数不同（KEGG 仅 12/28，GO:BP 23/28）、 背景不同（8,716 vs 20,972）、本体冗余结构不同。
- GO 的 BP / MF / CC 三支之间同样分开 —— 三个独立本体，条目体量差 6 倍。
- 每张图的图注必须写明：该源的词条数、背景基因数、 该源实际参与的基因数（KEGG 只有 12/28），以及 多重检验用的是 g:SCS 而非 BH-FDR。
★ 一条已更正的说法，别再沿用
早先版本写过「KEGG 与 GO 的 p 不可直接比较」。那是错的 —— g:SCS 返回的已是校正后 p，负担差异已被校正因子吸收（§5.6 有实测）。 分面的理由是覆盖与结构，不是 p 的可比性。
允许的做法：按 source 分面（facet）的多panel图，每个 panel 内部再按 p 排序； 或每个源单独一张图。
产物依据：by_source/ 下每个文件已自带该源的完整结果与 gscs_threshold； enrichment_source_metadata.csv 存各源词条数/背景/阈值。 作图脚本必须读产物取这些数，不许在图上手写。
★ 注意 gscs_threshold 是原始 p 的等效截断值，用于说明检验负担； 筛显著请用 p_adj < 0.05，不要拿 p_adj 去和 gscs_threshold 比。
为什么写在这里而不是等到第 7 步
「作图留到最后统一做」是已定原则，但规则要在有结论的时候就写死， 否则三个月后画图的人（包括我）只会看到一张混排的表，然后照着画。