简介ukbnmr 是面向 UK Biobank 与 Nightingale NMR 代谢组学数据分析者的 R 工具包重点解决原始 NMR 生物标志物字段提取、列名简化、基线与重复评估数据整理以及 Nightingale 比值计算等重复处理问题适合具备一定 R 基础、需要开展代谢组学数据清洗与统计的研究人员。压缩包共 19 个文件约 18KB包含 7 个 R 源码文件、4 个 Rd 帮助文档以及 README、DESCRIPTION、NAMESPACE 和内置 nmr_info.rda 数据等源码与说明文档配套便于查看函数实现与参数含义。已有 909 人学习下载。读者可参考 extract_biomarkers() 与 compute_nightingale_ratios() 的实现获取生物标志物信息表、字段组织方式和安装使用说明用于梳理 UK Biobank NMR 数据预处理流程并在此基础上完成自定义清洗与比值计算。1. 拿到 23xxx 字段就开跑的人八成会在半年后返工ukbnmr这个 R 包解决的是一个很具体的问题把 UK Biobank 里那两百多项 Nightingale NMR 生物标志物从一堆带点位后缀的原始字段变成一张能直接进模型的干净表格。我见过太多人拿到字段就fread进来算均值、跑回归结果看着漂亮半年后被人追问「你到底剔了哪些样本、板位效应校没校正」一句话答不上来。这个包的价值不在算法有多新而在于把「抽字段 → 查质控 → 去离群 → 重算衍生指标」这条链路固化成函数调用每一步都留日志。适合两类人手上已经有队列 NMR 数据、准备做代谢物关联分析的以及要给团队交付一版口径统一数据集、必须能被人复现的。2. 装包、读数和抽取把 extract_biomarkers() 的一次调用拆开看2.1 依赖与安装别在 CRAN 主库里白等先说结论这个包不在 CRAN 主库install.packages(ukbnmr)会直接报「package is not available」。常见做法是从项目仓库取源码包本地编译安装。它的计算几乎全压在data.table上底层没装好后面几十万行乘两百多列的操作会非常难受。# 1) 先把底层依赖装好ukbnmr 的抽取和质控都是 data.table 引用语义 install.packages(data.table) # 2) 从手上的源码包装。别写死版本号换机器时用通配符匹配更稳 pkg - list.files(pattern ^ukbnmr_.*\\.tar\\.gz$, full.names TRUE) stopifnot(length(pkg) 1) # 目录里有多个版本会装错先卡一下 install.packages(pkg, repos NULL, type source) # 3) 装完立刻验证能加载、能查到版本两个都得过 library(ukbnmr) packageDescription(ukbnmr)$Version这里有两个细节值得说。一是repos NULL, type source本地 tar.gz 必须走源码编译不能用二进制路径否则在 Linux 上会提示无法安装。二是stopifnot(length(pkg) 1)我踩过这个坑目录里同时放了新旧两个包install.packages按字母序装了旧的那个跑了两天才发现某个函数的参数对不上文档。版本号建议单独记进sessionInfo()的输出里和结果一起归档这是后面出问题唯一能查的证据。如果机器上没有编译环境先确认R CMD config CC能返回编译器路径返回空就说明缺 Rtools 或 build-essential这时候报的错通常是一堆看不懂的.c文件找不到头文件别去改包源码先把工具链补齐。2.2 读数据这一步就能决定后面所有结论抽取函数本身不难难在喂给它的东西对不对。整张主数据集动辄几十 GB全量读进来再挑列内存直接爆。我一般先用nrows 0只读表头把 NMR 字段筛出来再带着select二次读。library(data.table) library(ukbnmr) # 只读表头拿到全部列名这一步不占内存 all_cols - names(fread(ukb_main.csv, nrows 0)) # NMR 字段集中在 23xxx 区段且带 -0.0 这类评估点位后缀 nmr_cols - grep(^23[0-9]{3}-, all_cols, value TRUE) cat(命中 NMR 字段数, length(nmr_cols), \n) # 只把 eid 和 NMR 字段读进来 dt - fread(ukb_main.csv, select c(eid, nmr_cols)) # 抽取把字段名翻译成指标名统一类型和缺失编码 nmr - extract_biomarkers(dt)select参数让fread在读取阶段就跳过无关列几十 GB 的表能压到几 GB 以内。eid是样本主键必须带上后面所有质控日志、样本清单都靠它对回去。这里有个隐性问题同一批人可能有多次评估列名后缀是-0.0、-1.0、-2.0。上面的grep会把所有评估的字段都捞进来如果数据里确实存在重复评估extract_biomarkers()拿到的就是宽表套宽表结果列数翻倍。稳妥的做法是先用grep(^23[0-9]{3}-0\\.0$, all_cols, value TRUE)只取首次评估重复评估另存一份做纵向分析。另外 UK Biobank 的缺失编码用 -1、-3 表示「不知道」「不愿回答」这些值在数值列里会被当成合法数字读进来抽取阶段通常会被转成NA但保险起见读完先看一眼最小值# 有 -1/-3 混进来说明缺失编码没处理干净 summary(dt[[nmr_cols[1]]])2.3 extract_biomarkers() 的返回对象五个槽位逐个点名这个函数返回的不是一张表是一个 list。很多人拿到之后直接$biomarkers取完就走把另外几个槽位扔了等于把质控信息丢了一半。str(nmr, max.level 1) # 每个槽位的行列数一次看清样本有没有被悄悄丢掉这里最直观 data.table( slot names(nmr), nrow sapply(nmr, function(x) if (is.data.frame(x)) nrow(x) else NA_integer_), ncol sapply(nmr, function(x) if (is.data.frame(x)) ncol(x) else NA_integer_) )按经验五个槽位大致是这么分工的槽位内容用途biomarkers一行一个样本一列一个指标直接进模型的主表quality_metrics样本 × 指标的测量质量信息判断哪些值不可信、要不要置 NAsample_quality样本级的质量提示标记溶血、放置过久这类问题technical_summaries平台侧的汇总体检字段判断整次测量是否技术失败log每一步操作的记录出问题时的黑匣子具体列名以你装到的版本为准用names(nmr$biomarkers)[1:10]扫一眼再往下写代码别照抄别人的脚本。我吃过一次亏不同版本的指标列命名大小写和空格处理不一致脚本里写死Total Cholesterol换台机器就变成Total cholesterolNULL一路传下去最后整列全是NA还以为是数据缺失。log槽位特别值得单独看一次。它记录的是每次调用的动作和影响行数后面每做一步质控都会往里追加一条。养成print(nmr$log)的习惯比在代码里到处插cat()靠谱得多。3. 三层质控技术失败、低质量样本与离群值到底先删哪个质控的顺序是有讲究的顺序错了剔除记录就对不上。我的固定顺序是技术失败 → 样本质量 → 离群值。前两步是删行第三步通常是把异常值置为NA而不是删行性质完全不同。# 三步分开写每步都留一个对象方便后面生成样本清单 step1 - remove_technical_failures(nmr, data dt, verbose TRUE) step2 - remove_low_quality_samples(step1, verbose TRUE) step3 - remove_biomarker_outliers(step2, outlier_threshold 5, verbose TRUE) # 每一步删了多少全在 log 里 print(step3$log)3.1 技术失败样本先看平台侧的体检报告remove_technical_failures()有个容易忽略的点它需要你再传一次data。原因是判断技术失败的字段不在 23xxx 区段里而在主数据集的其它位置抽取阶段没被带出来。这个参数漏了会直接报错报错信息还算清楚但第一次用的人经常愣一下——「我不是刚传过数据吗」。判定依据是平台侧的汇总信息比如谱仪状态异常、信噪比不达标、样本量不足这类。这些属于「整次测量就没测好」不是指标层面的问题所以必须先做做完再谈别的。输出里technical_summaries槽位会同步更新把被剔除样本对应的行也去掉保证几个槽位的行数始终对得上。实操上建议先把技术失败的样本单独存一份eid列表。有些研究需要做「失败样本 vs 成功样本」的基线比较看有没有选择偏倚这时候原始列表就是唯一证据。3.2 低质量样本乳酸和葡萄糖这两条线最值得看第二步看的是样本层面的质量提示。NMR 测的是血浆里的小分子样本采集后如果放置时间过长红细胞继续消耗葡萄糖、产生乳酸谱图上就会出现一组特征性的偏移。所以质量提示里最值得盯的是葡萄糖偏低、乳酸偏高这类组合信号——它们往往指向同一个原因样本处理不及时。# sample_quality 通常是宽表每列一个提示标记先汇总触发次数 sq - nmr$sample_quality flag_cols - setdiff(names(sq), eid) sq[, lapply(.SD, function(x) sum(as.logical(x), na.rm TRUE)), .SDcols flag_cols]这段代码输出的是每个标记被触发的样本数。如果某个标记触发了上万例别急着全删——先想想是不是这批样本本来就来自某个采集中心全删等于把中心效应变成了选择偏倚。我的做法是触发比例低于 3% 的标记直接按提示剔除比例高的保留样本但在模型里加一个标记作为协变量同时做敏感性分析看主效应有没有翻转。这里要区分「质控」和「分析决策」。函数提供的是剔除的工具剔不剔是你的科学判断。写论文的方法部分时这两者要分开写清楚哪些是平台预设的失败判定哪些是你自己加的更严格标准。3.3 离群值与板位效应为什么删点之前要先校正第三步最容易出问题。remove_biomarker_outliers()做的事情不止是找离群点它同时处理板位效应——同一批样本在不同板、不同孔位上测孔与孔之间存在系统性的水平偏移尤其是边缘孔位蒸发和温度梯度更明显。顺序很重要先按板号和孔位分组估一个位置偏移量从原始值里减掉再在校正后的值上判断离群。如果反过来先判离群整板偏高会被判成「一整板都是异常值」一刀切下去样本量掉一大截而且掉得莫名其妙。# 校正前后分位数对比抽第一个指标看 bm - setdiff(names(step3$biomarkers), eid)[1] rbind( before quantile(step2$biomarkers[[bm]], c(.01, .5, .99), na.rm TRUE), after quantile(step3$biomarkers[[bm]], c(.01, .5, .99), na.rm TRUE) )正常情况下中位数应该几乎不动1% 和 99% 分位会被「收」进来一些——离群值被置成NA之后极端分位自然就收缩了。如果中位数动了超过 5%说明校正把系统性差异也当成了位置效应减掉这时候要么指标本身分布极度偏态要么分组粒度太细比如按板号 × 孔位 × 日期分组每组样本太少估出来的偏移就是噪声。outlier_threshold默认是 5 倍标准差。要不要放宽到 6 或收到 4取决于指标的测量精度。像脂肪酸比例这类由多项测量推导出来的指标误差会累积5 倍标准差可能偏严而像白蛋白这种直接测量的高丰度蛋白5 倍就很宽松了。我的习惯是先按默认跑一遍统计各指标被置 NA 的比例超过 1% 的指标单独看一眼分布再决定要不要为它单独设阈值。提示这一步把异常值置为NA而不是删行是因为删行会让同一个样本在不同指标上的有效 N 不一样。如果下游要做多指标联合模型最好在建模时显式报告每个指标的有效样本量而不是笼统写一个总 N。3.4 删完之后样本量怎么交代质控跑完必须有一张能放进补充材料的流程图用数字说话。我一般会生成这么一张表阶段剩余样本本步剔除剔除依据原始抽取全量—字段完整技术失败N1全量−N1平台侧判定低质量N2N1−N2样本质量提示离群处理N20置 NA5 SD 阈值注意最后一行离群处理不减少行数只减少有效值。这一点在写方法学时必须点明否则审稿人一看「最终 N」和「各指标 N」对不上就会追问。我见过有人把置 NA 的那部分样本悄悄删掉结果不同指标的分析样本量各不相同做出来的结果互相矛盾返工时全部重跑。4. 重算衍生指标为什么原厂给的比值对不上恒等式4.1 原厂衍生指标为什么会自相矛盾NMR 输出的两百多项指标里只有一部分是直接从谱图积分出来的剩下的是由这些直接测量推导出来的——各种脂蛋白亚组、脂肪酸比例、载脂蛋白比值。推导指标在原厂数据集里是分别算好的问题就出在这个「分别」上。第一个来源是舍入。直接测量值按固定小数位输出推导时用的是舍入前的内部值你拿输出值自己再算一遍末位一定对不上。第二个来源更麻烦质控步骤会把某些直接测量值置成NA或做板位校正但原来算好的衍生列不会跟着更新。于是出现这种局面——总胆固醇被校正过了HDL 胆固醇也被校正过了但「非 HDL 胆固醇」这一列还是老值两者之差对不上。第三个来源是缺失处理分母为 0 或低于检测限时原厂可能给一个特殊编码而不是NA直接拿去做除法会得到Inf。这些差异单看很小但你一旦用「非 HDL 胆固醇」做结局、用「总胆固醇 − HDL」做协变量两套口径混在一个模型里系数就没法解释了。4.2 recompute_biomarkers() 干了什么参数怎么给这一步的作用就是从直接测量列出发把全部衍生指标统一重算一遍保证内部自洽。# 在离群处理之后重算保证衍生指标用的是校正后的直接测量值 nmr_final - recompute_biomarkers(step3, verbose TRUE) # 先确认列名变化是原地覆盖还是新增列不同版本行为不一样 setdiff(names(nmr_final$biomarkers), names(step3$biomarkers))如果setdiff返回空说明是原地覆盖那就要在调之前把原始列备份一份否则回不去了。我的习惯是永远保留原始版本# 备份原厂值命名加后缀方便后续对照 original_bm - copy(step3$biomarkers) setnames(original_bm, setdiff(names(original_bm), eid), paste0(setdiff(names(original_bm), eid), _raw))参数上没什么玄学主要看verbose和是否需要保留原始列。真正需要你决策的是「重算之后用哪套」。做与文献对比的分析用原厂值更保险因为别人的结论都建立在那个口径上做内部一致性要求高的模型比如同时放总胆固醇、HDL、非 HDL 三个变量必须用重算值否则共线性会因为口径不一致被放大。4.3 用恒等式和相关性做冒烟测试重算完不要直接往下走跑几条硬性检查。恒等式是免费的验证非 HDL 胆固醇必须等于总胆固醇减 HDL 胆固醇这是定义不是近似。b - nmr_final$biomarkers # 列名匹配函数不同版本命名大小写不一别写死 find_col - function(nms, pat) { hit - grep(pat, nms, ignore.case TRUE, value TRUE) if (!length(hit)) stop(没匹配到列, pat) hit[1] } tc - find_col(names(b), ^Total Cholesterol$) hdl - find_col(names(b), ^HDL Cholesterol$) nhdl - find_col(names(b), ^Non-HDL Cholesterol$) # 恒等式误差正常应该是 0 或浮点级别 max(abs(b[[nhdl]] - (b[[tc]] - b[[hdl]])), na.rm TRUE)如果这个误差不是 0 而是一个可观的数比如 0.05 mmol/L 量级说明重算没有覆盖到那一列或者单位不一致——总胆固醇按 mmol/L 输出、某个衍生列按 mg/dL 输出这种混搭在跨批次合并的数据里真的会遇到。换算关系记牢胆固醇类 1 mmol/L ≈ 38.67 mg/dL甘油三酯 1 mmol/L ≈ 88.57 mg/dL。别拿文献里的 mg/dL 阈值直接往 mmol/L 的列上卡卡出来的患病率能差一倍。第二条检查是相关性冒烟测试。载脂蛋白 B 和 LDL 胆固醇在生理上高度相关如果这两列的 Spearman 相关系数掉到 0.8 以下基本可以判定某一列取错了或者单位错了。apob - find_col(names(b), ^Apolipoprotein B$) ldl - find_col(names(b), ^LDL Cholesterol$) cor(b[[apob]], b[[ldl]], method spearman, use pairwise.complete.obs)这个数正常应该在 0.9 上下。掉下来的常见原因有三个一是某一列被大量置了NA剩下的小样本里相关性不稳定二是单位不一致三是find_col匹配到了亚组列比如「LDL 中的胆固醇」而不是「LDL 胆固醇」名字很像但含义完全不同。第三种最阴因为代码不报错结果只是偏低一点点。5. 避坑与排查五个把结果带偏的血泪现场5.1 现象抽取后样本数比预期少了上万例原因主数据集里的缺失编码没处理干净。UK Biobank 用 -1、-3 表示「不知道」「不愿回答」这些值混在数值列里抽取阶段如果按「值不在合理生理范围内就置 NA」的规则过滤会连带着把一些合法但极端的真实值也干掉反过来如果没做范围检查-1 会被当成真实测量值参与计算把均值往下拉。解决读入之后、抽取之前先对这些编码做一次统一处理把 -1、-3、-9 这类哨兵值显式转成NA再核对样本数。写一个断言把预期样本数和实际样本数都打印出来差得多了立刻停下别往下跑。5.2 现象某个衍生指标整列是 NA且日志里没报错原因质控顺序错了。离群处理把某个直接测量值置成了NA而衍生指标是在这之前算好的或者重算时因为分母是NA直接得到NA。整列NA不会触发任何警告summary()里就是一片NA。解决固定顺序为「技术失败 → 样本质量 → 离群 → 重算」并且重算之后立刻统计每个指标的缺失率miss - sapply(nmr_final$biomarkers, function(x) mean(is.na(x))) sort(miss[miss 0.05], decreasing TRUE)缺失率超过 5% 的指标要么是质控太严要么是上游出了问题单独查。5.3 现象换台机器跑同一个脚本结果对不上原因包版本不同指标列命名规则变了大小写、空格、单位后缀脚本里写死的列名匹配不到某些步骤静默跳过。R 不会因为df[[Total Cholesterol]]返回NULL就报错。解决所有列名访问都走find_col()这类带stop()的包装函数匹配不到就炸不要用$直接取。同时把sessionInfo()和packageDescription(ukbnmr)$Version写进结果目录出问题时第一时间能定位。5.4 现象内存吃满、跑一整夜没结果原因全量读入 data.frame拷贝。nmr$biomarkers几十万行乘两百多列每次$取值、每次rbind都可能触发整表复制。加上quantile()默认不做原地计算频繁调用会把内存翻几倍。解决三件事——读数据用select只取需要的列中间对象用copy()显式控制不要隐式复制统计类操作尽量用data.table的分组语法一次算完别在循环里反复取子集。真跑不动就把样本分块每块独立跑完落盘最后合并。5.5 现象log 变换后出现 NaN 或 Inf原因某些指标存在 0 或极小值。NMR 指标里像某些脂肪酸比例、亚组浓度低于检测限时可能输出 0。直接log()得到-Inf后续任何回归都会失败。解决变换前先查最小值min(x, na.rm TRUE)小于等于 0 的指标单独处理——要么加一个基于该指标检测限的小常数再取对数要么直接保留原尺度改用稳健方法。别图省事统一加 1不同指标的量纲差几个数量级加同一个常数会彻底改变分布形状。6. 把流程钉死版本锁、样本清单与三行冒烟测试跑通一次不算数能被人复现才算。我现在的固定做法是每次分析都产出三样东西清洗后的指标表、样本清单、环境快照。library(arrow) # 1) 指标表落盘成 parquet比 CSV 小一个数量级类型也不丢 as.data.frame(nmr_final$biomarkers) | write_parquet(nmr_biomarkers.parquet) # 2) 样本清单每个样本走到哪一步、为什么被剔逐条可查 manifest - data.table(eid nmr$biomarkers$eid) manifest[, :( pass_technical eid %in% step1$biomarkers$eid, pass_quality eid %in% step2$biomarkers$eid, in_final eid %in% step3$biomarkers$eid )] manifest[, dropped_reason : fifelse(!pass_technical, technical_failure, fifelse(!pass_quality, low_quality, NA_character_))] fwrite(manifest, nmr_sample_manifest.csv) # 3) 一眼看清各环节损失 manifest[, .N, by dropped_reason]样本清单是我认为最被低估的交付物。有了它任何时候有人质疑「你为什么只剩这么点样本」直接把manifest按dropped_reason分组一拉数字摆在那里不用回忆。而且做敏感性分析时特别方便想试试「不剔低质量样本」的结果直接改pass_quality这一列重新筛不用重跑整条链路。环境快照更简单sessionInfo()输出重定向到文本文件和结果放一起。R 版本、包版本、依赖版本全在里面换机器复现时照着装。最后是我每次必跑的三行自检写成一个函数放在脚本末尾smoke_check - function(b) { stopifnot(nrow(b) 0) # 行数不为零 tc - find_col(names(b), ^Total Cholesterol$) hdl - find_col(names(b), ^HDL Cholesterol$) nhdl - find_col(names(b), ^Non-HDL Cholesterol$) stopifnot(max(abs(b[[nhdl]] - (b[[tc]] - b[[hdl]])), na.rm TRUE) 1e-6) apob - find_col(names(b), ^Apolipoprotein B$) ldl - find_col(names(b), ^LDL Cholesterol$) r - cor(b[[apob]], b[[ldl]], method spearman, use pairwise.complete.obs) stopifnot(r 0.8) # 生理相关性没崩 invisible(TRUE) } smoke_check(nmr_final$biomarkers)这三条检查花不了一秒但能挡住绝大部分「跑完了、结果不对、又不知道错在哪」的情况。恒等式管口径一致性相关性管列取错和单位错行数管上游丢样本。三条全过往下建模型才有意义。从那以后我每次接手一份 NMR 衍生数据不管是自己跑的还是一年前别人交付的都先把这三条跑一遍再谈分析。有两次就是因为恒等式差了 0.03顺藤摸瓜发现是两批数据的胆固醇单位不一致硬合并下去整篇结论都得推翻。希望帮到你。本文还有配套的精品资源点击获取
