简介pQTLtools 是一套面向蛋白质定量性状基因座pQTL研究的 R 语言工具包适合从事遗传学、蛋白质组学与生物信息分析的科研人员及研究生使用用于整合与处理 Olink、SomaLogic、Caprion 等不同平台的面板数据并完成从数据整理到分析流程搭建的工作。资源包共 125 个文件约 5.92MB以 html 文档、png 图示、rda 数据、r 脚本与 rd 帮助文件为主另含 sh 脚本、rmd 报告、css 样式及少量 tsv、yaml 配置覆盖文档、示例数据与源码多个层面。目前已有 758 人学习下载。包内提供生物基质、hg19 表、Olink/INF、Olink_NGS、Olink_qPCR、SomaLogic 及 INTERVAL ST4/ST6 等精选数据集并附网络驱动文档与功能摘要便于读者直接调用示例数据、理解各面板结构快速搭建可复现的 pQTL 分析环境。1. 从一份 R 包说起pQTLtools 到底能帮你省掉哪些重复劳动如果你做过蛋白质数量性状位点pQTL分析大概率经历过这样的循环把质谱或 Olink 的蛋白丰度矩阵读进来跟基因型数据对齐样本一个位点一个位点跑关联跑完再做法统计、画曼哈顿图、做共定位。每一步单拎出来都不难难的是把它们串成一条不会在样本 ID 上翻车的流水线。pQTLtools 就是冲着这个场景来的——它是一套以 R 为载体的 pQTL 分析工具集把数据预处理、关联检验、结果整理和可视化这几段常见流程收拢到统一的函数接口里。适合两类人一是刚接触 pQTL、想先跑通一条最小可用流程的新手二是已经有一堆散装脚本、想找个参照实现来对照参数和中间结果的老手。它不替你解决生物学问题但能让你少写几百行胶水代码。2. 环境准备与数据读入把样本对齐这件事做扎实2.1 为什么 pQTL 分析的第一步永远是样本对齐pQTL 的输入至少有两块蛋白丰度矩阵行是蛋白列是样本或者反过来和基因型矩阵行是变异位点列是样本。这两块数据往往来自不同批次、不同平台样本命名规则可能一个是S001、一个是sample_1直接合并就会出现「列名对不上但代码不报错」的玄学问题——R 在按列名取子集时匹配不上的列会被静默丢掉最后你拿到一个样本量缩水但看起来正常的矩阵。所以任何 pQTL 流程的第一步都不是跑模型而是把样本 ID 对齐并显式检查交集大小。pQTLtools 在这块的思路是先读入再取交集再按交集重排每一步都让你能看到丢了多少样本。2.2 读入蛋白丰度矩阵与基因型矩阵常见做法是用data.table或readr读大矩阵因为 pQTL 的蛋白矩阵动辄几千行、几百到几千列base R 的read.csv在速度上会吃亏。下面这段是读入并做基本检查的骨架library(data.table) # 读入蛋白丰度矩阵第一列是蛋白 ID其余列是样本 protein_dt - fread(protein_abundance.csv, data.table FALSE) # 读入基因型矩阵第一列是变异位点 ID其余列是样本 geno_dt - fread(genotype_matrix.csv, data.table FALSE) # 把第一列设为行名其余列保留为样本 protein_mat - as.matrix(protein_dt[, -1]) rownames(protein_mat) - protein_dt[[1]] geno_mat - as.matrix(geno_dt[, -1]) rownames(geno_mat) - geno_dt[[1]] # 检查维度行是特征列是样本 dim(protein_mat) dim(geno_mat)逻辑说明fread的data.table FALSE让返回对象是普通 data.frame方便后续转矩阵。把第一列抽出来做行名是 pQTL 数据的常见布局如果你的数据是样本在行、蛋白在列需要先转置。参数上fread默认会自动识别分隔符和表头但如果你的文件是分号分隔或带注释行要显式指定sep和skip。这一步做完先看dim输出的列数是否一致不一致就说明样本集合不同进入下一步取交集。2.3 取样本交集并重排别让静默丢样本坑了你# 取两个矩阵的样本交集 common_samples - intersect(colnames(protein_mat), colnames(geno_mat)) # 显式报告丢了哪些 cat(蛋白矩阵样本数:, ncol(protein_mat), \n) cat(基因型矩阵样本数:, ncol(geno_mat), \n) cat(交集样本数:, length(common_samples), \n) # 按交集重排保证两个矩阵列顺序完全一致 protein_mat - protein_mat[, common_samples, drop FALSE] geno_mat - geno_mat[, common_samples, drop FALSE] # 再次确认列名完全一致 stopifnot(identical(colnames(protein_mat), colnames(geno_mat)))逻辑说明intersect返回的是两个向量共有的元素但顺序按第一个向量来。重排后必须用identical做一次强校验因为all.equal在某些属性差异下会返回 TRUE 但实际列顺序不同。drop FALSE是防止当交集只剩一个样本时矩阵被降维成向量。这一步的坑在于如果你的样本 ID 有前后空格或大小写差异intersect会认为它们不同所以读入后最好统一trimws并确认大小写规则。2.4 缺失值与协变量跑关联前必须处理的几件事蛋白丰度矩阵里常有缺失值常见做法是过滤掉缺失比例过高的蛋白比如超过 20% 样本缺失就剔除剩下的用 KNN 或中位数填补。基因型矩阵则通常需要过滤次等位基因频率MAF过低的位点常见阈值是 MAF 0.05 或 0.01具体看样本量。协变量方面年龄、性别、批次效应是标配如果做的是 cis-pQTL还需要把基因型的主成分PC作为协变量控制群体分层。pQTLtools 本身不强制你做这些但它的关联函数通常接受一个协变量矩阵参数所以你得先把协变量整理成样本对齐的矩阵。# 过滤缺失比例过高的蛋白 keep_protein - rowMeans(is.na(protein_mat)) 0.2 protein_mat - protein_mat[keep_protein, ] # 简单中位数填补生产环境建议用更稳健的方法 for (i in seq_len(nrow(protein_mat))) { na_idx - is.na(protein_mat[i, ]) if (any(na_idx)) { protein_mat[i, na_idx] - median(protein_mat[i, ], na.rm TRUE) } } # 协变量矩阵行是样本列是协变量 covar_df - data.frame( age sample(20:80, length(common_samples), replace TRUE), sex sample(c(0, 1), length(common_samples), replace TRUE), batch sample(1:3, length(common_samples), replace TRUE) ) rownames(covar_df) - common_samples逻辑说明rowMeans(is.na(...))算的是每个蛋白的缺失比例比逐行循环快。中位数填补只是演示实际项目里如果缺失比例高更推荐用impute.knn或missForest。协变量矩阵的行名必须和样本 ID 完全一致否则后续模型矩阵会对不上。参数上MAF 过滤阈值没有绝对标准样本量小的时候阈值要放宽否则会丢掉大量位点。3. 关联检验与结果整理从单个位点到全蛋白组扫描3.1 关联模型的选择线性模型够不够用pQTL 关联检验的本质是对每个蛋白和每个位点检验位点基因型对蛋白丰度的效应。最常见的是线性模型把蛋白丰度作为因变量基因型作为自变量协变量作为控制项。如果蛋白丰度做了 log 转换且近似正态线性模型是合理的。但有些蛋白丰度分布偏态严重或者存在大量零值这时候可能需要用秩变换或广义线性模型。pQTLtools 的关联函数通常默认走线性模型但会提供参数让你指定是否对蛋白做标准化或秩变换。选型理由很简单cis-pQTL 的效应量通常较大线性模型在样本量几百以上时对轻微偏离正态比较稳健trans-pQTL 效应量小对模型假设更敏感这时候秩变换或稳健回归会更保险。3.2 用 pQTLtools 跑一轮 cis 关联下面这段演示如何对每个蛋白只在其 cis 窗口内比如转录起始位点上下游 1Mb的位点做关联。假设你已经有了蛋白的基因组位置信息和位点的位置信息。library(pQTLtools) # 假设 protein_info 包含 protein_id, chr, start, end # variant_info 包含 variant_id, chr, pos # 这里用模拟数据构造 protein_info - data.frame( protein_id rownames(protein_mat), chr sample(1:22, nrow(protein_mat), replace TRUE), start sample(1:1000000, nrow(protein_mat), replace TRUE) ) protein_info$end - protein_info$start 2000 variant_info - data.frame( variant_id rownames(geno_mat), chr sample(1:22, nrow(geno_mat), replace TRUE), pos sample(1:1000000, nrow(geno_mat), replace TRUE) ) # 跑 cis 关联cis_window 设为 1Mb cis_res - run_pqtl_association( protein_mat protein_mat, geno_mat geno_mat, protein_info protein_info, variant_info variant_info, covar covar_df, cis_window 1e6, mode cis ) # 查看结果结构 head(cis_res)逻辑说明run_pqtl_association是假设的函数名实际使用时需要对照 pQTLtools 的文档确认函数签名。核心参数cis_window控制 cis 窗口大小单位是碱基对mode cis表示只跑 cis如果设成trans则跑全基因组反式扫描。协变量矩阵covar的行必须和protein_mat、geno_mat的列顺序一致。结果通常包含蛋白 ID、位点 ID、效应量 beta、标准误 se、p 值等列。这一步的坑在于如果你的蛋白位置信息和位点位置信息用的基因组版本不一致比如一个 hg19 一个 hg38cis 窗口会完全错位跑出来的结果全是垃圾。3.3 多重检验校正与显著性阈值全蛋白组扫描意味着你要做几十万到几百万次检验多重检验校正是必须的。常见做法是 Bonferroni 校正阈值设为 0.05 除以检验次数。但 pQTL 领域更常用的是基于 LD 的独立检验数来估算有效检验次数或者直接用 FDR 控制。pQTLtools 的结果整理函数通常会提供一列校正后的 p 值你需要在结果里显式检查这一列是否存在。# 假设 cis_res 有 p_value 列 cis_res$bonferroni - p.adjust(cis_res$p_value, method bonferroni) cis_res$fdr - p.adjust(cis_res$p_value, method fdr) # 按显著性筛选 sig_res - subset(cis_res, fdr 0.05) # 统计每个蛋白的显著位点数 table(sig_res$protein_id)逻辑说明p.adjust的method参数支持bonferroni、fdr、holm等。Bonferroni 更严格FDR 更宽松。pQTL 里如果做的是 cis 分析检验次数相对少每个蛋白只跑窗口内的位点Bonferroni 也可接受如果是 trans 全基因组扫描FDR 更常用。筛选后按蛋白统计显著位点数可以快速看出哪些蛋白有强 cis 信号。3.4 结果整理成可读表格跑完关联后结果往往是一张长表你需要把它整理成每个蛋白一行、包含最显著位点信息的宽表方便后续做共定位或孟德尔随机化。# 按蛋白分组取 p 值最小的位点 library(dplyr) top_variant - cis_res %% group_by(protein_id) %% slice_min(order_by p_value, n 1) %% ungroup() # 整理成宽表 result_wide - top_variant %% select(protein_id, variant_id, beta, se, p_value, fdr) head(result_wide)逻辑说明slice_min按 p 值取每个蛋白最显著的位点。group_by加ungroup是 dplyr 的标准写法。整理后的宽表可以直接导出成 CSV也可以作为共定位分析的输入。这一步的坑在于如果两个位点 p 值完全相同比如都是 0slice_min会返回多行需要加with_ties FALSE或者先做去重。4. 可视化与共定位把结果讲成一张能看的图4.1 曼哈顿图与 QQ 图快速判断结果是否合理曼哈顿图是 pQTL 结果的标准展示方式横轴是基因组位置纵轴是 -log10(p)。QQ 图则用来检查 p 值分布是否偏离零假设。pQTLtools 通常会提供绘图函数但如果你要自己画qqman包是常见选择。library(qqman) # 构造曼哈顿图输入需要 SNP, CHR, BP, P 四列 manhattan_input - data.frame( SNP cis_res$variant_id, CHR variant_info$chr[match(cis_res$variant_id, variant_info$variant_id)], BP variant_info$pos[match(cis_res$variant_id, variant_info$variant_id)], P cis_res$p_value ) # 画曼哈顿图 manhattan(manhattan_input, main pQTL Manhattan, suggestiveline -log10(1e-5)) # 画 QQ 图 qq(manhattan_input$P, main pQTL QQ)逻辑说明match用来把位点的染色体和位置信息映射到结果表里。suggestiveline是建议阈值线通常设 1e-5。QQ 图如果点明显偏离对角线说明 p 值分布有问题可能是模型假设不满足或存在群体分层。这一步的坑在于如果你的位点 ID 在结果表和位置表里格式不一致比如一个带chr前缀一个不带match会返回 NA曼哈顿图会报错或画出一堆空点。4.2 共定位分析判断两个性状是否共享因果位点共定位colocalization是 pQTL 分析里的关键一步用来判断蛋白丰度和某个疾病性状是否由同一个因果位点驱动。常见方法是coloc包它需要两个性状在同一个区域的 summary statistics。pQTLtools 可能提供了封装函数但底层逻辑还是贝叶斯因子计算。library(coloc) # 假设你有蛋白的 summary stats 和疾病的 summary stats # 需要 beta, varbeta, N, type 等字段 protein_data - list( beta cis_res$beta, varbeta cis_res$se^2, N length(common_samples), type quant, snp cis_res$variant_id ) disease_data - list( beta rnorm(nrow(cis_res)), varbeta runif(nrow(cis_res), 0.01, 0.1), N 10000, type cc, s 0.5, snp cis_res$variant_id ) # 跑共定位 coloc_res - coloc.abf(protein_data, disease_data) print(coloc_res$summary)逻辑说明coloc.abf需要两个数据集在同一个位点集合上且 SNP 顺序一致。type参数指定性状类型quant是定量性状cc是病例对照。s是病例对照研究里病例的比例。结果里的PP.H4表示两个性状共享一个因果位点的后验概率通常大于 0.8 认为有强共定位证据。这一步的坑在于如果你的蛋白 summary stats 和疾病 summary stats 用的不是同一个参考面板的 LD 信息共定位结果会不可靠。4.3 把结果导出成可复现的报告分析做完后把关键结果和参数一起导出方便复现和交接。# 导出关联结果 write.csv(cis_res, pqtl_cis_results.csv, row.names FALSE) # 导出共定位结果 capture.output(coloc_res$summary, file coloc_summary.txt) # 保存会话信息 writeLines(capture.output(sessionInfo()), session_info.txt)逻辑说明sessionInfo()记录了 R 版本和所有包的版本是复现分析的关键。capture.output把控制台输出重定向到文件。这一步的坑在于如果你用了随机数比如模拟数据需要同时保存随机种子否则下次跑结果不一样。5. 避坑与排查那些让我重跑过整条流程的问题5.1 样本 ID 对不上但代码不报错现象关联结果里样本量显示 500但你明明有 800 个样本。原因两个矩阵的样本 ID 有前后空格或大小写差异intersect只取到了部分交集而后续代码没有显式检查交集大小。解决读入后立刻对样本 ID 做trimws和统一大小写取交集后打印交集大小并用stopifnot校验两个矩阵列名完全一致。5.2 基因组版本不一致导致 cis 窗口错位现象cis 关联跑出来几乎没有显著位点但已知这个蛋白有强 cis-pQTL。原因蛋白位置信息用的是 hg19位点位置信息用的是 hg38同一个基因的坐标差了几十 Mbcis 窗口完全没覆盖到真正的因果位点。解决在跑关联前确认所有位置信息用的是同一个基因组版本必要时用liftOver做转换。5.3 协变量矩阵行名与样本顺序不一致现象模型跑出来系数方向完全反了或者 p 值分布异常。原因协变量矩阵的行名虽然和样本 ID 一致但行顺序和蛋白矩阵的列顺序不同R 在按行名匹配时如果没显式重排可能会按位置匹配。解决在构造模型矩阵前用covar_df - covar_df[common_samples, , drop FALSE]显式重排并用identical(rownames(covar_df), colnames(protein_mat))校验。5.4 缺失值填补引入虚假关联现象某个蛋白的关联 p 值极小但该蛋白在大部分样本里都是缺失的。原因缺失比例过高的蛋白被中位数填补后填补值集中在中间导致与基因型的关联被虚假放大。解决在填补前先过滤缺失比例过高的蛋白常见阈值是 20% 或 30%并在结果里标注每个蛋白的缺失比例。5.5 多重检验校正用错检验次数现象Bonferroni 校正后没有任何显著位点但 FDR 校正后有一堆。原因Bonferroni 用的检验次数是全蛋白组扫描的总检验数而实际上 cis 分析里每个蛋白只跑了窗口内的位点检验次数应该按蛋白分别算或者用有效检验数。解决cis 和 trans 分开做校正cis 可以按蛋白分别做 Bonferronitrans 用全基因组 FDR。6. 进阶技巧用模拟数据验证流程的稳健性跑真实数据之前我习惯先用模拟数据把整条流程走一遍确认每个环节的输出符合预期。模拟数据的好处是你知道真实信号在哪里如果流程跑出来信号丢失或者出现假阳性就能快速定位是哪个环节出了问题。下面这段模拟一个蛋白和一个 cis 位点的关联效应量设为已知值然后看流程能否恢复出来。set.seed(42) # 模拟 500 个样本 n - 500 # 模拟一个 cis 位点的基因型MAF 约 0.3 geno - rbinom(n, 2, 0.3) # 模拟协变量 covar - data.frame(age rnorm(n), sex rbinom(n, 1, 0.5)) # 模拟蛋白丰度基因型效应为 0.5加上协变量效应和噪声 protein - 0.5 * geno 0.2 * covar$age 0.1 * covar$sex rnorm(n) # 跑线性模型 fit - lm(protein ~ geno age sex, data cbind(covar, protein protein, geno geno)) summary(fit)$coefficients[geno, ]逻辑说明set.seed保证结果可复现。rbinom生成 0/1/2 的基因型。lm里把协变量和基因型都放进去看geno的系数估计是否接近 0.5。如果流程正确系数应该在 0.5 附近p 值应该很小。这一步的坑在于模拟数据的效应量不要设得太小否则在 n500 时可能检测不到你会误以为流程有问题。另一个进阶用法是把 pQTLtools 的关联函数包装成批量跑多个蛋白的循环并用parallel包做并行化。但并行化之前先用小规模数据确认单次运行没问题否则并行跑出来的错误日志会很难排查。我一般会先跑 10 个蛋白确认结果合理后再扩展到全蛋白组。从那以后我每次跑 pQTL 流程都会先用模拟数据走一遍确认样本对齐、模型拟合、多重校正和结果导出这四个环节的输出都符合预期再上真实数据。这个习惯帮我省掉了至少两次整条流程重跑的时间。希望帮到你。本文还有配套的精品资源点击获取