企业官网建设流程全解析

简介：pQTLtools 是一套面向蛋白质定量性状基因座研究的 R 工具包，适合具备一定 R 语言与生物信息学基础、希望开展 pQTL 数据整理、注释、可视化和结果复现的研究人员。压缩包共 125 个文件，约 5.92MB，涵盖 html 教程与函数文档、png 图示、rda 数据对象、r 脚本、rd 帮助页、rmd 可复现报告以及 sh、yml 等配置和辅助文件，文档、示例数据与源码分层清晰，便于按需检索和对照学习。内容围绕 Caprion、Olink、SomaLogic 及 INTERVAL 等蛋白面板数据集展开，并附带 hg19ToHg38 转换链等资源，可帮助读者理解 pQTL 工具的数据组织方式、复现示例分析流程，并迁移到自己的蛋白质组学与遗传关联研究场景中。已有 759 人学习下载，适合作为 pQTL 工具入门、方法对照、脚本参考和排错思路积累的实用资料，对刚接触该工具包的用户也有较好的引导价值。

1. pQTLtools 是什么：一个 R 包把蛋白表和基因型接起来

手上有一批蛋白丰度数据，还有一套基因型，要把每个蛋白的 cis-pQTL 扫出来，多数人的第一反应是拼脚本：读表、取对数、按坐标开 ±1Mb 窗口、跑线性模型、算 FDR、画曼哈顿图。拼到第三个项目就会崩——列名对不上、协变量顺序错位、阈值说法年年变、换个人接手就复现不出来。pQTLtools 是一个 R 包，把这些重复动作收进统一接口：喂进蛋白丰度表、基因型三件套和协变量表，吐出一张可筛选的 pQTL 结果表、每个蛋白的独立信号，以及能直接进下游孟德尔随机化和共定位的格式。适合做蛋白组遗传学、需要批量跑几十到上千个蛋白的分析岗，也适合刚接手别人脚本、想找个能对照的实现把流程捋直的人。

2. 装包与三张表：列名对齐、ID 规范化和输入体检

装包这一步看着最没技术含量，实际上后面 80% 的报错都埋在这里：库路径混了、依赖版本对不上、样本 ID 尾部差一个字符。我的习惯是把它拆成三段——装、摆、检，每一步都留一个能打印出来的中间结果，出问题能立刻定位到是哪一段。

2.1 从源码包装起：依赖、编译和三个必查项

这类分析包一般不会在主干镜像里保持最新，源码包反而更可靠。拿到压缩包先别急着装，先解出 DESCRIPTION 看一眼，把依赖和 R 版本下限摸清楚，比装完再回头查缺什么快得多。

# 解包并查看包的元信息：依赖、R 版本下限、编译需求 tar -xzf pQTLtools_*.tar.gz sed -n '1,40p' pQTLtools/DESCRIPTION # 装到独立库目录，不污染系统库，也方便回滚到旧版本 mkdir -p ~/R/libs/pqtl R CMD INSTALL -l ~/R/libs/pqtl pQTLtools

-l指定安装目标库，这是我最推荐的一条参数：分析类包经常需要同时保留两三个版本做结果比对，全塞进系统库迟早互相覆盖。源码目录也别放在共享盘或同步盘上，编译过程会生成大量中间文件，网络盘上的文件锁能把安装卡死在 90%。

装完进 R，第一件事不是跑流程，是把导出函数列一遍。包内命名习惯（前缀、动词顺序）在不同构建里可能调整，硬套教程里的函数名是无效劳动。

# 把独立库放到搜索路径最前面，确保加载的是刚装的那份 .libPaths(c("~/R/libs/pqtl", .libPaths())) library(pQTLtools) # 1) 看导出了哪些函数，对一遍自己数据能用上几个 head(sort(ls("package:pQTLtools")), 50) # 2) 看关键函数的默认参数，尤其是窗口、MAF、缺失率这三个 args(pq_check) args(pq_cis)

三个必查项：R 版本是否达到 DESCRIPTION 里Depends的下限（低于下限往往不是报错，而是算出 quietly 错误的结果，更难查）；依赖包是否装齐（library()报 "there is no package called" 是最友好的一种错）；位数是否一致（64 位 R 配 32 位库，编译能过，加载会崩）。这三点确认完，再往下走。

2.2 蛋白表、基因型、协变量怎么摆

三张表的格式错误是新手最常翻车的地方，而且报错信息往往指向扫描函数内部，让人以为是算法出了问题。先把每张表的必需列列清楚。

蛋白表我一般保持「行是样本、列是蛋白」，因为下游取子集、算缺失率都是按列走。批次列如果存在蛋白表里，记得单独抽出来放进协变量表，不要让它跟着蛋白一起被当成表型。基因型这边，最常见的问题是 fam 文件的样本 ID 和你蛋白表里的 ID 长得像但不相等——差一个前缀、大小写不同、或者一个用了横线一个用了下划线。

library(data.table) prot <- fread("data/protein_abundance.csv") # 第一列必须是样本 ID geno <- fread("data/geno.fam") # bed/bim/fam 的 fam 表 cov <- fread("data/covariates.csv") # 显式取三张表的交集，别指望函数内部帮你对齐 sids <- Reduce(intersect, list(prot[[1]], geno[[1]], cov[[1]])) length(sids) # 这个数字如果明显小于单表行数，先停下查 ID setdiff(prot[[1]], geno[[1]])[1:10] # 打印前 10 个对不上的，肉眼比对格式

prot[[1]]取的是第一列，原因是 data.table 的[[按位置取列，不受列名影响，比写prot$SampleID更抗改名。Reduce(intersect, ...)一次拿到三表公共样本，比两两 intersect 少写一行也不容易漏。最后那行setdiff是关键：报错说不清哪里不一致的时候，把差集打出来看前十个，八成能一眼看出是前缀还是分隔符的问题。取完交集后，三张表都要按sids重排一遍，顺序不一致在某些实现里不会报错，但会把 A 的基因型配到 B 的蛋白上——这是最贵的一种错。

2.3 输入体检：把问题挡在扫描之前

扫描一次可能要几十分钟到几小时，等跑完才发现样本对不上、某个蛋白缺失率 60%，时间就白扔了。我一般强制先跑一遍输入检查，把每个蛋白的缺失率、方差、离群样本、基因型编码范围都打印出来。

chk <- pq_check( protein = prot[match(sids, prot[[1]]), ], genotype = "data/geno", # bed/bim/fam 的公共前缀 covariate = cov[match(sids, cov[[1]]), ], id_col = 1, # 三张表的样本 ID 都在第 1 列 log_transform = TRUE, # 强度型数据取 log2 missing_cutoff = 0.2, # 单蛋白缺失率上限 maf_cutoff = 0.05, sd_cutoff = 0.01 # 方差过小的蛋白直接剔除 ) chk$summary # 每个蛋白一行：缺失率、均值、标准差、是否保留 chk$dropped # 被剔掉的蛋白及原因，一定要看一眼

log_transform要看数据来源：原始强度值（抗体法给出的荧光计数、适配体法的 RFU）右偏严重，必须取对数；如果拿到的已经是平台归一化后的相对值，再取一次对数会把信号压平，这是很多人「跑出来一个显著位点都没有」的真实原因。missing_cutoff = 0.2是常用起点，缺失率超过 20% 的蛋白插补后基本不可信。sd_cutoff用来挡掉那些几乎恒定的蛋白，它们进线性模型只会产生数值警告。chk$dropped这张表不要跳过，看一眼被剔的是不是集中在某个批次或某个板上，这往往是前处理阶段出了问题，而不是蛋白本身质量差。

离群样本我的做法是先做 PCA，看前两个主成分有没有明显离群点，再决定是剔除还是加进协变量。直接把离群样本删掉又不在方法里写清楚，后面复现会很难看。协变量这边，性别编码统一成 0/1，批次转成因子，遗传主成分取前 5 到 10 个即可，取太多会把 cis 信号一起吸走。

3. cis 扫描：±1Mb 窗口、MAF 门槛和阈值口径

真正进入扫描环节，决定结果质量的其实只有三件事：窗口开多大、哪些位点允许进模型、显著性怎么定义。这三件事在文献里的口径并不统一，抄错一个，结论方向就可能变。

3.1 为什么 cis 窗口开 ±1Mb 而不是全基因组

cis-pQTL 的通行定义是位点落在蛋白编码基因转录起始位点上下 1Mb 以内。这个数字不是随便定的：绝大多数已报道的蛋白-遗传关联都落在这个范围内，窗口开到 1Mb 能覆盖绝大部分调控区，同时把多重检验的负担压到可接受量级。有人为了「不漏」直接把窗口开到 5Mb，结果是一个蛋白要检几十万个位点，Bonferroni 阈值被抬得极高，真正的信号反而检不出来。

trans 扫描是另一回事。全基因组 × 上千个蛋白，如果不做 LD 剪枝（把高度连锁的位点压成一个代表位点），位点数会到几百万，内存和时长都不可控。常见做法是先用 r² > 0.2 或 0.5 剪一遍，再按染色体分块跑，每跑完一条染色体就落盘，别指望一次跑完整个基因组。我的经验是：cis 优先，trans 作为补充，先把 cis 做干净，很多蛋白的遗传结构在 cis 层面就解释得差不多了，trans 的分析预算应该留给 cis 无信号的蛋白。

窗口边界还有一个细节：坐标体系要统一。蛋白编码基因的坐标如果用了一种版本，基因型文件的位点坐标用了另一种，窗口就整体偏移，扫描能跑完但结果是错的，这种错误几乎没有报错提示，属于典型的黑匣子。

3.2 一次可复现的 cis 扫描脚本

res_cis <- pq_cis( input = chk, window = 1e6, # 上下游各 1Mb maf = 0.05, # 只保留次等位基因频率 >= 5% 的位点 missing_rate = 0.05, # 单一位点缺失率上限 covariate_cols = c("age", "sex", "batch", "pc1", "pc2", "pc3", "pc4", "pc5"), model = "linear", # 加性模型 threads = 8, chunk_size = 200, # 每批处理多少个蛋白 out = "out/cis" # 结果落盘前缀 ) head(res_cis$table) res_cis$lead # 每个蛋白的 lead 位点，下游条件分析的输入

window按上面的理由取 1e6，如果研究对象是某个特定基因家族，可以按基因长度单独传一个向量，不必一刀切。maf = 0.05是平衡检出力与稳定性的折中，低于这个频率的位点在小样本里的效应估计极不稳定，容易制造假阳性。missing_rate卡的是位点层面而不是样本层面，两者别混。covariate_cols里写的名字必须和协变量表的列名完全一致，写错一个不会报错，只是那一列静默失效——这块我吃过亏，跑完发现前五个主成分根本没进模型。chunk_size是为了控内存：蛋白数上千的时候，一次读进内存的矩阵会撑爆 16G 的机器，按 200 一批跑，中间结果落盘，断了还能续。threads不要开到机器核数上限，留两核给系统，否则 IO 会成为瓶颈，加线程反而不提速。

res_cis$table的典型列包括蛋白名、位点 ID、染色体位置、效应等位基因、效应值、标准误、P 值、样本量。拿到之后先看两件事：每个蛋白的位点数是否和预期量级一致（对不上通常是窗口或 MAF 过滤把位点全滤掉了），以及效应等位基因那一列有没有空值。

3.3 阈值与 FDR 的取法

pQTL 的显著性阈值是个容易被含糊带过的地方。直接把 GWAS 那套 5e-8 搬过来不合适：蛋白组的样本量通常只有几百到几千，检出力达不到那个量级，硬套的结果是一个显著位点都没有。我一般同时报两套口径：蛋白内 Bonferroni（0.05 除以该蛋白实际参与检验的位点数）用来标记强信号，BH 校正后的 FDR < 0.05 用来做发现集的主口径。两套都写上，读者自己能判断。

还有一点常被忽略：不同蛋白的检验数量差异很大（基因长短不同、窗口内位点数不同），所以 FDR 要按蛋白分层校正，全表一把 BH 会让长基因的蛋白吃亏。结果表里留一列「该蛋白的检验位点数」，后面写方法时能直接引用，也方便别人复算阈值。复制集验证时，方向必须和发现集一致，只报 P 值不报方向是很常见的偷懒做法。

4. 下游三件事：条件分析、共定位与三张图

扫描出 lead 位点只是起点。一个蛋白区域里经常有多个独立信号，lead 位点可能只是其中最强的那个；而「这个位点和那个 GWAS 位点是不是同一个因果变异」这个问题，靠 P 值大小是回答不了的。下游这三步，决定了结果能不能写进文章。

4.1 条件分析：把连锁带出来的假信号剔掉

条件分析的逻辑很直白：把 lead 位点作为协变量加进模型再扫一遍，如果还有显著信号，说明该区域存在第二个独立信号；重复这个过程直到没有新信号出现。它对 LD 的依赖非常强，翻车也多半翻在这里。

cond <- pq_conditional( input = chk, lead = res_cis$lead, # 上一步的 lead 位点表 r2_cutoff = 0.1, # 与已选位点 LD 超过此值的直接跳过 max_iter = 10, # 迭代上限，防止弱信号来回震荡 p_cutoff = 1e-5, # 每一轮的入选门槛 ld_reference = "data/ld_panel" # LD 参考面板，必须和基因型同人群 ) cond$independent # 每个蛋白的独立信号列表 cond$iteration_log # 每轮加了哪个位点、P 值怎么变

r2_cutoff = 0.1卡的是「新位点和已入选位点之间的连锁程度」，设太松会把同一个信号拆成两条记录，设太紧会漏掉真正独立的弱信号。max_iter是后悔药：区域内信号一多，不加迭代上限可能一直在两个位点之间来回选。最关键的参数是ld_reference——参考面板的人群必须和你的基因型数据一致，拿一个人群的面板去算另一个人群的 LD，算出来的 r² 完全不能用，而这个过程不会报错。cond$iteration_log值得逐行看，如果某个蛋白迭代到第 8 轮还在出新位点，八成是 LD 阈值或者人群不匹配的问题，而不是真的有那么多个独立信号。

4.2 共定位：输入对齐比方法选择更重要

共定位回答的是「pQTL 信号和 GWAS 信号是否共享同一个因果变异」。方法和先验参数都有成熟套路，真正耗时间的是把两份汇总统计对齐：效应等位基因要一致，不一致的要么翻转效应值符号，要么直接剔除；位点 ID 的写法要统一（有的用 chr:pos，有的用 rsID）；等位基因频率的参考人群要说明。

coloc <- pq_coloc( region = "chr1:1000000-3000000", pqtl = res_cis$table, gwas = "data/gwas_summary.tsv", p12 = 1e-5, # 与 GWAS 共享因果变异的先验概率 p1 = 1e-4, # pQTL 侧存在信号的先验 p2 = 1e-4, # GWAS 侧存在信号的先验 align = "flip", # 等位基因不一致时按频率翻转 maf_cutoff = 0.05 ) coloc$summary # 每个区域的后验概率 PP.H0 ~ PP.H4

p12 = 1e-5是常用默认值，表示「一个区域里同时存在两个信号且共享因果变异」的先验相当低；如果你的候选区域是先验筛选过的（比如只在 GWAS 显著位点附近做），可以适当调高，但要在方法里写明。align = "flip"只处理能通过等位基因频率判断方向的情况，遇到回文位点（A/T、C/G）谨慎处理，这类位点在高频时无法确定链方向，我一般直接剔除而不是硬猜。coloc$summary里重点看 PP.H4（共享因果变异）和 PP.H3（两个独立信号）的比值，只看 H4 的绝对值容易被先验带偏。

4.3 三张图：曼哈顿、森林图与箱线图

出图这块包里有现成函数，参数不多但都值得调。曼哈顿图用来展示区域内的信号分布和 lead 位点位置；森林图展示独立信号的效应量和置信区间；箱线图按基因型分组展示蛋白丰度，是最直观的一张，也最容易被审稿人要求补。

# 区域曼哈顿图：只画候选区域，标注 lead 位点 pq_plot_manhattan(cond, region = "chr1:1000000-3000000", label_top = 3, point_size = 0.8, ylim = c(0, 12)) # 独立信号森林图：效应值 + 95% 置信区间 pq_plot_forest(cond$independent, xlab = "Effect on protein level (SD)") # 按基因型分组的箱线图，展示最強信号的剂量效应 pq_plot_box(chk, snp = res_cis$lead$snp[1], protein = res_cis$lead$protein[1])

label_top = 3控制标注几个位点，标太多图会糊；ylim手动设上限是为了跨图可比，几张图放在一起时坐标轴不一致会误导读者。森林图的 x 轴我统一用标准差单位，比原始强度值好解释。箱线图要注意基因型分组的样本量，最少一组不足 20 个样本时别画，容易得出反向结论。

5. 排查手册：五类报错的现场还原与处置

下面这五条是我在不同队列上真实撞过的，形态各不相同，但共同点是报错信息都不指向真正的病灶。

5.1 读数据阶段的三类

现象一：chk$summary里样本数 n 是 0 或者个位数。原因几乎永远是 ID 不匹配：fam 文件的样本 ID 带了0_或F_前缀，或者一边用横线一边用下划线，或者大小写不一致。解决方法是写一个规范化函数，把三张表的 ID 统一成「去前缀 + 转大写 + 分隔符统一为下划线」再取交集，并在脚本开头固定下来，不要每次临时改。

现象二：所有蛋白跑完一个显著位点都没有，P 值分布还偏保守。常见原因是表型方向搞反了——蛋白表行列颠倒，函数按列当样本处理，居然能算出结果。另一种是数据已经做过标准化，又开了一次log_transform，把方差压没了。处置：先dim()打印矩阵形状，再画一张 P 值 QQ 图，明显低于对角线就是模型设定有问题，不是「这个性状没有遗传基础」。

现象三：结果表里 P 值大面积 NA。检查基因型编码：有些流程输出的是 1/2/3，而模型期望 0/1/2，多出来的那个水平会被当成缺失。另外看一下 MAF 过滤后每个蛋白还剩多少位点，剩几十个的话，说明位点 ID 或坐标体系和注释对不上，窗口里根本没取到东西。

5.2 扫描与结果阶段的两类

现象四：条件分析迭代到max_iter还在出新位点。原因通常是 LD 参考面板和基因型数据人群不一致，或者r2_cutoff设得过松（比如 0.5），弱连锁的信号被反复当成新信号。处置：换同人群面板，把 r2 收到 0.1 以内，再对照cond$iteration_log看每一轮加进去的位点是不是都挤在同一个坐标区间。

现象五：跑到一半内存爆掉，或者一夜过去进度条没动。全基因组 trans 扫描没做 LD 剪枝是最常见的原因；其次是chunk_size设得太大，一批几百个蛋白的矩阵同时驻留内存。处置：按染色体分块、先剪枝、把chunk_size降到 50 到 100，每块跑完立刻落盘。还有一点容易忽略：threads开太大时每个线程各持一份数据副本，内存占用是乘上去的，不是加。

提示：这五条里有四条都不会抛出明确报错，属于「能跑完但结果不对」的类别。所以每次换数据集，先跑一遍输入体检并核对三个数字，比事后 debug 便宜得多。

6. 批量流水线：把参数抽出来，用三个数字验结果

单个队列跑通之后，真正的工作量在批量化——十几个队列、几百个蛋白、参数还要按队列微调。我现在的做法是把所有会变的参数抽成一张配置表，一行一个队列，脚本只读表不写死值。

# config.csv: cohort, protein_file, geno_prefix, cov_file, # window, maf, missing_cutoff, covariate_cols cfg <- fread("config/cohorts.csv") for (i in seq_len(nrow(cfg))) { p <- as.list(cfg[i]) chk <- pq_check(protein = fread(p$protein_file), genotype = p$geno_prefix, covariate = fread(p$cov_file), log_transform = TRUE, missing_cutoff = as.numeric(p$missing_cutoff), maf_cutoff = as.numeric(p$maf)) res <- pq_cis(input = chk, window = as.numeric(p$window), maf = as.numeric(p$maf), out = file.path("out", p$cohort)) fwrite(res$table, file.path("out", p$cohort, "cis_summary.tsv"), sep = "\t") }

参数抽出来后，配置表本身要能校验：window写成 1000 而不是 1e6，扫描能跑，但结果毫无意义；maf写成 5 而不是 0.05，位点会被全滤掉。所以我在循环开头加一道stopifnot(window >= 1e5, maf > 0 && maf < 0.5)，把明显的单位错误挡在跑批之前。

验证我只看三个数字，且每次都记在日志里：三表交集样本数占单表最大样本数的比例（低于七成要查 ID 和质控）、每个蛋白保留位点数的中位数（明显偏小说明窗口或 MAF 有问题）、lead 位点所在蛋白的 P 值中位数（应该在 1e-6 到 1e-15 量级，跑出来 1e-2 附近基本是模型没设对）。这三个数字对得上，结果才值得往下游送。从那以后我每接一个新队列，都强制先跑pq_check并把这三点写进日志，再开扫描——省下来的不是一次 debug 的时间，是整个结论不被推翻的底气。希望帮到你。

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