主题
10 · 下游分析脚本（代码详解）
对应脚本：analysis/03_phewas.R（PheWAS）、analysis/04_ldsc.R（遗传相关）、analysis/05_hyprcoloc.R（多性状共定位）
analysis/03_phewas.R/analysis/05_hyprcoloc.R把 MR 筛出的显著蛋白拿去做验证；analysis/04_ldsc.R不依赖 MR——它吃全量疾病 GWAS，是独立的前置分析。概念/结果/依赖关系见 下游分析；本页讲代码怎么实现的。
三个脚本对 R9/R13 的处理不一样（先看这个）
| 脚本 | 版本参数 | 为什么 |
|---|---|---|
analysis/03_phewas.R | 不分版本，R9+R13 合并跑 | PheWAS 查的是"这个 SNP 还关联多少表型"，是 SNP 自身属性，与结局用哪版无关 |
analysis/04_ldsc.R | R9 / R13 | 直接吃疾病 GWAS，版本不同结果不同 |
analysis/05_hyprcoloc.R | R9 / R13 | 同上，且目标蛋白取自对应版本的筛查结果 |
一、analysis/03_phewas.R（PheWAS 多效性扫描）
用途
把所有显著 cis-pQTL 工具 SNP 拿去 OpenGWAS 扫全表型，评估多效性。自动化原始实现 6_phewas 里脚本外的手动查询。
bash
Rscript analysis/03_phewas.R # 无版本参数：R9+R13 的显著 SNP 合并去重后一起查代码分段
r
# ① 收集显著 SNP + 蛋白注释（跨 R9+R13 的 *_significant.csv）
sig_files <- Sys.glob("results/screen_*/*_significant.csv")
snps <- 唯一 rs 开头的显著 SNP
# ② 分批查 OpenGWAS（每批20，失败重试3次）；有缓存则复用，加 force 重查
for (chunk in split(snps, 每20个)) {
r <- ieugwasr::phewas(variants=chunk, pval=QPVAL) # QPVAL=1e-5 只是"拉候选"的粗筛
}
# ③ ★ Bonferroni 校正：检验数 = SNP 数 × OpenGWAS 数据集数
n_db <- nrow(ieugwasr::gwasinfo()) # 取全库数据集数（失败则退化为返回结果里的数据集数）
n_tests <- length(snps) * n_db
BONF <- 0.05 / n_tests # 实测: 243 × 4687 → p < 4.39e-8（查询日期 2026-07-18）
ph[, bonferroni_sig := p < BONF]
ph[, p_bonf := pmin(1, p * n_tests)]
if (QPVAL < BONF) 告警("粗筛比 Bonferroni 还严，可能漏关联") # 自检
# ④ 汇总（以 Bonferroni 显著者为准）+ 脱靶明细
fwrite(phewas_all.csv) # 全部关联（含 bonferroni_sig / p_bonf）
fwrite(phewas_summary.csv) # 每 SNP: n_traits_bonf ★ / n_traits_raw / min_p
fwrite(phewas_offtarget_bonferroni.csv) # 脱靶明细，供阶段⑥方向矛盾判断粗筛阈值 ≠ 显著性阈值
QPVAL（默认 1e-5）只是去 OpenGWAS 拉候选用的服务端过滤。因为它比 Bonferroni 阈值 4.39e-8 宽松，Bonferroni 显著的关联必然都被拉回来了，不会漏。显著性一律以 bonferroni_sig 为准。
token 需先写入
~/.Renviron的OPENGWAS_JWT=...；ieugwasr::phewas自动读。
二、analysis/04_ldsc.R（疾病间遗传相关）
用途
用 GenomicSEM 的 munge+ldsc，算多个疾病两两的遗传相关 rg。还原原始实现 8_sensitivity 的 LDSC。
bash
Rscript analysis/04_ldsc.R R9 # → results/ldsc/（历史路径，保持兼容）
Rscript analysis/04_ldsc.R R13 # → results/ldsc_R13/纳入疾病（每版本的 id 在脚本 KEEP 里声明，文件名/病例数从 outcome_manifest.csv 读，不硬编码）：
| 版本 | 疾病数 | 纳入 |
|---|---|---|
| R9 | 8 | 糖网、黄斑、增殖DR、肾病、神经、AMD、T2D、T1D |
| R13 | 10 | 上述 + 糖网严格定义、湿性/干性 AMD 亚型（R9 没有） |
只纳入 FinnGen 自家的全量 sumstats——外部的 GCST（T1D）/ Mahajan（T2D）列格式不同，不进 LDSC。
代码分段
r
# ① 疾病清单：id 按版本声明，file/ncase/ncontrol 从 manifest 读
MAN <- fread("outcome_manifest.csv")[release == rel & id %in% KEEP[[rel]] & match_mode == "rsid"]
traits <- MAN[, .(id = sub("_finngen$", "", id), file, ncase, ncontrol)] # 短名好看
# ② 预处理：zcat + awk 取 SNP/A1/A2/beta/P（滤非 rs）
system("zcat 文件 | awk '... {print rsids, alt, ref, beta, pval, N}' > 输入.tsv")
# ③ munge 到 HapMap3
munge(files=输入.tsv, hm3="w_hm3.snplist", trait.names=..., N=...)
# ④ ldsc 遗传相关（观察尺度即可，rg 不受患病率影响）
res <- ldsc(traits=*.sumstats.gz, ld=eur_w_ld_chr, wld=eur_w_ld_chr, ...)
# ⑤ 提取 rg 矩阵：rg = S / (√diag %o% √diag)
fwrite(genetic_correlation_rg.csv); fwrite(h2_observed.csv)参考面板
eur_w_ld_chr+w_hm3.snplist放在/mnt/d/mrdata/ldsc/。
三、analysis/05_hyprcoloc.R（多性状共定位）
用途
对跨≥2 病的共享蛋白，用 UKB-PPP 区域数据 + FinnGen 疾病区域，跑 hyprcoloc 判断是否同一因果变异驱动多个疾病。还原原始实现 4_hyprcoloc。
bash
Rscript analysis/05_hyprcoloc.R R9 # → results/hyprcoloc_R9/
Rscript analysis/05_hyprcoloc.R R13 # → results/hyprcoloc_R13/
Rscript analysis/05_hyprcoloc.R R9 5 # 只跑前 5 个蛋白（测试用）
Rscript analysis/05_hyprcoloc.R R9 force # 忽略已完成的 part，全部重算纳入并发症结局（不含 T1D/T2D 基础病与外部数据集；文件名从 outcome_manifest.csv 读）：
| 版本 | 结局数 | 目标蛋白（跨≥2 病） |
|---|---|---|
| R9 | 12 | 24 个 |
| R13 | 9（含湿/干 AMD 亚型） | 31 个 |
数据来源
- 暴露区域：UKB-PPP Synapse 每蛋白
.tar（内含各染色体.gz，hg38） - 结局区域：全量 FinnGen 对应版本（hg38）
- 都 hg38 → 按
chr:pos位置匹配
断点续跑：每算完一个蛋白立即落盘
每个蛋白的结果当场写 _parts/<蛋白>.csv，最后再汇总成 hyprcoloc_results.csv。所以中途断了（关机/OOM/手滑 Ctrl+C）只丢当前这一个蛋白，重跑时自动跳过 _parts/ 里已有的，从断点接着算。
早期版本是全部跑完才一次性写文件——跑到一半被中断，已算好的结果全在内存里跟着没了。这个坑已经填了。
代码分段
r
# ① 蛋白→基因坐标(hg38) + tar 文件名；蛋白→其显著疾病（从对应版本 *_significant.csv）
sig <- Sys.glob(sprintf("results/screen_%s/*_significant.csv", rel))
prot_dis <- 跨>=2病的蛋白 及其疾病列表
todo <- setdiff(targets, 已完成的_parts) # ★ 断点续跑
# ② 暴露区域读取：从 tar 解压 cis 染色体文件 → 切 cis 窗口(±500kb)
read_prot_region <- function(prot) {
tar xf 蛋白.tar '内含的 chr<N> 文件' # 只解压需要的那条染色体
d[GENPOS 在 cis 窗口内]
返回 (key=chr:pos, ea=ALLELE1, oa=ALLELE0, beta=BETA, se=SE)
}
# ③ 结局区域读取：zcat + awk 按 chr:pos 窗口提取 FinnGen
read_dis_region <- function(疾病, chr, pos) {
zcat FinnGen | awk '$1==chr && $2 在窗口 {print chr:pos, alt, ref, beta, se}'
}
# ④ 逐蛋白：对齐等位（一致/翻转变号/模糊丢弃）→ 组 betas/ses 矩阵 → hyprcoloc
for (prot in todo) {
合并 [蛋白 + 各显著疾病] 在共同 SNP 上的 beta/se
hr <- hyprcoloc::hyprcoloc(betas矩阵, ses矩阵, trait.names=c(蛋白, 疾病...))
fwrite(hd, "_parts/<蛋白>.csv") # ★ 立即落盘，中断不丢
}
# ⑤ 汇总所有 part（含历史已算的）→ hyprcoloc_results.csv + hyprcoloc_shared.csv输出
results/hyprcoloc_R9/ （R13 同构，在 hyprcoloc_R13/）
├── _parts/<蛋白>.csv 每蛋白一份，断点续跑的依据
├── hyprcoloc_results.csv 全部结果（每行一个共定位簇/迭代）+ 判定列
├── hyprcoloc_shared.csv ★ 只留**真正支持**的簇（PP>0.7 且蛋白在簇内）
└── hyprcoloc_verdict.csv ★ **每蛋白一行的判决**（支持/证据不足/疑LD/无簇）——最省事的入口关键列：traits（同簇的性状，None=没检出共享）、posterior_prob（后验概率）、regional_prob、candidate_snp（候选因果变异）、posterior_explained_by_snp、protein、n_snp、n_traits、release，以及判定列 pp_pass / protein_in_cluster / coloc_support。
判定共定位「支持」需同时满足两条（都在 config.yaml 里配，不硬编码）：
posterior_prob > hyprcoloc_pp_threshold（本课题 0.7）- 该蛋白自己在簇里（
protein_in_cluster）——只有疾病成簇而蛋白不在，说明蛋白信号与疾病信号不是同一个变异，MR 关联多半是 LD 巧合
hyprcoloc()的reg.thresh/align.thresh也从 config 传入（各 0.5）。早期版本没传，用的是包默认值——恰好也是 0.5，结果一致纯属巧合；改 config 会以为生效其实没有。已修。
等位对齐（关键）
暴露和结局的效应等位可能相反。代码按 ea/oa 判断：一致→原样；翻转→ beta 变号；两边对不上（模糊）→丢弃该 SNP。保证多性状 beta 在同一等位方向上。
hyprcoloc 包安装
hyprcoloc 的 C++ 与新版 RcppEigen 不兼容，需先装旧版 RcppEigen 0.3.3.9.4 再 install_github("jrs95/hyprcoloc", build=FALSE)。
相关：阶段六·靶点网络 · 下游分析（概念+结果） · 批量筛查与出图脚本 · 分析方法通俗解释