R语言代谢组学数据分析全流程：从数据清洗到通路富集-北京尧图网络

简介这份资源面向具备一定R语言基础、希望进入代谢组学领域的科研人员与生物信息学初学者围绕代谢组学数据的处理、统计分析与可视化展开帮助读者在R环境中搭建从数据导入到结果报告的分析流程。压缩包共61个文件以52个R脚本为主体辅以3个Rmd报告文件、1个README说明文档及若干辅助文件整体约277KB脚本覆盖PCA、PLS-DA、OPLS-DA、通路分析、Bootstrap与热图相似性等典型分析环节Rmd文件则便于将代码与结果整合成可复现报告。目前已有3066人学习下载说明该资料在同类主题中具有一定参考价值。读者可借助这些脚本理解多元统计建模、差异分析与网络分析的实现思路并参考目录组织方式快速定位所需模块适合作为代谢组学分析入门的实操素材与代码模板。1. 代谢组学数据到手之后为什么 R 语言是绕不开的那条路一份从质谱平台导出的代谢组学定量表动辄几百个样本、上千个特征列名是样本编号行名是代谢物中间还夹着保留时间、m/z、加合物类型和一堆 NA。用 Excel 打开滚动两屏就卡想做个 PCA发现缺失值、批次效应、量纲差异全堆在一起。这时候多数人会把数据丢进 R因为代谢组学的整套分析链路——缺失值填补、归一化、多元统计、差异筛选、通路富集——在 R 里都有成熟且可复现的实现。标题里的「R语言代谢组学数据分析」说的不是某一个包而是一条从原始峰表到生物学解释的完整工作流。它适合两类人一类是刚拿到数据、只会跑现成脚本的实验人员另一类是做过转录组、想把这套统计思路迁移到代谢组的分析人员。下面按真实落地顺序拆开讲每一步都给能直接抄的命令和参数。2. 从峰表到分析矩阵数据清洗与预处理的落地步骤拿到峰表的第一件事不是建模而是判断这份数据能不能用。代谢组学数据的脏法很有特点缺失不是随机的低丰度代谢物在部分样本里直接测不到量纲跨度大同一个表里既有 1e3 的峰面积也有 1e8 的还有一批特征只在不到 20% 的样本里出现本质是噪声。这一章把清洗链路走完产出一个可以直接进统计模型的矩阵。2.1 读入峰表并做第一轮结构检查常见做法是把平台导出的 CSV 读进来先看维度和缺失分布再决定删哪些特征。不要一上来就填补先看清楚缺失的形态。# 读入峰表第一列通常是代谢物名称其余为样本 peak - read.csv(peak_table.csv, row.names 1, check.names FALSE) # 基本信息多少特征、多少样本 dim(peak) # 每个特征的缺失比例 na_ratio - rowMeans(is.na(peak)) summary(na_ratio) # 每个样本的缺失比例用来发现异常样本 sample_na - colMeans(is.na(peak)) hist(sample_na, breaks 30, main 样本缺失比例分布) # 先剔除缺失超过 50% 的特征这是最保守的阈值 keep_feat - na_ratio 0.5 peak_f - peak[keep_feat, ]这段代码的逻辑是先量化再动手。rowMeans(is.na(peak))得到每个特征的缺失率colMeans得到每个样本的缺失率。阈值 0.5 不是金标准样本量小、代谢物本身难测时可以放宽到 0.8但要在记录里写清楚。check.names FALSE是为了防止 R 把样本名里的连字符、空格自动改写后面做样本对齐时会用到原始名。2.2 缺失值填补三种策略怎么选代谢组学里缺失值填补没有万能解选错会直接污染下游的差异分析。常见三种策略用最小值的一半填补假设低于检出限、用该特征的中位数填补、用 KNN 基于样本间相似性填补。前两种简单但会压缩方差KNN 保留结构但计算量大且对异常值敏感。library(impute) # 策略一最小值一半填补适合缺失比例低且怀疑是检出限问题 min_half - function(x) { m - min(x, na.rm TRUE) / 2 x[is.na(x)] - m x } peak_min - t(apply(peak_f, 1, min_half)) # 策略二KNN 填补k 一般取 5 到 10 peak_knn - impute.knn(as.matrix(peak_f), k 10)$datamin_half用apply按行处理因为每一行的最小值不同。KNN 的k参数控制用多少个近邻样本估计缺失值样本量小于 30 时建议 k 取 5样本量大可以到 10。填补后一定要做一次分布对比把填补前后的密度图画在一起如果填补后出现明显的尖峰说明策略不合适。我一般会保留两份矩阵后续统计用 KNN 版本敏感性分析用最小值版本两者结论一致才敢往下写。2.3 归一化与标准化把量纲拉回同一尺度代谢组学数据必须做归一化否则高丰度代谢物会主导整个 PCA。常用的是概率商归一化PQN加对数变换或者中位数归一化。PQN 的思路是先用所有样本的均值算一个参考谱再对每个样本求它与参考谱的比值中位数作为校正因子。# PQN 归一化 pqn - function(x) { ref - apply(x, 1, function(v) mean(v, na.rm TRUE)) quot - sweep(x, 1, ref, /) factor - apply(quot, 2, median, na.rm TRUE) sweep(x, 2, factor, /) } peak_pqn - pqn(peak_knn) # 对数变换加 1 防止 log(0) peak_log - log1p(peak_pqn) # 自动缩放让每个特征均值为 0、标准差为 1 peak_scaled - t(scale(t(peak_log)))pqn函数里sweep两次第一次按行除以参考谱第二次按列除以校正因子。log1p是log(1x)比直接log安全。scale默认按列所以先转置再转回来保证是对每个特征做标准化。做完这三步矩阵就可以进 PCA 和 PLS-DA 了。注意标准化之后数据不再有绝对丰度含义后续做通路富集要回到未标准化的矩阵去算 fold change。3. 多元统计与差异筛选PCA、PLS-DA 和火山图的 R 实现清洗完的矩阵要回答两个问题样本之间有没有分组结构哪些代谢物在组间有差异。PCA 用来看整体分布和异常样本PLS-DA 用来找区分组别的方向差异分析用火山图把统计显著性和变化倍数一起呈现。这一章把这三步串起来参数怎么设、图怎么读都写清楚。3.1 PCA 的样本投影与异常点识别PCA 不需要分组信息适合在拿到数据后第一时间做质量检查。用prcomp做重点看前两个主成分的方差解释率和样本的聚集情况。pca_res - prcomp(t(peak_scaled), center TRUE, scale. FALSE) # 方差解释率 var_exp - pca_res$sdev^2 / sum(pca_res$sdev^2) var_exp[1:5] # 提取得分画图 scores - as.data.frame(pca_res$x[, 1:2]) scores$group - group_info # 外部读入的分组向量 library(ggplot2) ggplot(scores, aes(PC1, PC2, color group)) geom_point(size 3, alpha 0.8) stat_ellipse(level 0.95) labs(x paste0(PC1 (, round(var_exp[1] * 100, 1), %)), y paste0(PC2 (, round(var_exp[2] * 100, 1), %))) theme_bw()prcomp的输入是样本在前、特征在后所以对peak_scaled转置。center TRUE必须开scale. FALSE是因为前面已经做过自动缩放再开一次会重复。stat_ellipse画 95% 置信椭圆如果某个样本落在所有椭圆之外先别删回去查它的缺失率和总离子流大概率是进样问题。方差解释率低于 20% 时PCA 图只能看异常点不能用来下分组结论。3.2 PLS-DA 的有监督判别与置换检验PLS-DA 是有监督方法能放大组间差异但也容易过拟合。必须配置换检验来验证模型不是碰巧分对的。用mixOmics包实现。library(mixOmics) X - t(peak_scaled) Y - as.factor(group_info) plsda_res - plsda(X, Y, ncomp 2) # 置换检验nperm 一般取 100 到 1000 perm_res - perf(plsda_res, validation Mfold, folds 5, nrepeat 10, progressBar FALSE) # 看分类错误率 perm_res$error.ratencomp 2是先用两个成分看效果如果错误率高再加到 3。perf做的是交叉验证folds 5表示五折nrepeat 10表示重复十次取平均。错误率低于 0.2 才算模型可用。置换检验的零分布和真实模型的 Q2 对比如果真实 Q2 没有明显高于零分布说明模型不可信。这一步很多人跳过直接拿 PLS-DA 得分图写结论是代谢组学里最常见的翻车点之一。3.3 差异代谢物筛选fold change 与校正 p 值的联合阈值差异分析要同时看变化倍数和统计显著性。先用 t 检验算 p 值再做 BH 校正然后按 fold change 和校正 p 值双阈值筛选。# 假设 group_info 是二分类1 为对照组2 为处理组 idx1 - which(group_info 1) idx2 - which(group_info 2) pvals - apply(peak_log, 1, function(x) { t.test(x[idx1], x[idx2])$p.value }) padj - p.adjust(pvals, method BH) # fold change 用未标准化的数据算保留生物学含义 fc - rowMeans(peak_pqn[, idx2]) / rowMeans(peak_pqn[, idx1]) log2fc - log2(fc) # 双阈值筛选 sig - padj 0.05 abs(log2fc) 1 sum(sig)t.test默认不假设方差齐性代谢组学样本量小时更稳妥。p.adjust用 BH 控制 FDR比 Bonferroni 温和。fold change 一定要回到peak_pqn算因为peak_scaled已经中心化均值的比值没有意义。阈值padj 0.05和abs(log2fc) 1是常用组合样本量小于 6 时可以放宽到padj 0.1但要在图注里说明。火山图用ggplot2画横轴 log2fc、纵轴 -log10(padj)把显著的点标红。4. 通路富集与结果解释从代谢物列表到生物学故事差异代谢物列表本身不说明问题要映射到通路才能解释。代谢组学的通路富集和转录组不同因为代谢物是通路的底物和产物映射逻辑要用 KEGG 或 Reactome 的代谢物-通路关系。这一章讲怎么把差异列表变成可解释的通路结果。4.1 代谢物 ID 转换从名称到 KEGG 编号富集分析要求统一的 ID 体系KEGG 编号最通用。但峰表里的代谢物名往往是俗名需要先转换。library(KEGGREST) # 用 KEGG 的化合物数据库做名称匹配 metab_names - rownames(peak_log)[sig] kegg_ids - sapply(metab_names, function(nm) { res - tryCatch(keggFind(compound, nm), error function(e) NULL) if (is.null(res)) return(NA) names(res)[1] })keggFind按名称模糊匹配返回可能多个结果取第一个。匹配不上的代谢物不要硬凑单独列出来在附录里说明。常见问题是同义词比如「葡萄糖」和「D-葡萄糖」在 KEGG 里是同一个编号但名称匹配可能失败这时候要手动补一张映射表。ID 转换的命中率一般 60% 到 80%低于 50% 说明峰表的注释质量有问题要回去检查平台给的注释。4.2 基于超几何分布的代谢物集富集拿到 KEGG 编号后用超几何分布检验哪些通路富集了差异代谢物。这一步和 GO 富集逻辑一样只是背景集换成所有检测到的代谢物。# 背景集所有检测到的代谢物对应的 KEGG 编号 bg_ids - unique(na.omit(all_kegg_ids)) # 差异集 diff_ids - unique(na.omit(kegg_ids)) # 对每个通路算超几何检验 pathway_hits - sapply(pathway_list, function(pw) { length(intersect(diff_ids, pw)) }) pathway_size - sapply(pathway_list, length) N - length(bg_ids) n - length(diff_ids) pvals - phyper(pathway_hits - 1, pathway_size, N - pathway_size, n, lower.tail FALSE) padj - p.adjust(pvals, method BH)phyper的参数依次是命中数减一、通路大小、背景集减去通路大小、差异集大小。lower.tail FALSE算的是大于等于当前命中数的概率。pathway_list是从 KEGG 拉取的代谢物-通路关系可以用KEGGREST的keggLink构建。富集结果按 padj 排序取前 20 条画气泡图横轴富集因子、纵轴通路名、点大小表示命中数、颜色表示 padj。4.3 结果解释的边界相关不等于因果富集到某个通路只能说明这个通路里的代谢物在组间有差异不能直接推出通路被激活或抑制。代谢物是通路的中间产物一个代谢物升高可能是上游酶活性增强也可能是下游消耗减少。写结论时要用「与该通路相关的代谢物发生显著变化」这种表述不要写「该通路被激活」。如果要做因果推断需要结合酶活数据或同位素标记实验单靠代谢组学数据做不到。这一条是审稿人最常追问的点提前在讨论里写清楚能省很多返修时间。5. 避坑与排查代谢组学 R 分析里最容易翻车的五件事这一章列的五个问题每一个我都实际遇到过有的还反复出现。按「现象 → 原因 → 解决」写方便对照排查。5.1 样本名对不上导致分组错位现象PCA 图里分组完全混乱但检查数据发现组间确实有差异。原因读入分组表时样本名和峰表的列名顺序不一致R 按位置匹配而不是按名称匹配。解决永远用名称匹配不要用位置。读入分组表后先做all(colnames(peak) group_df$sample)检查返回 FALSE 就停下来用match重新对齐。group_df - read.csv(group.csv, stringsAsFactors FALSE) group_info - group_df$group[match(colnames(peak), group_df$sample)] stopifnot(!any(is.na(group_info)))match返回的是分组表里对应样本的位置stopifnot保证没有样本匹配失败。这一步花十秒能省掉后面几天的返工。5.2 缺失值填补后差异分析全变显著现象填补前差异代谢物只有几十个填补后变成几百个。原因用最小值一半填补时如果某个特征在对照组全缺失、处理组全检出填补后对照组全是同一个很小的值t 检验必然显著。解决填补前先过滤掉「在一组内缺失率超过 80%」的特征这类特征没有统计意义。# 按组算缺失率两组都要求缺失率低于 80% na_by_group - sapply(split(seq_len(ncol(peak_f)), group_info), function(idx) { rowMeans(is.na(peak_f[, idx])) }) keep - rowMeans(na_by_group 0.8) 1 peak_f - peak_f[keep, ]split按分组把列索引分开分别算缺失率要求两组都低于 0.8 才保留。这个过滤要在填补之前做顺序反了就没意义。5.3 PLS-DA 过拟合置换检验通不过还硬写结论现象PLS-DA 得分图两组分得很开但置换检验的 Q2 和零分布重叠。原因样本量小、特征多模型记住了噪声。解决先降维再建模用 PCA 前几个主成分或者先做差异筛选把特征数降到样本数的十分之一以下再跑 PLS-DA。如果置换检验还是通不过就放弃 PLS-DA改用单变量方法加多重校正。5.4 归一化方法选错导致假阳性现象换了归一化方法后差异代谢物列表几乎不重叠。原因不同归一化对批次效应和生物学变异的处理方式不同PQN 保留生物学变异中位数归一化会把它压掉。解决在方法部分写清楚用的哪种归一化并做敏感性分析至少用两种方法跑一遍结论一致才写进正文。如果两种方法结论矛盾说明数据里的批次效应没处理好要回去查质控样本。5.5 通路富集背景集选错现象富集结果里全是代谢物数量多的通路比如「代谢通路」这种大通路排第一。原因背景集用了整个 KEGG 数据库而不是本次实验检测到的代谢物。解决背景集必须限定在本次检测到的代谢物对应的 KEGG 编号里这样富集因子才有意义。phyper里的N要用length(bg_ids)不是 KEGG 的总化合物数。6. 把分析固化成可复现的 R 脚本参数管理与结果归档走到这一步分析已经跑通了但下次换一批数据还要重来一遍。真正省时间的做法是把整条链路写成一个可配置的脚本参数集中管理结果自动归档。我一般会把项目组织成四个文件config.R存所有阈值和路径01_preprocess.R做清洗02_stats.R做统计03_enrich.R做富集。每个脚本开头source(config.R)这样改阈值只改一个地方。参数管理的关键是把所有硬编码的数字抽出来。比如缺失率阈值、KNN 的 k、padj 阈值、log2fc 阈值全部放进config.R的列表里。# config.R params - list( na_feat_threshold 0.5, na_group_threshold 0.8, knn_k 10, padj_cutoff 0.05, log2fc_cutoff 1, ncomp 2, perm_folds 5, perm_repeat 10 )脚本里用params$padj_cutoff引用不要写0.05。这样换项目时只改config.R其他脚本不动。结果归档用saveRDS存中间对象用write.csv存最终表格图用ggsave统一存成 PDF 和 PNG 两份。PDF 给投稿用PNG 给快速查看用。验证分析是否可复现最直接的办法是清空环境重跑一遍。rm(list ls())之后source(01_preprocess.R)看能不能不报错地跑完。如果报错说某个变量找不到说明脚本之间有隐藏依赖要把依赖显式写进source顺序里。我还会在脚本末尾加一段sessionInfo()输出记录 R 版本和包版本投稿时附在补充材料里审稿人问起来直接给。最后一个习惯每次分析完把config.R里的参数和对应的结论写进一个README.md包括数据来源、样本数、用了哪种归一化、最终筛出多少个差异代谢物。隔三个月再回来看没有这份记录根本想不起来当时为什么把阈值设成 0.8 而不是 0.5。代谢组学分析里后悔药就是这份记录。希望帮到你。本文还有配套的精品资源点击获取

R语言代谢组学数据分析全流程：从数据清洗到通路富集

相关推荐

IEEE 754浮点数二进制解析：从内存位操作到C语言实现

从零了解LoRA微调：qwen-image-edit-plus-nsfw-lora入门前的完整知识指南

洗衣店服务管理系统开发实战：从需求拆解到答辩避坑

更多资讯

ABAP开发者必看：Eclipse ADT运行慢从14秒到0.8秒的优化全攻略

网络通讯基础一页纸：从IP、端口到DNS与TCP，看懂数据链路

云快充协议V1.8接入实战：从TCP长连接到订单对账的避坑指南

AI提示工程实战体系：结构化建模与可复现提示设计

2024全国铁路矢量数据：GIS空间分析与交通规划实战指南

一个 ComfyUI 适配版为何能冲上热度榜前列：节点生态才是放大器

今日更新

Cephadm实战指南：Linux系统工程师的生产级Ceph集群部署与排障

不想再被 Token 额度卡脖子？把免费 API 接入 WorkBuddy 的完整教程

ITK中SpatialObject几何体系解析：坐标系、变换与管道集成

本周精选

交通拥堵预测毕设：基于Python+LSTM的车辆流量预测系统全流程解析

MATLAB/Simulink搭建10机39节点电力系统暂态稳定仿真实战

无人机航拍三维重建全流程：从SfM到网格生成的避坑指南

本月热门

我发现了一个新思路：用 Remotion + Claude Code 像写代码一样自动化生成短视频

Windows下 Codex 中 Chrome 和 Computer Use 插件不可用问题排查及解决参考方式：TaoToken 统一 Key 配置与验证

2026 大模型集体涨价：用 Python 做企业 Token 成本测算与选型避坑（附配置）

读完文章，想聊聊您的网站？