简介pQTLtools 是面向蛋白质定量性状基因座pQTL研究的 R 工具包适合蛋白质组学、遗传关联和生物信息学分析人员使用主要解决多来源 pQTL 面板数据整理、注释查询与配套函数调用问题也可作为搭建分析环境和查阅参考数据的入口。压缩包共 125 个文件约 5.92MB内含 html 文档与教程页面、png 示意图、rda 数据对象、R 脚本、rd 帮助文件等覆盖安装说明、数据集与功能摘要、网络驱动文档和示例内容目录与文档组织较清晰便于按模块检索。数据层面涉及 Caprion、Olink、SomaLogic 等面板并包含 hg19 注释表与 hg38 链转换文件能帮助读者理解工具结构、复现分析流程并定位常用数据与函数也可辅助完成面板筛选、注释与结果核对。当前已有 759 人学习适合需要开展 pQTL 数据整理、面板资源查阅或 R 包二次开发的研究者参考。1. pQTL 分析在 R 里的接口断层从一份蛋白丰度矩阵说起手上有一份几千样本、上千个蛋白的丰度矩阵基因型躺在另一张 PLINK 表里协变量又在第三个 CSV 里——这是绝大多数 pQTL 项目开工时的真实状态。用 R 写关联检验本身不难几行lm就能出 beta 和 p 值真正耗时间的是把三张表的样本 ID、缺失模式、基因组坐标和批次结构对齐成一套能反复跑、换一批样本还能复用的东西。pQTLtools 这个包切的就是这段接口把蛋白表型侧和基因型侧捏到一起跑顺式与反式关联再产出能直接往下游走的结果表。它适合已经会用 R 读写表格、但对 pQTL 全流程还没形成固定套路的人也适合手上有一套祖传脚本、想找个参照实现来对照检查的老手。下面按装包、对齐、建模、拆结果、排错的顺序拆开讲。2. 装包、读表与三件套对齐把数据喂进 pQTLtools2.1 先核对包内导出函数再决定走哪条路径拿到一个不熟的 R 包我一般不会先去翻教程而是直接把导出对象列出来。原因很实际pQTL 这条线上不同团队的命名习惯差异极大有的把读入写成readPQTL有的叫loadProteomics靠猜函数名比翻文档还慢。先看清单再决定哪些步骤用包里的函数、哪些用自己写的代码兜底能省掉大量试错。# 加载包并列出全部导出对象函数 随包数据集 suppressPackageStartupMessages(library(pQTLtools)) exports - ls(package:pQTLtools) cat(导出对象数量:, length(exports), \n) cat(sort(exports), sep \n)ls(package:...)只给导出对象也就是作者希望外部调用的那部分。想连内部函数一起看换成ls(getNamespace(pQTLtools), all.names TRUE)遇到某个导出函数报错时顺着它调用的内部函数往下读比在网上搜报错信息快得多。# 对可疑函数确认参数签名避免位置参数传错 args(pQTLtools::你的目标函数) # 查看某个函数属于哪种 S3/S4 方法 methods(class 你的结果对象类名)args()的输出里带默认值的参数可以不动没有默认值的必须显式给。pQTL 场景下最容易传错的是「基因组版本」和「cis 窗口大小」这两类参数它们通常有默认值但不一定匹配你的数据靠位置传参很容易踩空一律用参数名 值的形式写。我习惯把导出函数按名字前缀粗分一下组读入类、关联类、cis 判定类、绘图类、工具类。分组之后你会发现包里通常只覆盖了其中两三块剩下的是留给你自己接的。这不是缺陷pQTL 项目的表型来源五花八门包作者不可能把 SomaScan、Olink、质谱三条线的前处理全写死。# 按关键词粗筛函数快速定位功能域 grep(read|load|import, exports, value TRUE, ignore.case TRUE) grep(assoc|lm|glm|test, exports, value TRUE, ignore.case TRUE) grep(cis|trans|window, exports, value TRUE, ignore.case TRUE) grep(plot|manhattan|qq, exports, value TRUE, ignore.case TRUE)这一步的真正价值在于建立预期哪一段流程包能接住、哪一段必须自己写。后面所有代码都按这个边界组织不硬套。2.2 表型、基因型、协变量三张表的 ID 对齐pQTL 项目里九成的「跑出来全是 NA」和「样本量对不上」根子都在 ID 对齐。三张表的 ID 来源天生不一致蛋白丰度表里可能是SampleId、样本编号、或者干脆是行名PLINK 的.fam里有 FID 和 IID 两列到底用哪个取决于建库时怎么写的协变量表经常是 Excel 导出的 CSV第一列没表头被读成V1。我先做一次硬对齐检查把三边的 ID 集合交出来看清楚损失多少再往下走。# 三张表的样本 ID 取交集先看清损失再决定补不补 pheno_ids - rownames(prot_mat) # 蛋白丰度矩阵行名 geno_ids - data.table::fread(geno.fam, header FALSE)[[V2]] # V2 IID cov_ids - rownames(cov_df) # 协变量表行名 common - Reduce(intersect, list(pheno_ids, geno_ids, cov_ids)) cat(表型:, length(pheno_ids), 基因型:, length(geno_ids), 协变量:, length(cov_ids), 交集:, length(common), \n)Reduce(intersect, ...)比intersect(intersect(a,b),c)好读样本表多于三张时直接往列表里加。跑完盯着交集数看如果比最小那张表还少一截说明 ID 格式有系统性差异常见的是大小写、前后空格、或者表型侧带了个前缀。# 找出「差一点就对不上」的 ID去空格、转大写后再交一次 norm - function(x) toupper(trimws(as.character(x))) pheno_n - norm(pheno_ids); geno_n - norm(geno_ids); cov_n - norm(cov_ids) common2 - Reduce(intersect, list(pheno_n, geno_n, cov_n)) cat(规范化后交集:, length(common2), \n) # 列出规范化才匹配上的那些人工扫一眼确认没有误伤 only_after - setdiff(common2, norm(common)) head(only_after, 20)规范化前后差出来的那批 ID 一定要人工看。我去过一个项目里样本号S-001和S001同时存在规范化后合并了实际是两个不同的样本这种误合并比丢样本更麻烦因为它不会报错只会让某个蛋白的关联结果悄悄偏移。对齐完成后把三张表统一裁到common并保持同一顺序。顺序这件事必须显式做不能指望读入顺序天然一致——merge默认会按 key 重排match不会两处混用就会出现「表型是第 3 个样本、基因型是第 5 个」的静默错位。# 统一裁切并强制同一顺序杜绝静默错位 prot_mat - prot_mat[common, , drop FALSE] cov_df - cov_df[common, , drop FALSE] fam - fam[match(common, fam$IID), , drop FALSE] stopifnot(identical(rownames(prot_mat), rownames(cov_df))) stopifnot(identical(rownames(prot_mat), fam$IID))stopifnot这两行看着啰嗦但它把「顺序一致」从约定变成了断言。后面无论怎么重构代码只要这两行还在错位就不可能悄悄发生。这是我从一次通宵排查里换来的习惯。2.3 一次能跑通的读入与预处理脚本蛋白丰度进关联模型前必须做变换这一步跳过去后面所有 p 值都不可信。SomaScan 这类 aptamer 平台的原始信号是右偏的Olink 的 NPX 已经是 log2 尺度但仍然存在样本间总体强度差异。通用做法是三步低检出蛋白过滤、log2 变换若尚未在 log 尺度、逆秩正态变换INT。# 蛋白丰度预处理过滤 - log2 - 逆秩正态变换 prep_prot - function(mat, min_detect 0.5, already_log2 FALSE) { # 1) 低检出蛋白过滤非缺失比例低于阈值的蛋白直接剔除 keep - colMeans(!is.na(mat)) min_detect mat - mat[, keep, drop FALSE] # 2) 未在 log 尺度的话先取 log2压右偏 if (!already_log2) mat - log2(mat 1) # 3) 逆秩正态变换把每个蛋白的秩映射到标准正态分位数 int_norm - function(v) { ok - !is.na(v) qnorm((rank(v[ok], ties.method average) - 0.5) / sum(ok)) } mat[] - lapply(seq_len(ncol(mat)), function(j) { out - rep(NA_real_, nrow(mat)) out[!is.na(mat[, j])] - int_norm(mat[, j]) out }) attr(mat, n_protein) - ncol(mat) mat } prot_norm - prep_prot(prot_mat, min_detect 0.5, already_log2 FALSE)min_detect 0.5是个起点不是铁律。蛋白数少、样本珍贵的时候可以放到 0.3代价是这批蛋白的检验效能很低后面得单独标注。ties.method average处理的是丰度值并列的情况质谱数据里并列不少用默认的first会让结果依赖行顺序换台机器跑就变。基因型侧我一般不在 R 里做质控直接用 PLINK 过一遍更省事也便于留下可复现的日志。# 基因型质控MAF、缺失率、HWE同时只保留常染色体 plink --bfile geno \ --autosome \ --maf 0.01 \ --geno 0.05 \ --hwe 1e-6 midp \ --make-bed --out geno_qc--maf 0.01对 pQTL 来说偏宽松rare variant 在蛋白关联里偶尔有惊喜但样本量不到几千的话检验效能撑不住会带来大量假阳性的顺式信号。--hwe 1e-6 midp用精确检验比渐近检验在小样本下稳。--autosome是为了跟后面 cis 区间的染色体命名对上避免chr1和1混用。跑完记得把.bim的第四列物理位置和你的基因组版本核对一次。这一列如果和蛋白注释用的版本不一致cis 判定会整体错位而且不会报任何错——这是我在第 5 章要重点讲的一个坑。3. 关联建模cis 窗口、协变量与多重检验的取舍3.1 cis 窗口定多大1Mb 只是默认值cis 和 trans 的划分是 pQTL 分析里最容易被当成「定死了的事」其实并不固定的一环。文献里 ±1Mb 最常见因为早期 eQTL 分析沿用了这个窗口但蛋白编码基因的调控结构跟基因表达不完全一样窗口大小会直接影响 cis 信号的召回率。窗口定小了落在 ±200kb 之外的远端调控 SNP 会被判成 trans等于把真 cis 信号扔进了多重检验更严的那一堆窗口定大了会把同一染色体上相邻基因的调控区卷进来cis 结果里混进假的近端 trans。我的做法是先按 ±1Mb 跑一版再对显著蛋白单独试 ±500kb 和 ±2Mb看信号边界落在哪儿把这个边界写进方法描述里而不是笼统写一句「±1Mb」。# 由蛋白编码基因坐标构造 cis 窗口供后续 SNP 归属判定 library(GenomicRanges) make_cis_window - function(gene_gr, flank 1e6) { # gene_gr 需含 seqnames / start / end / gene_id 四列 w - promoters(gene_gr, upstream flank, downstream flank) # promoters() 对负链基因的 upstream/downstream 语义会翻转 # 这里改成直接对 start/end 做算术避免链方向干扰 w - GRanges(seqnames(gene_gr), IRanges(pmax(1, start(gene_gr) - flank), end(gene_gr) flank)) mcols(w)$gene_id - mcols(gene_gr)$gene_id w } cis_win - make_cis_window(gene_gr, flank 1e6)上面刻意没用promoters()就是因为它对负链的upstream定义是「转录起始位点往上游」染色体坐标上反而更大混在批量循环里会出一半对一半错的结果。直接用start - flank/end flank加pmax(1, ...)兜住染色体起点越界行为可预测。窗口建好后把质控过的 SNP 按坐标注入用findOverlaps判定归属。这一步的输出是一张SNP × gene的布尔表后面关联时按它分流 cis 和 trans。# SNP 归属 cis / trans按坐标重叠判定不依赖 SNP 名 snp_gr - GRanges(snp_info$chr, IRanges(snp_info$pos, width 1)) seqlevels(snp_gr) - seqlevels(cis_win) # 统一染色体命名后再比 hits - findOverlaps(snp_gr, cis_win) cis_map - data.frame( snp snp_info$snp[queryHits(hits)], gene mcols(cis_win)$gene_id[subjectHits(hits)], stringsAsFactors FALSE )seqlevels(snp_gr) - seqlevels(cis_win)这行不能省。一个来自 PLINK、一个来自注释包1和chr1、MT和chrM的差异会让findOverlaps返回空结果而不报错。这也是血泪经验之一。3.2 线性模型还是混合模型绝大多数 pQTL 用lm就够了前提是协变量里已经放进了足够的前若干个遗传主成分。样本存在明显亲缘结构、或者同一个体贡献了多个时间点的样本时才需要上混合模型。判断标准很实际如果同一个蛋白在lm和lme4::lmer个体作随机截距下的 beta 差了两个数量级说明样本间相关性没被协变量吸收干净这时候上混合模型如果两者接近lm跑得快、结果好解释没必要为了「看起来更严谨」硬上。# 单蛋白全 SNP 扫描lm 版本协变量含年龄性别批次 前 10 个 PC scan_one_protein_lm - function(y, geno_mat, cov_df, snp_info) { ok - !is.na(y) y - y[ok] X - as.matrix(cov_df[ok, , drop FALSE]) res - lapply(seq_len(ncol(geno_mat)), function(j) { g - geno_mat[ok, j] if (sd(g, na.rm TRUE) 0) return(NULL) # 单态 SNP 跳过 fit - lm(y ~ g X) co - summary(fit)$coefficients if (!g %in% rownames(co)) return(NULL) data.frame(snp snp_info$snp[j], beta co[g, Estimate], se co[g, Std. Error], p co[g, Pr(|t|)]) }) do.call(rbind, res) }sd(g, na.rm TRUE) 0这个跳过条件看着不起眼但质控后的.bim里仍然可能残留单态位点比如只在某个亚群里有变异、整体 MAF 恰好过了阈值。不跳过的话lm会给NA系数rbind之后再筛会多出一堆无意义的行。lm(y ~ g X)里X是矩阵R 会自动展开成多个协变量列。这个写法比分步residualize再回归要稳因为后者在协变量有缺失时会悄悄丢样本而lm的na.action是统一的。混合模型版本结构一样只是把lm换成lmer并加(1 | individual)# 有重复样本时改用随机截距模型 library(lme4) scan_one_protein_lmm - function(y, geno_mat, cov_df, ind_id, snp_info) { ok - !is.na(y) dat0 - data.frame(y y[ok], cov_df[ok, , drop FALSE], ind factor(ind_id[ok])) res - lapply(seq_len(ncol(geno_mat)), function(j) { dat0$g - geno_mat[ok, j] if (sd(dat0$g, na.rm TRUE) 0) return(NULL) fit - lmer(y ~ g . - ind (1 | ind), data dat0, REML FALSE) co - summary(fit)$coefficients data.frame(snp snp_info$snp[j], beta co[g, Estimate], se co[g, Std. Error], p co[g, Pr(|t|)]) }) do.call(rbind, res) }REML FALSE是必须的混合模型做固定效应检验时要用 ML 估计REML 下的似然比不能直接比。这个参数忘了设lmerTest给的 p 值会偏保守显著位点会少一批。3.3 分层跑法与并行、内存控制规模一上来「一个蛋白一个蛋白跑」的循环很快就会撞上两个墙时间和内存。时间上上千蛋白 × 百万 SNP 的单线程lm是按天算的内存上如果每个蛋白的结果都rbind到一个大表里几百个蛋白之后 R 会话就开始吃满。我的做法是分层外层按蛋白分批内层按染色体分块每块跑完立刻落盘主进程只保留汇总后的显著结果。# 分批 分块每块跑完立刻写盘主进程不留全量结果 library(BiocParallel) run_pqtl_chunked - function(prot_norm, geno_mat, cov_df, snp_info, batch_size 50, out_dir pqtl_out) { dir.create(out_dir, showWarnings FALSE) prot_ids - colnames(prot_norm) for (b in seq(1, length(prot_ids), by batch_size)) { idx - b:min(b batch_size - 1, length(prot_ids)) chunk_res - bplapply(idx, function(k) { y - prot_norm[, k] r - scan_one_protein_lm(y, geno_mat, cov_df, snp_info) if (is.null(r)) return(NULL) r$protein - prot_ids[k] r }, BPPARAM MulticoreParam(workers 8)) chunk_df - do.call(rbind, chunk_res) saveRDS(chunk_df, file.path(out_dir, sprintf(batch_%04d.rds, b))) rm(chunk_res, chunk_df); gc() # 及时释放压住内存峰值 } }batch_size 50是平衡点批太小写盘和进程调度的开销占比过高批太大一批的结果在内存里堆着gc()也救不回来。workers 8要按机器的物理核数设超过物理核数只增加上下文切换实测没有收益。gc()在循环里显式调用某种程度上是玄学但配合rm()确实能把内存峰值压下来一截。更可靠的做法是在跑之前用ulimit或容器限额把内存卡死让它 OOM 报错而不是拖垮整台机器。结果汇总放在最后单独一步做只读显著的那部分# 汇总只把达到阈值的行拉进内存 files - list.files(pqtl_out, pattern ^batch_.*\\.rds$, full.names TRUE) sig - lapply(files, function(f) { d - readRDS(f) d[d$p 1e-5, , drop FALSE] # 粗筛阈值正式阈值在 4.3 定 }) sig - do.call(rbind, sig) saveRDS(sig, pqtl_suggestive.rds)1e-5只是个粗筛线目的是把全量结果压到能放进内存的体量。正式阈值要按 SNP 数和 LD 结构单独算见 4.3。4. 结果对象拆解从原始输出到可发表的顺式/反式显著表4.1 结果字段与二次整理关联跑完拿到的原始表通常长这样protein、snp、beta、se、p。这五列不足以支撑后面的分析至少还要补三块信息SNP 的基因组坐标、MAF、以及 cis/trans 标签。补齐的方式是 join而不是在关联循环里一并算——循环里加字段会把单次扫描的耗时抬高一个量级而这些信息本来就是常量表事后merge一次就够。# 补坐标、MAF 和 cis/trans 标签 enrich_res - function(res, snp_info, cis_map) { res - merge(res, snp_info[, c(snp, chr, pos, maf)], by snp, all.x TRUE, sort FALSE) # cis 标签用 snp gene 联合键匹配避免同名 SNP 误判 key_res - paste(res$snp, res$protein, sep |) key_cis - paste(cis_map$snp, cis_map$gene, sep |) res$is_cis - key_res %in% key_cis res[order(res$p), ] } res_all - readRDS(pqtl_suggestive.rds) res_all - enrich_res(res_all, snp_info, cis_map) table(res_all$is_cis)merge(..., sort FALSE)这个参数值得单独提一句默认的sort TRUE会按 key 重排行顺序变了不打紧但如果下游用行号索引过什么东西就会跟着错。all.x TRUE保证关联结果一行不少坐标补不上的行会留NA这正好是一个筛子——NA坐标的行说明 SNP 在注释表里找不到多半是染色体命名或版本问题要单独查。paste(snp, protein)联合键比单独按 SNP 匹配靠谱。一个 SNP 可以同时是基因 A 的 cis 和基因 B 的 trans只按 SNP 匹配会把 trans 也标成 cis。4.2 Manhattan 与 QQ 图怎么画才不误导pQTL 的 Manhattan 图和 GWAS 的看起来像但有一处结构性差异必须体现出来横轴不是按染色体排的 SNP而是按蛋白排的位点或者反过来。最常见的画法是「一个蛋白一张图、横轴按染色体排 SNP」密集的 cis 峰在图上非常显眼但多蛋白拼图时容易被误读成「这些蛋白共享同一个信号」。# 单蛋白 Manhattan横轴按染色体排cis 位点单独着色 plot_manhattan_one - function(res, protein_id, cis_flank_label TRUE) { d - res[res$protein protein_id !is.na(res$chr), ] d - d[order(d$chr, d$pos), ] d$chr - factor(d$chr, levels unique(d$chr[order(as.numeric( sub(^chr, , d$chr)))])) # 累积坐标跨染色体不留空隙 chr_len - tapply(d$pos, d$chr, max, na.rm TRUE) offset - c(0, cumsum(as.numeric(chr_len))[-length(chr_len)]) d$cumpos - d$pos offset[as.integer(d$chr)] plot(d$cumpos, -log10(d$p), pch 20, cex 0.4, col ifelse(d$is_cis, #B22222, #4D4D4D), xaxt n, xlab Chromosome, ylab expression(-log[10](italic(p))), main protein_id) axis(1, at offset as.numeric(chr_len) / 2, labels levels(d$chr), cex.axis 0.6, tick FALSE) abline(h -log10(5e-8), lty 2, col grey40) }sub(^chr, , ...)转数值是为了让染色体按 1、2、…、22、X 的顺序排直接按字符串排会得到 1、10、11、2 这种顺序图会乱。offset的算法是让每条染色体首尾相接cumsum少取最后一项是为了长度对齐。cis 用单独颜色标出来是因为它在 pQTL 图里的意义和 trans 完全不同cis 峰往往直接在基因座上位置本身就是证据trans 峰的位置几乎不提供信息只能靠 p 值判断。把两者画成一种颜色读图的人会下意识按位置解释 trans 信号。QQ 图要分层画全 SNP 一张、cis 单独一张。全 SNP 的 QQ 图在 pQTL 里几乎必然严重偏离对角线因为 cis 信号太强这不能说明模型有问题。# 分层 QQ全量 vs cis两条线叠在同一张图上 plot_qq_strat - function(res, protein_id) { d - res[res$protein protein_id, ] obs_all - -log10(sort(d$p)) obs_cis - -log10(sort(d$p[d$is_cis])) exp_all - -log10(ppoints(length(obs_all))) exp_cis - -log10(ppoints(length(obs_cis))) plot(exp_all, obs_all, pch 20, cex 0.4, col grey50, xlab expression(Expected~~-log[10](italic(p))), ylab expression(Observed~~-log[10](italic(p))), main paste0(protein_id, (QQ))) points(exp_cis, obs_cis, pch 20, cex 0.6, col #B22222) abline(0, 1, col black, lty 2) legend(topleft, c(all SNP, cis SNP), pch 20, col c(grey50, #B22222), bty n, cex 0.8) }ppoints()默认用(i - 0.5) / n在极端尾部比i / (n 1)更稳p 值很小的时候不会把期望值压到 0 导致-log10爆掉。4.3 阈值与 FDR分层校正的实际做法全基因组显著线5e-8是 GWAS 沿用的直接搬到 pQTL 上会偏严因为 pQTL 的 SNP 数通常比 GWAS 少一个量级。按 SNP 数做 Bonferroni 是个起点0.05 / n_snp。但更实际的做法是 cis 和 trans 分开定阈值。cis 检验的数比 trans 小两三个数量级用同一个阈值会两边都不合适trans 那边过严cis 这边过松。我一般按下面的口径# 分层阈值cis 用 Bonferronitrans 用 BH-FDR set_threshold - function(res, fdr_cut 0.05, alpha 0.05) { is_c - res$is_cis n_cis - sum(is_c) n_tr - sum(!is_c) res$thr_cis_bonf - alpha / max(n_cis, 1) res$fdr - NA_real_ res$fdr[!is_c] - p.adjust(res$p[!is_c], method BH) res$sig - ifelse(is_c, res$p res$thr_cis_bonf, !is.na(res$fdr) res$fdr fdr_cut) attr(res, n_cis) - n_cis attr(res, n_trans) - n_tr res } res_sig - set_threshold(res_all) table(cis res_sig$is_cis, significant res_sig$sig)p.adjust(..., method BH)只对 trans 那一部分做这是关键。如果把 cis 也混进去一起校正cis 那批极小的 p 值会把 FDR 的整体分布拉偏trans 的显著性会被系统性高估。分开算两边各自在自己的检验族里校正口径才干净。alpha / max(n_cis, 1)里的max(..., 1)是防零除某些蛋白在窗口内一个 SNP 都没有n_cis为 0 时阈值会变成Infp Inf恒真所有结果都被标成显著。这个坑不显眼但在蛋白覆盖度不均的数据里一定会撞上。5. 避坑与排查五个真会卡住的地方5.1 数据层读进来就已经错了现象三张表交集数看着正常跑出来的 cis 结果里某些蛋白的所有 SNP 都是NA系数summary(fit)里g那一行不存在。原因data.table::fread默认返回data.table取行名时rownames()拿到的是1、2这种自动编号不是真实样本 ID。后面的match全部按编号走看着对上了实际每个样本配的是另一份数据。另一种来源是 CSV 第一列没表头被读成V1行名彻底丢失。解决读入时显式指定行名列或者读成普通 data.frame 后手动设行名。加一道断言把「行名是否像自动编号」卡住。# 读入时锁定行名并断言不是自动编号 prot_df - data.table::fread(protein.csv, data.table FALSE) rownames(prot_df) - prot_df[[1]] prot_df[[1]] - NULL stopifnot(!all(grepl(^[0-9]$, rownames(prot_df))))现象所有蛋白的 cis 信号全部落在染色体两端中间几乎没有。原因.bim的物理位置列和注释包用的基因组版本不一致。比如基因型是 build 37、注释是 build 38同一个基因在两边差了几十兆碱基findOverlaps只在染色体靠边的地方偶然撞上。解决跑之前用少量已知位点核对坐标或者统一用 liftOver 转一遍。断言可以这样写抽查几个 MAF 接近 0.5 的常见位点比对两边位置差是否在千碱基量级以内超出就停下来查版本。5.2 建模层跑完结果不对现象cis 结果一批一批地空但 trans 里反而冒出大量集中在某条染色体上的「强信号」。原因基因型编码方向反了。PLINK 的.bed里 0/1/2 的含义取决于建库时的--recodeA或--keep-allele-order设置有些流程会把它转成 2/1/0。编码反了不影响 p 值但 beta 的符号全反而且 cis 判定用的坐标没错、信号却和预期基因方向对不上排查时会误以为是注释问题。解决在关联前加一步用table(geno_mat[, j])抽查几个位点确认 0 是纯合参考、2 是纯合替代。更稳的做法是统一从.bed用--recode A导出加性编码不要用管道中间产物。现象同一个蛋白lm和lmer的 beta 差了两个数量级或者加了某个协变量之后所有 cis 信号消失。原因协变量里混进了和蛋白本身强相关的变量。常见的是把「批次」和「样本处理时间」同时放进去两者共线更隐蔽的是把某个下游蛋白或临床指标当协变量而它本身就是被调控的产物等于把信号残差化掉了。解决进模型前先算一遍协变量之间的方差膨胀因子VIF 超过 10 的成对变量只留一个。同时给每个协变量做单变量回归看它和蛋白的关联强度关联过强的比如 |r| 0.7单独讨论要不要放。5.3 输出层结果对但用不了现象分批跑完之后汇总发现总行数比预期少一截某些蛋白完全没出现。原因bplapply返回的列表里NULL会被静默保留do.call(rbind, ...)遇到NULL会跳过它但不会报错。单态 SNP 全跳过的蛋白、或者整列都是NA的蛋白就这么消失了。解决汇总时先统计每批的NULL数量和预期蛋白数对一遍或者把NULL显式替换成带NA的占位行让缺口在结果里可见。# 汇总时保留缺口不静默丢蛋白 n_null - sum(vapply(chunk_res, is.null, logical(1))) if (n_null 0) warning(sprintf(本批有 %d 个蛋白无有效结果, n_null))现象重跑一次结果数值有细微差异显著位点名单对不上。原因并行 随机过程的组合没设种子。MulticoreParam默认的 RNG 流是LEcuyer-CMRG不显式设RNGseed时每次调度到哪个 worker 是不确定的用到了随机数的地方比如某些置换检验、bootstrap 置信区间结果就会漂。解决在MulticoreParam()里显式给RNGseed并且把结果按蛋白 ID 排序后再落盘把「顺序」这个自由度也锁死。顺序锁死之后即使数值有极小漂移也能立刻定位是哪个蛋白。6. 下游衔接共定位与 MR 的取数技巧pQTL 跑完只是拿到一张关联表真正支撑因果推断的是后面两步共定位和孟德尔随机化。这两步对输入格式的要求跟关联表不一样取数方式选错会在下游浪费掉大半天。共定位要的是「一个区域内的全 SNP 汇总统计」不是只挑显著位点。常见翻车是把p 5e-8的行筛出来喂给coloc结果 SNP 数不够、区域边界被截断PP.H4直接算成 0 或者NA。正确做法是按 cis 窗口把该区域内的所有 SNP 全取出来不管显著与否。# 共定位取数按区域取全 SNP不按显著性筛 extract_region - function(full_res, protein_id, chr, start, end, flank 2e5) { d - full_res[full_res$protein protein_id full_res$chr chr full_res$pos (start - flank) full_res$pos (end flank), ] # coloc 要求snpid / beta / varbeta / p / N缺一不可 data.frame(snp d$snp, beta d$beta, varbeta d$se^2, p d$p, N d$n_sample) }varbeta se^2这列很多人会漏coloc的type quant模式必须显式给缺了会直接报错。flank 2e5比 cis 窗口再放宽一截是为了让共定位有足够的背景 SNP 估计先验窗口卡太紧会让PP.H4偏乐观。孟德尔随机化那边取数逻辑反过来要的是「工具变量」也就是达到全基因组显著的独立位点而且必须做 LD 剪枝。常见翻车是把同一区域内的强 LD 位点全塞进去等于把同一个工具变量复制了好几遍IVW 的方差会被严重低估。# MR 工具变量取数显著 LD 剪枝 取每簇最强位点 library(data.table) pick_iv - function(full_res, protein_id, p_cut 5e-8, clump_kb 500) { d - full_res[full_res$protein protein_id full_res$p p_cut, ] if (nrow(d) 0) return(NULL) setDT(d); setorder(d, p) kept - integer(0) for (i in seq_len(nrow(d))) { if (length(kept) 0) { kept - i; next } prev - kept same_chr - d$chr[prev] d$chr[i] near - abs(d$pos[prev] - d$pos[i]) clump_kb * 1000 if (!any(same_chr near)) kept - c(kept, i) # 窗口内已有更强位点则丢弃 } d[kept, ] }clump_kb 500是个粗剪枝严格做法应该用参考面板算 r²常见阈值 0.001。用距离剪枝在样本量不大时够用但写方法部分时要如实标注用的是距离而非 r²否则审稿阶段会被追问。取完工具变量还要和 GWAS 侧的暴露数据对齐等位基因。这一步我用match按 rsID 匹配后显式检查效应等位基因是否一致不一致的取 beta 的相反数而不是简单丢弃——丢弃会让工具变量数减少MR 的效能跟着掉。我从这几个项目里养成的习惯是每次拿到的关联结果先不急着往下游走先花十分钟把「有多少 SNP、坐标版本是什么、beta 方向对不对」这三件事各抽查一遍。翻车的成本远高于这十分钟。希望帮到你。本文还有配套的精品资源点击获取