主题
换数据重新开始一个课题
已有的 config 配置详解 §E 讲的是
run_all.R单条精细分析怎么换课题。 但本课题的主力是analysis/批量筛查主线(01→16),那条线不吃config.yaml的 exposures/outcomes, 换数据要动的东西完全不同。这一页专门讲主线。
先分清你要换的是哪种
| 场景 | 看哪一页 |
|---|---|
| 同一个课题,换 FinnGen 新版本(R13→R14) | 本地环境 §四 — 只加 manifest 一行 |
| 同一批暴露,换一批疾病 | 本页 §二(改 manifest 就够) |
| 换暴露(不再是 UKB-PPP 蛋白) | 本页 §三(要动脚本) |
| 单条暴露×单条结局的精细分析 | config 配置详解 |
一、先搞清楚:哪些是"课题无关"的,哪些写死了本课题
这套代码分三层,换数据时受影响的程度完全不同:
| 层 | 内容 | 换课题要动吗 |
|---|---|---|
R/ 引擎 + run_all.R | 统计方法本体 | 永不动 |
analysis/_*.R _*.py 共用层 | 出图/区域读取/LD/单细胞载入 | 只改路径常量 |
analysis/01–16 | 本课题分析驱动 | 看情况,见下 |
按"课题绑定程度"给 16 个阶段分个类:
| 绑定程度 | 阶段 | 说明 |
|---|---|---|
| 🟢 通用 | 01 筛查、02/12 出图、09/16 汇总 | 只吃 manifest 和结果目录,换什么数据都能跑 |
| 🟡 半通用 | 03 PheWAS、04 LDSC、05 共定位、06 网络、13/14 出图 | 逻辑通用,但有写死的结局清单或路径,要改几行 |
| 🔴 本课题专用 | 07 可成药性、08 IAMDGC 复制、10 视网膜 eQTL、11/15 单细胞 | 绑死了 AMD 课题的特定外部数据,换课题基本要重写或直接不跑 |
二、场景 A:同一批暴露,换一批疾病(最常见)
这是最轻的,不用改任何 .R 文件。
1. 数据放进 D 盘
D:\mrdata\outcome\你的新GWAS.gz不用改路径——data/outcome 已软链到 /mnt/d/mrdata/outcome。
在 Windows 侧下载,别在 WSL 里下
WSL 是 NAT 网络,够不到 Windows 的 localhost 代理。Windows 下载 → WSL 经 /mnt/d 读,最省事。
2. 在 outcome_manifest.csv 加行
表头 12 列,逐列含义:
| 列 | 填什么 |
|---|---|
release | 版本标签,自己取(如 R13、MYSTUDY)。这个值会变成 results/screen_<release>/ 目录名 |
id | 结局短名,出现在所有结果文件名里 |
file | D:\mrdata\outcome\ 下的文件名 |
phenocode | 原始 endpoint 代码,仅备查 |
ncase / ncontrol | 病例/对照数,务必填官方 manifest 的真值(会进 F 统计量和发表表格) |
match_mode | rsid 或 position。没有 rsID 列就填 position,且必须与暴露同 build |
effect_type | beta / OR / HR |
pval_encoding | raw / neglog10 |
se_from_ci | 没有 SE 列、只有置信区间时填 TRUE |
group | 分组标签,出图分面用 |
notes | 中文备注 |
3. 跑
bash
cd ~/mr-pipeline
Rscript analysis/01_screen.R MYSTUDY # 断点续跑:已算过的结局自动跳过
Rscript analysis/12_figure_screen.R results/screen_MYSTUDY
python3 analysis/09_report.py results下游阶段有写死的结局清单
04_ldsc.R 和 05_hyprcoloc.R 里各有一份 COMPL / 清单常量,写死了本课题纳入的 endpoint id:
r
COMPL <- list(
R9 = c("Retinopathy","Maculopathy", ...),
R13 = c("Retinopathy","RetinopathyStrict", ...))换了 release 标签或结局集合,必须在这两个文件里加一条对应的清单,否则 stopifnot(rel %in% names(COMPL)) 会直接报错。
三、场景 B:换暴露(不再是 UKB-PPP 蛋白)
这一步重得多,因为多个阶段绑死了 UKB-PPP 的目录结构。
必须准备的两样东西
| 文件 | 现在是什么 | 换暴露后要提供什么 |
|---|---|---|
data/exposure/protein_info.csv | UKB-PPP ST9(每暴露一个 cis 哨兵,含 CHROM / GENPOS(hg38) / 复合 ID / beta / se / log10p / cis-trans 标记) | 同样的列结构,或改 01_screen.R 的读取段 |
D:\mrdata\ukbppp_regional\*.tar + ukbppp_european_index.csv | 每蛋白一个 tar,内含各染色体全区域统计量 | 共定位与 locus 图必需。没有区域级数据 → 05/13 跑不了 |
只有哨兵、没有区域数据会怎样
01 筛查能跑(MR 只要工具 SNP),但 05 共定位 和 13 locus 图 会全部跳过—— 共定位需要 cis 窗口内几千个 SNP,不是一个 top hit。 本课题的验证支柱就是共定位,没有区域数据等于失去主要验证手段,方案要重新想。
要改的硬编码路径
这些常量散在脚本顶部,换数据源时逐个改:
| 文件 | 常量 | 现值 |
|---|---|---|
05_hyprcoloc.R、08、10、13_figure_locus.R | TARDIR | /mnt/d/mrdata/ukbppp_regional |
| 同上 | 索引 csv | /mnt/d/mrdata/ukbppp_european_index.csv |
04、05、10、13 | FGDIR | /mnt/d/mrdata/outcome |
04_ldsc.R | LD scores | /mnt/d/mrdata/ldsc/eur_w_ld_chr |
_ld_io.R | LD_BFILE | /mnt/d/mrdata/ldpanel/EUR |
_singlecell_io.py | H5AD | /mnt/d/mrdata/singlecell/hrca_allcells.h5ad |
10_eqtl_coloc.R | eQTL 目录 | /mnt/d/mrdata/eqtl_retina_strunz |
08_replicate_iamdgc.R | 复制队列 | /mnt/d/mrdata/outcome/IAMDGC_...txt |
人群要对得上
_ld_io.R 用的是 1000G EUR 面板,04_ldsc.R 用的是 欧洲 LD scores。 换成非欧洲人群的 GWAS,这两处必须一起换(EAS 等人群面板在 Server1 /root/ldpanel/1kg.v3.tgz 里), 否则 r² 和遗传相关都是错的——而且不会报错。
四、场景 C:整个课题换掉(比如不做眼病了)
建议做法:复制整个仓库,删掉专用阶段,而不是在原地改。
bash
cp -r ~/mr-pipeline ~/newproject
cd ~/newproject
rm -rf results .git # 结果和历史都不要带过去然后处理这几类文件:
| 处理 | 文件 |
|---|---|
| 原样保留 | R/、run_all.R、analysis/_viz_common.R、analysis/_region_io.R、analysis/_ld_io.R、tests/ |
| 改路径常量 | 上表那些 |
| 重写清单 | outcome_manifest.csv(清空重填)、04/05 里的 COMPL |
| 删掉或搁置 | 07_druggability.py(AMD 靶点专用)、08_replicate_iamdgc.R(AMD 复制队列)、10_eqtl_coloc.R(视网膜 eQTL)、11/15(视网膜单细胞图谱) |
| 同步删阶段 | run_pipeline.sh 里对应的 run_stage 行 |
run_pipeline.sh 的断点标记在 results/pipeline_logs/*.done,删掉 results/ 就等于全部重来。
五、最小可用清单(只想跑 MR + 出图)
如果新课题只要"筛查 + 图 + 表",不做共定位那套,只需要 4 个阶段:
bash
Rscript analysis/01_screen.R MYSTUDY # MR 批量筛查
Rscript analysis/02_visualize.R results/screen_MYSTUDY # UpSet/森林/网络
Rscript analysis/12_figure_screen.R results/screen_MYSTUDY # 火山图/热图
python3 analysis/09_report.py results # Excel需要的输入只有两样:data/exposure/ 的工具表 + outcome_manifest.csv 里列出的结局文件。 这条最小路径不需要区域数据、不需要 LD 面板、不需要任何外部 API。
六、换完之后必须自查的四件事
ncase/ncontrol填的是真值吗——错了 F 统计量和发表表格全错,而且不报错- build 一致吗——
match_mode: position时暴露与结局必须同 build(本课题都是 GRCh38) - 人群一致吗——GWAS 人群、LD 面板、LDSC 的 LD scores 三者要对上
- 跑一遍 smoke test:
Rscript tests/smoke_test_m1.R,先确认引擎没被改坏
改了共用层,所有依赖阶段都要重验
_region_io.R / _viz_common.R / _ld_io.R 是被多个阶段共用的。改动后不要只看新阶段跑通了—— 本项目真实踩过:给区域读取器加一列,导致 fread 静默截断, 若不是顺手比对了共定位结果,05 会悄悄给出错误结论(详见 代码结构详解)。 验证办法:拿改动前的结果文件做 diff,逐字节一致才算通过。