主题
08 · ARIC 队列 cis-pQTL 系统性取数与接入
执行日期:2026-08-07 状态:✅ 完成(§1–§7 写于运行之前,§8 为运行后回填) 性质:数据取数步,不产生任何判定。它是 第 4X 步的前置条件。
为什么插在 4X 之前
4X 要区分「平台问题」还是「人群问题」,靠的是 deCODE vs ARIC: 二者同为 SomaScan 适配体平台、但人群不同(冰岛 vs 美国)。
只有 deCODE 的话,与 Olink 比较时平台和人群是一起变的, 方向不一致可以被推给人群差异,归因不了。
而现状是:ARIC 只有 3 个蛋白的手工结果、没有脚本, 5 对待查里只有 APOL1 这一维可用。不补齐,4X 有四对是瞎的。
★ 编号顺延
原计划 08 是 4X。因本步插入,4X 顺延为 09。
一、七要素 · 数据
| 项 | 内容 |
|---|---|
| 队列 | ARIC(Atherosclerosis Risk in Communities),欧洲裔子集 EA |
| 样本量 | n = 7,213 |
| 平台 | SomaScan v4 适配体 |
| 适配体数 | 4,657 |
| 文献 | Zhang J, Chatterjee N, et al. Nat Genet 2022;54:593–602(PMID 35501419) |
| 门户 | http://nilanjanchatterjeelab.org/pwas/ (Shiny 应用) |
| 实际取数地址 | https://jh-pwas.s3.amazonaws.com/results/EA.zip(静态 S3,不经 Shiny) |
| 对照表 | https://jh-pwas.s3.amazonaws.com/results/seqid.txt(SeqId → UniProt / 基因 / 染色体 / TSS) |
版本与口径(★ 必须写进 Methods)
| 项 | 值 | 为什么重要 |
|---|---|---|
| 生成日期 | 2021-09-28(S3 Last-Modified 2023-01-23) | 对应 bioRxiv v2 修订版,非初版 |
| cis 窗口 | TSS ±500 kb | ★ 与我们第 1 步的 cis 定义完全一致,无需重新划窗 |
| 基因组版本 | GRCh38 | ★ 与主线一致,不需要 liftover |
| 关联模型 | PLINK2 --glm linear,ADD | 列格式即 PLINK2 标准 |
| 体积 | EA.zip 352,021,558 B;解压 1,124,222,148 B(4,657 文件) |
门户说明里的一句话必须照抄
"Updated summary statistics … using +/-500Kb cis window … based on the latest revision"
即:当前 S3 上的就是 ±500 kb 版本。早期版本窗口不同,若日后有人重下要核对此说明。
二、七要素 · 方法
1. 服务器(美国机房)下载 EA.zip → unzip -t 完整性校验
2. seqid.txt 把我们的存活蛋白映射到 SeqId(MICB_MICA 拆成 MICB / MICA)
3. 只抽取需要的 .glm.linear(不解压全部 1.05 GB)
4. 回传本地 D:\mrdata\aric\
5. 写脚本做等位对齐 + 质控,产出统一查询表★★ 唯一的技术陷阱:PLINK2 的 A1 不是 ALT
实测 IFNAR1 哨兵:
#CHROM POS ID REF ALT A1 A1_FREQ BETA P
21 33353501 rs914142 G A G 0.262581 0.391707 8.04357e-99A1 = G = REF,BETA 与 A1_FREQ 都是针对 A1 而不是 ALT。
我们 Olink 侧的效应等位是 A,所以对齐后:
beta→ −0.391707(变号)EAF→ 1 − 0.262581 = 0.737419
谁假定 A1 == ALT,符号就整个反了
这与 T1D 源等位频率反转事故是同一类错误: 参照等位搞错不会报错,只会安静地把结论反过来。
→ 脚本里必须显式读 A1 列,并对每个位点断言 A1 ∈ {REF, ALT},不满足即 stop()。
三、七要素 · 预期(★ 本节写于运行之前)
3.1 覆盖率预期(已用 seqid.txt 实测,非猜测)
第 3 步存活 29 个蛋白,映射到 ARIC 适配体:
| 结果 | 数量 | 明细 |
|---|---|---|
| 有适配体 | 20 / 29 = 69% | 其中 APOE APOL1 NOTCH2 各 2 条适配体 |
| 无适配体 | 9 | ACRBP ATP6V1G2 BCL2L15 BTN3A2 HCG22 LACTB2 LTB TPPP3 TRIM40 |
★ MICB_MICA 值得单独说:Olink 是合并探针(一条探针测两个蛋白), 而 SomaScan 分开测——MICB = SeqId_5102_55,MICA = SeqId_2730_58。 这对「两个适配体互相矛盾」的调查本身就是证据,不是麻烦。
3.2 对 4X 五对的预期
| 对 | ARIC 可用? | 说明 |
|---|---|---|
APOL1 × Maculopathy | ✅ 2 条适配体 | 阳性对照 |
MICB_MICA × Maculopathy | ✅ MICB + MICA 分开 | |
MICB_MICA × Nephropathy | ✅ 同上 | |
AOC1 × Retinopathy | ✅ SeqId_15486_126 | |
TRIM40 × Retinopathy | ❌ 无 | ★ 但它走结局异质性路由,本就不需要 ARIC |
→ 预期结论:4a 路由的 4 对全部可评估,ARIC 这一维在 4X 上不留空白。
3.3 缺口结构的预期
9 个未覆盖里,ATP6V1G2 BTN3A2 HCG22 LTB TRIM40 都在 MHC 区。
→ 预期 ARIC 的缺口与 deCODE 同向、都集中在 MHC。 若成立,说明这是平台层面的系统性缺失,而非随机漏测,应写进 Limitations。
一个待验的巧合
deCODE 覆盖也是 20/29 = 69%。是不是同一批 20 个蛋白,本步要查。 若高度重叠 → 两个 SomaScan 队列共享同一套适配体清单, 则「换人群」的独立性比看起来弱,这会削弱 4X 的归因力,必须如实写。
3.4 ★★ 验收断言(抓不到即 stop())
本步的正确性有现成的检验对象:2026-08-01 手工算过 3 个蛋白。 系统性取数若正确,必须逐位复现它们。
| 蛋白 | 适配体 | 变异 | 对齐后 beta | p | 其他 |
|---|---|---|---|---|---|
IFNAR1 | SeqId_9183_7 | rs914142 | −0.391707 | 8.04357e-99 | EAF 0.737419;cis SNP 2,413 |
APOL1 | SeqId_11510_31 | rs136168 | +0.461581 | 7.81792e-111 | |
APOL1 | SeqId_9506_10 | rs136168 | +0.008252 | 0.694272 | 第二适配体无信号 |
ERMAP | SeqId_8631_13 | rs11210710 | −0.049519 | 0.00292733 |
这一条是 4a 教训的直接落地
4a 连翻四版,根因之一是「手里有阳性对照却没当验收标准」。
本步先把验收标准写死在这里,脚本里写成断言。 复现不出 = 新取数错了(或旧手工数字错了),两种都必须停下来查,不许"看着差不多"就过。
3.5 明确不做的事
- ❌ 不做任何 MR、不出任何判定 —— 本步只取数与对齐
- ❌ 不下 AA(非洲裔)子集 —— 我们主线是欧洲裔,混入会引入人群分层
- ❌ 不碰 PWAS 预测模型(
PWAS_EA.zip)—— 那是另一套方法,与本课题无关
四、七要素 · 文件位置(预期产物)
| 产物 | 路径 |
|---|---|
| 原始压缩包(服务器暂存) | myserver:/root/aric/EA.zip |
| SeqId 对照表 | D:\mrdata\aric\seqid.txt |
| 抽取的汇总统计子集 | D:\mrdata\aric\EA\SeqId_*.PHENO1.glm.linear |
| 取数脚本 | analysis/42_aric_fetch.R(编号待定) |
| 统一查询表 | results/external_replication/aric_cis_lookup.csv |
产物约定照旧三条(违反即 stop()): CSV 取值一律 ASCII / 不许出现内部流程编号 / 内部结局 ID 不上图。
五、七要素 · 风险与已知局限
| 风险 | 说明 | 处置 |
|---|---|---|
★ A1 ≠ ALT | 见 §2 | 显式读 A1 + 断言 |
| 功效低 | n=7,213,远小于 deCODE(≈3.5 万) | ★ 「无信号」≠「反驳」 |
| F 统计量不达标 | ERMAP 在 ARIC 上 p=0.0029、F=8.9 < 10 | 按既定规则写「无有效工具,此检验做不了」,不写阴性 |
| 适配体清单可能与 deCODE 高度重叠 | 见 §3.3 | 本步查明并如实写 |
| MHC 区缺失 | 见 §3.3 | Limitations |
ERMAP 那两格的措辞是已定的,别改
分母不合格时 Wald 比值会被放大成荒谬数值(本课题在 FABP6 与 PON1 上见过 −72 与 +82)。 所以正确表述是「另一平台没有有效工具,此检验做不了」—— 既不能说复制成功,也不能说被推翻。
六、七要素 · 正文还是补充
暂定:本步属方法与数据来源,进 Methods + 补充。 覆盖率与 MHC 缺口结构进 Limitations。 最终归属第 7 步统一复核。
七、执行前置检查(已完成)
| 检查 | 结果 |
|---|---|
| S3 三个文件可达 | ✅ 全部 HTTP 200,Accept-Ranges: bytes(可断点续传) |
| 服务器下载速度 | ✅ ≈13 MB/s(本机直连仅 ≈61 KB/s,慢 215 倍) |
| 服务器磁盘 | ⚠️ 仅 14 G 可用 —— 故只抽子集、不全量解压 |
| 本地磁盘 | ✅ D: 余 6.2 T |
EA.zip 完整性 | ✅ unzip -t 通过 |
| 结构确认 | ✅ EA/SeqId_<id>.PHENO1.glm.linear,4,657 个 |
| 列格式确认 | ✅ PLINK2 标准 14 列 |
| 与旧台账勾稽 | ✅ 三个手工蛋白的适配体号、cis SNP 数、beta、p 全部对上 |
八、运行结果
执行日期:2026-08-07 状态:✅ 完成 脚本:analysis/42_aric_lookup.R
8.1 取数
| 项 | 结果 |
|---|---|
| 抽取适配体 | 24 条 / 20 个蛋白(APOE APOL1 NOTCH2 各 2 条) |
| 体积 | 7.7 MB(打包 2.36 MB) |
| 传输校验 | ✅ 两端 md5sum 一致 308f5e04de5673406535e43a3913ea48 |
| 落地 | D:\mrdata\aric\EA\ |
| 产物 | results/external_replication/aric_cis_lookup.csv(33 行 × 29 列,非 ASCII 字节 0) |
查询状态:ok 23 · sentinel_absent_in_aric 1 · no_aptamer_on_platform 9 等位对齐:a1_equals_effect_allele 20 · a1_equals_other_allele_flipped 3(即 §2 那个陷阱真实发生了 3 次) 有效工具:F≥10 的 17 条、不足的 6 条 频率质控:eaf_diff 中位数 0.0065、最大 0.0535,无一条超 0.08
8.2 ★★ 验收断言:7/7 全部 PASS
| # | 断言 | 结果 |
|---|---|---|
| V1 | IFNAR1/SeqId_9183_7/rs914142 = −0.391707, p=8.04e-99, eaf=0.7374, cis=2413 | ✅ |
| V2 | APOL1/SeqId_11510_31/rs136168 = +0.461581, p=7.82e-111 | ✅ |
| V3 | APOL1/SeqId_9506_10/rs136168 = +0.008252, p=0.694 | ✅ |
| V4 | ERMAP/SeqId_8631_13/rs11210710 = −0.049519, p=0.00293 | ✅ |
| V5 | 所有位点 A1 ∈ {REF, ALT} | ✅ |
| V6 | 24 条适配体 / 20 个蛋白 | ✅ |
| V7 | 29 个存活蛋白全部有交代 | ✅ |
★ 2026-08-01 的三个手工数字被系统化流程逐位复现, 说明旧结论的数据来源确实是本页记录的这个源,且当时的手工对齐没做错。
8.3 预期 vs 实际(逐条对账)
| §3 的预期 | 实际 | 判定 |
|---|---|---|
| 覆盖 20/29 = 69% | 20/29 = 69% | ✅ 命中 |
| 缺 9 个,名单已列 | 完全一致 | ✅ 命中 |
| 4X 的 4 对全覆盖 | APOL1✅ MICB+MICA✅ AOC1✅ | ✅ 命中 |
TRIM40 无、但不需要 | 无适配体;走结局路由 | ✅ 命中 |
| 缺口偏 MHC | MHC 6/11 = 55% vs 非 MHC 14/18 = 78% | ✅ 命中 |
| deCODE 与 ARIC 是不是同一批 20 个? | ★★ 完全相同 | ❌ 落到最坏那一种 |
★★ §3.3 那个「待验的巧合」不是巧合
同口径对比(都按「平台上有没有这个蛋白」):
deCODE 20 个 · ARIC 20 个 · 交集 20 个 · 各自独有 0 个。
→ 换人群不带来任何覆盖增益,9 个缺口来自 SomaScan 面板本身。
这必须写进 Limitations:ARIC 与 deCODE 在「有没有这个蛋白」这一维上并不独立。 ARIC 的价值仅在于同平台、不同人群的效应量对照, 不能用它来补 deCODE 的覆盖缺口,也不能把两者当成两次独立抽样。
8.4 运行中查出并修掉的四个缺陷
| # | 缺陷 | 后果(若不修) | 处置 |
|---|---|---|---|
| 1 | 状态名 sentinel_not_in_cis_window 是错的 | CLPS 哨兵 rs2766588(6:35,785,450) 实际落在 ARIC 窗口 35,297,440–36,296,774 之内,只是 ARIC 没有这个变异。原标签会把「队列缺变异」误报成「窗口口径不一致」,把错误归因直接带进 4X | 拆成 sentinel_absent_in_aric / sentinel_outside_cis_window,并新增 in_cis_window 列 |
| 2 | fread 把 P 读成 character | 下游 log10(p) 静默出错(本次就是它让断言崩了) | 实测确认取值全部是合法数值字面量,改显式 as.numeric + 断言不产生 NA |
| 3 | 覆盖对比两个口径混用 | ARIC 用「有适配体」(20)、deCODE 用「拿到 cis 工具」(15),得出**「ARIC 覆盖更广」的假象** | 统一为「平台上有无该蛋白」,结果 20 vs 20 |
| 4 | p 下溢未标记 | MICA 的 p_aric = 0(PLINK2 对极显著位点输出 0),下游 -log10(p) 得 Inf 而不报错 | 新增 p_underflow 列 |
★ 第 1 条和第 3 条都属于**「结论看起来合理、但归因是错的」**—— 正是最难靠通读发现、只能靠对账发现的那类。
8.5 交付给 4X 的数据
| 对 | ARIC 结果 | 可用性 |
|---|---|---|
APOL1 × Macu | SeqId_11510_31 +0.4616(F=518) / SeqId_9506_10 +0.0083(F=0.15) | ✅ 主适配体可用 |
MICB_MICA × Macu/Neph | MICA SeqId_2730_58 −0.8468(F=3035) / MICB SeqId_5102_55 +0.2492(F=190) | ✅ 两条都可用 |
AOC1 × Reti | SeqId_15486_126 −0.4908(F=269) | ✅ |
TRIM40 × Reti | 无适配体 | ✅ 本就不需要(走结局异质性路由) |
★ 一个留给 4X 判、本步不下结论的现象
同一个变异 rs3132467,在同一队列、同一平台上: MICA 适配体给 −0.847、MICB 适配体给 +0.249——方向相反。
这与 4a 报出的「两适配体互相矛盾」是同一件事, 但归因是 4X 的职责,本步只交付数据,不判。
8.6 意外收获:PAV 注释可能不用外部工具
data/exposure/protein_info.csv 自带 Annotated gene consequence / CADD_phred / SIFT / PolyPhen 四列——正是 4X 判「哨兵是不是 PAV」要用的注释。
→ 4X 开工前先查这张表,可能省掉一次外部注释取数。待 4X 核实其完整性。
8.7 遗留
analysis/_diag_aric.R_diag_aric2.R_diag_cov.R是本次的诊断脚本,待确认后清理- ARIC AA(非洲裔)子集未下载 —— 主线是欧洲裔,按 §3.5 明确不做