簡介ukbnmr是一款專門面向R語言用戶的生物信息學工具包用于處理英國生物樣本庫UK Biobank中的Nightingale NMR代謝組學數(shù)據(jù)解決原始數(shù)據(jù)集里生物標志物字段提取與比率計算的問題。核心功能由extract_biomarkers()和compute_nightingale_ratios()兩個函數(shù)承載前者能從下載的原始數(shù)據(jù)中提取NMR生物標志物并給出簡潔列名且支持基線與首次重復評估等多行測量記錄后者則基于168種可下載的生物標志物計算Nightingale比率并附有nmr_info數(shù)據(jù)框方便查閱指標信息。資源壓縮后僅18KB包含19個文件主要有R函數(shù)源碼、Rd幫助文檔、Markdown說明、DESCRIPTION配置等結構緊湊既能通過remotes直接安裝也便于本地加載研究。目前已有909人學習下載特別適合熟悉R語言并從事UK Biobank代謝組學分析的研究人員參考學習。1. 拿到一批UKB導出的Nightingale NMR數(shù)據(jù)難的不是測量本身而是讀進來之后的那一步ukbnmr這個R工具解決的是一個很具體、也很容易被低估的問題UKB里Nightingale NMR生物標志物數(shù)據(jù)導出來是一張巨大的寬表列名全是f.23456_0_0這種帶數(shù)字的編碼單看列名你根本分不清哪一列是脂蛋白顆粒濃度、哪一列是氨基酸、哪一列是炎癥標志物。ukbnmr把字段映射、樣本訪問次數(shù)篩選和單位整理自動化讓原本要手工折騰一兩天的清洗工作變成幾行代碼的事。這篇筆記適合剛拿到UKB數(shù)據(jù)導出文件、準備做NMR代謝組學或遺傳關聯(lián)分析的R用戶。我會按讀取、清洗、質(zhì)控、接下游分析的順序拆開講中間包含具體代碼和參數(shù)說明最后記錄幾個我自己翻過車的地方。2. 先搞清楚Nightingale數(shù)據(jù)在UKB里到底長什么樣字段區(qū)間與命名規(guī)則2.1 導出文件里那一堆f.開頭的列是哪兒來的UKB的UKB數(shù)據(jù)導出不是按“一張表一個主題”給您的而是把所有測量項拼成一張超寬表每一行是一個樣本每一列是一個字段。Nightingale NMR的數(shù)據(jù)分布在以23400開頭的字段區(qū)間里導出來后列名形如f.23401_0_0、f.23402_0_0。這些字段并不是連續(xù)編號一一對應某個標志物中間有跳號也有語義分組所以沒有一份映射表的話光靠肉眼無法判斷哪列是誰。我拿到導出的第一件事是用data.table把表讀進來然后用grep把NMR相關列挑出來看數(shù)量。常見做法是這樣library(data.table) raw - fread(ukb_nmr_export.txt, sep \t, header TRUE) # 把所有以 f.234 開頭的列名篩出來 nmr_cols - grep(^f\\.234, names(raw), value TRUE) cat(NMR相關列數(shù)量:, length(nmr_cols), \n) head(nmr_cols, 20)這里fread是讀取大寬表最快的路徑幾百MB到幾GB的文件都能扛住。grep(^f\\.234, ...)匹配以f.234開頭的列名注意正則里點號要轉(zhuǎn)義。如果您的導出版本字段區(qū)間有變化把前綴換成實際區(qū)間即可比如f.235。這一步不依賴任何包先確認數(shù)據(jù)規(guī)模再決定后面怎么處理。2.2 把idp編號拆開看一個列名包含三層信息ukbnmr內(nèi)部把這些測量項統(tǒng)一叫idp列名格式是idp_字段號_訪問次數(shù)_重復索引。拆開看其實很好理解組成部分示例含義前綴idp表示這是一個測量項標識符字段號23401對應UKB導出里的原始字段編號訪問次數(shù)00初始評估1第一次重復訪問重復索引0同一次訪問內(nèi)的測量重復編號一般為0這個命名規(guī)則非常關鍵因為UKB里部分參與者不止一次來過采集中心。如果您的分析只應包含初始訪問就必須在篩選時把訪問次數(shù)鎖定為0如果不鎖同一個人的兩次測量會被當成兩個獨立樣本下游的關聯(lián)分析就會有重復樣本問題。這個我在第5章會專門講。另一個容易忽略的點是idp_字段號與Nightingale的可讀名稱之間不是簡單對應關系。比如某幾個字段號都指向“高密度脂蛋白中的膽固醇”只是亞類不同。ukbnmr包里內(nèi)置了一張元數(shù)據(jù)表記錄每個字段號的短名稱、全稱、所屬類別和單位這比您自己去對照平臺說明書要省很多事。2.3 寬表不是給分析用的從tidyverse視角看為什么要轉(zhuǎn)置UKB導出的寬表在存儲上很合理但對R里的統(tǒng)計建模來說卻不好用。您做回歸時通常需要“每行一個樣本每列一個變量”的形態(tài)寬表是“每行一個樣本每列一個字段編號”看似也是這個形態(tài)問題出在列名不可讀、列數(shù)太多、并且很多列之間高度相關。更麻煩的是字段編號不只包含指標信息還混著訪問次數(shù)直接拿寬表去merge臨床表型時特別容易串列。我一般會先把寬表轉(zhuǎn)成“長格式”每一行是樣本×標志物的組合然后按需再轉(zhuǎn)回寬矩陣。ukbnmr的raw_to_nightingale函數(shù)本質(zhì)上就是在做這個轉(zhuǎn)換同時把字段號映射成Nightingale可讀名稱、按訪問次數(shù)篩選樣本、把缺失編碼統(tǒng)一處理好。下一章我會把整個調(diào)用過程展開。3. ukbnmr核心工作流從原始寬表到整潔的生物標志物矩陣3.1 安裝與數(shù)據(jù)讀入先說清楚依賴關系ukbnmr是一個常規(guī)R包安裝方式與生態(tài)里其他包一致。建議用remotes安裝因為該包不一定在CRAN上而且依賴data.table直接裝能省去手動處理依賴的麻煩# 在R會話里執(zhí)行 install.packages(remotes) remotes::install_github(your_account/ukbnmr)這里的your_account是占位符實際安裝時以您獲取到的倉庫地址為準。裝好后加載同時加載data.tablelibrary(ukbnmr) library(data.table)數(shù)據(jù)讀入這一步仍然用fread原因很簡單UKB導出的文件常常有幾個GBread.csv會讀得人想棄坑。注意fread會自動識別分隔符但UKB的導出文件有時是空格分隔、有時是制表符分隔穩(wěn)妥起見在調(diào)用里明確指定sep。讀入之后不要急著分析先做兩個檢查一是列名里是否確實含有f.234二是行數(shù)是否與您申請時拿到的樣本數(shù)一致。不一致大概率是導出時勾選字段遺漏后面映射會缺一片變量。3.2 字段映射raw_to_nightingale做了哪三件事ukbnmr里最核心的函數(shù)是raw_to_nightingale它接收原始寬表返回一個包含多種對象的列表。這個函數(shù)實際做了三件事把f.開頭的列名轉(zhuǎn)成idp命名按訪問次數(shù)過濾樣本把字段號映射為Nightingale標志物的可讀名稱。# 只保留初始訪問生成處理后的對象 nmr_list - raw_to_nightingale( raw raw, visit 0, idp_meta ukbnmr::nmr_metadata ) # 看一眼返回對象里有什么 names(nmr_list)參數(shù)visit 0表示只保留初始評估數(shù)據(jù)。這在我的經(jīng)驗里幾乎是必填項因為默認行為不一定符合您的分析設計明確寫出來既防自己忘也讓腳本可讀性更好。idp_meta是元數(shù)據(jù)表一般用包內(nèi)置的那份即可如果您手里有更新版的UKB字段映射也可以自己傳入同格式的數(shù)據(jù)框。為什么要有這一步而不是直接按列名篩因為字段號和Nightingale名稱之間并非一一對應。比如某些標志物在平臺上有多個測量變體字段號不同但名稱近似又比如同一標志物可能在不同訪問次數(shù)下各有獨立字段。raw_to_nightingale在內(nèi)部處理這種一對多關系把寬表整理成“每個標志物一列”的整潔矩陣這是手工復制列名做不到的。3.3 輸出對象解讀samples、biomarkers、meta與flagsraw_to_nightingale返回的列表通常包含幾個固定對象樣本信息、標志物矩陣、元數(shù)據(jù)和標記信息。我習慣一拿到就先看維度確認沒有樣本丟失# 觀察各對象的維度 dim(nmr_list$samples) # 樣本信息 dim(nmr_list$biomarkers) # 標志物矩陣 # 查看元數(shù)據(jù)表前幾行 head(nmr_list$meta)samples里至少有eid和visit兩列eid是UKB的樣本主鍵后面所有表關聯(lián)都用它。biomarkers則是一個行數(shù)等于樣本數(shù)、列數(shù)等于標志物數(shù)的矩陣列名已經(jīng)是可讀的Nightingale名稱比如VLDL_P_L這種短代碼。meta記錄每個標志物的短名稱、全稱、單位、所屬類別脂蛋白、脂肪酸、氨基酸、糖酵解代謝物、酮體、炎癥標志物等。flags是包幫你做的越界標記比如測量值超出平臺線性范圍、缺失、以及訪問次數(shù)異常。這部分很多人不看但它在質(zhì)控里非常有用直接用它把不合格樣本篩掉比拿原始值自己設閾值要可靠。3.4 和dplyr管道一起用篩選樣本、拼接臨床表型整理好標志物矩陣后下一步通常是和臨床表型拼接。我一般先把samples轉(zhuǎn)成data.table再用merge關聯(lián)協(xié)變量表library(dplyr) pheno - fread(ukb_clinical.txt, sep \t) # 保留初始訪問樣本關聯(lián)年齡、性別等協(xié)變量 analysis_df - nmr_list$samples %% left_join(pheno, by eid) %% filter(visit 0) %% filter(!is.na(age), !is.na(sex)) # 與標志物矩陣合并成建模用寬表 dat - cbind(analysis_df, nmr_list$biomarkers)這里left_join按eid關聯(lián)能自動檢查兩邊樣本是否對得上。filter(visit 0)是雙保險因為前面轉(zhuǎn)換時已經(jīng)只保留初始訪問這里再篩一次并不會增加計算成本。拼好后的dat可以直接丟進lm、glm或者后續(xù)的GWAS工具。需要注意如果您的表型文件里有些樣本在NMR數(shù)據(jù)里不存在left_join會留下NA行需要在建模前再做一次缺失過濾。4. QC與單位換算Nightingale數(shù)據(jù)里最容易出錯的隱蔽區(qū)4.1 缺失編碼-1、-3和空值混在一起時的處理順序UKB導出的數(shù)據(jù)里缺失不是R的NA而是用特殊數(shù)字編碼的-1通常表示“不適用”-3表示“數(shù)據(jù)缺失”空值在fread讀入后會變成NA。如果您不做處理直接拿去算均值、做回歸-1和-3會被當成真實測量值參與運算結果會很離譜。正確做法是先統(tǒng)一缺失編碼再做任何統(tǒng)計。我一般這樣處理# 把常見的缺失編碼統(tǒng)一替換為 NA biomarker_mat - nmr_list$biomarkers biomarker_mat[biomarker_mat -1 | biomarker_mat -3] - NA # 檢查每個標志物的缺失率 missing_rate - colMeans(is.na(biomarker_mat)) summary(missing_rate)注意替換順序先替換-1和-3再做缺失率統(tǒng)計。colMeans(is.na(...))可以快速算出每一列的缺失比例如果某標志物缺失率超過20%我會認真考慮是否納入后續(xù)分析。這里不設固定閾值因為不同分析對缺失容忍度不同但至少心里要有數(shù)。4.2 單位換算表從平臺原始單位到臨床常用單位Nightingale NMR平臺輸出的單位與常規(guī)臨床檢驗單位經(jīng)常不一致。最典型的是脂蛋白顆粒濃度、膽固醇、甘油三酯這類指標平臺常用mmol/L但有些導出版本或某個字段會給出mg/dL還有脂肪酸組成平臺給的是摩爾百分比不是濃度。單位不統(tǒng)一建模系數(shù)解讀會出大問題。標志物類別平臺導出單位臨床常用單位換算關系膽固醇相關mmol/Lmg/dLmg/dL mmol/L × 38.67甘油三酯相關mmol/Lmg/dLmg/dL mmol/L × 88.57氨基酸mmol/Lμmol/Lμmol/L mmol/L × 1000脂肪酸組成mol%mol%無需換算直接使用我的習慣是先不改數(shù)據(jù)先把元數(shù)據(jù)表里每個標志物的單位列打印出來確認哪些需要換再批量處理。直接全部乘以一個系數(shù)是翻車重災區(qū)因為不同指標系數(shù)不同。# 以 meta 為準按單位類型批量轉(zhuǎn)換 meta - nmr_list$meta # 假設要對膽固醇類標志物換算 cholesterol_cols - tolower(meta$short_name) target - grep(chol, cholesterol_cols, value TRUE) for (col in target) { biomarker_mat[, col] - biomarker_mat[, col] * 38.67 }這里grep(chol, ...)匹配短名稱里帶chol的標志物把它們從mmol/L換算成mg/dL。實際使用時我是先在meta里確認哪些標志物單位確實是mmol/L再決定乘哪個系數(shù)。千萬不要閉著眼睛全表乘。4.3 樣本級標記極端值、重復樣本與visit選擇單位搞定后樣本級的質(zhì)控也很關鍵。Nightingale數(shù)據(jù)里有些標志物在極端情況下會超出定量范圍比如某些脂蛋白亞類在重度高脂血樣本里可能測不準。ukbnmr的flags對象一般會把這些樣本標記出來。我自己還會額外加一道IQR粗篩flag_outliers - function(x) { q - quantile(x, probs c(0.25, 0.75), na.rm TRUE) iqr_val - q[2] - q[1] x q[1] - 4 * iqr_val | x q[2] 4 * iqr_val } outlier_mat - apply(biomarker_mat, 2, flag_outliers) row_outlier_count - rowSums(outlier_mat, na.rm TRUE) # 標記出超過20個標志物離群的樣本 bad_samples - names(which(row_outlier_count 20)) bad_samples4倍IQR這種閾值比較寬松不會誤殺太多正常樣本但能抓住明顯的儀器異常或錄入錯誤。rowSums(outlier_mat)把每個樣本在所有標志物上的離群次數(shù)加起來超過20就值得進一步查原始記錄。這一步不是ukbnmr自帶的但我接每個項目都會多做這一道因為Nightingale數(shù)據(jù)幾十個指標相關性極強如果某樣本大量指標離群那大概率不是生物學變異而是測量問題。5. 處理UKB NMR數(shù)據(jù)時的五個典型翻車點與排查方法5.1 現(xiàn)象讀進R后全是NA一個真值都沒有這種情況我遇到過一次用fread讀檔后NMR相關列全是NA但其他臨床列正常。原因出在UKB導出文件的表頭里NMR字段列名帶有特殊前綴或空格fread自動識別時把列名弄臟了導致列名匹配不上。解決方法是讀入后立刻打印names(raw)前30個肉眼檢查列名是否干凈。如果列名帶\ufeff之類的不可見字符用gsub把特殊字符去掉。更穩(wěn)妥的做法是在fread調(diào)用里加check.names TRUE讓R自動清洗列名。5.2 現(xiàn)象同一行樣本在最終矩陣里出現(xiàn)兩次分組回歸結果對不上這個問題非常隱蔽。原因是有人做過重復訪問如果您在raw_to_nightingale里沒指定visit 0返回的samples會有同一個eid出現(xiàn)兩次。后續(xù)關聯(lián)分析時同一個人的兩條記錄被當成兩個獨立樣本標準誤被錯誤縮小。解決方法是轉(zhuǎn)換后立刻檢查重復any(duplicated(nmr_list$samples$eid))如果有重復回頭把visit參數(shù)加上。另一個可能是您手動合并了兩個版本的導出文件比如一次下載了初始訪問、另一次下載了重復訪問而您沒在merge前去重。5.3 現(xiàn)象單位換算后某些標志物的分布出現(xiàn)雙峰分布雙峰未必是錯的但如果是本來應該單峰的指標出現(xiàn)雙峰先別急著下結論大概率是單位不統(tǒng)一。比如某個標志物部分樣本以mmol/L為單位、部分樣本以mg/dL為單位合并后就會形成兩個峰。我一般先把該標志物按批次畫密度圖看雙峰是否與下載批次一致。如果一致就是單位混了如果雙峰與批次無關再考慮是否真實存在亞組差異。元數(shù)據(jù)表里單位列一定要逐項檢查不要只看前幾行。5.4 現(xiàn)象字段編號對不上舊筆記代碼里硬編碼失效UKB的字段編號在歷次數(shù)據(jù)版本中可能調(diào)整舊筆記里寫的f.23401到新版本里可能變成別的字段。我踩過的坑是按照舊教程寫的列名去篩結果篩出來的列全空?，F(xiàn)在我的做法是永遠不硬編碼字段號全部從nmr_list$meta里按標志物短名稱取。代碼里只出現(xiàn)Nightingale名稱不出現(xiàn)f.開頭字段號這樣無論字段編號怎么變腳本都能跑通。5.5 現(xiàn)象數(shù)據(jù)讀入沒問題一跑全量建模就內(nèi)存溢出UKB寬表本身幾個GB讀進來后R里的對象占用會翻幾倍。如果這時候還同時保留原始寬表、轉(zhuǎn)換后的列表、合并后的分析表內(nèi)存很容易爆。我的習慣是用fread讀入原始文件后完成raw_to_nightingale轉(zhuǎn)換立刻rm(raw)并調(diào)用gc()釋放內(nèi)存。后續(xù)只保留biomarker_mat和samples兩個對象。如果還不夠就把biomarkers按類別拆開處理先做脂蛋白類再做其他類別而不是一次性載入全部分析。6. 進階把ukbnmr輸出接進GWAS/PheWAS流程的最后一個關鍵動作6.1 生成plink表型文件NMR標志物做完QC后下游最常見的用途是GWAS或PheWAS。這類工具一般要求表型文件是純文本、每行一個樣本、包含F(xiàn)ID和IID列。ukbnmr給的是eid直接寫成兩列FIDIID即可pheno_out - data.frame( FID nmr_list$samples$eid, IID nmr_list$samples$eid, biomarker_mat ) write.table( pheno_out, file nmr_pheno.txt, sep \t, row.names FALSE, quote FALSE )注意quote FALSE必須寫否則plink會被引號搞亂。sep \t也是慣例空格的版本容易在標志物名稱里混入空格時解析出錯。6.2 殘差化之后再跑關聯(lián)而不是直接扔進去很多NMR標志物受年齡、性別、BMI影響很大直接拿原始值跑GWAS會出一堆假關聯(lián)。我一般先做殘差化dat$age - as.numeric(dat$age) dat$sex - as.factor(dat$sex) resid_mat - apply(biomarker_mat, 2, function(y) { m - lm(y ~ age sex BMI, data dat, na.action na.exclude) resid(m) }) out_pheno - data.frame(FID dat$eid, IID dat$eid, resid_mat)這里逐列回歸并取殘差相當于把所有標志物統(tǒng)一調(diào)整到相同協(xié)變量水平。注意na.action na.exclude保證殘差向量長度與輸入一致否則樣本錯位很難查。這一步做完再把殘差表輸給plink出來的關聯(lián)結果要干凈得多。從那以后我每次拿到新的UKB導出都會強制走一遍字段映射、訪問次數(shù)篩選、單位核對、缺失編碼替換、離群標記、殘差化這一整套流程中間任何一步偷懶后面返工的成本都更高。希望這份拆解能幫您在處理NMR數(shù)據(jù)時少走幾步彎路。本文還有配套的精品資源點擊獲取