ukbnmr实战：UK Biobank NMR生物标志物清洗与质控

简介：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_) )

按经验，五个槽位大致是这么分工的：

具体列名以你装到的版本为准，用names(nmr$biomarkers)[1:10]扫一眼再往下写代码，别照抄别人的脚本。我吃过一次亏：不同版本的指标列命名大小写和空格处理不一致，脚本里写死"Total Cholesterol"，换台机器就变成"Total cholesterol"，NULL一路传下去，最后整列全是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 删完之后，样本量怎么交代

质控跑完，必须有一张能放进补充材料的流程图，用数字说话。我一般会生成这么一张表：

注意最后一行：离群处理不减少行数，只减少有效值。这一点在写方法学时必须点明，否则审稿人一看「最终 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，顺藤摸瓜发现是两批数据的胆固醇单位不一致，硬合并下去整篇结论都得推翻。希望帮到你。

本文还有配套的精品资源，点击获取