09 · 批量筛查与出图脚本（代码详解）

对应脚本：analysis/01_screen.R、legacy/screen_t2d.R、legacy/batch_screen.R、analysis/02_visualize.R、analysis/09_report.py（均在项目根目录，是"入口驱动"，调用 R/ 引擎）

前面 01–08 讲的是 R/ 引擎模块；本页讲批量筛查这条线的驱动脚本——它们把引擎串起来，对"一个蛋白 × 几十个病"做大规模扫描。

---

一、analysis/01_screen.R（主力：manifest 驱动的全结局筛查）

用途

读 UKB-PPP cis-pQTL 工具 + outcome_manifest.csv 结局清单，逐结局跑单 SNP Wald，输出显著蛋白表、F 值、跨结局重叠、发文章级汇总。支持断点续跑。

代码分段

① 读配置 + 断点续跑开关

rel <- commandArgs(TRUE)[1]         # R9 或 R13
cfg <- load_config("config.yaml")
MAN <- fread("outcome_manifest.csv")[release == rel]   # 只取该版本的结局
RESUME <- isTRUE(cfg$project$resume) && !("force" %in% ...)  # force 参数强制重算

② 准备暴露（cis 工具）+ hg37 位置键 + 注释

e_std <- read_gwas(exp_unit, cfg, "exposure")     # 读标准化 cis-pQTL
exp_fmt <- to_tsmr_exposure(e_std, exp_unit)       # 转 TwoSampleMR 格式
raw[, hg37key := sub(":imp:v1.*", "", `Variant ID (...hg37...)`)]  # 供 position 匹配
ann <- unique(raw[, .(SNP=rsID, protein, gene, gene_consequence, biotype, ...)])

③ 每工具算 F/R²（相关性假设，发文章必备）

add_f <- function(h) {
  maf <- pmin(eaf, 1-eaf); num <- 2*beta^2*maf*(1-maf)
  R2 := num/(num + 2*N*maf*(1-maf)*se^2); Fstat := R2*(N-2)/(1-R2)
}

④ 主循环：逐结局（缓存复用 or 新算）

for (每个结局 ou in MAN) {
  if (RESUME && 缓存存在) ss <- fread(缓存)           # 断点续跑：直接读
  else {
    h <- (position匹配 ? read_pos_outcome : rsid读取+harmonise_pair)  # 两条路径
    ss <- mr_singlesnp(h) → 算FDR → 并入 add_f(h) + 两侧等位 + ann
    fwrite(ss, 缓存)
  }
  sig <- ss[FDR < 0.05]; 写 <结局>_significant.csv
  记录 sig_list（重叠用）+ summary_rows（病例/工具/minF/显著数）
}

⑤ 汇总：跨结局重叠矩阵 + SUMMARY_publication.csv

overlap_matrix.csv    # 每蛋白在几个病里显著
SUMMARY_publication.csv  # 每结局: ncase/ncontrol/n_instruments/minF/medianF/n_significant

position 匹配（Mahajan 无 rsID）由 read_pos_outcome() 用 chr:pos:等位 的 hg37 键对齐。

---

二、legacy/screen_t2d.R（T2D 位置匹配单跑）

还原原始实现 2_diabete_validation 的 Mahajan 处理：暴露 hg37 变异 ID 拼 chr:pos:A0:A1，与 Mahajan 的 chr:pos:NEA:EA 匹配（都 build37，无需 liftover），再 harmonise + Wald。逻辑已并入 analysis/01_screen.R 的 position 分支，此脚本保留作独立复现。

---

三、legacy/batch_screen.R（旧版，对照用）

analysis/01_screen.R 的前身：5 个结局写死、参数硬编码。用来和原始实现的结果逐行核对（已验证糖网 24/肾病 9/神经 5 等完全一致）。新分析用 analysis/01_screen.R。

---

四、analysis/02_visualize.R（出图：UpSet / 森林 / 网络）

用途

读某版本的 screen_*/ 结果，出三张图（同时 PNG 600dpi + 矢量 PDF）。

代码要点

DPI <- 600
PDFDEV <- if (capabilities("cairo")) cairo_pdf else pdf   # 矢量 PDF 设备
save_both <- function(g, name, w, h) {                     # 每图存 png+pdf
  ggsave(name.png, g, dpi=DPI); ggsave(name.pdf, g, device=PDFDEV)
}
# 1) UpSet: UpSetR::upset(overlap_matrix)          → 蛋白跨病重叠
# 2) 森林: 跨>=3病共享蛋白的 OR/CI (log 轴)
# 3) 网络: ggraph, 红=风险OR>1 / 蓝=保护OR<1

---

五、analysis/09_report.py（汇总成 Excel 多表）

用途

把 R9/R13 的 screen_*/ 结果汇总成发文章用的多表 Excel + CSV。

代码要点

# 1) 读所有 *_significant.csv 合并 → 主明细表(含 F/等位/OR/CI/P/FDR/注释)
# 2) 读 SUMMARY_publication.csv → 结局汇总
# 3) 读 overlap_matrix.csv → 重叠矩阵
# 4) 写多表 Excel (openpyxl) + 各 CSV 到 results/REPORT/

输出：MR_publication_tables.xlsx（多 sheet）+ table_S_*.csv。

---

相关：批量筛查（怎么用） · 下游分析脚本详解 · 代码结构详解