给生物信息学新手的实战指南：如何复现UK Biobank蛋白质组学中的pQTL分析（附代码思路）

生物信息学实战UK Biobank蛋白质组学pQTL分析全流程解析在生物信息学领域蛋白质数量性状位点(pQTL)分析已成为连接基因组学与蛋白质组学的重要桥梁。对于刚进入这一领域的研究者而言如何将复杂的统计遗传学方法转化为可操作的代码流程往往是最具挑战性的环节。本文将基于UK Biobank蛋白质组学数据(UKB-PPP)拆解pQTL分析的全流程从数据获取到结果可视化提供可直接复用的代码框架和避坑指南。1. 环境准备与数据获取1.1 软件依赖安装pQTL分析需要一系列生物信息学工具和统计包的支持。建议使用conda创建独立环境conda create -n pqtl python3.8 conda activate pqtl conda install -c bioconda plink regenie susie r-ggplot2 r-qqman pip install pandas numpy scipy statsmodels1.2 UKB-PPP数据下载与预处理UK Biobank蛋白质组学数据可通过官方申请获取。获批后关键文件包括蛋白表达数据ukb_proteomics.bed/bim/fam基因型数据ukb_cal_chr{1-22}_v2.bed协变量文件covariates.txt使用PLINK进行初步质控# 样本QC plink --bfile ukb_cal_chr1_v2 --mind 0.02 --maf 0.01 --hwe 1e-6 --make-bed --out ukb_sqc # 蛋白数据归一化 Rscript -e library(preprocessCore) protein - read.table(ukb_proteomics.bed) norm_protein - normalize.quantiles(as.matrix(protein)) write.table(norm_protein, ukb_proteomics_norm.txt) 注意UKB数据使用需遵守相关协议禁止非授权分发2. 关联分析核心流程2.1 REGENIE两步法实现REGENIE是处理大规模生物库数据的首选工具其两步分析能有效控制群体分层# 第一步构建预测模型 regenie \ --step 1 \ --bed ukb_sqc \ --phenoFile ukb_proteomics_norm.txt \ --covarFile covariates.txt \ --bsize 1000 \ --loocv \ --out step1_pred # 第二步关联测试 regenie \ --step 2 \ --bed ukb_sqc \ --phenoFile ukb_proteomics_norm.txt \ --covarFile covariates.txt \ --pred step1_pred_pred.list \ --bsize 400 \ --out pqtl_results关键参数说明参数作用推荐值--bsize块大小1000(步骤1)/400(步骤2)--loocv留一交叉验证必选--threads并行线程数根据服务器配置调整2.2 多重检验校正pQTL分析需考虑数千种蛋白质的多次测试问题。Bonferroni校正过于保守推荐使用FDRimport statsmodels.stats.multitest as multi pvals [...] # 从regenie结果中读取 _, fdr_pvals, _, _ multi.multipletests(pvals, methodfdr_bh) significant fdr_pvals 0.053. 精细定位与结果解析3.1 SuSiE精细定位SuSiE可识别独立信号并计算后验包含概率(PIP)library(susieR) # 准备LD矩阵和Z分数 z_scores - read.table(pqtl_zscore.txt) ld_matrix - read.table(ld_matrix.txt) # 运行SuSiE fit - susie_rss(z_scores, ld_matrix, L10) summary(fit)$cs # 输出可信集3.2 可视化分析曼哈顿图和QQ图是评估pQTL结果质量的核心工具library(qqman) # 曼哈顿图 manhattan(pqtl_results, chrCHR, bpBP, pP, snpSNP, suggestiveline-log10(1e-5), genomewideline-log10(5e-8)) # QQ图 qq(pqtl_results$P, mainQ-Q plot of pQTL p-values)常见问题排查QQ图膨胀λ值1.05提示群体分层未完全校正曼哈顿图异常某些染色体过度显著可能指示批次效应4. 高级分析与生物学解释4.1 共定位分析使用coloc评估pQTL与eQTL的共享信号library(coloc) # 准备pQTL和eQTL数据 pqtl_dataset - list(pvaluespqtl_p, Npqtl_n) eqtl_dataset - list(pvalueseqtl_p, Neqtl_n) # 运行共定位 result - coloc.abf(pqtl_dataset, eqtl_dataset) if(result$summary[6] 0.8) print(强共定位证据)4.2 通路富集分析对显著pQTL基因进行通路富集from gseapy import enrichr significant_genes [...] # 提取显著关联基因 enr_results enrichr(gene_listsignificant_genes, gene_sets[KEGG_2021_Human]) enr_results.results.head()典型富集通路包括补体激活通路炎症反应通路脂蛋白代谢通路5. 实战中的挑战与解决方案5.1 群体分层控制即使使用REGENIE的LOCO策略在混合群体中仍需额外注意# 计算前20个主成分 plink --bfile ukb_sqc --pca 20 --out pca_results # 在协变量文件中加入PCs paste covariates.txt pca_results.eigenvec | awk {print $1,$2,$4,$5,$6,$7,$8,$9,$10,$11,$12,$13,$14,$15,$16,$17,$18,$19,$20,$21,$22,$23} covariates_pcs.txt5.2 连锁不平衡处理不同人群需要不同的LD参考面板人群推荐LD面板获取方式欧洲1000G EURPLINK官网东亚1000G EASPLINK官网非洲HapMap3 AFRHapMap项目5.3 计算资源优化大规模分析时的实用技巧分染色体分析并行处理22条染色体for chr in {1..22}; do plink --chr $chr --bfile ukb_sqc --make-bed --out ukb_chr$chr done wait内存管理使用--memory参数限制REGENIE内存使用临时文件设置--tmp指向大容量临时存储6. 扩展应用与前沿方法6.1 跨祖先pQTL分析非欧洲人群分析需特别注意# 使用METAL进行跨祖先meta分析 metal_cmd - MARKER SNP ALLELE A1 A2 EFFECT BETA STDERR SE PVALUE P WEIGHT N PROCESS european.txt PROCESS asian.txt ANALYZE writeLines(metal_cmd, metal_script.txt) system(metal metal_script.txt)6.2 孟德尔随机化应用以PCSK9为例的代码框架library(TwoSampleMR) # 选择工具变量 exposure_dat - extract_instruments(prot-a-1645) # PCSK9蛋白ID outcome_dat - extract_outcome_data(exposure_dat$SNP, ieu-a-7) # CAD数据 # 运行MR harmonised_dat - harmonise_data(exposure_dat, outcome_dat) res - mr(harmonised_dat) mr_forest_plot(res)关键检查点F统计量10避免弱工具偏差异质性检验Cochrans Q检验P0.05多效性检验MR-Egger截距P0.05在实际项目中我们常遇到REGENIE结果文件过大的问题。一个实用的解决方案是使用Tabix建立索引bgzip pqtl_results.txt tabix -s 1 -b 2 -e 2 pqtl_results.txt.gz这样即可快速查询特定区域的关联结果而无需加载整个文件。对于包含2,923种蛋白质的全基因组分析这种方法可将查询时间从分钟级降至秒级。