Skip to content

14 · 第 6 步 · 描述性注释

日期:2026-08-08 状态:6-1 跑前记录(尚未运行

★★ 本步的铁律:只加分,不否决

v4 明文:6-1~6-4 只加分不否决全流程的否决点只有两个 (第 3 步反向 MR、4X 表位/PAV),第 5 步只升降措辞。

第 6 步无权改变任何候选的去留。 → 断言 N2 强制校验 candidate_status.csvretained 集合跑前跑后完全相同


〇、输入集(用户 2026-08-08 拍板)

只做共定位支持的 11 个蛋白PP.H4 ≥ 0.8,全部非 MHC)。 MHC 且 PP.H4 < 0.8 的不进入本步,改走「MHC 专用核验」。

四个子步骤 6-1 / 6-2 / 6-3 / 6-4 全部只做这 11 个

已核对的输入(实测,非推断)

蛋白哨兵 SNPMR 结局数最高 PP.H4
APOErs4293583(★ 其中 Reti 在第 3 步被否决)0.998
NUDT5rs1050843810.977
ERMAPrs1121071010.967
PAMrs14980297810.959
WARSrs227380410.956
GALNT3rs211654610.948
IFNAR1rs91414210.940
NOTCH2rs264134820.897
SIGLEC5rs110647620.887
ACRBPrs795965810.867
LACTB2rs19158809910.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 文献底本

#来源做法
C1RmedRxiv 2026.07.14.26358022 · Methods 第 407–414 行GWAS ATLAS,4,756 个性状 / 28 个表型大类,阈值 P < 5×10⁻⁸
C1RFigure 4 图注方向感知:"associations of risk alleles … triangles up representing positive associations and down representing negative"
C1RDiscussion 第 296–298 行★ 自陈局限:"PheWAS based on a cis-pQTL cannot fully recapitulate pharmacological inhibition"
本课题analysis/03_phewas.R + 37_drug_target_phewas.ROpenGWAS,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_siggws_sig)供读者自行取用。

理由:Bonferroni 分母应与实际做的检验数一致(这是它的定义), 但一条 SNP–性状关联若连全基因组显著都达不到,本身就不可信。两条都要过。

1.3 方法与代码

脚本本轮须改
查询analysis/03_phewas.R输入集由「所有显著 SNP」改为这 11 个;缓存指纹随之失效重查
判读analysis/37_drug_target_phewas.R★ 现读 results/targets/target_decision_table.csvtier 表,用户已定删除)→ 改读 results/candidate_status.csvcoloc/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 条关联
E3Bonferroni 分母 == 实际 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 库有变,须记录
E6APOE远超其余 10 个,量级差 1–2 个数量级rs429358 是已知极端多效位点
E7疾病类脱靶:ERMAPIFNAR1 接近 02026-07-31 实测 ERMAP 0 条、IFNAR1 1 条(红细胞分布宽度,属定量性状不参与判读)
E8LACTB2 结果未知首次查询;★ 无论结果如何都必须与「没查」区分开

预期与实际不符怎么办

一律先按缺陷处理,不许现场解释。查清是数据变了、代码错了、还是预期本身写错了, 判完再写进本页。


三、验收断言(抓不到即 stop()

3.1 全步通用(N 系列)

#断言
N1输入集 == coloc_verdict.csvcoloc_tier=="strong" & in_mhc==0 的蛋白集合,且恰为 11 个11 个唯一哨兵 SNP
N2★★ 跑前跑后 candidate_status.csvretained 集合完全相同(第 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 消歧,消不掉的一律剔除并记数
P3Bonferroni 分母 == 实际查询 SNP 数 × 运行时数据集数(写入 provenance)
P4跨仓对账:10 个旧仓覆盖的 SNP,在同一阈值 1.06×10⁻⁸ 下关联数与旧仓比对;不一致必须逐条列出并说明(可能是 OpenGWAS 库变动,不是自动通过)
P5性状三分类穷尽且互斥:每条关联恰属疾病 / 定量 / 分子之一,未分类数 == 0
P6implied_drug_action 逐蛋白与 coloc_verdict.csv 完全一致(抑制 5 / 激动 6)

四、★ 已知坑(本轮必须防住)

来自 2026-07-31 的实测,见记忆 project_drug_target_phewas_20260731

#本轮对策
1in_MHC 静默丢数据 —— 直接 merge target_decision_table.csvin_MHC,该表只覆盖 31 个蛋白,其余是 NA,in_MHC == FALSE 把它们整批悄悄滤掉且不报错改用 candidate_status.csvin_mhc(覆盖全部 57 行);断言无 NA
2「没测」≠「测了没有」 —— LACTB2 等三个蛋白的 SNP 根本不在 phewas_all.csv 里,被记成「0 条脱靶」断言 P1LACTB2 正是这次的当事者
3只按 ID 前缀滤分子性状会漏 —— ebi-a-GCST900102xx 那批「Galectin-4 levels」不带 prot- 前缀补性状名判据(以 " levels" 结尾);断言 P5
4data.tablej 里变量名不能叫 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 互不为输入、都不否决),但有两个实在好处:

  1. 成药性给 PheWAS 提供对照基准 —— 若靶点已有上市药,PheWAS 的脱靶谱可与该药已知不良反应谱互相印证。C1R 就是这个写法:先查出 C1R 有 conestat alfa / C1-酯酶抑制剂(遗传性血管性水肿),再说 PheWAS 无显著关联。
  2. 不依赖会过期的 token —— Open Targets 公共 API 无需鉴权;OpenGWAS 的 JWT 有效期只有 14 天(本轮就是因过期而中断)。

★ 澄清:PheWAS 的方向感知不依赖 6-2。用药方向已在 coloc_verdict.csvimplied_drug_actioninhibit 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
inhibitINHIBITOR / ANTAGONIST / NEGATIVE MODULATOR / BLOCKER / DEGRADER
activate(= agonize_or_supplementAGONIST / 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 ★ 预期(写于运行之前)

机械可验

#预期
D111 个蛋白全部解析到 Ensembl ID(NOT_FOUND 与「查到但无药」必须分开)
D2WARS 须经别名表命中 WARS1(HGNC 现行符号)—— 与富集那次栽的是同一个改名问题
D3用药方向与 coloc_verdict.csv 一致:inhibit 5 / activate 6

内容预期(★ 最要紧的一条)

#预期依据
D4★★ IFNAR1 会查到已上市药(anifrolumab,抗 IFNAR1 单抗,SLE,Phase 4),但 dir_match 应为 direction_opposite_or_unknown我们推得 IFNAR1激动/补充(风险降低方向),而现有药是拮抗——方向正好相反
D5NOTCH2 可能查到抑制类药物(γ-分泌酶抑制剂 / 抗 NOTCH 抗体),方向同样相反我们推得需激动
D6其余多数蛋白(ACRBP ERMAP GALNT3 LACTB2 NUDT5 PAM SIGLEC5 WARS)预期 no_known_drug均非经典药靶
D7tractability 多数有 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

蛋白用药方向小分子抗体药物数已上市最高期方向判定
IFNAR1activatenoyes12124direction_consistent
NOTCH2activateyesyes102direction_opposite_or_unknown
APOEinhibityesyes000no_known_drug
PAMactivateyesyes000no_known_drug
NUDT5inhibityesno000no_known_drug
WARSinhibityesno000no_known_drug
ACRBPinhibitnoyes000no_known_drug
ERMAPactivatenoyes000no_known_drug
GALNT3activatenoyes000no_known_drug
SIGLEC5inhibitnoyes000no_known_drug
LACTB2activatenono000no_known_drug

可及性(tractability,预测分桶,不等于有药):抗体可及 8/11、小分子可及 5/11、 两者皆无 1/11LACTB2)。

★★ 预期对账: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利好,且指向药物重定位的可能性。

但不能就此说「可以重定位干扰素」

三条必须同时写:

  1. cis-pQTL 推得的方向 ≠ 药理干预净效应(C1R 自陈的同一条局限)
  2. 干扰素类有明确且严重的不良反应谱,IFNAR1 通路双向都有临床用途(anifrolumab 治 SLE)
  3. 我们的 IFNAR1 证据来自血浆蛋白水平,与受体激动的组织效应不是一回事

★ 抓到的缺陷:APPROVAL 被静默算成 0 期

第一次运行时 IFNAR1 显示 drugs=12maxphase=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 curlrc=28 超时 / 返回非 JSON
myserver curlhttp=000 / 返回 EBI 的 HTML 页面
WebFetchHTTP 500

三条路径都不通 → 是 EBI 端的问题,不是我们的网络。 (间歇性:GALNT3CHEMBL4523291NUDT5CHEMBL4105713 曾成功解析。)

最终结果(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)
APOErs4293580.1561262951896
PAMrs1498029780.054884234
NUDT5rs105084380.336683834
WARSrs22738040.260533227
SIGLEC5rs11064760.113593229
NOTCH2rs26413480.111643025
ACRBPrs79596580.361412019
ERMAPrs112107100.476666
IFNAR1rs9141420.271944
GALNT3rs21165460.298922
LACTB2rs1915880990.007100★ 旧仓从未查过

7.2 ★ 预期对账

#预期实际
E111 个 SNP 全查,含 rs191588099
E2每 SNP 都有「已查询」记录,哪怕 0 关联✅ 覆盖表 11 行,LACTB2 = queried + 0 条
E3Bonferroni 分母 = SNP 数 × 数据集数✅ 11 × 50,164 = 551,804
E5各 SNP 关联数 ≥ 旧仓(本轮阈值更松)10/10 全部 ≥(见上表右两列)
E6APOE 远超其余 1–2 个数量级✅ 951 vs 次高 42,约 23 倍
E7疾病类脱靶 ERMAP/IFNAR1 接近 0⏸ 待方向感知判读(需性状分类)
E8LACTB2 结果未知✅ 返回 0 条 —— 但见下

其他:别名 rsID 按坐标归一 3 → 0 条未能归一;蛋白归属缺失 0 条。

7.3 ★★ LACTB2 的 0 条不能读成「干净靶点」

rs191588099MAF = 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.241APOE MAF 0.156 却有 1262 条), 所以不是 MAF 普遍驱动关联数LACTB2 是单点的稀有性问题,不是全局趋势。

7.4 ⏸ 方向感知判读:待决

余下一半(获益/风险判读、疾病 vs 定量 vs 分子的三分类,即 E4/E7) 由 37_drug_target_phewas.R 完成,但它有两处依赖不可用:

依赖问题
results/targets/target_decision_table.csvtier 表,用户已定删除,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 条全 PASST1T7 + 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'sT2DEarly AMDCataracts 全部正确判为 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,改为新脚本,原因两条:

  1. 37 读 results/targets/target_decision_table.csvtier 表,已定删除,且只覆盖 31 个蛋白 —— 直接 merge 会让其余蛋白的 in_MHC 变 NA 后被整批静默滤掉,第一版就踩过)
  2. ★★ 37 内含一条已被本轮推翻的性状再分类规则:grepl("\\blevels$", trait) => molecular

移植了 37 最有价值的部分:等位翻转矩阵、回文位点用 EAF 消歧、随机抽样自检。

等位对齐:1,641 条可对齐 / 18 条无法对齐(回文不可消歧或碱基不匹配,已剔除)。 对齐旁证:命中自身结局 1 条,方向与 MR 同号 1/1

断言 12 条全 PASSN1N3 + O1O9),其中 N2 校验 跑前跑后 retained 集合完全相同 —— 第 6 步无否决权,已在代码里落实。

★ 判读前抓到的两处「疾病类」渗漏

#现象根因修法
1SIGLEC5/NOTCH2/APOE 的「疾病脱靶」里出现 basophil / monocyte / lymphocyte cell count这三个的 categoryContinuous,但 subcategoryImmune system;而我把 subcategory 关键词规则排在 Continuous 之前immune 命中了。「Immune system」是身体系统标签,不是疾病标志连续性状规则前移;关键词由裸 immune 改为 autoimmune;立 T8
2PAM 出现 3 条「脱靶获益」:metformin 用药码ICD10 E11.9on_target_axis 只查了 diabet 字样,漏掉 UKB 的 ICD 编码型与用药编码型糖尿病代理E10/E11 ICD 正则与 16 种降糖药名;立 O8/O9

⚠️ 反向也要守住:APOE他汀/依折麦布用药码是真脱靶(反映脂质轴),不能一起剔除。

修复前后:PAM 疾病脱靶 3 → 0SIGLEC5 2 → 0NUDT5 2 → 1;NOTCH2 4 → 3。


八、★ 第 6 步当前结论:11 个蛋白

蛋白用药方向索引结局OR显著脱靶疾病类定量分子获益风险已上市药
APOEinhibit肾病1.151926131469326251060
NUDT5inhibit视网膜3.080291244100
SIGLEC5inhibit黄斑1.105290236000
WARSinhibit视网膜1.385280208000
NOTCH2activate视网膜0.598243156030
ACRBPinhibit视网膜1.450200812000
PAMactivate视网膜0.888180135000
ERMAPactivate黄斑0.3296006000
IFNAR1activate黄斑0.78240130012
GALNT3activate视网膜0.8462002000
LACTB2activate视网膜0.674NANANANANANA0

8.1 三条主要发现

① 8/10 可评估的蛋白疾病类脱靶 = 0。 SIGLEC5 WARS ACRBP PAM ERMAP IFNAR1 GALNT3 全为 0,NUDT5 仅 1 条(FEV1/FVC < 0.7,判为获益)。 ★ 其中 ERMAPIFNAR1 的 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: ButterIllnesses of mother:Major dietary changes 等)。

⚠️ 而我为量化它写的关键词筛查本身也过度捕获 —— 它把 Nonalcoholic fatty liver diseaseDelirium, 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 文献底本

#来源内容
R1Advani 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)。★ 本地有 PDFF:\project\s41467-024-46063-8.pdf
R2Single-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 开放
R3C1R(medRxiv 2026.07.14.26358022)Methods "Gene expression analyses"FUMA GENE2FUNC(基于 GTEx 的通用组织)

我们比 C1R 强的地方要写进 Discussion:他们用 GTEx 通用组织做描述性表达定位; 我们用视网膜专属 eQTLMR + 共定位(因果层),再叠单细胞定位(描述层)。 ⚠️ 两层性质不同,不能混着说 —— 见记忆 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 个基因:

基因结局方法bpPP.H4
ERMAP黄斑Wald ratio+0.09292.52e-050.9488
ERMAP视网膜Wald ratio+0.00580.6610.0081
IFNAR1黄斑IVW+0.00310.7900.0189
IFNAR1视网膜IVW−0.00010.9750.0098
APOL10 条记录(非 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.RR9 <- "/home/research/mr-pipeline-r9"旧仓);靶点读 common/targets_r9dm.csv(只有 3 个)指向 mr-pipeline-r9v3;靶点改由 coloc_verdict.csv 生成 11 个
mr-pipeline-r9v3/analysis/_singlecell_io.pyload_targets()results/network_R9_dm/target_list.csvtier 表,v3 不存在改读 coloc_verdict.csv

9.6 ★ 预期(写于运行之前)

#预期依据
S1Advani 只覆盖 9,408 个基因(非全转录组)→ 11 个里会有若干个不是 eGene(0 条记录)记忆 reference_retina_eqtl 明载
S2ERMAP 黄斑 PP.H4 复现 ≈ 0.949、b 为IFNAR1 两个结局均阴性9.4 实测
S3APOL1 不在本轮 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.csvretained 集合完全相同
N3产物全 ASCII、无内部流程编号与内部结局 ID
N4★ 输入文件缺失 / 读取失败必须 stop()不得当成「该基因无记录」
R1★ ENSG 号一律取 druggability.csvensembl,并用坐标与哨兵 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
GALNT3ENSG00000115339924782
NUDT5ENSG00000165609278241
IFNAR1ENSG00000142166172167
LACTB2ENSG00000147592147147
ERMAPENSG00000164010136136
WARSENSG000001401059696
ACRBPENSG000001116442525
APOEENSG0000013020344
NOTCH2 PAM SIGLEC500✘ 非 eGene

8/11 是视网膜 eGeneS1 预期命中)。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 行)

基因结局nsnpbpFDRPP.H4vs 血浆
ERMAP黄斑1+0.09292.5e-051.5e-040.9488相反
APOE黄斑1−0.00795.3e-030.016NA相反
WARS黄斑1−0.00680.0230.0450.1735血浆不显著
NUDT5黄斑7−0.00810.200.300.0131
GALNT3黄斑6+0.00770.570.690.0161
IFNAR1黄斑3+0.00310.790.790.0189相反
WARS视网膜1−0.00772.0e-051.2e-04★★ 0.8086相反
NUDT5视网膜7−0.02018.7e-042.6e-030.0346相反
APOE视网膜1−0.00300.0830.17NA相反
GALNT3视网膜6+0.01130.330.490.0089
ERMAP视网膜1+0.00580.660.790.0081
IFNAR1视网膜3−0.00010.970.970.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,557full n=54,219;2025-10 Olink+SomaScan 整合图谱(>90,000)⚠️ 故意不用最大的(discovery 有留出复制;full 原文自称 putative),已在定稿决定里
视网膜 eQTLAdvani 2024,403 眼同一个已是最大的开放数据。EyeGEx 406 眼受控访问且 312/406 是病例;mega-analysis 2020 仅 311 眼
单细胞scRNA 265,767 / 18 类snRNA 3,177,310 / 31 类差 12 倍 —— 已换,见 11.1
血液 eQTLeQTLGen Phase 1 · 31,684Phase 2 · 43,301(medRxiv 2026-02 · PMID 42051578)落后 37%,属 6-3 范围
sQTL / 多组织 eQTLGTEx v8GTEx v10(2024-11,eQTL 样本 +23%★ 落后一个大版本,属 6-3 范围
代谢物 · 血液 mQTLNightingale 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,767182.54 GB
snRNA-seq all cells(已下载)3,177,3103137.67 GB

逐类型差集核对结论

  • RPE、小胶质、Müller、星形胶质,两个版本都有 —— 关键大类没丢
  • 多出的 15 类几乎全是双极细胞与节细胞的精细亚型diffuse bipolar 1/2/3a/3b/4/6ON/OFF midget ganglion…),对本课题不关键
  • 真正收益在细胞数(12 倍) —— 小胶质这类稀有细胞在 265k 里可能只有几百个, 在 3.18M 里是几千,对检出功效与 AUC 富集计算差别是实的

换数据集也解决不了的硬伤

两个版本都没有内皮细胞与周细胞。 糖尿病视网膜病变是血管病, 这一层缺失必须写进 Limitations。

下载:wget -c 直连 11.1 MB/s(并发 4 路仅 8 MB/s,瓶颈在上行带宽,不值得分块)。


十二、当前范围审计(★ 逐产物点过,非印象)

步骤蛋白数
第 2 步 MR 显著31
第 3 步 反向 MR95(含未显著者)
4b 换队列66
4c 两侧都换12
第 5 步 共定位28
6-1 PheWAS 查询 / 判读11
6-2 成药性11
6-4 视网膜 eQTL11

| 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.pdfpdftotext -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
RPERPE 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
2MAF ≥ 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 WARS31/31✅ 全覆盖
GALNT329/31✅ 全覆盖
ACRBP15/31❌ 无
SIGLEC51/31❌ 无

ACRBP / SIGLEC5 的缺席是**「未进入检验」不是「不表达」**(未过表达过滤), 用的时候必须三态记账,不得写成阴性。

★★ WARS 在这份数据里就叫 WARS,不是 WARS1 —— 与 HRCA 相反。 别把别名规则写死,要按数据集分别处理。

血管细胞里的表达比例(pct_cluster_exp

基因脉络膜毛细血管动脉静脉周细胞
PAM0.3690.4050.4310.278
WARS0.1700.2150.2070.074
NOTCH20.0650.0810.0860.173
IFNAR10.0750.0980.0790.077
APOE0.0690.0780.1010.090
NUDT50.0360.0500.0440.059
GALNT30.0180.0160.0420.013
ERMAP0.0170.0190.0180.020
LACTB20.0160.0160.0220.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 已实测核实版本)

用什么版本判定
血液 eQTLeQTLGen 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
多组织 sQTLGTEx v10 GTEx_Analysis_v10_sQTL.tar(1.96 GB)可公开下载,非 requester-paysgs://adult-gtex/bulk-qtl/v10/single-tissue-cis-qtl/。★ 应从现用的 v8 升级
(备选)多组织 eQTLGTEx_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 ★ 预期(写于运行之前,与实际不符一律先按缺陷处理)

#预期依据
E111 个基因里 8–11 个在 eQTLGen 有 cis-eQTL旧仓实测 IFNAR1 799 条、ERMAP 1354 条;eQTLGen 覆盖 88.6% 基因
E2LACTB2 很可能不可评估其哨兵 MAF=0.0071,eQTLGen 用的是常见变异
E3WARS 必须用 WARS1HGNC 已改名;这是本轮踩过的坑
E4跨层方向至少有 2–3 个相反旧仓 ERMAP 血浆 vs 视网膜已相反;跨层非同一变异
E5IFNAR1 血液 eQTL 层有信号(与视网膜阴性形成对照)旧仓 07_r9dm 已跑过 IFNAR1,可直接对账
E6GTEx 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参照等位必须与血浆层对齐后才比较方向 —— 未对齐的方向比较一律作废
X4GTEx 版本号写进产物元数据(v10),且断言不是 v8
X5Steiger filtering 的 unitsprevalence 必须齐全,缺失时显式跳过并记录,不得静默通过
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
★★ 2SNPPosGRCh37/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 个并发症(视网膜/黄斑/肾病/神经病变),不限于索引结局全血是系统性组织,没有理由只对眼病测;且四个一起跑才有跨结局对照
3cis 窗口 ±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 = ...) —— keydata.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 后)
PAMENSG00000145730ensembl7,0985,03250
SIGLEC5ENSG00000105501symbol9,2351,05864
LACTB2ENSG00000147592ensembl6,0431,44421
WARSENSG00000140105ensembl7,0461,39955
ERMAPENSG00000164010ensembl5,6981,35438
IFNAR1ENSG00000142166ensembl6,04379939
NUDT5ENSG00000165609ensembl7,61937218
NOTCH2ENSG00000134250ensembl3,6712129
GALNT3ENSG00000115339ensembl5,8982102
ACRBPENSG00000111644ensembl6,667954
APOEENSG00000130203ensembl7,0760

流程量:抽 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.9773OR 3.080risk_increasing
全血转录本(本层)0.9821b = +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/seZscore + 等位频率换算(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.h5ad3,177,310 细胞 × 35,475 基因 / 31 类, 非零 8,248,015,746(每行均 2,596)。用的是 X(归一化值),不是原始计数。

16.1 ★ SIGLEC5 第二次踩同一个坑,被双键拦下

图谱 var没有 ENSG00000268500druggability.csv 给的), 有的是 ENSG00000105501。若只用主 ENSG,它会第二次被静默报成「不表达」。

resolved_by = ensembl_alternate 已逐行记进 singlecell_gene_resolution.csv。 → 跨数据集的基因标识必须多路回退 + 逐个留痕,这已是同一天第二次。

16.2 ★★ AUC 首跑全错:平局按行序破了

首跑用 argsort 给 1..n 的唯一秩,等于对并列的 0 值按行号排先后。 而本图谱的细胞按类型排序存放

细胞类型n行号相对位置
retinal pigment epithelial cell863[1.000, 1.000](全在文件最末)
microglial cell4,894[0.998, 1.000]
Mueller cell221,612[0.924, 0.994]
retinal rod cell1,066,056[0.523, 0.859]

于是排在末尾的细胞类型凭空拿到高秩。露馅的数字SIGLEC5 × RPE 表达比例 0.0% 却得到 AUC 0.996

  • 受影响:只有 auc_vs_restpct_expressing / mean_expression 不用秩,不受影响
  • 修法:scipy.stats.rankdata(method="average")
  • 旧产物改名 .badauc-171100 留档,未删

新增断言 SC5

一个细胞都不表达的类型,AUC 不可能 > 0.5 (平均秩下全零组 AUC = 0.5 × 对照组中同为 0 的比例 ≤ 0.5)

★ 这条若首跑就写,当场就会被拦下。「不可能的数字」值得写成断言,而不是靠眼看。

16.3 结果(修正后)

基因AUC 最高的细胞类型AUC表达比例次高
PAMGABA 能无长突细胞0.86689.4%OFF 伞状神经节 0.837(97.8%)
APOEMüller 细胞0.79466.9%小胶质 0.697(51.5%)
NOTCH2Müller 细胞0.76761.2%星形胶质 0.746(60.7%)
NUDT5OFF 伞状神经节0.67259.7%ON 中央凹侏儒神经节 0.623
ERMAP视锥细胞0.66135.5%S 视锥 0.575
GALNT3Müller 细胞0.65131.1%星形胶质 0.539
WARSOFF 伞状神经节0.59133.5%视锥 0.558
IFNAR1OFF 伞状神经节0.56036.2%H1 水平细胞 0.553
LACTB2OFF 伞状神经节0.53213.8%S 视锥 0.523
SIGLEC5小胶质0.5163.7%侵入型侏儒双极 0.504
ACRBP弥散双极 3b0.5051.6%— 无富集

16.4 与旧版(scRNA 子集 265,767 细胞 / 18 类)的差异

基因旧值新值说明
IFNAR1伞状神经节 74.4% / AUC 0.75736.2% / 0.560细胞数 12 倍、snRNA vs scRNA,数值变化在预期内
ERMAPAUC 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 GB

FTP 根:https://ftp.ebi.ac.uk/pub/databases/spot/eQTL/sumstats/QTS000015/<QTD*>/ 索引表:github.com/eQTL-Catalogue/eQTL-Catalogue-resourcestabix/tabix_ftp_paths.tsv

★★ 下面这段结论已被 §19.1 部分推翻,保留原文以便对照

当时只看了变异层的 p 分布就下了「不是显著性过滤」的判断。正确说法是: 变异层没筛(零分布完整,故共定位有效),但性状层就是按显著性筛的.cc 只含 FDR 显著的分子性状,保留比例 0.58%–9.9%)。以 §19.1 的表为准。

★ 实测澄清:.cc 不是显著性过滤。 取头部 3 MB 解压看 pvalue 分布:

文件nminmedianmax
QTD000261.all(肾 ge)19,9983.95e-060.4890.99999
QTD000261.cc(肾 ge)19,9987.15e-120.4380.99998
QTD000290.cc(神经 sQTL)19,9983.05e-180.4660.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非索引结局信号需独立层判真伪
E5GTEx 全血(n=670)是流水线正对照:应能复现 eQTLGen(n=31,684)已确认信号中的至少 1 个,但功效低得多同组织不同队列;全 0 即接线错误而非生物学
E6SIGLEC5 的 ENSG 又会不一致druggability.csvENSG00000268500同一坑第三次(eQTLGen …105501、HRCA …105501
E7rsID 命中率 >95%本层 GRCh38 且 1000G 30x,与 FinnGen R9 同 build(eQTLGen 因跨 build 只有 90.8%)
E8sQTL 层部分基因缺席 .cc.cc 按性状筛;须靠 .permuted 定三态
E9GTEx 无眼组织 —— 眼病结局在本层只能用血/神经替代,这是结构性限制不是阴性结果与 6-3 的 E6
E10LACTB2 大概率仍不可评估哨兵 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},由 E3APOE 对照决定取值;不得单写 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
K8tabix 切片行数 >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 最小组织之一E3APOE 功效对照 + K3negative_reason;Limitations 必写
BsQTL 只有 .cc 没有 .all.permuted 定三态(K2/K8
CGTEx 无眼组织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 个基因在三组织的被测与显著情况(.permutedp_beta

基因胫神经 n=532肾皮质 n=73全血 n=670
ACRBP★ 7.15e-040.2520.0527
APOE5.46e-170.3970.567
ERMAP0.05680.2730.705
GALNT3★ 5.32e-410.1080.410
IFNAR11.34e-110.200★ 0.0392
LACTB2★ 1.32e-040.622★ 1.54e-04
NOTCH2★ 1.96e-020.7960.967
NUDT5★ 3.26e-090.931★ 0.0492
PAM0.2670.984★ 9.45e-63
SIGLEC5★ 3.05e-050.657★ 2.85e-27
WARS0.2640.543★ 2.81e-26
p_beta<0.058/11★★ 0/116/11
被测性状总数23,78824,31016,701

跑前预期对账(★ 逐条,含我预测错的)

#预期实际判定
E1胫神经 6–9 个显著8/11✅ 命中
E2肾皮质显著者 ≤30/11✅ 命中(比预期更极端)
E3APOE 作肾层功效对照p_beta=0.397 不显著;同基因胫神经 5.46e-17★★ 触发:判定肾皮质 n=73 完全无功效
E5全血能复现 eQTLGen 已确认信号 ≥16/11(含 IFNAR1/NUDT5/WARS/PAM/SIGLEC5/LACTB2✅ 正对照成立
E6SIGLEC5 的 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 → ★ E4IFNAR1 × 神经病变那个 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
PAM1654.471e-16rs385827
WARS1311.669e-16rs941926
LACTB2962.501e-05rs201659904
IFNAR1936.881e-12rs2040109
GALNT338.001e-06rs1432275
NUDT521.842e-05rs10906083
SIGLEC512.188e-05rs8104955
ACRBP0无显著 cis-eQTL
APOE★★ 0无显著 cis-eQTL
ERMAP0无显著 cis-eQTL
NOTCH20无显著 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=686GTEx 肾 n=73 (p_beta)读数
PAM165 条,p=4.47e-160.984落空
WARS131 条,p=1.67e-160.543落空
IFNAR193 条,p=6.88e-120.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 ★ 另两个写脚本前查实的列语义坑

#事实(实测)后果
★★ 1maf 不是效应等位频率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 / rsidrsid 缺失 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_nameWARS1,再次印证改名。

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,用它诊断功效是范畴错误
2negative_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★★ eafac/an 而非 maf§17.9:21% 的行两者不等,用 maf 会静默翻转频率
5基因坐标写死(Ensembl GRCh38),加断言 K12b:TSS 必须落在该基因变异跨度内且染色体一致目录窗为 TSS±1 Mb,故 TSS 必在跨度内;写死使脚本不依赖网络可复现
6K4 只设单阈值,不像 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 全部 diseaseMaculopathy, Nephropathy, Neuropathy, Retinopathy ← 6-3 用这个
coloc_tier=="strong" 过滤Maculopathy, Nephropathy, Retinopathy ← 我错用了这个

Neuropathy 只出现在非 strong 的候选对里,按 strong 过滤会把它整个丢掉 —— 于是 E4IFNAR1 × 神经病变的 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「断点续跑:文件存在=已完成」。 修法:先写 .buildingfile.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-21243120:42只有 3 个结局结局集漂移,丢了 Neuropathy(见 §17.11)
.commasplit-22034621:444 结局,但 rsID 拆逗号对接★ 自创无先例的做法,已撤回(见 §18.6)
无后缀(定版)22:194 结局 + 精确 rsID 匹配以此为准

18.1 ★★ PP.H4 ≥ 0.5 的全部结果(只有 3 条,全在全血)

组织基因结局PP.H4nsnpbpvs 血浆
bloodWARSRetinopathy0.92528−0.12222.59e-06opposite
bloodIFNAR1Neuropathy0.84212+0.41870.262不可比
bloodWARSNeuropathy0.66018−0.25074.25e-10不可比

**胫神经与肾皮质无一条 ≥0.5。**各自最高:

组织Top 3
nerve_tibialACRBP×Reti 0.318 · GALNT3×Reti 0.318 · SIGLEC5×Neuro 0.272
kidney_cortexNOTCH2×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,6840.946原始发现
6-5 全血GTEx n=670(独立队列)0.8421复现了
6-5 胫神经GTEx n=5320.0584不外推到神经组织
6-5 肾皮质GTEx n=730.0366该层整体功效不足,不作解读

★★ 两条都重要:

  1. 不是 eQTLGen 的偶然 —— 换独立队列仍 0.84
  2. 但它是「血液层」现象,不是「神经组织」现象 —— 且这是有功效的阴性IFNAR1 在胫神经是强 eGene(p_beta=1.34e-11),不是测不出

18.3 ★★★ NUDT5 的头号结论被这一层挑战

计划里写着「NUDT5 是唯一两个独立分子层、共定位均 >0.97、方向一致的候选」。

数据NUDT5 × Retinopathy PP.H4
血浆 pQTLUKB-PPP0.9773
6-3 全血eQTLGen n=31,6840.9821
6-5 全血GTEx n=670★★ 0.0328
6-5 胫神经GTEx n=5320.2479
6-5 肾皮质GTEx n=730.0418

同为全血,两个独立数据集给出 0.982 与 0.033。

功效解释成立吗?部分成立,但不够干净:

GTEx 全血 p_beta全血 cis 变异 / p<0.05GTEx 全血 PP.H4
WARS2.81e-26(强 eGene)3,609 / 7980.925
IFNAR10.0392(弱 eGene)3,099 / 2570.842
NUDT50.0492(弱 eGene)4,821 / 3580.033

→ 样本量差 48 倍(670 vs 31,684),NUDT5 在 GTEx 全血只是勉强的 eGene, 功效不足是合理解释;但 IFNAR1 的 eGene 强度相当(0.0392 vs 0.0492)却出了 0.840, 所以不能把功效当成完整解释

写作必须处理这一条,不得回避NUDT5 的「两层共定位」证据现在有一个 独立血液队列的不复现。要么补更有功效的血液 eQTL 队列,要么在正文明确写出这个不一致。

18.4 ★ APOE × 肾病:三层全阴,且原因已判定

组织PP.H4p_betanegative_reasonSusztak n=686 显著对
kidney_cortex0.04990.397no_cis_eqtl_even_at_n6860 条
blood0.02950.5670
nerve_tibial0.00505.46e-170

→ 血浆层 APOE×肾病 coloc 0.998,但肾组织层在 n=686 都检不出 cis-eQTL。 措辞用「未检出显著 cis-eQTL」,不得写「无 eQTL」(Susztak 只发布显著对)。

★ 注意胫神经那一行:APOE 在胫神经是极强 eGenep_beta=5.46e-17)却仍得 0.0050 —— 这是有功效的阴性,比肾层那条更有信息量。

18.5 MR FDR<0.05 的 4 条(★ 与共定位分开读)

有 MR 估计的行 112 / 132(其余为该基因在该组织无合格工具)。

组织基因结局nsnpbpFDRPP.H4
bloodWARSNeuropathy8−0.25074.25e-103.83e-090.660
bloodWARSRetinopathy8−0.12222.59e-062.33e-050.925
nerveNUDT5Retinopathy−0.09835.77e-056.35e-040.248
kidneyACRBPRetinopathy−0.05291.45e-031.16e-020.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:

形态占比
单个 rsNNN94.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 / rs397806081beta 相同)→ 必须靠一个任意的 break 防止同一行被当成两个独立工具
4共定位侧实测零收益:132 行中 PP.H4 最大只动 0.0071WARS×Reti 0.9181→0.9252),无一条跨过 0.8 判定线

⚠️ 我上次那句「实测零收益」范围说小了,在此更正:零收益只覆盖共定位 PP.H4MR 侧不是零 —— 它决定了 ERMAP 在全血有没有工具(见 §18.5 的 tip)。 事实与读法要分开写,不能拿「结果差不多」当撤回的理由。

★ 现在全流程的口径(已逐脚本实查,不是凭印象)

取 FinnGen 的实现匹配方式
第 1–5 步主线R/adapters/read_tsv.Rgrep -Fw -fR/05_harmonise.Rout_std[SNP %in% snps]精确
6-3 全血 eQTL 51_crossomics_eqtl.Rsource(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)
WARS0.9560.8090.7160.9250.0550.064
IFNAR10.940阴性0.946(神经)0.842(神经)0.0580.057
NUDT50.9770.0350.9820.0330.2480.042
ERMAP0.967★ 0.949(黄斑,反向)0.0100.060
APOE0.998(肾)0.0320.0140.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=53247,4824,716(9.9%1.02e-031.07e-3138
全血 n=67036,0382,878(8.0%8.13e-045.60e-1720
肾皮质 n=7353,122308(0.58%5.76e-058.38e-111

正确的说法

  • :被保留的性状里,所有变异都在(含大量 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 里争议中心的那个蛋白,这条假阳性会直接写进正文。

定下的规则

  1. 归因只用 gene_id(数据生产方的注释),坐标法一律不得用于归因
  2. 坐标法只用于覆盖账(判「该区间有没有簇被检验过」),且产出列必须带 _positional 后缀,并由断言禁止其进入分析列
  3. ★ 这条限制要写进 Limitations:重叠基因的剪接归因依赖 leafcutter 的注释选择

19.3 ★★ 本层的可分析面比计划预想的小得多(实查,非估计)

gene_id 直连,11 个基因在三个组织 .cc 中的命中:

组织可分析的基因行数 / 簇 / 内含子
胫神经 n=532WARS · ERMAP · PAM23,436/1/3 · 6,223/1/1 · 23,158/2/3
全血 n=670WARS · ERMAP · PAM · ACRBP7,787/1/1 · 12,508/1/2 · 7,781/1/1 · 7,685/1/1
肾皮质 n=73WARS5,228/1/1

可分析的 基因×组织 组合共 8 个(满格是 33 个)。

★★★ 必须提前说明:6-6 回答不了 §18.3 的 NUDT5 争议

NUDT5IFNAR1 在任何组织的 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
S2NUDT5/IFNAR1 在三组织全部为「测了但不显著」已实查确定,写成断言 Q7
S3WARS 至少在一个组织出 PP.H4 ≥ 0.5预测。理由:它 eQTL 层全血 0.925,且 sQTL p_beta=5.93e-92(全血)/1.28e-95(神经)极强
S4ERMAP 在神经或全血出 PP.H4 ≥ 0.5预测。其 sQTL p_beta 8.42e-17 / 6.69e-18,若成立会是本层第二条结论
S5肾层(只有 WARS出不了 ≥0.5预测,理由同 6-5:n=73 功效不足
S6rsID 命中率落在 70–85%预测。6-5 实测 78.7%,本层变异集不同故给区间不给点值
S7会出现「MR 显著但 PP.H4 低」的行预测。可分析组合只有 8 个 → FDR 分母极小,须按 LD 混杂读,不得当阳性

19.5 验收断言(抓不到即 stop();★ 未全过只准写 .partial 且非零退出)

沿用 6-5 的 N1N4 / O1O2 / K12aK12b(同一批写死的 GRCh38 坐标),新增:

断言内容
Q1★★ 归因只用 gene_id:任何来自坐标法的列必须带 _positional 后缀,且不得出现在分析列中
Q2覆盖表 33 行(11 基因 × 3 组织),四态齐全无 not_queried
Q3coloc_input 全部 = cc_subset不是 6-5 的 all_nominal_cis),元数据同步记
Q4★ 分析单元留痕:每行必须带 n_clusters / n_introns / lead_intron
Q5★★ 硬断言:可分析的 基因×组织 组合恰为 §19.3 那 8 个,接线错立刻暴露
Q6rsID 命中率 ≥ 70%(沿用 6-5 的实测基线阈值,不再凭先验拍
Q7★★ NUDT5/IFNAR1 全部组织的 coverage_status 必须是 no_significant_sclusternegative_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.9400.7670.714
Neuropathy0.7220.7140.710
Nephropathy0.3570.3430.349
Maculopathy0.2740.1350.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.50.767WARS×Reti)★★ 我把 6-5 得到的「肾 n=73 功效不足」笼统外推到了本层。实为:n 小是平均意义上的功效不足,单个极强的分子性状照样能共定位WARS 肾 sQTL p_beta=1.3e-11)。「该层功效不足」不能用来预判该层里最强的那个信号

命中的三条:S3WARS 至少一组织 ≥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 ★ IFNAR1NUDT5有功效的阴性,不是没数据

基因胫神经全血肾皮质
IFNAR13 簇被检验,最佳 p_beta 0.5471 簇,0.4453 簇,0.495
NUDT55 簇,2.92e-07(★ 见下)4 簇,0.234 簇,0.0496

三个组织全部记 coverage_status = no_significant_sclusternegative_reason = tested_not_significant(断言 Q7 守住,不得写成「未测」)。

可写的结论IFNAR1 在 6-3/6-5 有明确的表达层共定位(全血 eQTLGen 0.946 / GTEx 0.842),而剪接层三个组织全阴其遗传信号定位在表达调控而非剪接。 这是对前两层的机制层面界定,不是重复。

⚠️ NUDT5 胫神经那个 2.92e-07 是 _positional(坐标重叠),不是它的 —— 该簇的 gene_idSEC61A2。详见 §19.2。产物里这个数只出现在覆盖表, 断言 Q1 保证它进不了分析表。

NUDT5 的 0.982 vs 0.033 之争,本层给不出任何新证据(跑前 §19.3 已声明)。

20.4 MR FDR<0.05 的 14 条(★ 与共定位分开读)

前 6 条(按 FDR):

组织基因结局nsnpbFDRPP.H4
kidneyWARSRetinopathy1+0.08991.81e-050.767
nerveWARSNeuropathy15−0.09981.03e-040.722
nervePAMRetinopathy8−0.06455.29e-040.198
kidneyWARSNeuropathy1+0.13487.25e-040.714
nerveWARSRetinopathy5+0.12422.44e-030.940
kidneyWARSNephropathy1+0.09733.11e-030.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 格)

状态格数含义
analyzable8有 FDR 显著的剪接簇 → 做了 MR + 共定位
no_significant_scluster20测了但不显著tested_not_significant
not_tested5该基因无任何内含子簇被注释/检验SIGLEC5 三组织全无、LACTB2 神经与肾无

SIGLEC5 在三个组织都没有内含子簇 —— 与它在 6-4a/6-5 屡次「非 eGene」一致, 写 Limitations 时应合并成一句:该蛋白在多个分子层均缺可用工具

20.6 ★ 实现上值得记的两点

  1. 同一基因在不同结局可能选中不同内含子lead_intron 列留痕)。 如 WARS 在胫神经:Reti/Macu 选 14:100361921:100375283,Neuro/Neph 选 14:100369258:100374136 —— 同簇不同内含子。这是设计使然(按 PP.H4 取最大), 但报告时必须带 lead_intron,否则读者无法复现是哪条。
  2. .cc.all 对应物,故 coloc_inputcc_subset(断言 Q3)。 与 6-5 的 all_nominal_cis 不是一回事,跨层比较时不能混。

20.7 候选格局(六层之后)

蛋白血浆视网膜全血 eQTLGen全血 GTEx胫神经 eQTL肾 eQTL★ sQTL(最高)
WARS0.9560.8090.7160.9250.0550.064★★ 0.940(神经)三组织全 ≥0.71
IFNAR10.940阴性0.9460.8420.0580.057三组织全阴(有功效)
NUDT50.9770.0350.9820.0330.2480.042无显著剪接簇
ERMAP0.9670.9490.0100.0600.374
APOE0.998(肾)0.0320.0140.063无显著剪接簇

★★ WARS 的领先进一步扩大:血浆 + 视网膜 + 两个全血 eQTL 队列 + 三个组织的剪接层, 是唯一在两种分子机制(表达与剪接)上都有共定位的候选。


相关:共定位方法选型 · 第 5 步实录 · 分析流程定稿 v4 · 第 7 步定稿与作图计划

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