主题
流水线代码逐文件详解(入门向)
面向基础较弱的读者,从零讲清楚
mr-pipeline里每个代码文件是做什么的、它们怎么串起来。 对应目录:~/mr-pipeline(本地 WSL / 腾讯云 RStudio 同构)。
一、这套代码整体是干嘛的?
一句话:它是一条自动化流水线,用来回答"某个蛋白的水平高低,会不会导致某种病"这类因果问题——这就是孟德尔随机化(Mendelian Randomization, MR)。
你给它三样东西:
- 暴露(exposure)——比如某个蛋白,用它的基因变异 SNP 当"代理"
- 结局(outcome)——比如糖尿病视网膜病变的 GWAS 数据
- 配置(config)——告诉它文件在哪、参数怎么设
它自动跑完一整套统计分析,最后吐出结果表 + 图 + 报告。
一个比喻
把它想成一条工厂流水线:原材料(GWAS 数据)从一头进去,经过十几道工序(清洗、筛选、对齐、计算、检验、出报告),成品(因果结论)从另一头出来。
二、一次运行是怎么走的(10 道工序)
当你敲 Rscript run_all.R,代码依次做这些事:
① 读配置 看你要分析哪些暴露/结局、参数是啥
② 读+标准化数据 把各种格式的原始数据洗成统一格式
③ 选工具变量 从暴露里挑出"合格的" SNP 当工具
④ 谐化对齐 把暴露和结局的等位基因对齐(防方向搞反)
⑤ MR 估计 算"蛋白→疾病"的因果效应(OR / beta)
⑥ 敏感性检验 各种方法互相验证,排除假阳性
⑦ 共定位 确认是"同一个变异"在起作用,不是巧合
⑧ 多重检验校正 跑了几百个,校正假阳性(FDR)
⑨ 证据整合 把 MR 和共定位的结论合起来判断
⑩ 出报告 生成 HTML 报告 + 图表每道工序 = R/ 文件夹里的一个编号文件。编号 01→12 就是工序顺序。
三、逐个文件详解
A. 入口脚本(run_all.R 在根 + analysis/ 编号)
| 文件 | 说明 |
|---|---|
| run_all.R(根,框架层) | 总开关。跑完整流水线:一个暴露 × 一个结局 → 出完整报告。开头把 R/ 里所有工序"装配"进来,再按 ①~⑩ 顺序指挥 |
| analysis/01_screen.R | 批量筛查版(主力)。一次把几十个结局都跑一遍,做"一个蛋白扫一堆病"。读 outcome_manifest.csv 决定跑哪些 |
| legacy/batch_screen.R | 上面那个的旧版(只有 5 个病、参数写死)。留作与原始实现结果对照 |
| legacy/screen_t2d.R | 专门跑 T2D(Mahajan 数据无 rsID,需按染色体位置匹配)的单独脚本 |
| analysis/02_visualize.R | 把结果画成图(森林图、网络图等) |
| validate_real.R | 用真实数据做验证测试 |
下游验证脚本(MR 筛完之后跑,见 下游分析):
| 文件 | 说明 | 用法 |
|---|---|---|
| analysis/03_phewas.R | 把显著工具 SNP 拿去 OpenGWAS 扫全表型,评估多效性 | Rscript analysis/03_phewas.R(不分版本,R9+R13 合并) |
| analysis/04_ldsc.R | 算疾病两两的遗传相关 rg | Rscript analysis/04_ldsc.R R9(或 R13) |
| analysis/05_hyprcoloc.R | 多性状共定位:一个蛋白 + 多个病是不是同一因果变异 | Rscript analysis/05_hyprcoloc.R R9(或 R13) |
| analysis/09_report.py | 把结果汇总成发文章用的多表 Excel(筛查层 + 12 张下游表) | python analysis/09_report.py results |
| analysis/10_eqtl_coloc.R | 血浆 pQTL × 视网膜 eQTL × AMD 三性状共定位 | Rscript analysis/10_eqtl_coloc.R |
| analysis/11_singlecell.py | 候选基因在视网膜各细胞类型的表达 | python3 analysis/11_singlecell.py |
| analysis/12_figure_screen.R | 火山图 + 蛋白×结局效应热图 | Rscript analysis/12_figure_screen.R results/screen_R13 |
| analysis/13_figure_locus.R | 区域共定位多轨图(每个候选一张) | Rscript analysis/13_figure_locus.R R13 |
| analysis/14_figure_downstream.R | PheWAS 负担图 + LDSC 热图/遗传力 | Rscript analysis/14_figure_downstream.R R13 |
| analysis/15_figure_singlecell.py | UMAP + feature 图 + 细胞类型富集检验 | python3 analysis/15_figure_singlecell.py |
各图怎么读、有哪些不能踩的坑:见 图表清单与解读。
A2. 共用层(下划线开头,不单独运行)
这三个文件不是"阶段",是被别的脚本 source / import 的函数库。抽出来的目的很实际:同一段逻辑只写一遍,改的时候不会漏改。
| 文件 | 提供什么 | 谁在用 |
|---|---|---|
| analysis/_viz_common.R | save_both()(一次落 PNG 600dpi + 矢量 PDF)、统一主题 theme_mr()、统一配色 | 12 / 13 / 14 |
| analysis/_region_io.R | 区域数据读取器(从 tar 解 pQTL、从 gz 切 FinnGen 窗口),结果缓存到 results/_region_cache/ | 05 / 13 |
| analysis/_singlecell_io.py | 图谱 backed 载入、候选基因符号匹配(含别名表) | 11 / 15 |
| analysis/_ld_io.R | 调 plink 算 LD r²(1000G EUR 面板)+ LocusZoom 分档配色 | 13 |
为什么 _region_io.R 值得单独存在
原来这两个读取器写在 05_hyprcoloc.R 里,区域数据算完就丢。13 要画区域图得再读一遍——如果复制一份代码,以后改提取逻辑就得改两处。 抽出来之后顺带加了缓存:区域图第一次跑约 20 s/基因,第二次秒出。
一个真实踩过的坑:awk 空字段会让 fread 静默截断
_region_io.R 从 FinnGen 切区域时要带上 rsids 列(第 5 列)供 LD 匹配。但部分行 rsids 是空的,awk 输出那行就只有 5 个字段,而 fread(col.names=6个) 遇到列数不齐会 "Stopped early" 直接截断整个区域——区域数据被砍到几十行,图几乎空白,而且不报错只给 warning。 修法:awk 里对空值补占位符 .,fread 再加 fill=TRUE 兜底,同时把缓存 key 升版(旧缓存里存的是截断数据,必须作废)。 教训:共用层一改,所有依赖它的阶段都要重新验证——这个 bug 如果没发现,重跑 05_hyprcoloc.R 会算出错误的共定位结果。
B. 配置文件(不是代码,是"设置")
| 文件 | 说明 |
|---|---|
| config.yaml | 总设置表:暴露/结局文件路径、显著性阈值、cis 窗口、LD 面板、人群祖先等。换课题只改这个 |
| outcome_manifest.csv | 结局清单:每个结局的文件名、病例数、对照数、匹配方式。analysis/01_screen.R 照它跑 |
| QUICKSTART.md / USAGE.md | 使用说明 |
C. R/ = 引擎(被入口脚本调用,不单独运行)
两个"基础设施"(通用工具,不算工序):
| 文件 | 说明 |
|---|---|
| utils.R | 工具箱:日志、存/读文件、断点续跑等小工具,所有模块共用 |
| 00_setup.R | 开机准备:自动装/加载所需 R 包(TwoSampleMR、coloc 等)、设随机种子 |
12 道工序(按顺序):
| 文件 | 工序 | 说明 |
|---|---|---|
| 01_config.R | 读配置 | 读 config.yaml 并检查有没有填错(build 一致吗?祖先匹配吗?),错了立即报错;还负责算 cis 窗口 |
| 02_standardize.R | 洗数据 | 把五花八门的原始格式洗成统一表:-log10P 还原成 P、OR 转成 beta、去坏行/去重 |
| 04_iv_select.R | 选工具 | 从暴露挑合格 SNP:够显著吗?在 cis 窗口内吗?去连锁冗余(clump)?算 F 值剔除弱工具 |
| 05_harmonise.R | 对齐 | 把暴露和结局的等位基因对齐(A/T、正负链),防效应方向搞反;Steiger 检验方向 |
| 06_mr_core.R | 算因果 | 核心计算:按工具数量自适应——1 个用 Wald ratio、多个用 IVW / Egger 等 |
| 07_sensitivity.R | 稳健性 | 用多种方法互相验证:Cochran's Q、Egger 截距、MR-PRESSO、留一法(LOO),排除多效性/离群点假象 |
| 08_coloc_abf.R | 共定位 | 确认暴露和结局是同一个因果变异在起作用(PP.H4>0.8),而非两个不同变异碰巧挨着 |
| 09_hyprcoloc.R | 多性状共定位 | 升级版:一个蛋白 + 多个疾病,在同一位点一起做共定位 |
| 10_multiple_testing.R | 校正 | 几百个组合做 FDR 校正;把 MR 与共定位证据整合成一个判断 |
| 11_report.R | 出报告 | 生成出版级 HTML 报告(散点/森林/漏斗图,遵循 STROBE-MR) |
| 12_bidirectional.R | 反向验证 | 反过来问"疾病会不会影响蛋白"(结局→暴露),确认因果方向没反 |
为什么没有 03 号?
早期编号遗留——谐化那步后来合并进了 05,跳号属正常。
D. R/adapters/ = 读取适配层(专门读不同格式的数据)
| 文件 | 说明 |
|---|---|
| read_dispatch.R | 调度员:看数据是什么格式,分派给对应读取器 |
| read_tsv.R | 读普通表格式 GWAS(FinnGen 这种) |
| read_instrument_table.R | 读工具表(UKB-PPP 那种一行一个 pQTL) |
| read_vcf.R | 读 VCF 格式(GWAS-VCF) |
为什么单独抽一层:不同数据库格式不一样,把"读取"独立出来,主流水线就不用管格式差异——换数据源只要加个读取器,引擎不动。
四、它们怎么串起来的
你运行: run_all.R (或 analysis/01_screen.R)
│ 开头 source() 装配下面所有引擎
▼
┌──────── R/ 引擎 ────────┐
utils + 00_setup (工具箱 / 开机)
01→02→04→05→06→07→08→09→10→11→12 (12 道工序)
adapters/ (读数据时调用)
│
▼
results/ (结果表 + 图 + HTML 报告)核心思想:
run_all.R/analysis/01_screen.R是指挥R/是干活的工序adapters/是搬运原料的config.yaml/outcome_manifest.csv是任务单
改任务只动任务单,工序引擎不用碰——这就是为什么它能"换课题只改配置"。
五、三层分工(常见疑问)
三处都是代码,区别不是语言,而是角色:
R/= 引擎库(零件):只被source()加载的函数模块,自己不单独运行。run_all.R= 通用引擎入口(框架层):留在根目录,跑"单暴露深挖"。analysis/= 本课题分析入口(编号 01→16):你实际Rscript运行的驱动脚本,做批量筛查 + 下游。
💡 入口脚本里写的是
source("R/…")、fread("data/…")这类相对路径,默认"站在项目根目录运行"。所以哪怕脚本挪进analysis/子目录,只要从仓库根运行(cd ~/mr-pipeline && Rscript analysis/01_screen.R),相对路径照常解析——Rscript 的工作目录取自调用位置,不取自脚本所在目录。另有
legacy/(被取代脚本,留档对照)。跑全链用根目录的run_pipeline.sh。
六、每个 R 文件在哪份文档里讲(完整对照表)
你会发现 R/ 里有十几个文件,而「方法引擎」只有 8 个页面——这不是漏了,而是方法引擎按"工序"分页,有的一页对应多个文件;基础设施类文件则在本页讲。下表让每个文件都有据可查:
| R 文件 | 角色 | 详细文档 |
|---|---|---|
R/utils.R | 工具箱(日志/IO/续跑) | 本页 §三-C |
R/00_setup.R | 依赖加载/开机 | 本页 §三-C |
R/01_config.R | 读配置+校验+cis窗口+祖先检查 | config 配置详解 |
R/02_standardize.R | 数据标准化 | 方法引擎 01-数据读取与标准化 |
R/adapters/*.R | 读取适配(tsv/工具表/vcf/调度) | 方法引擎 01-数据读取与标准化 |
R/04_iv_select.R | 工具变量选择 | 方法引擎 02-工具变量选择 |
R/05_harmonise.R | 谐化与定向过滤 | 方法引擎 03-谐化与定向过滤 |
R/06_mr_core.R | MR 估计 | 方法引擎 04-MR估计方法 |
R/07_sensitivity.R | 敏感性分析 | 方法引擎 05-敏感性分析 |
R/08_coloc_abf.R + R/09_hyprcoloc.R | 共定位 | 方法引擎 06-共定位 |
R/10_multiple_testing.R | 多重检验+证据整合 | 方法引擎 07-多重检验与证据整合 |
R/11_report.R + R/12_bidirectional.R | 报告 + 反向MR | 方法引擎 08-报告与双向MR |
⇒
R/里每个文件都有文档覆盖:方法工序在「方法引擎」8 页,基础设施在本页,配置读取在「config 配置详解」。
七、我到底要改哪个文件?(一句话)
- 跑单条分析 → 改
config.yaml(Rscript run_all.R) - 批量扫多个病 → 改
outcome_manifest.csv(Rscript analysis/01_screen.R R9) - 改方法/参数 → 也在
config.yaml里(iv/ld/coloc段),R/引擎不用动
延伸阅读:结果怎么看 · config 配置详解 · 流水线架构 · 使用教程 · 方法引擎详解