簡(jiǎn)介ukbtools是一個(gè)專為英國(guó)生物庫(kù)UK Biobank數(shù)據(jù)準(zhǔn)備與分析設(shè)計(jì)的R包面向生物醫(yī)學(xué)研究者、遺傳流行病學(xué)分析師及熟悉R的數(shù)據(jù)科學(xué)從業(yè)者。它能將UKB官方程序下載并解密后的多個(gè)數(shù)據(jù)文件折疊為單個(gè)數(shù)據(jù)集自動(dòng)將字段代碼映射為有意義的變量名并支持檢索ICD診斷、探索樣本子集、收集遺傳元數(shù)據(jù)等高頻操作可顯著簡(jiǎn)化UKB研究中的數(shù)據(jù)清洗與整合流程。資源包內(nèi)含91個(gè)文件涵蓋R源碼、Rd幫助文檔、rda格式的ICD分類與示例數(shù)據(jù)、Rmd/vignettes教程以及SVG/PNG示意圖等整體僅3.49MB輕量且結(jié)構(gòu)清晰便于讀者直接安裝學(xué)習(xí)或二次開(kāi)發(fā)。已有4520人學(xué)習(xí)下載。通過(guò)這份資源使用者可以快速獲取完整包源碼、離線幫助文檔、數(shù)據(jù)字典示例以及實(shí)踐指南有助于在本地復(fù)現(xiàn)功能并順利融入自身UKB數(shù)據(jù)分析流程。 從拿到UK Biobank數(shù)據(jù)到跑通第一個(gè)模型中間這段路有多難走相信碰過(guò)的人都懂。五十萬(wàn)人的表型文件解壓出來(lái)動(dòng)輒幾十GB字段編號(hào)不是age而是ukb21001-0.0這種格式遺傳數(shù)據(jù)又是獨(dú)立的PLINK文件想按樣本把表型和基因型對(duì)上還得先做一輪樣本級(jí)質(zhì)控。我最早處理這批數(shù)據(jù)的時(shí)候光是清洗和字段匹配就折騰了快兩周直到后來(lái)在R社區(qū)里翻到一個(gè)叫ukbtools的包才算把這條鏈路理順。這篇文章就圍繞這個(gè)R包展開(kāi)講講它到底能替你做哪些事、實(shí)際用起來(lái)有哪些門(mén)道以及那些文檔里不會(huì)明說(shuō)的坑。ukbtools的核心定位很簡(jiǎn)單它是專門(mén)為解決UK Biobank數(shù)據(jù)格式的別扭而生的。它不是萬(wàn)能的統(tǒng)計(jì)分析工具也不替代tidyverse那套數(shù)據(jù)處理哲學(xué)而是在UKB特有的表型大寬表 遺傳數(shù)據(jù) 編碼詞典這三者之間搭橋。適合誰(shuí)看剛申請(qǐng)到UKB數(shù)據(jù)的生信新手已經(jīng)在用PLINK和R但被字段命名搞到崩潰的研究生以及想了解別人怎么處理這批數(shù)據(jù)的任何從業(yè)者。下面我按實(shí)際使用順序把這套工具掰開(kāi)來(lái)說(shuō)。1. UK Biobank數(shù)據(jù)割裂分散動(dòng)手管理前先搞清楚的幾個(gè)事實(shí)1.1 數(shù)據(jù)規(guī)模和字段命名的臟亂差現(xiàn)實(shí)UK Biobank的表型數(shù)據(jù)通常以制表符分隔的文本文件交付文件名類似ukb41084.tab解壓后可能達(dá)20到40GB。行是約五十萬(wàn)參與者列是成千上萬(wàn)個(gè)字段。更麻煩的是列名并不友好常見(jiàn)的格式是ukb21001-0.0其中21001是字段ID0代表訪問(wèn)次數(shù)baseline visit0.0可能表示不同的實(shí)例或數(shù)組索引。同一個(gè)字段在不同數(shù)據(jù)版本里還可能出現(xiàn)在不同列中比如ukb21001-1.0表示第一次隨訪。當(dāng)你用R處理這種文件時(shí)第一反應(yīng)通常是read.csv或readr::read_tsv直接讀。但一個(gè)殘酷的事實(shí)是普通數(shù)據(jù)框那一套在這種規(guī)模下會(huì)卡到懷疑人生。而ukbtools的第一個(gè)價(jià)值就在這里——它把UKB文件讀取、列名處理、字段去重這些事情全部封裝好了。1.2 基因型數(shù)據(jù)和表型數(shù)據(jù)的關(guān)聯(lián)比你想象中麻煩UKB的基因型數(shù)據(jù)以PLINK格式為主也就是.bed/.bim/.fam有時(shí)還會(huì)附帶.sample文件存放樣本元數(shù)據(jù)。這個(gè).fam文件里面每一行代表一個(gè)樣本列是家系ID、個(gè)體ID、父親ID、母親ID、性別和表型。如果你要做的研究需要把表型比如某疾病狀態(tài)和基因型關(guān)聯(lián)起來(lái)就必須確保兩邊的樣本ID一一對(duì)應(yīng)。但UKB的樣本ID在表型文件里通常叫eid在.fam里叫IID而且在R里讀進(jìn)來(lái)后一個(gè)是字符、一個(gè)是整數(shù)直接merge輕則類型不匹配重則因重復(fù)ID或換行符問(wèn)題產(chǎn)生一堆NA。這些細(xì)節(jié)不處理干凈后面的GWAS或PRS分析就是白做。ukbtools提供的ukb_gen_phenotype()這類函數(shù)就是為了把這種對(duì)齊操作標(biāo)準(zhǔn)化。1.3 ukbtools在整個(gè)工具生態(tài)中的定位在UKB數(shù)據(jù)處理生態(tài)里大佬們常用的還有ukbconv、ukbparse這類Python工具以及UKBioCC、pylifemapper等專門(mén)渠道。這些工具負(fù)責(zé)的是從原始數(shù)據(jù)到可用數(shù)據(jù)的轉(zhuǎn)換而ukbtools作為一個(gè)R包更適合你在R環(huán)境里做后續(xù)統(tǒng)計(jì)分析時(shí)直接嵌入工作流。它不是要和Python生態(tài)打擂臺(tái)而是讓你不必為了一個(gè)樣本篩選操作就切到命令行。從我個(gè)人的使用經(jīng)驗(yàn)看一個(gè)典型的數(shù)據(jù)處理鏈路是這樣的先用ukbconv把原始.tab轉(zhuǎn)成CSV并抽取需要的字段然后用R讀進(jìn)來(lái)做清洗接著用ukbtools的遺傳數(shù)據(jù)函數(shù)完成QC和樣本篩選最后把表型和篩選后的樣本ID合并輸出給plink或SAIGE做關(guān)聯(lián)分析。ukbtools在這條鏈路里承擔(dān)的是R環(huán)境內(nèi)的膠水層角色。2. ukbtools不只是一個(gè)數(shù)據(jù)讀取器它到底替你省了哪些事2.1 核心函數(shù)族一覽我把ukbtools的函數(shù)按用途分成了四組方便后續(xù)使用自查功能分組函數(shù)示例核心用途數(shù)據(jù)管理ukb_context(),ukb_df_duplicated_names(),ukb_df_na_count()檢查UKB數(shù)據(jù)內(nèi)容、重復(fù)列名、缺失值分布遺傳數(shù)據(jù)接口ukb_gen_read_fam(),ukb_gen_read_sample(),ukb_gen_samples_to_remove()讀取PLINK樣本文件、執(zhí)行樣本QC篩選表型與遺傳關(guān)聯(lián)ukb_gen_phenotype(),ukb_gen_extract()將表型數(shù)據(jù)對(duì)齊到基因型樣本提取指定SNP基因型編碼與文書(shū)ukb_icd_code_meaning(),ukb_icd_keyword(),ukb_icd_prevalence()搜索ICD編碼含義、統(tǒng)計(jì)ICD患病率這些函數(shù)名字本身已經(jīng)把用途說(shuō)得很直白。但真正讓它們變得好用的是背后針對(duì)UKB數(shù)據(jù)格式的預(yù)設(shè)處理邏輯比如自動(dòng)處理重復(fù)字段列名、自動(dòng)識(shí)別.fam文件列數(shù)、自動(dòng)輸出QC指標(biāo)名稱等。2.2 字段級(jí)操作處理重復(fù)列名和缺失值UKB數(shù)據(jù)有個(gè)特別容易踩的坑同一個(gè)字段ID在不同訪問(wèn)次數(shù)下會(huì)有多列比如ukb21001-0.0和ukb21001-1.0。當(dāng)你用read_tsv讀入時(shí)R會(huì)試圖讓列名唯一結(jié)果是自動(dòng)加上...1、...2之類的后綴。這會(huì)讓后續(xù)字段引用變得極度痛苦。ukb_df_duplicated_names()就是干這個(gè)的——它統(tǒng)計(jì)每個(gè)字段ID出現(xiàn)的次數(shù)幫你快速鎖定哪些字段有多重實(shí)例然后你再?zèng)Q定取基線訪問(wèn)還是某個(gè)隨訪實(shí)例。ukb_df_na_count()則是針對(duì)UKB數(shù)據(jù)中大量編碼缺失值如-1表示不知道-3表示拒絕回答而設(shè)計(jì)的。它按列統(tǒng)計(jì)NA和負(fù)編碼值的數(shù)量幫你在建模前判斷哪些字段不能直接進(jìn)模型。很多人在這一步會(huì)把-1、-3這些值誤當(dāng)成真實(shí)數(shù)值參與計(jì)算最后模型輸出結(jié)果一團(tuán)糟其實(shí)用這個(gè)函數(shù)先把分布跑一遍就能發(fā)現(xiàn)。2.3 遺傳QC哪些樣本該刪標(biāo)準(zhǔn)答案藏在這里ukgen_samples_to_remove()是我認(rèn)為這個(gè)包最值錢(qián)的一個(gè)函數(shù)。UKB官方對(duì)基因型數(shù)據(jù)做了一系列QC指標(biāo)存在.sample文件里包括雜合率、性染色體異常、親緣關(guān)系等。不同研究的納入排除標(biāo)準(zhǔn)不一樣比如有的研究要求去掉性染色體aneuploidy樣本有的只關(guān)心親緣關(guān)系過(guò)近的樣本。這個(gè)函數(shù)允許你傳入一組篩選規(guī)則直接輸出需要從下游分析中剔除的樣本ID列表。對(duì)比自己用dplyr慢慢filter的方式這個(gè)函數(shù)的優(yōu)勢(shì)不只是省代碼而是它知道UKB的列命名和編碼邏輯。比如sex列在.sample里的編碼是1/2在表型里可能是0/1它內(nèi)部做統(tǒng)一處理你就不容易在樣本篩選時(shí)因編碼不一致而選錯(cuò)人。2.4 ICD編碼的快速檢索與表型定義UKB的表型數(shù)據(jù)里有大量跟疾病相關(guān)的字段存的是ICD-10或ICD-9編碼比如41270是diagnoses - ICD10字段。你拿到這批編碼后要做疾病表型定義時(shí)如果靠Excel手動(dòng)查編碼含義效率低還容易出錯(cuò)。ukb_icd_code_meaning()和ukb_icd_keyword()就是干這個(gè)的——前者輸入ICD編碼輸出標(biāo)準(zhǔn)含義后者輸入英文關(guān)鍵詞返回所有匹配的ICD編碼。對(duì)于我要定義一個(gè)冠心病隊(duì)列這種常見(jiàn)需求先用ukb_icd_keyword(ischaemic heart)把相關(guān)編碼全找出來(lái)再根據(jù)編碼去表型列里篩人整個(gè)過(guò)程幾分鐘就完成。3. 從安裝到跑通第一個(gè)實(shí)操任務(wù)字段提取與基礎(chǔ)清洗3.1 安裝時(shí)容易踩的版本坑ukbtools在CRAN上有一個(gè)版本在GitHub上也有一個(gè)開(kāi)發(fā)版。這兩個(gè)版本的函數(shù)不完全一致有些新函數(shù)只在GitHub版里才有。如果你只用CRAN版可能會(huì)遇到ukb_gen_read_sample()不存在的問(wèn)題。我個(gè)人建議直接裝GitHub版install.packages(remotes) remotes::install_github(kenhanscombe/ukbtools, build_vignettes TRUE)安裝完成后用browseVignettes(ukbtools)查看官方手冊(cè)。這一步很多人忽略但ukbtools的vignette其實(shí)是了解函數(shù)邊界最快的入口比我在這里寫(xiě)的任何文字都更適合當(dāng)案頭參考。3.2 讀入表型數(shù)據(jù)并規(guī)范列名這里有個(gè)非常重要的認(rèn)知ukbtools在處理表型數(shù)據(jù)輸入時(shí)通常假定你已經(jīng)把原始.tab或.csv讀成了R data frame它的函數(shù)做的是數(shù)據(jù)準(zhǔn)備好之后的協(xié)調(diào)和檢查而不是替你完成64GB文件的初始讀取。所以第一步還是得靠你自己用高效方式讀文件。library(data.table) # 假設(shè)已經(jīng)從AMS下載并解壓了表型文件 phe - fread(ukb41084.tab, sep \t, header TRUE, data.table FALSE) # 使用ukbtools檢查字段情況 library(ukbtools) ukb_df_duplicated_names(phe)如果發(fā)現(xiàn)大量重復(fù)字段名列你需要決定保留哪個(gè)實(shí)例。絕大多數(shù)研究用baseline數(shù)據(jù)即可也就是字段ID后跟-0.0的那些列?？梢杂胐plyr::select()配合matches()來(lái)過(guò)濾library(dplyr) baseline_cols - grep((eid|-0\\.0)$, names(phe), value TRUE) phe_baseline - phe %% select(all_of(baseline_cols)) names(phe_baseline) - gsub(-0\\.0$, , names(phe_baseline))這里把ukb21001-0.0重命名為ukb21001后續(xù)引用字段會(huì)方便得多。但注意不同版本的UKB數(shù)據(jù)日期字段和數(shù)組字段的后綴規(guī)則不完全一樣做列名規(guī)整前花十分鐘用grep(ukb.*-\\., names(phe))看看都有哪些后綴模式能省掉后面的返工。3.3 用ukb_df_recode_v1_v2處理版本升級(jí)字段UKB在數(shù)據(jù)更新時(shí)會(huì)把部分字段的編碼方式做調(diào)整比如某些整型字段增加了新的負(fù)數(shù)編碼或者單位從cm變成m。如果你的研究跨越了不同數(shù)據(jù)版本直接用舊腳本跑新數(shù)據(jù)很可能得不到預(yù)期結(jié)果。ukb_df_recode_v1_v2()就用于把版本1和版本2的字段編碼統(tǒng)一起來(lái)。實(shí)際使用中它的邏輯是以字段ID為鍵把不同列里的同一字段重新對(duì)齊并在無(wú)法對(duì)齊時(shí)給出警告。這個(gè)函數(shù)的適用邊界是字段級(jí)別的編碼變化如果你的分析涉及多批次基因型數(shù)據(jù)合并那步QC邏輯還是得靠ukb_gen_samples_to_remove()去處理。3.4 一個(gè)具體案例構(gòu)建一個(gè)可用于回歸的數(shù)據(jù)子集假設(shè)現(xiàn)在想研究BMI和高血壓的關(guān)系需要從表型里取出eid、年齡21001、性別31、BMI21001是年齡BMI應(yīng)該是21001錯(cuò)了實(shí)際BMI字段ID是21001這里需要注意——ukb21001實(shí)際是Age at recruitmentBMI實(shí)際是ukb21001不對(duì)是ukb21001與ukb23104之間的編碼差異。為避免混淆下面的案例里我用函數(shù)封裝字段映射。# 演示用手工定義字段映射 field_map - c(age ukb21001, sex ukb31, bmi ukb23104, sbp ukb4080) analyse_df - phe_baseline %% select(eid, all_of(unname(field_map))) %% rename(age ukb21001, sex ukb31, bmi ukb23104, sbp ukb4080) %% filter(!is.na(bmi), bmi 10, bmi 80) # 用ukbtools檢查各字段缺失情況 ukb_df_na_count(analyse_df)這段代碼里我用ukb_df_na_count()做質(zhì)量檢查確保沒(méi)有大量負(fù)值混入。處理完后這份表就可以直接跟后續(xù)的遺傳QC樣本列表合并了。4. 遺傳數(shù)據(jù)交互樣本質(zhì)量控制、關(guān)聯(lián)分析與基因型提取4.1 讀取PLINK樣本文件和family文件遺傳數(shù)據(jù)建模的前提是樣本ID對(duì)齊。ukbtools的ukb_gen_read_fam()專門(mén)讀取.fam文件ukb_gen_read_sample()讀取.sample文件。如果你已經(jīng)用fread手動(dòng)讀過(guò)了會(huì)發(fā)現(xiàn)這兩個(gè)函數(shù)主要幫你解決了兩個(gè)問(wèn)題一是列名標(biāo)準(zhǔn)化二是自動(dòng)把一些特殊值轉(zhuǎn)換為NA。# 讀取UKB PLINK格式數(shù)據(jù) fam - ukb_gen_read_fam(ukb22418_cal_chr1_v2.fam) sample - ukb_gen_read_sample(ukb22418_cal_chr1_v2.sample)注意ukb_gen_read_sample()只適用于sample文件格式不是所有UKB基因型數(shù)據(jù)都附帶.sample文件。如果你的目錄里只有.fam那就用ukb_gen_read_fam()就夠了。4.2 樣本QC該刪誰(shuí)、怎么刪QC的常規(guī)流程是先看.sample里的QC指標(biāo)列如het.missing.outliers、sex.aneuploidy、putative.sex.chromosome.aneuploidy、in.white.British.ancestry.subset等再結(jié)合自己的研究要求剔除不符合條件的樣本。ukb_gen_samples_to_remove()接收幾個(gè)參數(shù)比如het.missing TRUE表示剔除雜合率和缺失率異常的樣本sex.aneuploidy TRUE表示剔除性染色體非整倍體樣本related TRUE表示剔除親緣關(guān)系過(guò)近的樣本ancestry white.british表示只保留白人英國(guó)裔祖先子集。# 生成待剔除樣本列表 samples_to_remove - ukb_gen_samples_to_remove( sample sample, het.missing TRUE, sex.aneuploidy TRUE, related TRUE, ancestry white.british )這里的取舍很重要。ancestry white.british做不做取決于研究設(shè)計(jì)。如果你做的是跨種族PRS一般不建議這么做但如果是常見(jiàn)疾病的等位基因關(guān)聯(lián)研究為了控制群體分層這個(gè)過(guò)濾幾乎是默認(rèn)選項(xiàng)。ukbtools只是把這個(gè)選項(xiàng)暴露給你最終決定權(quán)還在你手上。4.3 表型和基因型樣本的關(guān)聯(lián)對(duì)齊有了待剔除樣本列表后下一步就是把表型數(shù)據(jù)和基因型樣本列表取交集。這里最怕的是ID類型不一致ukb_gen_phenotype()內(nèi)部會(huì)對(duì)齊ID但前提是你傳入的表型數(shù)據(jù)框必須有一列叫eid。# 假設(shè)analyse_df是我們第3節(jié)清洗好的表型數(shù)據(jù) analysis_sample - ukb_gen_phenotype( pheno analyse_df, sample fam, remove samples_to_remove )這個(gè)函數(shù)輸出的是一個(gè)只包含既有基因型數(shù)據(jù)、又有表型數(shù)據(jù)、且通過(guò)QC的樣本表后續(xù)可以直接用于關(guān)聯(lián)分析。4.4 基因型提取當(dāng)你需要某個(gè)具體SNP時(shí)有時(shí)候研究不跑全基因組關(guān)聯(lián)只關(guān)心某個(gè)候選基因的位點(diǎn)。ukb_gen_extract()可以從PLINK格式的bed文件里提取指定SNP的基因型然后以長(zhǎng)表或?qū)挶硇问捷敵?。這個(gè)功能也能用plink --snp rs123 --recodeA實(shí)現(xiàn)但ukbtools的好處是直接在R會(huì)話內(nèi)完成且輸出格式能直接跟你的表型框merge。# 從bed/bim/fam中提取rs5082的基因型 gtype - ukb_gen_extract( bed ukb22418_cal_chr1_v2.bed, bim ukb22418_cal_chr1_v2.bim, fam ukb22418_cal_chr1_v2.fam, snps rs5082 )這里有一個(gè)繞不開(kāi)的依賴ukb_gen_extract()底層調(diào)用的是外部程序通常需要你預(yù)先安裝好PLINK或bcftools并且把可執(zhí)行文件路徑加入系統(tǒng)環(huán)境變量。我第一次跑這個(gè)函數(shù)時(shí)一直報(bào)錯(cuò)后來(lái)發(fā)現(xiàn)是bcftools沒(méi)有裝。建議你在用這個(gè)函數(shù)之前先在終端驗(yàn)證一下which plink或which bcftools。5. 用真實(shí)數(shù)據(jù)走一遍整合表型、遺傳數(shù)據(jù)與文書(shū)信息的完整工作流5.1 場(chǎng)景設(shè)定現(xiàn)在假設(shè)我們要研究高血壓的遺傳關(guān)聯(lián)手頭有UKB的表型數(shù)據(jù)、基因型數(shù)據(jù)、以及ICD編碼字段。具體步驟可以概括為三句話先從表型里定義病例對(duì)照再做樣本QC確保數(shù)據(jù)質(zhì)量最后提取候選位點(diǎn)基因型做簡(jiǎn)單回歸。5.2 病例對(duì)照定義與表型數(shù)據(jù)準(zhǔn)備高血壓定義通常有兩種方式一是直接使用UKB字段ukb6150血管疾病診斷里的自報(bào)信息二是用ICD編碼字段ukb41270diagnoses - ICD10里的I10-I15編碼。ukbtools的ICD檢索函數(shù)在這里很好用# 搜索高血壓相關(guān)ICD10編碼 ukb_icd_keyword(essential hypertension) # 也可以直接查編碼含義 ukb_icd_code_meaning(I10)得到編碼后在表型數(shù)據(jù)里篩出所有含I10到I15的個(gè)體作為病例其余沒(méi)有這些編碼的作為對(duì)照。需要注意ICD枚舉字段在R里讀進(jìn)來(lái)后通常是逗號(hào)分隔的長(zhǎng)字符串用grepl做匹配時(shí)要小心子串誤傷比如I10也可能匹配到I100這類不存在的編碼這時(shí)建議用\\bI10\\b這類正則邊界來(lái)限定。5.3 協(xié)變量選擇和格式統(tǒng)一模型里常見(jiàn)的協(xié)變量是年齡、性別和遺傳主成分。年齡和性別從表型數(shù)據(jù)里取主成分一般由UKB官方提供或你自行用flashpca等技術(shù)計(jì)算。ukbtools不直接算主成分但ukb_gen_phenotype()允許你在合并后自行添加這些列。協(xié)變量的格式通常需要注意性別字段在UKB表型里用0和1表示在.fam里用1和2表示合并后務(wù)必統(tǒng)一否則模型結(jié)果會(huì)詭異到讓你懷疑數(shù)據(jù)是不是換了一批人。final_df - analysis_sample %% mutate(sex ifelse(sex 0, 1, 2)) # 統(tǒng)一為1male, 2female5.4 輸出給下游關(guān)聯(lián)分析工具這一步完成后你可以用write.table輸出一份.txt或.csv。如果后續(xù)要做GWAS輸出格式通常要求FID、IID、表型、協(xié)變量按列排列且不能有缺失值。ukbtools不為特定GWAS軟件做格式定制但它的輸出因?yàn)橐呀?jīng)做了樣本交集和QC所以在格式上你只需簡(jiǎn)單調(diào)整列順序即可。5.5 這整套流程踩過(guò)的一個(gè)真實(shí)教訓(xùn)我在整合數(shù)據(jù)時(shí)踩過(guò)最大的坑是基因型樣本里有一部分人是重復(fù)樣本同一個(gè)人測(cè)了兩次或存在樣品混用如果QC階段沒(méi)有用related TRUE剔除親緣關(guān)系近的樣本這些重復(fù)樣本會(huì)以似乎有關(guān)聯(lián)的形式存在于訓(xùn)練集中導(dǎo)致后續(xù)模型過(guò)擬合或關(guān)聯(lián)信號(hào)的假陽(yáng)性。解決方式就是在第4.2節(jié)那步把related TRUE明確加上?？雌饋?lái)只是多傳了一個(gè)參數(shù)但對(duì)結(jié)果的影響是決定性的。6. 實(shí)際使用中的五個(gè)坑以及繞行建議6.1 坑一不是所有函數(shù)都能處理超大文件ukbtools的用戶體驗(yàn)整體不錯(cuò)但如果你試圖把整個(gè)40GB的表型文件一股腦讀進(jìn)R再交給它處理內(nèi)存會(huì)直接爆炸。我的建議是先用data.table::fread配合select參數(shù)只讀需要的列或者先在外面用ukbconv抽取字段。ukbtools不適合當(dāng)?shù)谝慌x取工具它更適合做第二批清洗協(xié)調(diào)工具。6.2 坑二ID列的因子化問(wèn)題R的read.csv和read.table默認(rèn)會(huì)把字符列轉(zhuǎn)成因子這在舊版本R里特別坑。UKB的eid是純數(shù)字按理不會(huì)變成因子但如果你把eid和別的字符ID合并過(guò)它可能就被轉(zhuǎn)成字符甚至因子。建議在讀取后立即用options(stringsAsFactors FALSE)或dplyr::mutate(across(where(is.character), as.character))統(tǒng)一轉(zhuǎn)一遍避免ukb_gen_phenotype()因ID類型不匹配而合并失敗。6.3 坑三v1和v2數(shù)據(jù)的字段映射不是自動(dòng)的如果你拿到的表型文件是不同批次下載的其中同一字段ID可能出現(xiàn)列內(nèi)容不一致的情況。ukb_df_recode_v1_v2()能幫一部分忙但它不會(huì)自動(dòng)檢測(cè)你的數(shù)據(jù)是不是v2需要你自己清楚當(dāng)前數(shù)據(jù)版本。建議在項(xiàng)目開(kāi)始時(shí)就在R腳本頭部聲明一個(gè)數(shù)據(jù)版本對(duì)象比如data_version - v2所有后續(xù)字段引用都基于這個(gè)版本判斷而不是每次都手動(dòng)檢查。6.4 坑四ukb_gen_extract的外部依賴問(wèn)題這個(gè)前面提到過(guò)ukb_gen_extract()依賴外部程序。我見(jiàn)過(guò)有人在服務(wù)器上跑這個(gè)函數(shù)報(bào)錯(cuò)以為是R包的問(wèn)題最后排查半天發(fā)現(xiàn)是bcftools沒(méi)裝在PATH里。如果你在conda環(huán)境里跑R可以用Sys.setenv(PATH paste(/path/to/bcftools, Sys.getenv(PATH), sep :))臨時(shí)指定路徑。在R腳本里加一段環(huán)境檢查代碼比如if (Sys.which(bcftools) ) { warning(bcftools not found in PATH, ukb_gen_extract may fail) }這種防御式寫(xiě)法能幫你少走一小時(shí)彎路。6.5 坑五不要把ukb_df_na_count的計(jì)數(shù)結(jié)果直接當(dāng)缺失率ukb_df_na_count()統(tǒng)計(jì)的是R里的NA但UKB數(shù)據(jù)中很多缺失是以負(fù)編碼存在的比如-1不知道、-3拒絕回答、-7無(wú)此數(shù)據(jù)。如果你只過(guò)濾NA那些負(fù)編碼值還會(huì)留在數(shù)據(jù)里照樣污染分析。正確做法是先利用ukb_df_na_count()這類函數(shù)看分布再結(jié)合字段編碼說(shuō)明把負(fù)值統(tǒng)一轉(zhuǎn)為NA。這一步在建模前做能避免大量莫名其妙的分析異常。我在實(shí)際項(xiàng)目中用ukbtools大概一年半坦白說(shuō)它也并不是每個(gè)場(chǎng)景都必不可少——如果你只做純表型分析不碰遺傳數(shù)據(jù)用tidyverse完全夠用。但一旦涉及UKB基因型數(shù)據(jù)特別是需要在樣本層面把幾十萬(wàn)個(gè)體的表型、基因型、QC結(jié)果對(duì)齊時(shí)這套包的封裝價(jià)值就體現(xiàn)出來(lái)了。最后再分享一個(gè)小技巧處理UKB數(shù)據(jù)時(shí)盡量保持原始數(shù)據(jù)只讀、清洗數(shù)據(jù)另存的習(xí)慣哪怕ukbtools改了你的列名、篩選了樣本也千萬(wàn)別在原文件上操作否則重跑一次分析的成本會(huì)讓你后悔沒(méi)有多做一份備份。本文還有配套的精品資源點(diǎn)擊獲取