简介：pQTLtools 是面向蛋白质组学与遗传学研究的 R 语言工具包，聚焦蛋白质定量性状基因座（pQTL）分析流程，适合具备一定 R 基础、需要处理 Olink、SomaLogic、Caprion 等亲和力面板数据的科研人员与生信分析者，可缓解不同平台蛋白标记物来源分散、基因组坐标不统一带来的整合困难。压缩包共 125 个文件、约 5.92MB，其中 21 个 html 与 9 个 rd 构成 pkgdown 文档站点和函数帮助，10 个 R 源码与 4 个 Rmd 呈现具体实现，11 个 rda 存放精选数据集，另有 15 个 png 图示以及 css、js、yml 等站点与构建配置。包内按 Caprion、Olink 的 NGS 与 qPCR 版本、SomaLogic、INTERVAL 论文 ST4 与 ST6 等面板分类整理数据，并附带基于 UCSC 基因组浏览器的 hg19 注释表及 hg19 到 hg38 的链转换文件，便于统一基因坐标、跨平台比对蛋白标记物。目前已有 759 人学习下载，适合作为 pQTL 数据处理与流程搭建的参考起点。

1. 别急着跑关联：pQTLtools 把 pQTL 的五步流水线收进一个 R 包

手上有一批蛋白丰度矩阵，再加一批 PLINK 格式的基因型，被问一句「这些蛋白有没有 cis-pQTL」，很多人的第一反应是打开表格手工对 rsID，再写两百行 for 循环套lm()。跑是能跑，跑完发现样本对不上、效应等位基因翻错、批次没进协变量，前面两周白干。pQTLtools 是一个用 R 写的 pQTL 分析工具，把蛋白丰度、基因型剂量、协变量三者的对齐、关联检验、cis/trans 注释和结果可视化串成一条可重跑的流水线，中间那些「不报错但结果全错」的环节被显式做成了检查项。它适合手上有几十到上千个蛋白、样本量几百到几千、需要产出一份能进补充材料的 pQTL 结果表的人；如果你只做单个位点的验证，直接写lm()更快，不必上这套。

2. 三件套对齐：把蛋白丰度、基因型、协变量压成同一张表

pQTL 分析里 80% 的翻车不在模型，在数据准备。三份输入分别来自蛋白检测平台、基因型检测平台和临床/实验记录，样本 ID 命名习惯各不相同，缺失模式也各不相同。这一章的目标很朴素：让三个矩阵的行列索引严格一致，任何一处对不上就当场报错，而不是等到关联结果出来才发现某个蛋白的 n 少了 40。

2.1 输入的最小可用集与格式约定

先把四类文件摆齐，缺一个都别往下走。包里的函数名以ls("package:pQTLtools")的实际输出为准，下面按主流程顺序给出调用，参数名建议先用args()对一遍再抄，不同版本的默认值会有出入。

基因型一定要先做一轮质控，否则后面所有的 MAF 过滤都是在垃圾上做筛选：

# SNP 层面与样本层面各做一轮，阈值按阵列密度调整 plink --bfile raw \ --maf 0.01 \ --geno 0.05 \ --hwe 1e-6 midp \ --mind 0.05 \ --make-bed --out qc # 导出加性编码剂量（0/1/2），列名形如 rs1234_A，A 即被计数的效应等位基因 plink --bfile qc --recode A --out qc_dosage

--maf 0.01对 pQTL 来说是底线，稀有变异在几百到几千样本里几乎没有检验效能，留着只会把多重检验的分母撑大；--hwe 1e-6 midp用精确检验，是因为 pQTL 研究样本量普遍不大，渐近检验在小样本、低 MAF 位点上偏差明显；--mind 0.05是样本层面的缺失率，比 SNP 层面的--geno更要紧——一个样本在某个蛋白上缺失，可能整行都进不了模型。--recode A输出的列名后半段就是你后续做等位基因对齐的依据，务必保留。

读入 R 之后先做一次维度自检，别急着合并：

library(pQTLtools) prot <- read_protein("npx_matrix.tsv", id_col = "Protein", log2 = TRUE) geno <- read_dosage("qc_dosage.raw", id_col = "FID") covar <- read_covar("covariates.tsv", id_col = "IID") # 三个矩阵的样本集合必须能取到交集，交不出来就说明 ID 规则不一致 stopifnot(length(intersect(colnames(prot), colnames(geno))) > 100) # 蛋白缺失率：整行全缺失的蛋白直接剔除，别指望插补救回来 keep_prot <- rowMeans(is.na(prot)) < 0.2 prot <- prot[keep_prot, , drop = FALSE]

log2 = TRUE针对的是已经做过平台内标准化的 NPX 值，如果你的输入是原始荧光强度，先做 log2 再做中心化；rowMeans(is.na(prot)) < 0.2这条阈值不是硬标准，蛋白缺失率超过 20% 的位点即使插补出来，关联结果的可信度也很低，我一般直接剔。

2.2 样本 ID 与等位基因方向的对齐

样本 ID 的问题通常很具体：蛋白矩阵里叫S-001，基因型 fam 里叫S001，协变量表里叫sample_1。写一个显式的映射函数，别用gsub一把梭，因为一旦某条规则误伤（比如把样本号里的-全删掉，结果A-1-2和A-12撞车），后果是样本错配且不报错。

norm_id <- function(x) { x <- toupper(trimws(x)) x <- gsub("[_\\-\\.]", "", x) # 统一去掉分隔符 x <- sub("^SAMPLE", "S", x) # 前缀归一 x } colnames(prot) <- norm_id(colnames(prot)) colnames(geno) <- norm_id(colnames(geno)) covar$IID <- norm_id(covar$IID) common <- Reduce(intersect, list(colnames(prot), colnames(geno), covar$IID)) # 交集之后必须严格按同一顺序重排，否则后面的矩阵运算会静默错位 prot <- prot[, common, drop = FALSE] geno <- geno[, common, drop = FALSE] covar <- covar[match(common, covar$IID), , drop = FALSE]

等位基因方向是更隐蔽的一类错误。--recode A给出的计数等位基因取自 bim 文件的 A1，不同批次、不同来源的 bim 里 A1/A2 顺序可能相反。方向翻错的直接后果是 beta 符号整体反号，p 值一模一样——它不会报错，只会让你在写讨论的时候把「升高」写成「降低」。

# 从列名里取回被计数的效应等位基因，统一成大写 ea <- toupper(sub("^.*_", "", colnames(geno))) # 用样本内部等位基因频率做一致性检查：频率应落在 (0,1) 且与剂量均值吻合 eaf <- rowMeans(geno, na.rm = TRUE) / 2 bad <- which(eaf < 0.01 | eaf > 0.99) if (length(bad)) geno <- geno[-bad, , drop = FALSE] # 结果表里必须同时输出 ea / oa / eaf 三列，方便下游按参考面板统一方向

我一般会在结果表里强制保留ea（效应等位基因）、oa（其他等位基因）、eaf三列，后续不管是和外部数据集做 meta 分析，还是和公共参考面板比对方向，都靠这三列翻盘。少了它们，半年后回看结果表只能靠猜。

2.3 协变量矩阵：批次、年龄、性别与前十个主成分

协变量放什么，直接决定 QQ 图的形状。pQTL 里最少要放：年龄、性别、检测批次（plate/run）、前 5～10 个基因型主成分。蛋白丰度本身受技术批次影响极大，跨批次合并而不校正，几乎必然产生基因组膨胀。

# 基因型主成分：用 LD 剪枝后的位点算，避免高 LD 区域主导 system("plink --bfile qc --indep-pairwise 50 5 0.2 --out pruned") system("plink --bfile qc --extract pruned.prune.in --pca 10 --out pcs") pcs <- read.table("pcs.eigenvec", header = FALSE) colnames(pcs) <- c("FID", "IID", paste0("PC", 1:10)) pcs$IID <- norm_id(pcs$IID) covar <- merge(covar, pcs[, c("IID", paste0("PC", 1:10))], by = "IID") covar <- covar[match(common, covar$IID), , drop = FALSE] # 共线性自检：设计矩阵的秩必须等于列数，否则某些协变量是其他列的线性组合 X <- as.matrix(covar[, c("age", "sex", "plate", paste0("PC", 1:10))]) X <- cbind(1, X) if (qr(X)$rank < ncol(X)) warning("协变量存在共线，检查 plate 与批次分组是否重复编码")

--indep-pairwise 50 5 0.2的窗口和 r² 阈值是常规起点，位点密度高的时候把窗口调小。共线性检查这一步看着多余，实际很值得：plate如果和另一个「批次」变量是一一对应的，模型会把方差拆到两个变量上，系数估计不稳定，条件分析阶段尤其容易出问题。蛋白侧的 PEER 这类概率性因子推断方法也可以用，但对样本量有要求，样本少于 100 时因子估计本身就不稳，不如老老实实放几个主成分。

3. 关联检验：加性模型、cis 窗口与全基因组阈值

数据对齐之后，模型本身反而不复杂：每个蛋白对每个位点跑一次加性线性回归。真正需要想清楚的是三件事——协变量怎么写进公式、cis 窗口取多大、以及多重检验的分母怎么定。这三件事定错了，结果表的解读会整体跑偏。

3.1 线性模型怎么写：NPX ~ 剂量 + 协变量

加性模型是 pQTL 的默认选择，因为绝大多数已报道的蛋白定量性状位点效应可以用剂量线性解释。乘性模型、显性模型只在有明显证据时才试，全基因组扫一遍 dominance 模型的代价是多一倍的计算量和一堆难以解释的显著点。

res <- pqtl_assoc( protein = prot, # 蛋白 × 样本，已 log2 并中心化 geno = geno, # 位点 × 样本，0/1/2 加性剂量 covar = covar[, c("age", "sex", "plate", paste0("PC", 1:10))], model = "additive", # 加性模型；"dominant"/"recessive" 按需 maf = 0.01, # 与上游 QC 保持一致，不要在这里放宽 missing = 0.05, # 位点级缺失率上限 min_n = 200, # 有效样本数下限，低于此不出结果 threads = 8 )

参数里最容易被忽略的是min_n。默认不设的话，一个缺失严重的蛋白可能在 80 个样本上「显著」，p 值漂亮得离谱，实则是小样本噪声。missing = 0.05是位点级过滤，蛋白级的缺失在 2.1 已经处理过，两层过滤不要混。threads建议设成物理核数的一半，pQTL 关联是内存带宽敏感型任务，开满核往往不升反降。

如果不用包内的封装，等价的手写形式是逐个蛋白建公式，注意协变量必须作为整体一次性放进lm()，不要先对表型做残差化再回归——两步法在样本有缺失时会让残差化和基因型的样本集不一致。

fit_one <- function(y, g, cv) { ok <- !is.na(y) & !is.na(g) d <- data.frame(y = y[ok], g = g[ok], cv[ok, , drop = FALSE]) fit <- lm(y ~ g + ., data = d) cf <- summary(fit)$coefficients["g", ] c(beta = cf[1], se = cf[2], p = cf[4], n = sum(ok)) }

3.2 cis 窗口取多大：±500kb 还是 ±1Mb

cis 的定义直接影响结果表里 cis 位点的数量。常见做法是以蛋白对应基因的 TSS 为中心，取 ±500kb 或 ±1Mb。窗口越大，捕获远端调控变异的概率越高，但同时每个蛋白要检验的位点数增加，Bonferroni 阈值变严，弱效应容易被淹没。

一个蛋白对应多个基因的情况必须处理：蛋白检测探针可能映射到多个 ENSG，或者蛋白是某个基因的产物但注释有两个转录本 TSS。我的做法是把所有候选 TSS 都保留，一个位点只要落在任意一个 TSS 的窗口内就记为 cis，同时在结果表里记下最近的 TSS 距离，方便后续人工判断。

anno <- read_gene_anno("gene_tss.tsv") # 列：ENSG, chr, tss map <- read_probe_map("probe2gene.tsv") # 列：Protein, ENSG res <- annotate_cis(res, map = map, anno = anno, window = 500000) # cis 与 trans 分开标记，后续阈值不同 table(res$cis)

窗口参数window单位是碱基，写成 500000 而不是 500。这一步还应该顺手把「位点-蛋白」的距离算出来，因为后续做共定位和条件分析时，距离近的信号优先级更高。

3.3 结果表的列含义与 FDR 控制

一张合格的 pQTL 结果表至少有这些列，缺任何一列后面都会卡住：

多重检验的分母要分两层定。cis 检验里每个蛋白的位点数在几千量级，用蛋白内 Bonferroni（0.05 / 位点数）或者对 cis 子集单独做 BH 都可以；trans 检验的位点数是全基因组量级，用 5e-8 这类常规全基因组阈值，或者对 trans 子集做 BH。把 cis 和 trans 混在一起做一次 BH，会导致 cis 信号被 trans 的巨大分母稀释掉。

res$fdr <- NA is_cis <- res$cis res$fdr[is_cis] <- p.adjust(res$p[is_cis], method = "BH") res$fdr[!is_cis] <- p.adjust(res$p[!is_cis], method = "BH")

4. 从显著位点到可信信号：条件分析、共定位与出图

一张有几百行显著位点的结果表，离「可信信号」还有距离。同一个 cis 区域里往往有多个强 LD 的位点同时在跑分，看起来像多个独立信号，实际是一个；一个 pQTL 和某个疾病的 GWAS 位点位置重合，也不代表它们共享同一个因果变异。这一章处理的就是这两类判断。

4.1 条件分析：同一个 cis 区域里的第二个信号

逐步条件分析的做法很直白：把当前最显著的位点作为协变量加进模型，重新扫这个区域，如果还有位点显著，说明存在第二个独立信号。关键是加进去的位点要先做 LD 剪枝，否则两个 r² > 0.9 的位点互为协变量，模型直接崩。

cond <- pqtl_conditional( protein = prot, geno = geno, covar = covar[, c("age", "sex", "plate", paste0("PC", 1:10))], region = list(chr = "6", start = 25000000, end = 35000000), ld_r2 = 0.1, # 加入协变量的位点之间 LD 上限 max_iter = 5 # 最多找 5 个独立信号 )

ld_r2 = 0.1是常用起点，样本量小的时候可以放宽到 0.2，但不要超过 0.2，否则条件模型会出现 VIF 爆表。max_iter不必设太大，超过 3 个独立 cis 信号的区域已经很少见，设到 5 是给自己留个观察窗口。条件分析的结果要和新模型的样本量对照看——每加一个协变量位点，缺失会导致 n 略微下降，如果 n 掉得厉害，说明位点缺失模式有问题。

4.2 共定位：pQTL 与 eQTL/疾病 GWAS 是不是同一个信号

位置重合不等于共定位。判断两个性状是否共享因果变异，常规做法是接coloc做近似贝叶斯因子计算，输入需要两个性状在同一批位点上的 beta 和 varbeta，以及各自的 MAF。

library(coloc) d1 <- list(beta = res_pqtl$beta, varbeta = res_pqtl$se^2, snp = res_pqtl$snp, MAF = res_pqtl$eaf, N = median(res_pqtl$n), type = "quant") d2 <- list(beta = res_gwas$beta, varbeta = res_gwas$se^2, snp = res_gwas$snp, MAF = res_gwas$eaf, N = 50000, type = "cc", s = 0.4) # s 为病例比例 coloc_res <- coloc.abf(d1, d2) coloc_res$summary["PP.H4.abf"] # 共享因果变异的后验概率

type = "cc"时s必须给对，病例比例填错会让后验概率整体偏移。样本重叠是这里最大的陷阱：如果 pQTL 的样本和 GWAS 样本有大量重叠，两个性状的效应估计不独立，PP.H4 会被系统性高估。公开数据做共定位时，先查两边的样本来源，能查到重叠就做敏感性分析，查不到就在方法里明确写出这是未校正的结果。

4.3 出图：曼哈顿图、QQ 图与顺式区域图

出图不只是为了好看，它是排错的主要手段。QQ 图看整体膨胀，曼哈顿图看有没有染色体级别的异常堆积，顺式区域图看信号是不是集中在 TSS 附近。

# 全基因组曼哈顿图，只标 cis 且 FDR < 0.05 的位点 plot_manhattan(res, chr_col = "chr", pos_col = "pos", p_col = "p", highlight = res$cis & res$fdr < 0.05, label_top = 10, threshold = 5e-8) # QQ 图：lambda 超过 1.1 就该回头查协变量 plot_qq(res, p_col = "p", ci = 0.95) # 单个蛋白的 cis 区域图，x 轴用相对 TSS 的距离更好读 plot_cis(res, protein = "PROT_0012", anno = anno, window = 500000)

label_top = 10只标前 10 个点，标太多图就没法看了。threshold = 5e-8画一条参考线，实际判断以结果表里的fdr为准。QQ 图上的置信带在样本量小时会很宽，别看到偏离就慌，先看 lambda 的数值。

5. 排查：pQTL 跑出来全是「假阳性」时的五条血泪经验

下面五条都是我自己或者身边人实际踩过的，共同点是：不报错，但结果全错或者不可信。放在这里当检查清单用。

5.1 结果层面的三个异常

现象一：所有蛋白的 QQ 图整体上扬，lambda 到 1.2 以上。原因通常不是「有真实信号」，而是群体分层或者批次效应没进模型——蛋白丰度对技术批次极其敏感，跨 plate 合并而不校正，炎症相关蛋白会集体偏移。解决：先确认协变量里有plate，再补基因型主成分；如果仍然膨胀，按 plate 做分层分析或者对蛋白做批次内中心化，看 lambda 是否回落到 1.05 附近。

现象二：cis 位点的 beta 符号和已发表结果相反，p 值却完全一致。典型的效应等位基因方向问题。p 值对符号不敏感，所以这类错误永远不会触发报错。解决：把结果表里的ea、eaf和基因型原始列名逐条核对，统一以--recode A输出的计数等位基因为准，或者按参考面板的等位基因频率做一次方向统一。

现象三：结果表里同一个蛋白的n逐行跳动，最大最小差出几十。原因是个位点独立删除缺失样本，不同位点的有效样本集不一样。这在位点缺失率不高时影响有限，但如果某个蛋白的缺失和基因型平台批次相关，会引入结构性偏差。解决：先做样本级 QC，把缺失率高的样本整体剔除，把min_n设成一个硬门槛，同时在结果表保留n列供审稿人核查。

5.2 运行层面的两个异常

现象四：条件分析加到第二个位点后模型系数爆炸，或者直接报奇异矩阵。两个位点强 LD，互为协变量后设计矩阵近似共线。解决：把ld_r2压到 0.1 以下再跑，或者先对候选位点做 LD 剪枝，只把剪枝后的位点集喂给条件分析；同时用car::vif()检查一遍，VIF 超过 5 就说明还有共线。

现象五：全基因组 × 全蛋白的扫描跑到一半内存爆掉。把基因型一次性展开成稠密矩阵是最常见的原因，位点数上万、样本上千的时候内存占用是平方级增长的。解决：按染色体分块读入、算完就释放，threads控制在物理核数一半，蛋白侧也分批处理；如果条件允许，把剂量矩阵换成稀疏存储，只保留非零项。

6. 进阶与验收：用阳性对照位点给自己的流程打分

流程能跑通只是第一步，能证明它跑对了才是交付标准。我习惯在上生产之前做三层验收，成本不高，但能挡住绝大多数低级错误。

第一层是阳性对照。挑几个已有充分文献报道的顺式作用位点作为锚点，这些位点在不同人群、不同平台上的效应方向应该稳定。跑完流程后检查它们是否出现在结果表的 cis 显著位点里，方向是否一致。这一步不通过，后面的所有新发现都不用看。

第二层是置换检验。把蛋白丰度矩阵的样本顺序整体打乱，重跑一遍关联，统计 p < 0.05 的比例。理论上这个比例应该接近 5%，如果跑到 20% 以上，说明模型里还有未校正的结构。

set.seed(20240101) n_perm <- 200 perm_rate <- numeric(n_perm) for (i in seq_len(n_perm)) { # 只打乱样本维度，位点与协变量的对应关系保持不变 prot_perm <- prot[, sample(ncol(prot)), drop = FALSE] colnames(prot_perm) <- colnames(prot) r <- pqtl_assoc(protein = prot_perm, geno = geno, covar = covar[, c("age", "sex", "plate", paste0("PC", 1:10))], maf = 0.01, missing = 0.05, threads = 4) perm_rate[i] <- mean(r$p < 0.05, na.rm = TRUE) } summary(perm_rate) # 中位数应接近 0.05，明显偏高说明模型仍有系统性偏差

n_perm = 200是精度和耗时的折中，200 次足够看出中位数是否偏离 0.05 一个量级以上；如果中位数在 0.05 附近但尾部很宽，说明某些蛋白的方差估计不稳，需要回头检查这些蛋白的缺失率和离群值。置换时只打乱样本维度，位点与协变量的对应关系保持不动，否则检验的就不是「表型与基因型无关」这个零假设了。

第三层是参数敏感性。同一份数据，用不同的窗口和阈值各跑一遍，看核心结论是否稳定：

如果某个 cis 位点只在 ±1Mb 窗口下出现、MAF 又接近 0.01，它大概率不是稳健信号，写进正文要谨慎。反过来，核心位点在四组参数下都稳定，才值得往下做共定位和 MR。

从那以后我每次跑完 pQTL，都强制走一遍「阳性对照 + 200 次置换」这两步，再去看新发现——顺序不能反，先看新结果容易被漂亮的 p 值说服。这套验收习惯救过我至少两次，一次是发现某个批次的蛋白样本 ID 全被错误归一化，另一次是揪出协变量表里plate被重复编码。希望帮到你。

本文还有配套的精品资源，点击获取