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