Skip to content

共定位证据分档与 SuSiE 降级

两项改动,同一天做的,起因都是与 baseline 文章 Yuan 2023 对标: ① 共定位加中间档(抄它的,我们原来缺) ② SuSiE 降级为敏感性检查(不抄它的,我们有数据说明为什么不能那么用)


一、共定位加中间档

为什么要加

原来是二值:PP.H4 ≥ 0.8 算共定位,否则算阴性。这会把两件完全不同的事写成同一个"阴性":

情形后验分布特征真实含义
APP.H3 高两个性状各有各的因果变异 —— 确实不共定位
BPP.H0 / H1 / H2 高谁都没有明确信号 —— 功率不足,无法判断

审稿人问"这 76 对阴性里有多少是真阴性",二值表答不出来。

实测分布(coloc_abf_pairs.csv,142 对)

定义对数占比
strongPP.H4 ≥ 0.83223%
medium0.5 < PP.H4 < 0.82820%
weak0.3 < PP.H4 ≤ 0.564%
nonePP.H4 ≤ 0.37654%

中间档这 28 对的 PP.H3 全部在 0.03–0.47,多数 < 0.15 —— 全属上表的情形 B(功率不足),不是情形 A。二值阈值把它们错误地表述成了"不共定位"。

实现

分档函数只有一份实现,放在 R/utils.R

r
coloc_tier <- function(pp_h4, strong = 0.8, medium = 0.5, weak = 0.3) {
  stopifnot(weak < medium, medium < strong)
  data.table::fifelse(is.na(pp_h4), NA_character_,
    data.table::fifelse(pp_h4 >= strong, "strong",
      data.table::fifelse(pp_h4 > medium, "medium",
        data.table::fifelse(pp_h4 > weak, "weak", "none"))))
}

阈值走 config,不硬编码:

yaml
coloc:
  pp_h4_threshold: 0.8    # 强支持下界 [Giambartolomei 2014]
  pp_h4_medium: 0.5       # 中等提示下界 [Yuan 2023 同口径]

调用点两处:

  • analysis/24_coloc_abf.R —— 写 coloc_abf_pairs.csv 时加 coloc_tier 列(下次全量跑 24 时生效
  • analysis/25_coloc_compare.R —— 写 coloc_compare_hypr_vs_abf.csv 时加 coloc_tier 列,并另出中间档补充表

新产物:results/hyprcoloc_R9_dm/coloc_abf_medium_supplement.csv(28 行)

25 只读 CSV、几秒跑完,已执行并核对:strong 32 / medium 28 / weak 6 / none 76,与手工计数一致。

中间档只进补充表,不进候选层

中间档里 ITGB7×Retinopathy (0.756)SIGLEC5×Retinopathy (0.731)TIGIT×Retinopathy (0.704) 正是反向 MR并发症方向反向显著的那批。 把它们抬进候选层会与反向 MR 的结论自相矛盾。 coloc_tier == "medium" 只作为补充材料如实披露,不参与 A/B 级分层。


二、SuSiE 降级为敏感性检查

结论

正文只放 HyPrColoc + coloc.abf;SuSiE 连同其一致性参数与收敛率写进补充材料的方法学局限,不作判据。

依据一:出结果的比例太低,且冲突无法裁决

尝试对数32
两侧均收敛12
给出 PP.H48(25%)
其中与 coloc.abf 改判4(50%)

改判的四对:

蛋白 × 结局coloc.abfSuSiE
APOE × T2D0.9990
PAM × Retinopathy0.9590
PAM × T2D0.9240
ACRBP × Retinopathy0.8670

对三个 A 级候选的实际影响是:ERMAP、IFNAR1 两侧都不收敛;APOL1 收敛且与 abf 一致(0.926 → 1)。

依据二:失配是定量的,且救不回来

susieR::estimate_s_rss() 给出 z 值与 LD 矩阵的一致性参数 s(0 = 完全一致,越大越失配)。在四个位点上实测:

位点s(蛋白 / 疾病)现状修法②修法①①+②
IFNAR1 × 黄斑0.925 / 0.717✗✗✗✗✗✗✗✗
ERMAP × 黄斑0.877 / 0.735✗✗✗✗✗✓✓✓(各 2 个 CS)
NOTCH2 × 糖网0.799 / 0.591✗✗✗✗✗✗✗✗
APOL1 × 黄斑(对照)0.676 / 0.393✓✓✓✓✗✓✓✓
  • 修法①=按 kriging_rss() 剔除离群 SNP(logLR > 2 且 |z| > 2
  • 修法②estimate_residual_variance = FALSE

三条结论:

  1. 修法①单用有害 —— 它把本来收敛的 APOL1 蛋白侧弄成不收敛。不可单独使用。
  2. 修法②单用无害,但救不回任何一个失败位点。
  3. ①+② 合用只救回四分之一(ERMAP)。IFNAR1、NOTCH2 依然全败。

s 不能单独用来预测能否救回

NOTCH2 的 s(0.799)低于 ERMAP(0.877),却没被救回来。 准确的说法是:s 大致排序了失配严重程度,但不决定能否收敛。

根因

1000 Genomes Phase 3 EUR 面板只有 503 个样本,而各区域 SNP 数在 1,361–4,477 之间,全部远超面板秩上限。蛋白侧与疾病侧在同一个蛋白 cis 区共用同一个 LD 矩阵,所以区域一旦失配,两侧同时崩 —— 这解释了为什么 12 对是"两侧都不收敛"。

换更大的参考面板(如 UKB ~337k)理论上可解,但成本高,且不影响任何现有结论,不做

写作口径

Methods 里可以照 Yuan 2023 的句式说明 SuSiE 的用途(那句话本身准确):

"传统共定位方法无法处理暴露与结局在同一区域共享一个以上因果变异的情形,故补充 SuSiE 方法。"

必须补上它没写的三样

  1. 尝试对数与收敛率(32 对 / 12 对两侧收敛 / 8 对给出 PP.H4)
  2. 一致性参数 s 的实测范围(中位 0.7–0.9)
  3. 与 coloc.abf 冲突时以 coloc.abf 为准,及其理由

不要抄 Yuan 2023 对 SuSiE 的用法

它全文 SuSiE 只产出 4 个数字(ARG1、THG1L 升为强支持;MANSC4 × 低血糖 0.67→0.96、× 糖网 0.63→0.85), 4 个全是把 abf 没到线的结果往上抬,一次都没报告过 SuSiE 把 abf 强结果打下来的情况。 而且它把这两个升档直接算进了摘要的头号计数(abf 强支持 9 个 + SuSiE 2 个 = 摘要的"11")。 我们的实测里降档占给出结论的一半,藏不住也不该藏。


三、改动清单

文件改动
config.yaml新增 coloc.pp_h4_medium: 0.5
R/utils.R新增 coloc_tier()(唯一实现)
analysis/24_coloc_abf.R写盘前加 coloc_tier 列(下次全量跑生效)
analysis/25_coloc_compare.Rcoloc_tier 列 + 输出中间档补充表 + 分档计数
results/hyprcoloc_R9_dm/coloc_abf_medium_supplement.csv新产物(28 行)

已回灌 F 盘副本并逐文件 md5 核对一致。

顺带修掉的一个静默失败

sync_f.sh 的代码同步段原用 cp -a。NTFS 不支持保留属主/时间戳,cp -a 整体返回非零、 又被 2>/dev/null 吞掉,结果是一条 ok 都不打印、代码根本没同步。已改为 cp -r 并显式报错中止。 本次正是靠"第 1 节零输出"发现的。


四、成对 UpSet 图(2026-07-31 补)

对标导师的 hypco_new.tiff,出**严格版(正文)+ 宽松版(补充)**两张, 脚本 analysis/38_figure_coloc_upset_pair.R,图在 public/r9/coloc_upset_{strict,loose}.png

★ 先纠正一个口径不一致

待办计划里那张集合大小表的「PP>0.7」列,与同页"做法"写的严格版定义不是一回事。 三个口径实测如下(results/hyprcoloc_R9_dm/hyprcoloc_results.csv,108 行):

口径定义蛋白对数各结局蛋白数
A完全不过滤4174糖网 19 / 黄斑 17 / T1D 13 / T2D 12 / 糖肾 9 / 神经 4
B仅 PP>0.7,不要求蛋白在簇内2756黄斑 15 / 糖网 14 / 糖肾 9 / T2D 7 / T1D 7 / 神经 4
CPP>0.7 且蛋白在簇内1319黄斑 5 / T2D 5 / 糖网 4 / T1D 4 / 糖肾 1 / 神经 0

计划书表里的「PP>0.7」列是 B,但"做法"一节写的严格版是 C

B 不能用作蛋白–疾病共定位图:蛋白不在簇内时,那是几个疾病彼此共定位、与该蛋白无关, 画进「蛋白 × 疾病」的 UpSet 会被读成蛋白与疾病共定位。正文图用 C

另一个必须处理的细节

traits == "None" 表示该位点没找到任何簇(posterior_probNA),共 66 行。 展开时若不剔除,会凭空多出一个叫 None 的"结局"、挂着 63 个蛋白。已剔除。

两张图

口径用途图注要点
coloc_upset_strictC(19 对 / 13 蛋白)正文蓝色;标题注明「PP>0.7 且蛋白在簇内」
coloc_upset_looseA(74 对 / 41 蛋白)补充灰色;标题直接写「未过滤,仅示分布,不构成共定位证据」

脚本内含断言:严格版展开结果必须与流水线既有的 hyprcoloc_shared.csvcoloc_support == TRUE 口径完全一致(实测 19 对全等), 否则两张图会与正文其他数字互相打架。

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