简介ukbnmr是一款专门面向R语言用户的生物信息学工具包用于处理英国生物样本库UK Biobank中的Nightingale NMR代谢组学数据解决原始数据集里生物标志物字段提取与比率计算的问题。核心功能由extract_biomarkers()和compute_nightingale_ratios()两个函数承载前者能从下载的原始数据中提取NMR生物标志物并给出简洁列名且支持基线与首次重复评估等多行测量记录后者则基于168种可下载的生物标志物计算Nightingale比率并附有nmr_info数据框方便查阅指标信息。资源压缩后仅18KB包含19个文件主要有R函数源码、Rd帮助文档、Markdown说明、DESCRIPTION配置等结构紧凑既能通过remotes直接安装也便于本地加载研究。目前已有909人学习下载特别适合熟悉R语言并从事UK Biobank代谢组学分析的研究人员参考学习。1. 拿到一批UKB导出的Nightingale NMR数据难的不是测量本身而是读进来之后的那一步ukbnmr这个R工具解决的是一个很具体、也很容易被低估的问题UKB里Nightingale NMR生物标志物数据导出来是一张巨大的宽表列名全是f.23456_0_0这种带数字的编码单看列名你根本分不清哪一列是脂蛋白颗粒浓度、哪一列是氨基酸、哪一列是炎症标志物。ukbnmr把字段映射、样本访问次数筛选和单位整理自动化让原本要手工折腾一两天的清洗工作变成几行代码的事。这篇笔记适合刚拿到UKB数据导出文件、准备做NMR代谢组学或遗传关联分析的R用户。我会按读取、清洗、质控、接下游分析的顺序拆开讲中间包含具体代码和参数说明最后记录几个我自己翻过车的地方。2. 先搞清楚Nightingale数据在UKB里到底长什么样字段区间与命名规则2.1 导出文件里那一堆f.开头的列是哪儿来的UKB的UKB数据导出不是按“一张表一个主题”给您的而是把所有测量项拼成一张超宽表每一行是一个样本每一列是一个字段。Nightingale NMR的数据分布在以23400开头的字段区间里导出来后列名形如f.23401_0_0、f.23402_0_0。这些字段并不是连续编号一一对应某个标志物中间有跳号也有语义分组所以没有一份映射表的话光靠肉眼无法判断哪列是谁。我拿到导出的第一件事是用data.table把表读进来然后用grep把NMR相关列挑出来看数量。常见做法是这样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相关列数量:, length(nmr_cols), \n) head(nmr_cols, 20)这里fread是读取大宽表最快的路径几百MB到几GB的文件都能扛住。grep(^f\\.234, ...)匹配以f.234开头的列名注意正则里点号要转义。如果您的导出版本字段区间有变化把前缀换成实际区间即可比如f.235。这一步不依赖任何包先确认数据规模再决定后面怎么处理。2.2 把idp编号拆开看一个列名包含三层信息ukbnmr内部把这些测量项统一叫idp列名格式是idp_字段号_访问次数_重复索引。拆开看其实很好理解组成部分示例含义前缀idp表示这是一个测量项标识符字段号23401对应UKB导出里的原始字段编号访问次数00初始评估1第一次重复访问重复索引0同一次访问内的测量重复编号一般为0这个命名规则非常关键因为UKB里部分参与者不止一次来过采集中心。如果您的分析只应包含初始访问就必须在筛选时把访问次数锁定为0如果不锁同一个人的两次测量会被当成两个独立样本下游的关联分析就会有重复样本问题。这个我在第5章会专门讲。另一个容易忽略的点是idp_字段号与Nightingale的可读名称之间不是简单对应关系。比如某几个字段号都指向“高密度脂蛋白中的胆固醇”只是亚类不同。ukbnmr包里内置了一张元数据表记录每个字段号的短名称、全称、所属类别和单位这比您自己去对照平台说明书要省很多事。2.3 宽表不是给分析用的从tidyverse视角看为什么要转置UKB导出的宽表在存储上很合理但对R里的统计建模来说却不好用。您做回归时通常需要“每行一个样本每列一个变量”的形态宽表是“每行一个样本每列一个字段编号”看似也是这个形态问题出在列名不可读、列数太多、并且很多列之间高度相关。更麻烦的是字段编号不只包含指标信息还混着访问次数直接拿宽表去merge临床表型时特别容易串列。我一般会先把宽表转成“长格式”每一行是样本×标志物的组合然后按需再转回宽矩阵。ukbnmr的raw_to_nightingale函数本质上就是在做这个转换同时把字段号映射成Nightingale可读名称、按访问次数筛选样本、把缺失编码统一处理好。下一章我会把整个调用过程展开。3. ukbnmr核心工作流从原始宽表到整洁的生物标志物矩阵3.1 安装与数据读入先说清楚依赖关系ukbnmr是一个常规R包安装方式与生态里其他包一致。建议用remotes安装因为该包不一定在CRAN上而且依赖data.table直接装能省去手动处理依赖的麻烦# 在R会话里执行 install.packages(remotes) remotes::install_github(your_account/ukbnmr)这里的your_account是占位符实际安装时以您获取到的仓库地址为准。装好后加载同时加载data.tablelibrary(ukbnmr) library(data.table)数据读入这一步仍然用fread原因很简单UKB导出的文件常常有几个GBread.csv会读得人想弃坑。注意fread会自动识别分隔符但UKB的导出文件有时是空格分隔、有时是制表符分隔稳妥起见在调用里明确指定sep。读入之后不要急着分析先做两个检查一是列名里是否确实含有f.234二是行数是否与您申请时拿到的样本数一致。不一致大概率是导出时勾选字段遗漏后面映射会缺一片变量。3.2 字段映射raw_to_nightingale做了哪三件事ukbnmr里最核心的函数是raw_to_nightingale它接收原始宽表返回一个包含多种对象的列表。这个函数实际做了三件事把f.开头的列名转成idp命名按访问次数过滤样本把字段号映射为Nightingale标志物的可读名称。# 只保留初始访问生成处理后的对象 nmr_list - raw_to_nightingale( raw raw, visit 0, idp_meta ukbnmr::nmr_metadata ) # 看一眼返回对象里有什么 names(nmr_list)参数visit 0表示只保留初始评估数据。这在我的经验里几乎是必填项因为默认行为不一定符合您的分析设计明确写出来既防自己忘也让脚本可读性更好。idp_meta是元数据表一般用包内置的那份即可如果您手里有更新版的UKB字段映射也可以自己传入同格式的数据框。为什么要有这一步而不是直接按列名筛因为字段号和Nightingale名称之间并非一一对应。比如某些标志物在平台上有多个测量变体字段号不同但名称近似又比如同一标志物可能在不同访问次数下各有独立字段。raw_to_nightingale在内部处理这种一对多关系把宽表整理成“每个标志物一列”的整洁矩阵这是手工复制列名做不到的。3.3 输出对象解读samples、biomarkers、meta与flagsraw_to_nightingale返回的列表通常包含几个固定对象样本信息、标志物矩阵、元数据和标记信息。我习惯一拿到就先看维度确认没有样本丢失# 观察各对象的维度 dim(nmr_list$samples) # 样本信息 dim(nmr_list$biomarkers) # 标志物矩阵 # 查看元数据表前几行 head(nmr_list$meta)samples里至少有eid和visit两列eid是UKB的样本主键后面所有表关联都用它。biomarkers则是一个行数等于样本数、列数等于标志物数的矩阵列名已经是可读的Nightingale名称比如VLDL_P_L这种短代码。meta记录每个标志物的短名称、全称、单位、所属类别脂蛋白、脂肪酸、氨基酸、糖酵解代谢物、酮体、炎症标志物等。flags是包帮你做的越界标记比如测量值超出平台线性范围、缺失、以及访问次数异常。这部分很多人不看但它在质控里非常有用直接用它把不合格样本筛掉比拿原始值自己设阈值要可靠。3.4 和dplyr管道一起用筛选样本、拼接临床表型整理好标志物矩阵后下一步通常是和临床表型拼接。我一般先把samples转成data.table再用merge关联协变量表library(dplyr) pheno - fread(ukb_clinical.txt, sep \t) # 保留初始访问样本关联年龄、性别等协变量 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关联能自动检查两边样本是否对得上。filter(visit 0)是双保险因为前面转换时已经只保留初始访问这里再筛一次并不会增加计算成本。拼好后的dat可以直接丢进lm、glm或者后续的GWAS工具。需要注意如果您的表型文件里有些样本在NMR数据里不存在left_join会留下NA行需要在建模前再做一次缺失过滤。4. QC与单位换算Nightingale数据里最容易出错的隐蔽区4.1 缺失编码-1、-3和空值混在一起时的处理顺序UKB导出的数据里缺失不是R的NA而是用特殊数字编码的-1通常表示“不适用”-3表示“数据缺失”空值在fread读入后会变成NA。如果您不做处理直接拿去算均值、做回归-1和-3会被当成真实测量值参与运算结果会很离谱。正确做法是先统一缺失编码再做任何统计。我一般这样处理# 把常见的缺失编码统一替换为 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再做缺失率统计。colMeans(is.na(...))可以快速算出每一列的缺失比例如果某标志物缺失率超过20%我会认真考虑是否纳入后续分析。这里不设固定阈值因为不同分析对缺失容忍度不同但至少心里要有数。4.2 单位换算表从平台原始单位到临床常用单位Nightingale NMR平台输出的单位与常规临床检验单位经常不一致。最典型的是脂蛋白颗粒浓度、胆固醇、甘油三酯这类指标平台常用mmol/L但有些导出版本或某个字段会给出mg/dL还有脂肪酸组成平台给的是摩尔百分比不是浓度。单位不统一建模系数解读会出大问题。标志物类别平台导出单位临床常用单位换算关系胆固醇相关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%无需换算直接使用我的习惯是先不改数据先把元数据表里每个标志物的单位列打印出来确认哪些需要换再批量处理。直接全部乘以一个系数是翻车重灾区因为不同指标系数不同。# 以 meta 为准按单位类型批量转换 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再决定乘哪个系数。千万不要闭着眼睛全表乘。4.3 样本级标记极端值、重复样本与visit选择单位搞定后样本级的质控也很关键。Nightingale数据里有些标志物在极端情况下会超出定量范围比如某些脂蛋白亚类在重度高脂血样本里可能测不准。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)把每个样本在所有标志物上的离群次数加起来超过20就值得进一步查原始记录。这一步不是ukbnmr自带的但我接每个项目都会多做这一道因为Nightingale数据几十个指标相关性极强如果某样本大量指标离群那大概率不是生物学变异而是测量问题。5. 处理UKB NMR数据时的五个典型翻车点与排查方法5.1 现象读进R后全是NA一个真值都没有这种情况我遇到过一次用fread读档后NMR相关列全是NA但其他临床列正常。原因出在UKB导出文件的表头里NMR字段列名带有特殊前缀或空格fread自动识别时把列名弄脏了导致列名匹配不上。解决方法是读入后立刻打印names(raw)前30个肉眼检查列名是否干净。如果列名带\ufeff之类的不可见字符用gsub把特殊字符去掉。更稳妥的做法是在fread调用里加check.names TRUE让R自动清洗列名。5.2 现象同一行样本在最终矩阵里出现两次分组回归结果对不上这个问题非常隐蔽。原因是有人做过重复访问如果您在raw_to_nightingale里没指定visit 0返回的samples会有同一个eid出现两次。后续关联分析时同一个人的两条记录被当成两个独立样本标准误被错误缩小。解决方法是转换后立刻检查重复any(duplicated(nmr_list$samples$eid))如果有重复回头把visit参数加上。另一个可能是您手动合并了两个版本的导出文件比如一次下载了初始访问、另一次下载了重复访问而您没在merge前去重。5.3 现象单位换算后某些标志物的分布出现双峰分布双峰未必是错的但如果是本来应该单峰的指标出现双峰先别急着下结论大概率是单位不统一。比如某个标志物部分样本以mmol/L为单位、部分样本以mg/dL为单位合并后就会形成两个峰。我一般先把该标志物按批次画密度图看双峰是否与下载批次一致。如果一致就是单位混了如果双峰与批次无关再考虑是否真实存在亚组差异。元数据表里单位列一定要逐项检查不要只看前几行。5.4 现象字段编号对不上旧笔记代码里硬编码失效UKB的字段编号在历次数据版本中可能调整旧笔记里写的f.23401到新版本里可能变成别的字段。我踩过的坑是按照旧教程写的列名去筛结果筛出来的列全空。现在我的做法是永远不硬编码字段号全部从nmr_list$meta里按标志物短名称取。代码里只出现Nightingale名称不出现f.开头字段号这样无论字段编号怎么变脚本都能跑通。5.5 现象数据读入没问题一跑全量建模就内存溢出UKB宽表本身几个GB读进来后R里的对象占用会翻几倍。如果这时候还同时保留原始宽表、转换后的列表、合并后的分析表内存很容易爆。我的习惯是用fread读入原始文件后完成raw_to_nightingale转换立刻rm(raw)并调用gc()释放内存。后续只保留biomarker_mat和samples两个对象。如果还不够就把biomarkers按类别拆开处理先做脂蛋白类再做其他类别而不是一次性载入全部分析。6. 进阶把ukbnmr输出接进GWAS/PheWAS流程的最后一个关键动作6.1 生成plink表型文件NMR标志物做完QC后下游最常见的用途是GWAS或PheWAS。这类工具一般要求表型文件是纯文本、每行一个样本、包含FID和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 残差化之后再跑关联而不是直接扔进去很多NMR标志物受年龄、性别、BMI影响很大直接拿原始值跑GWAS会出一堆假关联。我一般先做残差化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)这里逐列回归并取残差相当于把所有标志物统一调整到相同协变量水平。注意na.action na.exclude保证残差向量长度与输入一致否则样本错位很难查。这一步做完再把残差表输给plink出来的关联结果要干净得多。从那以后我每次拿到新的UKB导出都会强制走一遍字段映射、访问次数筛选、单位核对、缺失编码替换、离群标记、残差化这一整套流程中间任何一步偷懒后面返工的成本都更高。希望这份拆解能帮您在处理NMR数据时少走几步弯路。本文还有配套的精品资源点击获取