通用两样本 MR 流水线 · 使用教程

一句话：改一个 config.yaml，就能对任意暴露、任意结局跑标准化的两样本孟德尔随机化（MR）+ 共定位，严格遵循 STROBE-MR。方法引擎不用动，换课题只改配置。

本页讲的是 run_all.R 单条精细分析

换数据要跑的若是 analysis/ 批量筛查主线（01→16），那条线不吃 config.yaml 的 exposures/outcomes，要改的是 outcome_manifest.csv 和几处写死的清单——见 换数据重新开始。

📖 本页讲怎么用。想了解引擎内部的编排见流水线架构，想了解每一步的统计方法与公式见方法引擎详解，真实落地范例见案例研究。

当前进度（2026-07-11）：M1–M4 全部完成并验证通过，整套流水线可一键从原始数据跑到出版级 HTML 报告。

M1 读取/标准化/校验 · M2 IV选择/谐化/自适应MR/敏感性 · M3 共定位(coloc.abf + hyprcoloc)+证据整合 · M4 多重检验(FDR)独立模块 + STROBE-MR HTML 报告(散点/森林/漏斗/LOO 图 600dpi 内嵌 + 另存矢量 PDF)

1. 项目在哪、怎么进

部署在腾讯云 RStudio（就是 MR 分析主力机）：

网页版 RStudio：https://rstudio.stlog.cn　账号 research / 密码 〈密码·见Vaultwarden〉
命令行：本机 Git Bash 里 ssh tencent
项目路径：~/mr-pipeline（即 /home/research/mr-pipeline，登录 RStudio 后即在此账号下）

LD 参考面板（1000G EUR，clumping 用）在 /mnt/d/mrdata/ldpanel/EUR，plink 用系统 PATH 上的 plink1.9，都不用你管。默认用欧洲人面板（ld$ancestry: EUR）；换非欧暴露要同步换面板，见 §4.4。

2. 30 秒跑通（验证环境）

RStudio 的 Terminal 或 SSH 里：

bash

cd ~/mr-pipeline
Rscript run_all.R

会读 config.yaml 里配的暴露/结局，跑完在 results/ 下出结果。第一次用自带的示例配置即可跑通。

3. 目录结构（心里有个图）

mr-pipeline/
├── config.yaml            ← 唯一需要你改的文件
├── run_all.R              ← 一键入口
├── R/                     ← 方法引擎（一般不用动）
│   ├── 01_config.R  02_standardize.R  04_iv_select.R
│   ├── 05_harmonise.R  06_mr_core.R  07_sensitivity.R
│   └── adapters/          ← 读取适配层（tsv / instrument_table / vcf）
├── reference/
│   ├── ld/EUR.*           ← 1000G 欧洲人群 LD 面板（clumping 用）
│   └── bin/plink          ← plink 二进制
├── data/
│   ├── exposure/          ← 你的暴露文件放这
│   └── outcome/           ← 你的结局文件放这
└── results/
    ├── standardized/      ← 每个数据集标准化后的缓存
    ├── iv/                ← 每个暴露选出的工具变量
    ├── pairs/<暴露>__<结局>/   ← 每个暴露-结局对的独立结果 ★
    │   ├── harmonised.rds/.csv     谐化后的 IV（含 F/R²）
    │   ├── mr.rds/.csv             全部方法的估计
    │   ├── sensitivity.rds         Q/Egger/PRESSO/LOO/方向性
    │   └── pair_summary.rds/.csv   ← 标准化汇总（下游分析读这个）
    └── mr_all_pairs_primary.csv    ← 所有对的主结果 + FDR 总表

4. 核心：怎么改 `config.yaml`

换一个新课题，你只做三件事：① 把文件放进 data/ → ② 在 config 里加一条暴露、一条结局 → ③ Rscript run_all.R。

4.1 加一个「结局」（标准 GWAS summary statistics，`tsv` 格式）

以 FinnGen 为例，在 config.yaml 的 outcomes: 下加一块：

yaml

outcomes:
  - id: "DR_finngenR13"                 # 自己起的短名，会用于结果目录名
    format: "tsv"                       # 扁平表（tsv/csv/.gz 都行）
    file: "data/outcome/finngen_R13_DR.gz"
    build: "GRCh38"                     # 基因组版本（必填，用于一致性校验）
    ancestry: "EUR"                     # 样本祖先（STROBE-MR 需报告；与暴露不同=跨人群MR，报告会警示）
    trait: "Diabetic retinopathy (FinnGen R13)"
    type: "binary"                      # binary→出 OR；continuous→出 beta
    sample_size: 312000
    ncase: 12000                        # binary 才需要
    ncontrol: 300000
    prevalence: 0.038                   # binary 共定位用
    pval_encoding: "raw"                # p 值原样；若列是 -log10 则填 neglog10
    effect_type: "beta"                 # beta | OR | HR（OR/HR 会自动取 log）
    column_map:                         # 列名映射（留空则自动探测）
      SNP: rsids
      chr: "#chrom"
      pos: pos
      effect_allele: alt               # FinnGen: alt 是效应等位
      other_allele: ref                # ref 是其他等位
      eaf: af_alt
      beta: beta
      se: sebeta
      pval: pval

列名不确定？ 先不写 column_map，让它自动探测（内置了 FinnGen/UKB/GWAS-SSF 常见别名）。探测不到会明确报错告诉你缺哪列，再补映射即可。

4.2 加一个「暴露」

情况 A：合并工具表（如 UKB-PPP 的 protein_info.csv） —— 一张表里多个蛋白、一行一个已选 pQTL：

yaml

exposures:
  - id: "IL6R"
    format: "instrument_table"          # ← 工具表模式
    file: "data/exposure/ukbppp_instruments.csv"
    build: "GRCh38"
    ancestry: "EUR"                     # 样本祖先（必填）；须与 ld$ancestry 一致，否则 clumping 报错拦截
    trait: "IL6R protein"
    type: "continuous"
    gene: "IL6R"                        # cis 模式用它定位窗口
    chr: 1                              # 或直接给坐标
    tss: 154377669
    tes: 154441926
    sample_size: 34557               # UKB-PPP discovery(欧裔)；须与所取的 discovery 效应量列一致
    pval_encoding: "neglog10"           # UKB-PPP 的 log10p / pvallog 是 -log10
    effect_type: "beta"
    variant_id_split:                   # 复合变异 ID 拆列（没有就删掉这块）
      column: "Variant.ID"
      pattern: "chr:pos:OA:EA"          # 声明冒号各段含义：A0=OA, A1=EA
    filter: { "cis.trans": "cis" }      # 只取 cis 行

情况 B：全量 GWAS sumstats（一个风险因子一个文件，如 LDL） —— 用 format: "tsv"，和结局写法一样（type: continuous），流水线会按阈值自己选显著 SNP。

4.3 关键字段速查

字段	作用	常见取值
`format`	输入类型	`tsv`（扁平表）/ `instrument_table`（工具表）/ `vcf`（M3 支持）
`build`	基因组版本	`GRCh38` / `GRCh37`；暴露与结局必须一致（否则报错）
`ancestry`	样本祖先	`EUR`/`EAS`/`AFR`/`AMR`/`SAS`；暴露须与 LD 面板一致（否则 clumping 报错），供 STROBE-MR 报告
`pval_encoding`	p 值编码	`raw`（原始 p）/ `neglog10`（-log10p，如 FinnGen mlogp、UKB log10p）
`effect_type`	效应量类型	`beta` / `OR` / `HR`（后两者自动取 log 转成 beta）
`type`	结局性质	`binary`（出 OR）/ `continuous`（出 beta）
`column_map`	列名映射	留空=自动探测；探测不到再显式写
`filter`	行过滤（工具表）	如 `{ "cis.trans": "cis" }`
`variant_id_split`	复合 ID 拆列	`pattern` 用 `chr:pos:OA:EA` 声明各段

4.4 换非欧人群暴露（必读）

用非欧人群 pQTL/GWAS（如东亚 pQTL 图谱、非裔队列）当暴露时，LD clumping 的参考面板必须换成对应人群，否则连锁结构错配会让选出的工具变量出现偏差 —— 这是 STROBE-MR 关注的人群分层问题。

流水线已内置防呆：暴露 ancestry 与 LD 面板祖先不一致时直接报错拦截（加载 config 时就 fail-fast，不会跑出错误结果）。换人群要同步改三处：

yaml

iv:
  clump_pop: "EAS"              # api 模式：改成对应人群
ld:
  bfile: "/mnt/d/mrdata/ldpanel/EAS"   # local 模式：换成对应人群面板
  ancestry: "EAS"              # 面板祖先，须与下面暴露一致
exposures:
  - id: "IL6R_EAS"
    ancestry: "EAS"            # 暴露样本祖先
    # ...其余字段同前

单一 config 只能服务同一祖先

一份 config 只有一套全局 LD 配置，不能同时跑不同祖先的暴露。要跑多祖先（如欧洲 + 东亚交叉验证），请每个祖先各写一份 config 分别 Rscript run_all.R。这也符合分祖先分析、避免跨人群混杂的最佳实践。

暴露与结局祖先不一致

若暴露与结局祖先不同（跨人群两样本 MR），流水线不拦截（LD 面板只跟暴露走），但报告的「人群祖先与分层」段会自动打出 ⚠️ 警示，提示等位频率/LD 差异可能引入偏差，结论需谨慎。

5. 两种分析模式

在 config 顶部 analysis_mode: 切换：

cis（分子暴露：cis-pQTL / cis-eQTL）：按基因 TSS/TES ±窗口取 SNP，阈值 5e-6，强制做共定位（M3）。单 IV 时用 Wald ratio。
genome_wide（普通风险因子）：全基因组显著 5e-8，跑全套敏感性。

阈值、窗口都在 iv: 段可调：

yaml

iv:
  pval_cis: 5.0e-6
  pval_genome_wide: 5.0e-8
  cis_window_kb: 500        # ±500kb，可改 1000 = ±1Mb
  clump_r2: 0.001
  clump_kb: 10000
  f_stat_min: 10            # 弱工具阈值
  steiger_filtering: true

6. 运行与断点续跑

bash

Rscript run_all.R                 # 用默认 config.yaml
Rscript run_all.R my_config.yaml  # 用别的配置

断点续跑：config.yaml 里 project.resume: true 时，已算好的中间结果（.rds）直接复用，改一处重跑不会全部重来。想全新跑就把某个中间目录删掉，或设 resume: false。当前仓库即为 resume: true；2026-07-22 分家那次曾临时设 false 以做「从零验证有效性」，投稿定稿前的可复现性证明也应设 false 跑一次。
阶段裁剪：run.stages 控制跑哪些步骤，如只想重出 MR 不重读数据。

7. 结果怎么读

📖 逐列/逐图的完整讲解见结果怎么看（含真实数字实例、判读检查清单）。本节只列产物清单。

最常看的三个产物：

results/report/MR_report.html —— 出版级 STROBE-MR 报告：方法学（含全部文献）、主结果表、逐对散点/森林/漏斗/leave-one-out 图（600dpi 内嵌，另存矢量 PDF 于 results/figures/）、共定位结论。用浏览器打开这个就够看全貌。
results/mr_all_pairs_primary.csv —— 所有暴露-结局对的主结果一张表，含 OR/beta、95%CI、P、FDR、F 统计量、evidence。
results/pairs/<暴露>__<结局>/ —— 单对的全部细节：
- mr.csv：IVW/Egger/加权中位数/众数所有方法
- sensitivity.rds：Cochran's Q、Egger 截距、MR-PRESSO、LOO、方向性
- pair_summary.csv：标准化汇总，下游 LDSC/PheWAS/网络分析直接读这个（字段稳定，coloc 列由 M3 填）

pair_summary.csv 字段：exposure_id, outcome_id, method, n_iv, estimate_type, point, ci_low, ci_high, pval, overall_f, isq_gx, evalue, evalue_ci, egger_intercept_p, presso_global_p, coloc_status, coloc_pp_h4, coloc_shared, hyprcoloc_status, hyprcoloc_coloc, hyprcoloc_pp, rev_status, rev_pval, rev_significant。主表 mr_all_pairs_primary.csv 另有 p_fdr / significant_fdr / evidence（证据整合：MR 显著且 coloc 支持=一致；MR 显著但 coloc 不支持=证据不一致）。

STROBE-MR 合规（M5）

报告已按 STROBE-MR 20 条清单逐条覆盖分析类条目：

item 4/10b 数据源表（含各数据集 ancestry 样本祖先列）、item 5/7 三大 IV 假设显式声明+评估映射、item 6/10d 人群祖先与分层（LD 面板祖先 + 暴露/结局祖先，跨人群自动警示）、item 9a 软件与包版本、item 10a IV 筛选流程表（原始→显著→cis窗→clump→F 剩余数）、item 10d-ii 样本重叠说明、item 12b I²_GX + E-value、item 13c 双向 MR、item 19 数据/代码可及性。
人群/LD 面板匹配：每个暴露声明 ancestry，clumping 前校验其与 LD 面板（ld$ancestry/clump_pop）一致，错配 fail-fast 拦截，避免连锁结构错配的工具变量偏差。
新增方法：MR-RAPS（弱工具/多效性/样本重叠稳健）已并入 ≥3 IV 的方法集。
双向 MR：config 里 bidirectional.enabled: true 开启（需暴露具备全量/区域 sumstats 以匹配反向工具，否则守卫跳过）。
样本重叠 / 数据可及性 / 代码可及性：在 config 的 sample_overlap.note 和 report.data_availability / code_availability 里填，报告自动引用。

共定位怎么用（M3）

config 的 coloc.engines: ["coloc_abf","hyprcoloc"] 控制跑哪些引擎；min_region_snps（默认 50）是守卫，区域重叠 SNP 不足就跳过并标注"数据不足"，不会硬跑出误导结果。
coloc 需要"区域级"数据（该基因位点附近成百上千个 SNP），而不是只有几个入选 IV：
- 结局（tsv 全量 sumstats）：自动按 cis 窗口切出区域，无需额外操作。
- 暴露若是 instrument_table（只有工具 SNP）：请在该暴露下填 region_file: "data/exposure/xxx_region.tsv" 指向该基因位点的全量 summary stats；不填则 coloc 因 SNP 不足而跳过。
按 rsID 对齐；hyprcoloc 会自动把各性状 β 翻转到统一等位方向。

8. 换疾病 / 换 FinnGen 版本 = 改一条 config

场景：把结局从 FinnGen R13 的 DR 换成 R13 的 AMD

把 finngen_R13_AMD.gz 放到 data/outcome/
config 的 outcomes: 复制一份、改 id/file/trait/ncase/ncontrol
Rscript run_all.R

换 FinnGen 版本（R9 → R13）：只改 file 路径和 sample_size/ncase/ncontrol。列名一致就不用动 column_map。

加更多蛋白：instrument_table 一张表里已有多个蛋白，config 里给每个想分析的蛋白各写一条 exposure（filter 会自动只取该蛋白的 cis 行）；或后续做成批量循环。

8.5 批量筛查与下游分析的命令速查

本页到这里讲的都是 run_all.R（单条精细分析）。本课题实际主力是批量筛查 + 下游三件套，命令一并列在这：

bash

# 批量筛查（结局清单在 outcome_manifest.csv）
Rscript analysis/01_screen.R R9            # 或 R13；加 force 忽略缓存全部重算
Rscript analysis/02_visualize.R results/screen_R9
python analysis/09_report.py results

# 下游三件套
Rscript analysis/03_phewas.R                # 多效性扫描
Rscript analysis/04_ldsc.R R9               # 遗传相关（或 R13）
Rscript analysis/05_hyprcoloc.R R9       # 多性状共定位（或 R13）
Rscript analysis/06_network.R R9          # 阶段六：整合出靶点网络（或 R13）

哪些分 R9/R13，哪些不分

脚本	分版本吗	输出
`analysis/01_screen.R`	✅ 分	`results/screen_R9/`、`results/screen_R13/`
`analysis/03_phewas.R`	❌ 不分——多效性是 SNP 自身属性，跟结局用哪版无关，合并查一次即可	`results/phewas/`
`analysis/04_ldsc.R`	✅ 分	R9→`results/ldsc/`；R13→`results/ldsc_R13/`
`analysis/05_hyprcoloc.R`	✅ 分	`results/hyprcoloc_R9/`、`results/hyprcoloc_R13/`

详见批量筛查与下游分析。

9. 常见报错与排查

报错/现象	原因	解决
`require_rsid=TRUE 但仅 X% 是 rsID`	SNP 列不是 rsID（是 chr:pos）	换带 rsID 的文件，或 `column_map` 指到 rsID 列
`build 不一致`	暴露/结局基因组版本不同	统一到同一 build（都 GRCh38 最省事），或开 liftover（可选模块）
`暴露祖先(X) 与 LD 面板(Y) 不一致`	暴露 `ancestry` 与 `ld$ancestry`/`clump_pop` 不匹配	换匹配人群的 LD 面板（见 §4.4），或纠正暴露 `ancestry`；多祖先分 config 跑
`暴露未声明 ancestry` 告警	暴露条目缺 `ancestry`	补上 `ancestry`（供 STROBE-MR 报告 + 面板匹配校验）；不补则跳过校验、祖先匹配性不保证
`列自动探测失败，缺 X`	列名太非标准	在该数据集的 `column_map` 里显式写 `X: 真实列名`
clumping 很慢	每次加载 850万 SNP 面板（~8s）	正常；多暴露时是一次性开销
内存吃紧	读全基因组文件	reader 已只读必要列；真不够就分染色体或升配
p 值全变 1 或极端	`pval_encoding` 填错	FinnGen mlogp / UKB log10p 要填 `neglog10`
OR 当成 beta	`effect_type` 没填 OR	效应量是 OR/HR 时务必填对，会自动取 log

10. 当前进度与后续里程碑

✅ M1：读取（tsv + instrument_table）/ 标准化 / rsID·build 校验
✅ M2：IV 选择（cis窗口/clump/F&R²/弱工具）/ 谐化+回文+Steiger / 自适应 MR（Wald·IVW·Egger·加权中位数/众数）/ 敏感性（Q·Egger截距·MR-PRESSO·LOO·方向性）
✅ M3：共定位 —— coloc.abf（两两，PP.H4>0.8 判共享）+ hyprcoloc（多性状），按 rsID 对齐 + 等位定向对齐，区域 SNP 守卫，MR×coloc 证据整合
✅ M4：多重检验独立模块（by_outcome/global FDR）+ 出版级 STROBE-MR HTML 报告（散点/森林/漏斗/LOO 图，600dpi 内嵌 + 矢量 PDF）
⏸ 可选下游（不在本轮，但 pair_summary 已为其留好接口）：LDSC 遗传相关、PheWAS、网络图
✅ M5：STROBE-MR 20 条合规增强 —— MR-RAPS、I²_GX、E-value、双向 MR、IV 流程表、数据源表、软件版本、数据/代码可及性
⏸ 可选未来：GWAS-VCF 适配（read_vcf stub）、cisMR-cML / cis-MRBEE（2024 稳健 cis 法，需额外包）、多变量 MVMR、下游 LDSC/PheWAS/网络

整套流水线到此完整。 换课题只改 config.yaml，Rscript run_all.R 一键出符合 STROBE-MR 的报告。

附录：常见数据源列名映射速查

源	SNP	效应等位	其他等位	eaf	beta	se	pval	备注
FinnGen	`rsids`	`alt`	`ref`	`af_alt`	`beta`	`sebeta`	`pval`(或 `mlogp`=neglog10)	GRCh38
UKB-PPP	`rsid`	`ALLELE1`/A1	`ALLELE0`/A0	`A1FREQ`	`BETA`	`SE`	`LOG10P`=neglog10	GRCh38，工具表
GWAS-SSF	`rsid`	`effect_allele`	`other_allele`	`effect_allele_frequency`	`beta`	`standard_error`	`p_value`	GWAS Catalog 标准
deCODE	`rsids`	`effectAllele`	`otherAllele`	`effectAlleleFreq`	`Beta`	`SE`	`Pval`	常为 GRCh38

这些别名大多已内置自动探测，通常不写 column_map 也能认出来；认不出再按上表补映射。

通用两样本 MR 流水线 · 使用教程 ​

1. 项目在哪、怎么进 ​

2. 30 秒跑通（验证环境） ​

3. 目录结构（心里有个图） ​

4. 核心：怎么改 config.yaml ​

4.1 加一个「结局」（标准 GWAS summary statistics，tsv 格式） ​

4.2 加一个「暴露」 ​

4.3 关键字段速查 ​

4.4 换非欧人群暴露（必读） ​

5. 两种分析模式 ​

6. 运行与断点续跑 ​

7. 结果怎么读 ​

STROBE-MR 合规（M5） ​

共定位怎么用（M3） ​

8. 换疾病 / 换 FinnGen 版本 = 改一条 config ​

8.5 批量筛查与下游分析的命令速查 ​

9. 常见报错与排查 ​

10. 当前进度与后续里程碑 ​

附录：常见数据源列名映射速查 ​

通用两样本 MR 流水线 · 使用教程

1. 项目在哪、怎么进

2. 30 秒跑通（验证环境）

3. 目录结构（心里有个图）

4. 核心：怎么改 `config.yaml`

4.1 加一个「结局」（标准 GWAS summary statistics，`tsv` 格式）

4.2 加一个「暴露」

4.3 关键字段速查

4.4 换非欧人群暴露（必读）

5. 两种分析模式

6. 运行与断点续跑

7. 结果怎么读

STROBE-MR 合规（M5）

共定位怎么用（M3）

8. 换疾病 / 换 FinnGen 版本 = 改一条 config

8.5 批量筛查与下游分析的命令速查

9. 常见报错与排查

10. 当前进度与后续里程碑

附录：常见数据源列名映射速查