尧图网站建设 
简介此份资源是面向蛋白质组学与遗传学研究的实用工具包聚焦于蛋白质定量性状基因座简称pQTL的识别与注释适合从事生物信息学数据挖掘的科研人员。资源包内共一百二十五个文件主要包含R语言源码、数据对象、说明文档与配置文件整体大小约为六兆字节涵盖多个蛋白质组学平台的注释信息。目前已有七百五十八人学习过该资源可作为判断其实用性的参考。通过这份工具包使用者能够直接获取整理好的基因型注释、蛋白丰度数据、示例工作流及可视化脚本减少前期数据清洗与格式转换的步骤快速开展pQTL关联分析。资源内附有多平台数据面板和文档教程尤其适合需要标准化流程的初学者快速上手。1. pQTLtools蛋白质数量性状位点分析为什么绕不开这套工具链pQTLtools这个名称对应的是蛋白质数量性状位点protein quantitative trait locus, pQTL分析中一整条工具链而不是某个单一脚本。做过血浆蛋白组学与基因型数据关联分析的人都有体会拿到上千例样本的蛋白定量矩阵和数百万级变异位点目标是找“哪个遗传变异影响了哪种蛋白的丰度”但实际时间几乎全耗在格式转换、样本匹配、循环跑回归、处理异常值这些杂活上。pQTLtools要解决的正是这件事把样本质控、蛋白数据预处理、cis/trans关联扫描、共定位和孟德尔随机化串成标准化流程。它适合三类人刚拿到蛋白组数据准备做遗传关联分析的读pQTL文献想复现结果的以及想用蛋白做工具变量开展因果推断的。2. pQTL和eQTL差在哪先搞懂数据逻辑再谈工具选型2.1 三种蛋白定量平台给pQTL带来的数据性格蛋白定量数据的来源差异决定了pQTL分析的第一步不能直接套用eQTL的标准化流程。目前大规模pQTL研究里最常见的是Olink、SomaScan和质谱三种平台它们产出的数据结构完全不同。Olink平台基于邻位延伸技术产物是NPX值做过log变换后近似正态分布直接跑线性回归通常问题不大。但它的检测下限LOD是硬边界低于LOD的值会被截断分布左尾被切除。若一个蛋白有超过30%的样本都压在检测下限附近这个蛋白在遗传关联分析里几乎没有信息量跑出来的显著位点也大概率是假阳性。SomaScan则是适配体技术覆盖的蛋白数量多、动态范围宽但适配体交叉结合会让部分蛋白的信号不纯某些“蛋白”本质上是多个蛋白的混合信号。质谱平台无偏覆盖全、蛋白数最多但缺失率动辄20%以上且缺失不是随机的经常和批次、样本质量绑定。这些差异落入参数选择就是Olink数据不需要再做什么log变换但要按LOD以下比例过滤蛋白SomaScan要额外做批次校正质谱数据要先看缺失模式再决定填充还是舍弃。pQTLtools这类工具链的价值恰恰是把这些平台特定的QC逻辑收拢成几个明确的参数而不是让每个用户从零写一套。很多人在这上面翻车是因为拿着一份Olink的数据套用了某个转录组流程做了重复的标准化把真实生物信号削掉了一截。2.2 cis-pQTL与trans-pQTL效应量不同分析策略完全不同pQTL分析分两个层次cis-pQTL看的是变异位点附近通常是转录起始位点上下游1Mb范围内对蛋白丰度的影响trans-pQTL看的是染色体上其他位置甚至其他染色体上的远距离调控。这两个层次的效应量和所需样本量差异悬殊。cis-pQTL的效应通常很大。实际项目里一个cis-pQTL解释10%到30%的蛋白丰度方差并不罕见这和eQTL的典型效应量分布完全不是一个量级。也正因为如此cis-pQTL在几百例样本里就可能有不错检出率。trans-pQTL则是另一回事效应小、需要检验的变异-蛋白对数量巨大动辄上亿次检验样本量低于2000基本不用指望稳定复现。这个差异直接影响分析策略样本量只有几百的时候老老实实只做cis样本量上了几千才值得考虑trans。窗口大小的设定也会改变检验次数cis窗口默认取1Mb有人习惯放到2Mb但每扩大一档需要校正的检验数量就成倍增加。我一般建议窗口参数保持1Mb如果某条通路需要看远端调控单独跑trans而不是硬把cis窗口拉大。2.3 为什么不能拿eQTL工具直接跑pQTL转录组和蛋白组的遗传调控并非一一对应。mRNA丰度受转录调控蛋白丰度还叠加了翻译效率、降解速率和分泌过程的影响所以同一个变异对mRNA和蛋白的效应经常不一致。直接用eQTL pipeline跑pQTL最容易出问题的是协变量设计。eQTL分析里的协变量通常是年龄、性别、遗传主成分加隐性批次变量而蛋白定量数据的批次效应显著更强不同板plate、不同上机时间、不同样本制备批次都能制造系统性偏差。如果这些批次变量不进入模型pQTL的曼哈顿图会出现整条染色体抬高的假象。反过来协变量加太多又会把真实信号稀释掉尤其是在样本量不够大的时候。这也是pQTLtools这类工具链存在的理由它不是只做关联扫描而是把数据准备、平台QC、协变量构建、关联检验、下游因果推断作为一个整体来设计。单看每一步用PLINK或Matrix eQTL也能实现但把中间文件的格式定义清楚把共定位和MR的输入要求提前考虑进去才是整套方案真正的价值。3. 用pQTLtools跑通最小cis-pQTL流程三条命令与五个必调参数3.1 输入文件四个文件一个都不能少一个标准的最小化cis-pQTL分析需要基因型、蛋白定量矩阵、协变量文件、蛋白位置注释四个输入。缺了任何一个流程都会在某个不起眼的环节停下来。文件常见格式关键列说明基因型PLINK bed/bim/fam变异ID、chr、pos、A1/A2建议先做MAF、HWE、缺失过滤蛋白定量矩阵tab分隔文本首列为样本ID后续每列一个蛋白Olink输出NPX值SomaScan输出RFU协变量文件tab分隔文本样本ID、年龄、性别、批次、PC1-10样本ID必须和基因型完全一致蛋白位置文件BED格式chr、start、end、protein_id基因组版本必须与基因型一致蛋白位置文件是最容易被忽略的。cis窗口要判断“某个变异是否落在蛋白编码基因附近”就必须知道每个蛋白的编码基因在染色体上的位置。如果这个文件的基因组版本和基因型VCF不一致后面所有cis映射都白做。常见做法是从Ensembl或UCSC下载对应版本的注释提取蛋白编码转录本的外显子边界合并成BED。3.2 第一步prepare所有格式问题在这里终结拿到原始数据后先用PLINK对基因型做一次常规过滤再调用prepare模块把蛋白矩阵转成统一格式并做样本匹配。这一步的核心目标是把所有ID格式问题暴露在分析之前而不是让它们在关联扫描中途跳出来。# 基因型QCMAF0.01、HWE p1e-6、位点缺失率5% plink --bfile geno_raw --maf 0.01 --hwe 1e-6 --geno 0.05 \ --make-bed --out geno_qc # 调用prepare模块匹配样本、生成BED格式蛋白矩阵、审计ID pqtools prepare --geno geno_qc \ --protein protein_matrix.txt \ --protein-bed protein_pos.bed \ --covar covariates.txt \ --genome-version GRCh38 \ --out workdir--geno 0.05表示过滤掉缺失率超过5%的位点这是eQTL分析常用的阈值pQTL分析中如果样本量只有几百可以放宽到0.1避免位点数量被砍得太狠。--genome-version建议写死不要相信默认值GRCh37和GRCh38混用是这类分析最常见的隐性问题。prepare执行完毕后会生成一个样本ID匹配报告重点看两件事一是蛋白矩阵里有多少样本在基因型文件中找不到二是重复样本ID有多少。如果匹配率低于95%先回去检查样本命名来源不要在流程后续阶段处理。3.3 第二步cis扫描窗口、阈值与协变量一次设对prepare做完之后cis扫描本身就是一个命令的事。但参数要在运行前想清楚因为重跑整个扫描的成本不小。# cis-pQTL关联扫描线性加性模型窗口1Mb pqtools cis --bfile workdir/geno \ --protein-bed workdir/protein.bed \ --covar workdir/covar.txt \ --window 1000000 \ --maf 0.01 \ --model linear \ --out workdir/cis_result.txt--window决定cis窗口的半径1000000表示以蛋白编码基因边界为中心上下游各扩1Mb。--maf 0.01是低频变异的过滤线样本量小于500时建议收紧到0.05否则个别样本携带的稀有变异可能制造出极具迷惑性的信号。--model linear指定的是加性线性模型这也是pQTL分析的事实标准基因型编码为0/1/2对应效应等位基因的个数。输出文件的核心字段包括variant_id、protein_id、chr、pos、ref、alt、beta、se、tstat、pvalue、maf。beta的解读是每增加一个效应等位基因蛋白丰度变化多少单位。Olink数据的NPX值本身在log2尺度beta就能直接解读为相对变化翻倍或减半这一点在后续写论文时非常方便。注意cis窗口的判断依据是变异位点到蛋白编码基因边界的最短距离不是到转录起始位点的距离。基因很大时用TSS算窗口会漏掉落在基因下游远端内含子里的真实调控位点。3.4 结果初筛与QQ图检查扫描跑完后不要急着看显著位点先做一步整体质控看p值分布是否正常。最有效的方式是画QQ图看尾部有没有异常抬升。import gzip import numpy as np import pandas as pd import scipy.stats as stats import matplotlib matplotlib.use(Agg) import matplotlib.pyplot as plt # 读取cis关联结果只保留有效p值 cis pd.read_csv(workdir/cis_result.txt, sep\t) pvals cis[pvalue].dropna() pvals pvals[(pvals 0) (pvals 1)] # 计算期望p值分布排序后均匀分布的分位数 obs -np.log10(np.sort(pvals)) exp -np.log10(np.linspace(1 / len(pvals), 1, len(pvals))) # 画出QQ图并计算基因组膨胀因子lambda slope np.polyfit(exp[:int(len(exp) * 0.9)], obs[:int(len(exp) * 0.9)], 1)[0] lam np.median(obs) / np.median(exp) plt.figure(figsize(6, 6)) plt.scatter(exp, obs, s4, alpha0.5) plt.plot([0, exp.max()], [0, exp.max()], colorred, linestyle--) plt.xlabel(Expected -log10(p)) plt.ylabel(Observed -log10(p)) plt.title(fcis-pQTL QQ plot (lambda{lam:.3f})) plt.savefig(workdir/qq_plot.png, dpi150)lambda值在1.0到1.2之间属于正常范围超过1.3说明存在系统性偏差最常见的来源是批次效应没有完全校正。此外看QQ图尾部如果出现一条与对角线平行但整体抬高的曲线通常是协变量缺失如果尾部剧烈上翘而中部贴合很好说明有少量真实信号这是理想形态。这一步能拦住绝大多数返工场景属于跑完pQTL之后的第一道质量门。4. 从关联到因果共定位与孟德尔随机化怎么接进pQTLtools4.1 共定位输入的坑beta/se比p值好使拿到cis-pQTL显著位点后下一个自然的问题是这个区域里影响蛋白丰度的变异和影响疾病风险的变异是不是同一个共定位分析colocalization就是回答这个问题的。实际做共定位时输入格式比模型本身更容易出错。import pandas as pd # cis关联结果 cis pd.read_csv(workdir/cis_result.txt, sep\t) # MAF文件从PLINK的freq输出提取 freq pd.read_csv(workdir/freq.frq, sep\t, delim_whitespaceTrue) # coloc要求的核心列snp、beta、se、MAF、N coloc_input cis[[variant_id, beta, se]].merge( freq[[SNP, MAF]], left_onvariant_id, right_onSNP, howleft ) coloc_input.rename(columns{variant_id: snp}, inplaceTrue) coloc_input[N] 3000 # 填入实际样本量不是越大越好 coloc_input coloc_input.dropna(subset[MAF, se]) coloc_input.to_csv(workdir/coloc_input.txt, sep\t, indexFalse)代码里的N是实际参与关联分析的样本量填错了会直接影响后验概率的计算精度。有人图省事填一个很大的数结果PP.H4被系统性抬高做完共定位以为找到了共享因果变异换个数据集一验证就露馅。另一个隐蔽问题是MAF的来源如果是从千人基因组参考panel里取的频率要和自己基因型数据的实际MAF做一遍相关性检查相关系数低于0.99就说明参考panel和实际人群不匹配这会扭曲共定位的后验概率。4.2 工具变量筛选pQTL到MR的必经筛选把蛋白丰度作为暴露、疾病作为结局做孟德尔随机化时pQTL显著位点就是工具变量的候选池。但显著不等于合格工具变量要满足三条和暴露强相关和混杂因素独立只通过暴露影响结局。第一条做起来相对机械后两条要靠下游敏感性分析。# LD clumping保留每个区域内的代表变异r20.1 plink --bfile workdir/geno \ --clump workdir/cis_result.txt \ --clump-p1 5e-8 \ --clump-kb 1000 \ --clump-r2 0.1 \ --out workdir/clumpedclump之后再用F统计量过滤弱工具变量。F的计算方式是(beta/se)^2F大于10才认为工具变量足够强。import pandas as pd iv pd.read_csv(workdir/clumped.clumped, sep\t, delim_whitespaceTrue) cis pd.read_csv(workdir/cis_result.txt, sep\t) iv iv.merge(cis[[variant_id, beta, se]], left_onSNP, right_onvariant_id) iv[F] (iv[beta] / iv[se]) ** 2 iv iv[iv[F] 10] iv.to_csv(workdir/iv_final.txt, sep\t, indexFalse)--clump-p1 5e-8是全基因组显著性的经典阈值cis-pQTL分析的校正阈值如果更宽比如1e-5那做MR时的工具变量筛选就要重新设定P阈值否则一堆弱工具变量混进来会让MR的因果估计偏向零。--clump-r2 0.1是比较严格的标准如果工具变量数量不足可以放宽到0.3但必须同步做MR-Egger截距检验确认没有水平多效性的痕迹。4.3 可视化曼哈顿图和区域关联图的轻量做法关联结果的可视化不必一上来就上重型工具先把显著位点标到一张曼哈顿图上足够用于日常质检。常见做法是用R的qqman包manhattan()一行就能出图如果要画某个区域的关联图LocusZoom风格的图可以用R的locuszoomR或python的gwaslab替代。pQTLtools如果集成了可视化模块核心输出一般就两类全基因组的曼哈顿图和一个区域内所有变异的关联信号叠加图。区域图的横坐标是基因组位置点按LD着色能把共定位信号直观地展示出来。这里我习惯将关联最强位点标为紫色LD r2大于0.8的标为红色0.4到0.8标为橙色小于0.4标为灰色会议报告和论文配图都能直接使用。5. pQTLtools实战避坑五个让结果翻车的常见问题5.1 蛋白矩阵样本ID与基因型样本ID不一致现象prepare步骤报错提示样本匹配率低于60%强行忽略后继续跑最终结果里大量蛋白的关联p值全是NaN。原因蛋白定量矩阵里的样本ID用的是“样本编号日期后缀”基因型数据用的是纯编号少数样本的ID里还有不可见空格。这类问题靠肉眼看不出来程序只做精确匹配匹配不上就是匹配不上。解决在prepare之前先写一个小脚本对两份样本ID做归一化处理统一去掉空格和多余分隔符再抽20个样本人工核对一遍。准备阶段多花半小时能省掉后面整轮重跑。5.2 基因组版本不一致导致cis映射全部错位现象cis-pQTL结果里某个蛋白的显著信号出现在完全不同的染色体上或者窗口内明明包含编码基因输出里却没有该基因的任何变异。原因基因型VCF是GRCh38坐标蛋白位置注释文件却来自GRCh37两个坐标系差几十到几百个碱基大片段倒位区域错得更离谱。解决拿一个已知的cis-pQTL阳性对照位点去验证。比如文献里明确报道过某蛋白的cis信号位置看看自己的结果能否映射到同一区域。这个检查应该在正式扫描前做而不是等结果出来再排查。5.3 协变量矩阵过度拟合导致信号被系统性稀释现象QQ图lambda值显著小于1p值分布整体右移之前能看到的一些信号全部消失。原因协变量里同时放入了批次、plate、上机时间、年龄、性别、PC前20自由度消耗过多。蛋白组学数据的批次效应确实要校正但协变量之间存在共线性时模型会把真实的遗传效应一并吸收掉。解决先检查协变量之间的方差膨胀因子VIFVIF大于10的变量只保留其中一个。通常批次、plate、上机时间是高度相关的保留plate就够。遗传主成分PC前10足够不要加到20。5.4 trans-pQTL显著位点经不起置换检验现象trans-pQTL扫描用5e-8阈值筛出上百个显著位点看起来成果丰富换个数据集一个都复现不了。原因trans-pQTL的变异-蛋白检验对数量比cis大几个数量级5e-8这个经典全基因组阈值是针对单性状GWAS设计的用在trans上仍然不够严格。检验次数上亿次时5e-8的期望假阳性数量已经相当可观。解决trans结果必须做两遍过滤。第一遍用5e-8初筛第二遍对每个候选位点做置换检验或者用FDR控制在1%以内。更稳妥的做法是只在cis结果基础上做下游共定位和MRtrans位点一律打上“待验证”标签。5.5 蛋白ID映射到基因时丢失大量条目现象蛋白数量是几百个但prepare结束时报告只有六成蛋白成功映射到基因组位置四分之一的蛋白静默失踪。原因平台输出的蛋白ID是UniProt访问号注释文件用的是Ensembl基因ID中间映射表不完整还有些分泌型蛋白的编码基因虽然存在但注释数据库里没有收录最新的转录本。解决不要用平台自带的映射表将就。从UniProt的ID mapping API拉最新映射跑之前把映射成功的比例打印出来看一眼。如果丢的是目标通路里的关键蛋白宁可停下来补全映射也不要带着缺失往下跑。这个检查做到位后面多少能避开一半以上的无效分析。6. 质量校验的进阶技巧方向一致性验证与保留LD结构的置换检验6.1 结果方向一致性上线前必做的快速体检pQTL分析做完最该做的不是急着看显著性列表而是评估结果的稳健性。我习惯把样本按批次一分为二各自独立跑一遍cis-pQTL然后只取两批分析中都达到名义显著水平的位点比较beta的方向一致率。如果一致率低于80%说明结果受批次影响严重回去检查协变量如果在85%到90%之间可以接受但需要说明超过90%说明关联信号比较稳固。这个检查的妙处在于成本极低只需两次并行扫描却能在论文投稿前挡住一批经不起复现的结果。6.2 置换检验要保住LD循环位移比随机打乱更可靠cis-pQTL的置换检验常犯一个错误随机打乱基因型或随机打乱蛋白表达值这样做会彻底破坏变异位点之间的连锁不平衡结构置换出来的零分布失真。正确做法是对蛋白表达向量做循环位移circular shift也就是把整列样本的蛋白值整体平移一个随机偏移量再重新跑关联。这样保留了样本间的相关结构只破坏了基因型和表型之间的配对置换p值才可信。用来算经验p值时1000次置换是下限5000次更稳虽然慢一点但算出来的显著性有说服力。做pQTL分析这几年我在置换检验上栽过跟头。有一次拿到一份cis-pQTL结果置换后几乎全军覆没排查到最后发现是随机打乱了基因型矩阵把LD结构全毁了。后来固定用循环位移方案再也没出过这种问题。现在拿到任何pQTL分析结果我第一件事就是拆半验方向一致性和循环位移置换这俩没过关后面的共定位和MR都白搭。希望这些细节能帮你在pQTL分析的路上少走几步弯路。本文还有配套的精品资源点击获取 pQTLtools实战：蛋白质数量性状位点分析全流程与避坑指南
简介此份资源是面向蛋白质组学与遗传学研究的实用工具包聚焦于蛋白质定量性状基因座简称pQTL的识别与注释适合从事生物信息学数据挖掘的科研人员。资源包内共一百二十五个文件主要包含R语言源码、数据对象、说明文档与配置文件整体大小约为六兆字节涵盖多个蛋白质组学平台的注释信息。目前已有七百五十八人学习过该资源可作为判断其实用性的参考。通过这份工具包使用者能够直接获取整理好的基因型注释、蛋白丰度数据、示例工作流及可视化脚本减少前期数据清洗与格式转换的步骤快速开展pQTL关联分析。资源内附有多平台数据面板和文档教程尤其适合需要标准化流程的初学者快速上手。1. pQTLtools蛋白质数量性状位点分析为什么绕不开这套工具链pQTLtools这个名称对应的是蛋白质数量性状位点protein quantitative trait locus, pQTL分析中一整条工具链而不是某个单一脚本。做过血浆蛋白组学与基因型数据关联分析的人都有体会拿到上千例样本的蛋白定量矩阵和数百万级变异位点目标是找“哪个遗传变异影响了哪种蛋白的丰度”但实际时间几乎全耗在格式转换、样本匹配、循环跑回归、处理异常值这些杂活上。pQTLtools要解决的正是这件事把样本质控、蛋白数据预处理、cis/trans关联扫描、共定位和孟德尔随机化串成标准化流程。它适合三类人刚拿到蛋白组数据准备做遗传关联分析的读pQTL文献想复现结果的以及想用蛋白做工具变量开展因果推断的。2. pQTL和eQTL差在哪先搞懂数据逻辑再谈工具选型2.1 三种蛋白定量平台给pQTL带来的数据性格蛋白定量数据的来源差异决定了pQTL分析的第一步不能直接套用eQTL的标准化流程。目前大规模pQTL研究里最常见的是Olink、SomaScan和质谱三种平台它们产出的数据结构完全不同。Olink平台基于邻位延伸技术产物是NPX值做过log变换后近似正态分布直接跑线性回归通常问题不大。但它的检测下限LOD是硬边界低于LOD的值会被截断分布左尾被切除。若一个蛋白有超过30%的样本都压在检测下限附近这个蛋白在遗传关联分析里几乎没有信息量跑出来的显著位点也大概率是假阳性。SomaScan则是适配体技术覆盖的蛋白数量多、动态范围宽但适配体交叉结合会让部分蛋白的信号不纯某些“蛋白”本质上是多个蛋白的混合信号。质谱平台无偏覆盖全、蛋白数最多但缺失率动辄20%以上且缺失不是随机的经常和批次、样本质量绑定。这些差异落入参数选择就是Olink数据不需要再做什么log变换但要按LOD以下比例过滤蛋白SomaScan要额外做批次校正质谱数据要先看缺失模式再决定填充还是舍弃。pQTLtools这类工具链的价值恰恰是把这些平台特定的QC逻辑收拢成几个明确的参数而不是让每个用户从零写一套。很多人在这上面翻车是因为拿着一份Olink的数据套用了某个转录组流程做了重复的标准化把真实生物信号削掉了一截。2.2 cis-pQTL与trans-pQTL效应量不同分析策略完全不同pQTL分析分两个层次cis-pQTL看的是变异位点附近通常是转录起始位点上下游1Mb范围内对蛋白丰度的影响trans-pQTL看的是染色体上其他位置甚至其他染色体上的远距离调控。这两个层次的效应量和所需样本量差异悬殊。cis-pQTL的效应通常很大。实际项目里一个cis-pQTL解释10%到30%的蛋白丰度方差并不罕见这和eQTL的典型效应量分布完全不是一个量级。也正因为如此cis-pQTL在几百例样本里就可能有不错检出率。trans-pQTL则是另一回事效应小、需要检验的变异-蛋白对数量巨大动辄上亿次检验样本量低于2000基本不用指望稳定复现。这个差异直接影响分析策略样本量只有几百的时候老老实实只做cis样本量上了几千才值得考虑trans。窗口大小的设定也会改变检验次数cis窗口默认取1Mb有人习惯放到2Mb但每扩大一档需要校正的检验数量就成倍增加。我一般建议窗口参数保持1Mb如果某条通路需要看远端调控单独跑trans而不是硬把cis窗口拉大。2.3 为什么不能拿eQTL工具直接跑pQTL转录组和蛋白组的遗传调控并非一一对应。mRNA丰度受转录调控蛋白丰度还叠加了翻译效率、降解速率和分泌过程的影响所以同一个变异对mRNA和蛋白的效应经常不一致。直接用eQTL pipeline跑pQTL最容易出问题的是协变量设计。eQTL分析里的协变量通常是年龄、性别、遗传主成分加隐性批次变量而蛋白定量数据的批次效应显著更强不同板plate、不同上机时间、不同样本制备批次都能制造系统性偏差。如果这些批次变量不进入模型pQTL的曼哈顿图会出现整条染色体抬高的假象。反过来协变量加太多又会把真实信号稀释掉尤其是在样本量不够大的时候。这也是pQTLtools这类工具链存在的理由它不是只做关联扫描而是把数据准备、平台QC、协变量构建、关联检验、下游因果推断作为一个整体来设计。单看每一步用PLINK或Matrix eQTL也能实现但把中间文件的格式定义清楚把共定位和MR的输入要求提前考虑进去才是整套方案真正的价值。3. 用pQTLtools跑通最小cis-pQTL流程三条命令与五个必调参数3.1 输入文件四个文件一个都不能少一个标准的最小化cis-pQTL分析需要基因型、蛋白定量矩阵、协变量文件、蛋白位置注释四个输入。缺了任何一个流程都会在某个不起眼的环节停下来。文件常见格式关键列说明基因型PLINK bed/bim/fam变异ID、chr、pos、A1/A2建议先做MAF、HWE、缺失过滤蛋白定量矩阵tab分隔文本首列为样本ID后续每列一个蛋白Olink输出NPX值SomaScan输出RFU协变量文件tab分隔文本样本ID、年龄、性别、批次、PC1-10样本ID必须和基因型完全一致蛋白位置文件BED格式chr、start、end、protein_id基因组版本必须与基因型一致蛋白位置文件是最容易被忽略的。cis窗口要判断“某个变异是否落在蛋白编码基因附近”就必须知道每个蛋白的编码基因在染色体上的位置。如果这个文件的基因组版本和基因型VCF不一致后面所有cis映射都白做。常见做法是从Ensembl或UCSC下载对应版本的注释提取蛋白编码转录本的外显子边界合并成BED。3.2 第一步prepare所有格式问题在这里终结拿到原始数据后先用PLINK对基因型做一次常规过滤再调用prepare模块把蛋白矩阵转成统一格式并做样本匹配。这一步的核心目标是把所有ID格式问题暴露在分析之前而不是让它们在关联扫描中途跳出来。# 基因型QCMAF0.01、HWE p1e-6、位点缺失率5% plink --bfile geno_raw --maf 0.01 --hwe 1e-6 --geno 0.05 \ --make-bed --out geno_qc # 调用prepare模块匹配样本、生成BED格式蛋白矩阵、审计ID pqtools prepare --geno geno_qc \ --protein protein_matrix.txt \ --protein-bed protein_pos.bed \ --covar covariates.txt \ --genome-version GRCh38 \ --out workdir--geno 0.05表示过滤掉缺失率超过5%的位点这是eQTL分析常用的阈值pQTL分析中如果样本量只有几百可以放宽到0.1避免位点数量被砍得太狠。--genome-version建议写死不要相信默认值GRCh37和GRCh38混用是这类分析最常见的隐性问题。prepare执行完毕后会生成一个样本ID匹配报告重点看两件事一是蛋白矩阵里有多少样本在基因型文件中找不到二是重复样本ID有多少。如果匹配率低于95%先回去检查样本命名来源不要在流程后续阶段处理。3.3 第二步cis扫描窗口、阈值与协变量一次设对prepare做完之后cis扫描本身就是一个命令的事。但参数要在运行前想清楚因为重跑整个扫描的成本不小。# cis-pQTL关联扫描线性加性模型窗口1Mb pqtools cis --bfile workdir/geno \ --protein-bed workdir/protein.bed \ --covar workdir/covar.txt \ --window 1000000 \ --maf 0.01 \ --model linear \ --out workdir/cis_result.txt--window决定cis窗口的半径1000000表示以蛋白编码基因边界为中心上下游各扩1Mb。--maf 0.01是低频变异的过滤线样本量小于500时建议收紧到0.05否则个别样本携带的稀有变异可能制造出极具迷惑性的信号。--model linear指定的是加性线性模型这也是pQTL分析的事实标准基因型编码为0/1/2对应效应等位基因的个数。输出文件的核心字段包括variant_id、protein_id、chr、pos、ref、alt、beta、se、tstat、pvalue、maf。beta的解读是每增加一个效应等位基因蛋白丰度变化多少单位。Olink数据的NPX值本身在log2尺度beta就能直接解读为相对变化翻倍或减半这一点在后续写论文时非常方便。注意cis窗口的判断依据是变异位点到蛋白编码基因边界的最短距离不是到转录起始位点的距离。基因很大时用TSS算窗口会漏掉落在基因下游远端内含子里的真实调控位点。3.4 结果初筛与QQ图检查扫描跑完后不要急着看显著位点先做一步整体质控看p值分布是否正常。最有效的方式是画QQ图看尾部有没有异常抬升。import gzip import numpy as np import pandas as pd import scipy.stats as stats import matplotlib matplotlib.use(Agg) import matplotlib.pyplot as plt # 读取cis关联结果只保留有效p值 cis pd.read_csv(workdir/cis_result.txt, sep\t) pvals cis[pvalue].dropna() pvals pvals[(pvals 0) (pvals 1)] # 计算期望p值分布排序后均匀分布的分位数 obs -np.log10(np.sort(pvals)) exp -np.log10(np.linspace(1 / len(pvals), 1, len(pvals))) # 画出QQ图并计算基因组膨胀因子lambda slope np.polyfit(exp[:int(len(exp) * 0.9)], obs[:int(len(exp) * 0.9)], 1)[0] lam np.median(obs) / np.median(exp) plt.figure(figsize(6, 6)) plt.scatter(exp, obs, s4, alpha0.5) plt.plot([0, exp.max()], [0, exp.max()], colorred, linestyle--) plt.xlabel(Expected -log10(p)) plt.ylabel(Observed -log10(p)) plt.title(fcis-pQTL QQ plot (lambda{lam:.3f})) plt.savefig(workdir/qq_plot.png, dpi150)lambda值在1.0到1.2之间属于正常范围超过1.3说明存在系统性偏差最常见的来源是批次效应没有完全校正。此外看QQ图尾部如果出现一条与对角线平行但整体抬高的曲线通常是协变量缺失如果尾部剧烈上翘而中部贴合很好说明有少量真实信号这是理想形态。这一步能拦住绝大多数返工场景属于跑完pQTL之后的第一道质量门。4. 从关联到因果共定位与孟德尔随机化怎么接进pQTLtools4.1 共定位输入的坑beta/se比p值好使拿到cis-pQTL显著位点后下一个自然的问题是这个区域里影响蛋白丰度的变异和影响疾病风险的变异是不是同一个共定位分析colocalization就是回答这个问题的。实际做共定位时输入格式比模型本身更容易出错。import pandas as pd # cis关联结果 cis pd.read_csv(workdir/cis_result.txt, sep\t) # MAF文件从PLINK的freq输出提取 freq pd.read_csv(workdir/freq.frq, sep\t, delim_whitespaceTrue) # coloc要求的核心列snp、beta、se、MAF、N coloc_input cis[[variant_id, beta, se]].merge( freq[[SNP, MAF]], left_onvariant_id, right_onSNP, howleft ) coloc_input.rename(columns{variant_id: snp}, inplaceTrue) coloc_input[N] 3000 # 填入实际样本量不是越大越好 coloc_input coloc_input.dropna(subset[MAF, se]) coloc_input.to_csv(workdir/coloc_input.txt, sep\t, indexFalse)代码里的N是实际参与关联分析的样本量填错了会直接影响后验概率的计算精度。有人图省事填一个很大的数结果PP.H4被系统性抬高做完共定位以为找到了共享因果变异换个数据集一验证就露馅。另一个隐蔽问题是MAF的来源如果是从千人基因组参考panel里取的频率要和自己基因型数据的实际MAF做一遍相关性检查相关系数低于0.99就说明参考panel和实际人群不匹配这会扭曲共定位的后验概率。4.2 工具变量筛选pQTL到MR的必经筛选把蛋白丰度作为暴露、疾病作为结局做孟德尔随机化时pQTL显著位点就是工具变量的候选池。但显著不等于合格工具变量要满足三条和暴露强相关和混杂因素独立只通过暴露影响结局。第一条做起来相对机械后两条要靠下游敏感性分析。# LD clumping保留每个区域内的代表变异r20.1 plink --bfile workdir/geno \ --clump workdir/cis_result.txt \ --clump-p1 5e-8 \ --clump-kb 1000 \ --clump-r2 0.1 \ --out workdir/clumpedclump之后再用F统计量过滤弱工具变量。F的计算方式是(beta/se)^2F大于10才认为工具变量足够强。import pandas as pd iv pd.read_csv(workdir/clumped.clumped, sep\t, delim_whitespaceTrue) cis pd.read_csv(workdir/cis_result.txt, sep\t) iv iv.merge(cis[[variant_id, beta, se]], left_onSNP, right_onvariant_id) iv[F] (iv[beta] / iv[se]) ** 2 iv iv[iv[F] 10] iv.to_csv(workdir/iv_final.txt, sep\t, indexFalse)--clump-p1 5e-8是全基因组显著性的经典阈值cis-pQTL分析的校正阈值如果更宽比如1e-5那做MR时的工具变量筛选就要重新设定P阈值否则一堆弱工具变量混进来会让MR的因果估计偏向零。--clump-r2 0.1是比较严格的标准如果工具变量数量不足可以放宽到0.3但必须同步做MR-Egger截距检验确认没有水平多效性的痕迹。4.3 可视化曼哈顿图和区域关联图的轻量做法关联结果的可视化不必一上来就上重型工具先把显著位点标到一张曼哈顿图上足够用于日常质检。常见做法是用R的qqman包manhattan()一行就能出图如果要画某个区域的关联图LocusZoom风格的图可以用R的locuszoomR或python的gwaslab替代。pQTLtools如果集成了可视化模块核心输出一般就两类全基因组的曼哈顿图和一个区域内所有变异的关联信号叠加图。区域图的横坐标是基因组位置点按LD着色能把共定位信号直观地展示出来。这里我习惯将关联最强位点标为紫色LD r2大于0.8的标为红色0.4到0.8标为橙色小于0.4标为灰色会议报告和论文配图都能直接使用。5. pQTLtools实战避坑五个让结果翻车的常见问题5.1 蛋白矩阵样本ID与基因型样本ID不一致现象prepare步骤报错提示样本匹配率低于60%强行忽略后继续跑最终结果里大量蛋白的关联p值全是NaN。原因蛋白定量矩阵里的样本ID用的是“样本编号日期后缀”基因型数据用的是纯编号少数样本的ID里还有不可见空格。这类问题靠肉眼看不出来程序只做精确匹配匹配不上就是匹配不上。解决在prepare之前先写一个小脚本对两份样本ID做归一化处理统一去掉空格和多余分隔符再抽20个样本人工核对一遍。准备阶段多花半小时能省掉后面整轮重跑。5.2 基因组版本不一致导致cis映射全部错位现象cis-pQTL结果里某个蛋白的显著信号出现在完全不同的染色体上或者窗口内明明包含编码基因输出里却没有该基因的任何变异。原因基因型VCF是GRCh38坐标蛋白位置注释文件却来自GRCh37两个坐标系差几十到几百个碱基大片段倒位区域错得更离谱。解决拿一个已知的cis-pQTL阳性对照位点去验证。比如文献里明确报道过某蛋白的cis信号位置看看自己的结果能否映射到同一区域。这个检查应该在正式扫描前做而不是等结果出来再排查。5.3 协变量矩阵过度拟合导致信号被系统性稀释现象QQ图lambda值显著小于1p值分布整体右移之前能看到的一些信号全部消失。原因协变量里同时放入了批次、plate、上机时间、年龄、性别、PC前20自由度消耗过多。蛋白组学数据的批次效应确实要校正但协变量之间存在共线性时模型会把真实的遗传效应一并吸收掉。解决先检查协变量之间的方差膨胀因子VIFVIF大于10的变量只保留其中一个。通常批次、plate、上机时间是高度相关的保留plate就够。遗传主成分PC前10足够不要加到20。5.4 trans-pQTL显著位点经不起置换检验现象trans-pQTL扫描用5e-8阈值筛出上百个显著位点看起来成果丰富换个数据集一个都复现不了。原因trans-pQTL的变异-蛋白检验对数量比cis大几个数量级5e-8这个经典全基因组阈值是针对单性状GWAS设计的用在trans上仍然不够严格。检验次数上亿次时5e-8的期望假阳性数量已经相当可观。解决trans结果必须做两遍过滤。第一遍用5e-8初筛第二遍对每个候选位点做置换检验或者用FDR控制在1%以内。更稳妥的做法是只在cis结果基础上做下游共定位和MRtrans位点一律打上“待验证”标签。5.5 蛋白ID映射到基因时丢失大量条目现象蛋白数量是几百个但prepare结束时报告只有六成蛋白成功映射到基因组位置四分之一的蛋白静默失踪。原因平台输出的蛋白ID是UniProt访问号注释文件用的是Ensembl基因ID中间映射表不完整还有些分泌型蛋白的编码基因虽然存在但注释数据库里没有收录最新的转录本。解决不要用平台自带的映射表将就。从UniProt的ID mapping API拉最新映射跑之前把映射成功的比例打印出来看一眼。如果丢的是目标通路里的关键蛋白宁可停下来补全映射也不要带着缺失往下跑。这个检查做到位后面多少能避开一半以上的无效分析。6. 质量校验的进阶技巧方向一致性验证与保留LD结构的置换检验6.1 结果方向一致性上线前必做的快速体检pQTL分析做完最该做的不是急着看显著性列表而是评估结果的稳健性。我习惯把样本按批次一分为二各自独立跑一遍cis-pQTL然后只取两批分析中都达到名义显著水平的位点比较beta的方向一致率。如果一致率低于80%说明结果受批次影响严重回去检查协变量如果在85%到90%之间可以接受但需要说明超过90%说明关联信号比较稳固。这个检查的妙处在于成本极低只需两次并行扫描却能在论文投稿前挡住一批经不起复现的结果。6.2 置换检验要保住LD循环位移比随机打乱更可靠cis-pQTL的置换检验常犯一个错误随机打乱基因型或随机打乱蛋白表达值这样做会彻底破坏变异位点之间的连锁不平衡结构置换出来的零分布失真。正确做法是对蛋白表达向量做循环位移circular shift也就是把整列样本的蛋白值整体平移一个随机偏移量再重新跑关联。这样保留了样本间的相关结构只破坏了基因型和表型之间的配对置换p值才可信。用来算经验p值时1000次置换是下限5000次更稳虽然慢一点但算出来的显著性有说服力。做pQTL分析这几年我在置换检验上栽过跟头。有一次拿到一份cis-pQTL结果置换后几乎全军覆没排查到最后发现是随机打乱了基因型矩阵把LD结构全毁了。后来固定用循环位移方案再也没出过这种问题。现在拿到任何pQTL分析结果我第一件事就是拆半验方向一致性和循环位移置换这俩没过关后面的共定位和MR都白搭。希望这些细节能帮你在pQTL分析的路上少走几步弯路。本文还有配套的精品资源点击获取