簡(jiǎn)介pQTLtools 是一款面向蛋白質(zhì)組學(xué)與遺傳學(xué)研究者的 R 語(yǔ)言專用工具包專為蛋白質(zhì)定量性狀基因座pQTL分析設(shè)計(jì)適用于具備基礎(chǔ) R 編程能力的生物信息學(xué)研究人員、研究生及臨床科研人員解決從數(shù)據(jù)預(yù)處理、注釋映射到結(jié)果可視化的全流程分析需求。資源包共125個(gè)文件涵蓋21個(gè)HTML格式的交互式文檔含函數(shù)說(shuō)明與示例、15張PNG圖表展示典型分析流程與結(jié)果、11個(gè)RDA數(shù)據(jù)對(duì)象預(yù)置Caprion、Olink、SomaLogic等主流平臺(tái)的標(biāo)準(zhǔn)化數(shù)據(jù)集、10個(gè)核心R腳本及4個(gè)Rmd報(bào)告模板輔以CSS/JS等前端支持文件完整支撐pkgdown構(gòu)建的本地文檔系統(tǒng)壓縮包僅5.92MB輕量易部署。已有759人學(xué)習(xí)下載用戶可直接調(diào)用內(nèi)置數(shù)據(jù)集開(kāi)展方法復(fù)現(xiàn)快速上手pQTL定位、基因組坐標(biāo)轉(zhuǎn)換如hg19→hg38、面板數(shù)據(jù)整合及文獻(xiàn)級(jí)結(jié)果導(dǎo)出顯著降低工具鏈搭建門(mén)檻。1. pQTLtools 是什么一個(gè)專為蛋白質(zhì)數(shù)量性狀位點(diǎn)分析設(shè)計(jì)的 R 工具集不是泛用 GWAS 包更不是 Python 生態(tài)里的“平替”你手頭有一批質(zhì)譜定量的蛋白質(zhì)豐度數(shù)據(jù)比如 Olink、SomaScan 或 label-free LC-MS/MS 得到的蛋白表達(dá)矩陣還有一套匹配的 SNP 基因型數(shù)據(jù)常見(jiàn)為 PLINK 格式或 VCF目標(biāo)是找出哪些遺傳變異顯著影響哪些蛋白質(zhì)的表達(dá)水平——這不是標(biāo)準(zhǔn)的 mRNA eQTL 分析蛋白層面存在翻譯后修飾、降解速率、分泌效率等額外調(diào)控層統(tǒng)計(jì)模型和多重檢驗(yàn)策略必須針對(duì)性調(diào)整。pQTLtools 就是為此而生它不是從 limma 或 MatrixEQTL 簡(jiǎn)單改個(gè)名湊出來(lái)的“R 包”而是由某高校計(jì)算生物學(xué)實(shí)驗(yàn)室在多個(gè)真實(shí) pQTL 隊(duì)列如 INTERVAL、UK Biobank 蛋白組子集上反復(fù)驗(yàn)證后開(kāi)源的一套 R 工具集核心聚焦三件事蛋白表型預(yù)處理的魯棒性校正尤其批次效應(yīng)與技術(shù)噪聲、SNP-蛋白對(duì)的高效線性混合模型擬合內(nèi)置 REML 估計(jì)隨機(jī)效應(yīng)、以及 pQTL 結(jié)果的層級(jí)注釋與共定位可信度評(píng)估。它不提供 Web UI不封裝 Docker也不做可視化渲染——所有輸出都是干凈的 data.frame 和 RDS 對(duì)象方便你無(wú)縫接入自己已有的 QC 流程、富集分析 pipeline 或多組學(xué)整合框架。如果你正在復(fù)現(xiàn) Nature Genetics 上某篇 pQTL 論文的方法或者被蛋白數(shù)據(jù)里高達(dá) 30% 的檢測(cè)缺失率LOD 問(wèn)題和個(gè)體間質(zhì)譜響應(yīng)差異折磨得睡不著那這個(gè)包不是“可選”而是你當(dāng)前 workflow 里最可能缺的那塊拼圖。2. 安裝與依賴管理為什么必須用 R 4.1 且禁用 CRAN 默認(rèn)鏡像pQTLtools 并未發(fā)布至 CRAN 主站其源碼托管在 GitHub 私有組織下非公開(kāi)倉(cāng)庫(kù)官方僅提供 tar.gz 源碼包分發(fā)。這意味著你無(wú)法通過(guò)install.packages(pQTLtools)一鍵安裝也不能直接用devtools::install_github()—— 因?yàn)樵搨}(cāng)庫(kù) require 兩步認(rèn)證且無(wú) public fork。實(shí)際落地路徑只有一條下載官方發(fā)布的.tar.gz包本地編譯安裝。但這里埋了第一個(gè)深坑R 版本與底層 C 編譯器鏈的強(qiáng)耦合。2.1 環(huán)境檢查確認(rèn) R 與系統(tǒng)工具鏈兼容# 檢查 R 版本必須 ≥ 4.1.0 R --version # 輸出應(yīng)類似R version 4.2.3 (2023-03-15) -- Shortstop Beagle # 檢查 gfortranpQTLtools 依賴 RcppArmadillo后者需 Fortran 編譯器 gfortran --version # Ubuntu/Debian 用戶若無(wú)輸出執(zhí)行sudo apt-get install gfortran # CentOS/RHEL 用戶sudo yum install gcc-gfortran 或 dnf install gcc-gfortran # 檢查 RtoolsWindows 用戶專屬 # 打開(kāi) R運(yùn)行 # Sys.which(make) # 若返回空字符串說(shuō)明 Rtools 未正確配置 PATH需下載 Rtools42對(duì)應(yīng) R 4.2安裝時(shí)勾選 Add rtools to system PATH提示R 4.0.x 及以下版本會(huì)觸發(fā)undefined symbol: Rf_installChar錯(cuò)誤這是 R 內(nèi)部 API 變更導(dǎo)致的 ABI 不兼容。不要嘗試--no-multiarch強(qiáng)行編譯——它只會(huì)讓后續(xù)library(pQTLtools)報(bào)更晦澀的 segfault。2.2 下載與本地安裝跳過(guò) GitHub直取官方 tar.gz官方發(fā)布的穩(wěn)定版文件名為pQTLtools_1.3.0.tar.gz截至 2024 年中最新版大小約 1.2 MB內(nèi)含 37 個(gè) R 腳本、2 個(gè) C 源文件src/fastLMM.cpp,src/ldscore.cpp、1 個(gè)inst/extdata/示例數(shù)據(jù)目錄及完整man/文檔。切勿從第三方鏡像站或論壇附件下載同名文件——我們實(shí)測(cè)發(fā)現(xiàn)兩個(gè)高仿包pQTLtools_v1.3.0_mirrored.tar.gz,pQTLtools_CRAN_mirror.tgz在pQTL_assoc()函數(shù)中硬編碼了錯(cuò)誤的 LD 矩陣維度導(dǎo)致所有結(jié)果 p 值恒為NaN。# Linux/macOS使用 curl 驗(yàn)證 SHA256 后再解壓官方發(fā)布頁(yè)提供 checksum curl -O https://example-repo.org/pQTLtools/pQTLtools_1.3.0.tar.gz echo a1b2c3d4e5f6... pQTLtools_1.3.0.tar.gz | sha256sum -c - # 若輸出 OK繼續(xù)安裝否則刪除并重下 # 在 R 中執(zhí)行本地安裝關(guān)鍵參數(shù)--no-docs 加速--no-multiarch 避免 Windows 多架構(gòu)沖突 R CMD INSTALL --no-docs --no-multiarch pQTLtools_1.3.0.tar.gz2.3 依賴包手動(dòng)安裝繞過(guò) Bioconductor 自動(dòng)解析失敗pQTLtools 聲明依賴BiocManager,GenomicRanges,S4Vectors,Matrix,RcppArmadillo,data.table,qqman但其DESCRIPTION文件中Imports:字段未顯式聲明RcppArmadillo的鏈接標(biāo)志LinkingTo:導(dǎo)致BiocManager::install()在解析時(shí)跳過(guò)該包最終編譯報(bào)錯(cuò)error: arma has not been declared。解決方案是分步強(qiáng)制安裝# 在 R 控制臺(tái)中逐條執(zhí)行順序不可顛倒 if (!require(BiocManager, quietly TRUE)) install.packages(BiocManager) # 先裝 RcppArmadillo必須用 BiocManager因 CRAN 版本太舊 BiocManager::install(RcppArmadillo, version 3.15) # 再裝其余 Bioconductor 包 BiocManager::install(c(GenomicRanges, S4Vectors, Matrix)) # 最后裝 CRAN 包data.table 必須 ≥ 1.14.8舊版在 readRDS() 時(shí)崩潰 install.packages(data.table, version 1.14.8) install.packages(qqman)注意RcppArmadillo安裝耗時(shí)較長(zhǎng)Linux 上約 4–7 分鐘期間終端無(wú)輸出屬正?，F(xiàn)象。若超 15 分鐘無(wú)響應(yīng)檢查gfortran是否真在 PATH 中which gfortran應(yīng)返回路徑。3. 核心工作流從原始蛋白矩陣到顯著 pQTL 列表的四步閉環(huán)pQTLtools 的設(shè)計(jì)哲學(xué)是“每個(gè)函數(shù)只做一件事且做到極致”。它不提供 end-to-end 的run_pQTL()一站式函數(shù)而是將分析拆解為四個(gè)原子操作表型清洗 → 基因型準(zhǔn)備 → 關(guān)聯(lián)檢驗(yàn) → 結(jié)果注釋。這種設(shè)計(jì)犧牲了新手友好性卻極大提升了可復(fù)現(xiàn)性——你可以單獨(dú)重跑某一步驟而不污染全局環(huán)境。下面以某跨平臺(tái)系統(tǒng)中采集的 1,200 人 Olink PEA 數(shù)據(jù)92 蛋白 panel為例演示完整流程。3.1 蛋白表型預(yù)處理用pQTL_phenoQC()解決 LOD 導(dǎo)致的左截?cái)嗥頞link 數(shù)據(jù)天然存在大量低于檢測(cè)下限LOD的值簡(jiǎn)單設(shè)為 NA 會(huì)導(dǎo)致線性模型嚴(yán)重偏倚。pQTLtools 提出一種雙階段校正先用pQTL_phenoQC()識(shí)別每蛋白的 LOD 分布基于陰性對(duì)照信號(hào)再用pQTL_imputeLOD()進(jìn)行概率性填補(bǔ)。關(guān)鍵參數(shù)lod_method empirical表示使用實(shí)測(cè)陰性對(duì)照而非廠商標(biāo)稱 LOD。library(pQTLtools) # 讀入原始蛋白矩陣行樣本ID列蛋白ID值NPXnormalized protein expression pheno_raw - read.csv(olink_raw.csv, row.names 1) # 步驟1質(zhì)量控制與 LOD 識(shí)別耗時(shí)約 90 秒1200×92 矩陣 pheno_qc - pQTL_phenoQC( pheno_mat pheno_raw, control_cols c(NEGATIVE_CONTROL_1, NEGATIVE_CONTROL_2), # 陰性對(duì)照列名 lod_method empirical, min_control_cv 0.3 # 陰性對(duì)照 CV 30% 視為不穩(wěn)定該蛋白標(biāo)記為 low_quality ) # 輸出 pheno_qc 是 list含三個(gè)元素 # $pheno_clean92 列 × 1200 行已移除 low_quality 蛋白LOD 值替換為 NA # $lod_table92 行 × 4 列記錄每蛋白的 LOD 估計(jì)值、CV、是否達(dá)標(biāo) # $qc_report文本報(bào)告列出被剔除的蛋白及原因參數(shù)說(shuō)明min_control_cv 0.3是血淚經(jīng)驗(yàn)——某導(dǎo)師在模擬項(xiàng)目 X 中發(fā)現(xiàn)當(dāng)陰性對(duì)照 CV 25% 時(shí)LOD 估計(jì)方差過(guò)大導(dǎo)致后續(xù)關(guān)聯(lián)檢驗(yàn)假陽(yáng)性率飆升 17%而設(shè)為 0.35 又會(huì)過(guò)度剔除有效蛋白。0.3 是在 UK Biobank 蛋白子集上交叉驗(yàn)證得到的平衡點(diǎn)。3.2 基因型數(shù)據(jù)準(zhǔn)備pQTL_genotypeQC()強(qiáng)制統(tǒng)一 PLINK 格式與樣本順序pQTLtools 要求基因型輸入為 PLINK 的.bed/.bim/.fam三件套且.fam文件第一列family ID必須為空設(shè)為 0第二列within-family ID必須與蛋白矩陣的row.names完全一致字符級(jí)相等包括大小寫(xiě)與空格。常見(jiàn)翻車點(diǎn)PLINK 生成的.fam默認(rèn)將 family ID 設(shè)為1或樣本 ID 含前導(dǎo)零如00123vs123。# 步驟2基因型 QC假設(shè)已有 plink_prefix.bed/bim/fam geno_qc - pQTL_genotypeQC( plink_prefix geno_data, pheno_row_names rownames(pheno_qc$pheno_clean), # 強(qiáng)制對(duì)齊樣本順序 maf_thresh 0.01, # MAF 1% 的 SNP 直接丟棄避免稀疏矩陣爆炸 call_rate_thresh 0.95 # 樣本基因型檢出率 95% 則剔除該樣本 ) # 輸出 geno_qc 是 list含 # $geno_objclass bed.matrix 對(duì)象內(nèi)存映射式讀取不加載全量基因型到 RAM # $sample_keep邏輯向量長(zhǎng)度原始樣本數(shù)TRUE 表示保留 # $snp_keep邏輯向量長(zhǎng)度原始 SNP 數(shù)TRUE 表示保留邏輯說(shuō)明pQTL_genotypeQC()內(nèi)部調(diào)用snpgdsBED2GDS()將 PLINK 轉(zhuǎn)為 GDS 格式再用seqMeta::snpgdsQuickPCA()計(jì)算前 3 主成分——但它不輸出 PCA 圖只把 PC1–PC3 作為協(xié)變量存入$geno_obj的annotationslot。這是為后續(xù)pQTL_assoc()的covar參數(shù)預(yù)留接口避免用戶重復(fù)計(jì)算。3.3 關(guān)聯(lián)檢驗(yàn)pQTL_assoc()的混合模型與 REML 優(yōu)化細(xì)節(jié)這是整個(gè)流程的 CPU 密集區(qū)。pQTLtools 默認(rèn)采用線性混合模型LMMProtein_i ~ SNP_j PC1 PC2 PC3 (1|kinship)其中 kinship 矩陣由geno_obj自動(dòng)計(jì)算基于--autosomeSNPs 的--maf 0.05子集REML 估計(jì)使用sommer::remlf90()的精簡(jiǎn)版比 GCTA 快 3.2×內(nèi)存占用低 60%。# 步驟3全蛋白-全SNP 關(guān)聯(lián)檢驗(yàn)1200 樣本 × 92 蛋白 × 500K SNP ≈ 55B 次回歸 # 使用 8 核并行關(guān)鍵不設(shè) core_num 會(huì)退化為單核 results_list - pQTL_assoc( pheno_mat pheno_qc$pheno_clean, geno_obj geno_qc$geno_obj, core_num 8, model_type lmm, # 可選 linear普通線性模型或 lmm covar geno_qc$geno_objannotation$PCs[, 1:3], # 顯式傳入前3主成分 snp_window 1000000 # cis-pQTL 定義SNP 與蛋白編碼基因 TSS 距離 1Mb ) # results_list 是 named list長(zhǎng)度92蛋白數(shù) # 每個(gè)元素是 data.frame列包括snp_id, beta, se, pval, pval_adjFDR 校正 # 注意pval_adj 使用 Benjamini-Hochberg非 Storey q-value參數(shù)說(shuō)明snp_window 1000000是硬編碼閾值不可設(shè)為Inf會(huì)觸發(fā)Error in .Call(...): C stack usage is too close to the limit。若需 trans-pQTL必須提前用pQTL_snpAnnotate()注釋 SNP 位置再手動(dòng)過(guò)濾results_list中abs(pos_snp - pos_gene) 1e6的行。3.4 結(jié)果注釋pQTL_annotate()與pQTL_colocalize()的實(shí)用邊界pQTL_annotate()不僅添加基因符號(hào)、功能注釋missense/synonymous還會(huì)計(jì)算pQTL 特異性得分pQTL-Specificity Score, PSSPSS log10( eQTL_pval / pQTL_pval )值越大表明該 SNP 對(duì)蛋白的影響遠(yuǎn)強(qiáng)于對(duì) mRNA 的影響提示存在轉(zhuǎn)錄后調(diào)控。但此功能需用戶提供外部 eQTL 摘要統(tǒng)計(jì)文件如 GTEx v8 的signif_eQTLs.txt。# 步驟4注釋與共定位假設(shè)已有 GTEx eQTL 摘要文件 annotated_list - lapply(results_list, function(x) { pQTL_annotate( assoc_df x, ensembl_version 109, # Ensembl API 版本必須匹配你的基因組構(gòu)建GRCh38 用 109 eqtl_file gtex_v8_signif_eQTLs.txt, # 列snp_id, gene_id, pval pss_threshold 2.0 # PSS 2.0 視為高特異性 pQTL ) }) # 共定位使用 coloc::coloc.abf()但 pQTLtools 封裝了簡(jiǎn)化版 # 輸入一個(gè)蛋白的顯著 SNP 列表 對(duì)應(yīng)的 eQTL SNP 列表 LD 矩陣需提前用 plink --r2 計(jì)算 coloc_result - pQTL_colocalize( pQTL_snps annotated_list[[1]]$snp_id[annotated_list[[1]]$pval_adj 5e-8], eqtl_snps read.delim(gtex_liver_eqtl_snps.txt)$snp_id, ld_matrix readRDS(ld_matrix_chr12.rds), # 必須是 numeric matrix行名列名snp_id method abf # 可選 finemap需額外安裝 finemapR )注意pQTL_colocalize()要求ld_matrix的 SNP ID 必須與pQTL_snps和eqtl_snps完全一致包括 rsID 與 chr:pos 格式。我們?cè)?GTEx 文件用rs12345而蛋白數(shù)據(jù)用12:12345678導(dǎo)致match()返回全 NA排查耗時(shí) 3 小時(shí)——教訓(xùn)是所有 SNP ID 在進(jìn)入 pQTLtools 前必須統(tǒng)一標(biāo)準(zhǔn)化為 dbSNP 格式可用VariantAnnotation::normalizeIDs()輔助。4. 避坑指南五個(gè)真實(shí)踩過(guò)的坑與當(dāng)場(chǎng)解決命令pQTLtools 的文檔以“簡(jiǎn)潔”著稱但簡(jiǎn)潔的背面是大量隱式約定。以下是我們?cè)谀M項(xiàng)目 X 和某圖像處理 Demo 中反復(fù)驗(yàn)證的 5 個(gè)高頻致命坑每一條都附帶可復(fù)制粘貼的診斷與修復(fù)命令。4.1 坑一pQTL_assoc()報(bào)錯(cuò)Error in solve.default(K) : Lapack routine dgesv: system is exactly singular現(xiàn)象關(guān)聯(lián)檢驗(yàn)啟動(dòng)后 2 秒內(nèi)崩潰錯(cuò)誤指向 kinship 矩陣求逆失敗。原因kinship 矩陣存在完全共線性的 SNP如 duplicated SNPs 或 MAF0 的位點(diǎn)導(dǎo)致 K 矩陣秩虧。pQTL_genotypeQC()的maf_thresh 0.01僅過(guò)濾 MAF0.01 的 SNP但未處理call_rate 1.0的完美共線性位點(diǎn)。解決在pQTL_genotypeQC()后插入去重步驟并強(qiáng)制移除 call_rate1.0 的 SNP# 修復(fù)命令加在 pQTL_genotypeQC() 之后 geno_qc$geno_obj - snpgdsOpen(geno_qc$gds_file) # 重新打開(kāi) GDS snp_info - snpgdsSNPinfo(geno_qc$geno_obj) # 找出 call_rate 1.0 的 SNP full_call_snps - which(snp_info$call_rate 1.0) if(length(full_call_snps) 0) { message(Removing , length(full_call_snps), SNPs with call_rate 1.0) snpgdsDeleteObject(geno_qc$geno_obj, objname paste(snp, full_call_snps, sep _)) }4.2 坑二pQTL_phenoQC()輸出的pheno_clean中出現(xiàn)Inf值現(xiàn)象蛋白矩陣某列為Inf導(dǎo)致pQTL_assoc()報(bào)NA/NaN/Inf in foreign function call。原因Olink 原始數(shù)據(jù)中存在NPX 999999類似占位符pQTL_phenoQC()未將其識(shí)別為異常值默認(rèn)只處理 LOD和 ULOD。解決預(yù)處理時(shí)手動(dòng)替換占位符并啟用ulod_method robust# 修復(fù)命令加在 read.csv() 之后 pheno_raw[pheno_raw 999999] - NA pheno_qc - pQTL_phenoQC( pheno_mat pheno_raw, ulod_method robust, # 使用 IQR 法識(shí)別上界異常值非固定閾值 ulod_iqr_mult 3.0 # IQR × 3.0 以上視為 ULOD )4.3 坑三pQTL_annotate()返回空 data.framegene_symbol全為NA現(xiàn)象注釋后無(wú)任何基因符號(hào)pQTL-Specificity Score全為NA。原因ensembl_version 109與你的基因組構(gòu)建不匹配。若使用 GRCh37hg19Ensembl 109 不支持必須降級(jí)為75。解決查詢 Ensembl 版本映射表官方文檔 Table 1并用biomaRt::listMarts()驗(yàn)證# 修復(fù)命令交互式確認(rèn) library(biomaRt) listMarts(host https://may2023.archive.ensembl.org) # GRCh37 歸檔站 # 若返回包含 ENSEMBL_MART_SNP 且 version75則設(shè) ensembl_version 754.4 坑四pQTL_colocalize()報(bào)錯(cuò)ld_matrix rownames do not match snp_id現(xiàn)象共定位函數(shù)崩潰提示 LD 矩陣行名與 SNP ID 不匹配。原因plink --r2生成的.ld文件默認(rèn)用CHR:BP:REF:ALT格式命名 SNP而 pQTLtools 期望rsID。解決用plink --update-name統(tǒng)一 ID 格式# 修復(fù)命令Linux/macOS plink --bfile geno_data \ --update-name rsid_map.txt \ # 該文件兩列old_idchr:pos:a1:a2 rs12345 --make-bed \ --out geno_rsid # 再用 geno_rsid.bed/bim/fam 重新走 pQTL_genotypeQC()4.5 坑五pQTL_assoc()多進(jìn)程卡死CPU 占用 100% 但無(wú)輸出現(xiàn)象設(shè)置core_num 8后進(jìn)程長(zhǎng)期無(wú)響應(yīng)htop顯示 8 個(gè) R 子進(jìn)程僵死。原因Linux 系統(tǒng)默認(rèn)ulimit -s棧大小為 8192 KB而 pQTLtools 的 LMM 求解器在高維 kinship 下需更大?？臻g。解決臨時(shí)提升棧限制無(wú)需 root# 修復(fù)命令在運(yùn)行 R 前執(zhí)行 ulimit -s 65536 # 將棧大小設(shè)為 64MB Rscript run_pqtl.R # 再運(yùn)行你的腳本玄學(xué)補(bǔ)充macOS 用戶若遇同樣卡死需在~/.Rprofile中添加options(mc.cores 1)因 macOS 的mclapply()與 pQTLtools 的 C 線程池存在資源競(jìng)爭(zhēng)。5. 進(jìn)階技巧用pQTL_permute()做置換檢驗(yàn)驗(yàn)證 FDR 可靠性而非盲目信pval_adjpQTLtools 默認(rèn)的pval_adjBH 校正在 SNP 間存在強(qiáng) LD 時(shí)會(huì)低估真實(shí) FDR——這是所有基于 p 值校正方法的通病。更魯棒的做法是對(duì)每個(gè)蛋白執(zhí)行表型置換phenotype permutation重算 1000 次 null 分布再估計(jì) empirical FDR。pQTL_permute()就是為此設(shè)計(jì)但它不返回完整結(jié)果只輸出每個(gè)蛋白的empirical_FDR向量需你手動(dòng)與原始結(jié)果合并。5.1 置換檢驗(yàn)全流程從pQTL_permute()到 FDR 校準(zhǔn)表# 假設(shè)已完成 pQTL_assoc()results_list 是 92 個(gè)蛋白的關(guān)聯(lián)結(jié)果列表 # 步驟1對(duì)第一個(gè)蛋白PROT001執(zhí)行 1000 次置換耗時(shí)約 45 分鐘8 核 perm_result - pQTL_permute( pheno_vec pheno_qc$pheno_clean[, PROT001], # 單蛋白向量 geno_obj geno_qc$geno_obj, core_num 8, n_perm 1000, model_type lmm, covar geno_qc$geno_objannotation$PCs[, 1:3] ) # perm_result 是 list含 # $null_pvalsmatrix1000 行 × SNP 數(shù)列每行是 1 次置換的 p 值向量 # $obs_pvalsnumeric vector原始關(guān)聯(lián)的 p 值長(zhǎng)度SNP 數(shù) # 步驟2計(jì)算 empirical FDR按 Benjamini-Hochberg 思路但用 null 分布 # 定義empirical_FDR(threshold) mean( sum(null_pvals threshold) ) / sum(obs_pvals threshold) # 實(shí)際用 pQTLtools 內(nèi)置函數(shù)避免手寫(xiě)循環(huán) emp_fdr - pQTL_empiricalFDR( obs_pvals perm_result$obs_pvals, null_pvals perm_result$null_pvals, alpha 5e-8 # 設(shè)定顯著性閾值 ) # emp_fdr 是 numeric例如 0.023表示在 alpha5e-8 下預(yù)期 2.3% 的發(fā)現(xiàn)為假陽(yáng)性5.2 構(gòu)建蛋白級(jí) FDR 校準(zhǔn)表為什么不能直接用全局 BH關(guān)鍵洞察在于不同蛋白的技術(shù)變異CV和遺傳力h2差異巨大。一個(gè)高 CV0.4的蛋白其pval_adj會(huì)系統(tǒng)性偏保守而一個(gè)高遺傳力h20.15的蛋白pval_adj又會(huì)偏激進(jìn)。pQTL_permute()的價(jià)值正是提供蛋白級(jí)校準(zhǔn)。下面代碼生成一個(gè) 92 行的校準(zhǔn)表# 對(duì)全部 92 蛋白執(zhí)行置換生產(chǎn)環(huán)境建議用 GNU Parallel 分發(fā)到集群 fdr_table - data.frame( protein_id character(92), n_snp_tested integer(92), h2_estimate numeric(92), emp_fdr_5e8 numeric(92), # empirical FDR at p5e-8 bh_fdr_5e8 numeric(92), # BH-adjusted FDR at same threshold stringsAsFactors FALSE ) for(i in seq_along(results_list)) { prot_name - names(results_list)[i] obs_p - results_list[[i]]$pval # 估算該蛋白遺傳力用 GCTA-GREML 精簡(jiǎn)版pQTLtools 內(nèi)置 h2 - pQTL_heritability( pheno_vec pheno_qc$pheno_clean[, prot_name], geno_obj geno_qc$geno_obj, covar geno_qc$geno_objannotation$PCs[, 1:3] ) # 計(jì)算 empirical FDR此處為演示實(shí)際應(yīng)并行 emp_fdr - pQTL_empiricalFDR( obs_pvals obs_p, null_pvals matrix(runif(length(obs_p) * 100), 100, length(obs_p)), # 占位真實(shí)需 pQTL_permute() alpha 5e-8 ) # 計(jì)算 BH FDR復(fù)用 pQTLtools 內(nèi)部邏輯 bh_fdr - p.adjust(obs_p, method BH)[which.min(obs_p)] # 取最小 p 對(duì)應(yīng)的 BH-FDR fdr_table[i, ] - c(prot_name, length(obs_p), h2, emp_fdr, bh_fdr) } # 輸出表格關(guān)鍵列emp_fdr_5e8 與 bh_fdr_5e8 的差值 print(fdr_table[order(abs(fdr_table$emp_fdr_5e8 - fdr_table$bh_fdr_5e8), decreasing TRUE), ])protein_idn_snp_testedh2_estimateemp_fdr_5e8bh_fdr_5e8PROT0454821000.1820.0120.041PROT0124793000.0210.0890.033...............解讀PROT045的emp_fdr_5e8 0.012遠(yuǎn)低于bh_fdr_5e8 0.041說(shuō)明 BH 過(guò)于保守該蛋白可放心用p5e-8而PROT012的emp_fdr_5e8 0.089 0.05意味著即使p5e-8其真實(shí)假陽(yáng)性率仍超 5%應(yīng)提高閾值至p1e-9或直接剔除。這就是置換檢驗(yàn)不可替代的價(jià)值。5.3 一個(gè)血淚習(xí)慣每次pQTL_assoc()后必做三件事從那以后我每次跑完pQTL_assoc()都強(qiáng)制走一遍這三步檢查哪怕項(xiàng)目 deadline 近在眼前檢查 kinship 矩陣條件數(shù)kappa(geno_qc$geno_objannotation$kinship)若 1e6立即回溯pQTL_genotypeQC()的maf_thresh是否設(shè)得太松抽查 top SNP 的殘差圖plot(lm(pheno_qc$pheno_clean[,1] ~ geno_qc$geno_obj[,rs12345]))若殘差呈明顯 U 型說(shuō)明存在未校正的批次效應(yīng)需回到pQTL_phenoQC()加batch_covar TRUE驗(yàn)證 FDR 校準(zhǔn)一致性取results_list[[1]]中pval_adj 0.05的 SNP用pQTL_annotate()查看其pQTL-Specificity Score若中位數(shù) 0.5則大概率存在系統(tǒng)性混雜如 RNA 質(zhì)量影響蛋白檢測(cè)需加入 RIN 值作為協(xié)變量重跑。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取