简介这份资源面向具备一定R语言基础、希望进入代谢组学领域的科研人员与生物信息学初学者系统整理了从数据导入到结果报告的全流程分析思路。压缩包共61个文件以52个R脚本为核心辅以3个Rmd报告文档、1个README说明及若干辅助文件整体约277KB体积轻便却覆盖PCA、PLS-DA、OPLS-DA、通路分析、bootstrap与热图相似性等典型分析模块。内容围绕数据预处理、多元统计降维、差异代谢物筛选、共表达网络构建及KEGG、HMDB注释等关键环节展开并涉及MetaboAnalystR等专用包的调用方式。目前已有3066人学习下载适合作为代谢组学课题的脚本模板与流程参考帮助读者快速搭建可复用的分析框架理解各步骤的代码实现与结果解读逻辑。1. 代谢组学数据到手之后为什么 R 语言是绕不开的那条路拿到一份代谢组学数据很多人的第一反应是打开 Excel 做差异倍数和 T 检验几十个样本还能忍一旦样本量上百、代谢物几千个手工操作立刻变成灾难。R 语言代谢组学数据分析之所以成为这个领域的主流选择核心原因在于它把数据清洗、缺失值填补、归一化、多元统计、差异筛选、通路富集这一整条链路都变成了可复现的脚本。你不需要记住每一步点了哪个菜单只需要保存一份能重跑的代码。这套流程适合三类人刚拿到 LC-MS 或 GC-MS 峰表的研究生、需要批量处理多组学数据的分析人员、以及想把分析过程标准化沉淀下来的实验室成员。标题里的“R语言代谢组学数据分析”本质上是一套从原始峰表到生物学解释的完整工作流而不是某一个单独的包或函数。下面我按实际做项目的顺序把这条链路拆开讲清楚。2. 从峰表到分析就绪矩阵数据清洗与预处理的完整链路代谢组学分析翻车最多的地方不在统计模型而在数据进入模型之前。原始峰表通常包含大量缺失值、批次效应、量纲差异直接丢进 PCA 或 PLS-DA出来的图看着漂亮但结论站不住。这一章把预处理拆成可执行的步骤每一步都给出 R 代码和参数解释。2.1 读入峰表与样本元数据的对齐常见做法是把峰表存成 CSV行是代谢物含 m/z、保留时间、名称列是样本另一份元数据表记录样本分组、批次、性别等。两份表靠样本名对齐这一步出错后面全错。# 读入峰表第一列作为行名代谢物标识 peak_table - read.csv(peak_table.csv, row.names 1, check.names FALSE) # 读入样本元数据 sample_meta - read.csv(sample_meta.csv, row.names 1, check.names FALSE) # 检查样本名是否完全一致顺序也要对齐 stopifnot(all(colnames(peak_table) rownames(sample_meta))) # 如果顺序不一致按元数据顺序重排峰表列 peak_table - peak_table[, rownames(sample_meta), drop FALSE]row.names 1表示把第一列当作行名check.names FALSE防止 R 把代谢物名称里的特殊字符自动改写。stopifnot是一道保险样本名对不上直接报错比后面发现分组错位再回头查要省事得多。元数据里分组列建议用因子类型并显式指定水平顺序否则后续画图时对照组和实验组的排列可能和你预期相反。2.2 缺失值处理删、补还是留缺失值在代谢组学里分两类真缺失低于检出限和随机缺失仪器波动。处理策略取决于缺失比例。我一般按代谢物缺失率过滤缺失超过 20% 的代谢物直接剔除剩余缺失用最小值一半或 KNN 填补。# 计算每个代谢物的缺失比例 na_ratio - rowSums(is.na(peak_table)) / ncol(peak_table) # 剔除缺失率大于 20% 的代谢物 peak_table - peak_table[na_ratio 0.2, ] # 剩余缺失用每行最小值的一半填补 min_half - apply(peak_table, 1, function(x) min(x, na.rm TRUE) / 2) for (i in seq_len(nrow(peak_table))) { peak_table[i, is.na(peak_table[i, ])] - min_half[i] }20% 这个阈值不是铁律样本量小的时候可以放宽到 30%但要在方法里写清楚。用最小值一半填补是代谢组学里比较保守的做法避免引入过大的假信号。如果数据来自多批次填补前最好先做批次校正否则填补值会带着批次偏差进入下游。2.3 归一化与转换让样本之间可比峰强度受样本上样量、仪器灵敏度影响必须归一化。常见方法有总峰面积归一化、内标归一化、概率商归一化PQN。PQN 对整体差异的校正更稳健但计算稍复杂。# 总峰面积归一化 total_area - colSums(peak_table) peak_table_norm - sweep(peak_table, 2, total_area, /) * median(total_area) # 对数转换压缩量纲差异 peak_table_log - log2(peak_table_norm 1) # 帕累托缩放常用于 PCA 前处理 pareto_scale - function(x) { (x - mean(x)) / sqrt(sd(x)) } peak_table_scaled - t(apply(peak_table_log, 1, pareto_scale))sweep按列除以总峰面积再乘中位数让所有样本回到同一量级。log2(x1)里的加一是为了防止零值取对数出错。帕累托缩放除以标准差的平方根相比自动缩放除以标准差不会过度放大低丰度代谢物的噪声这是代谢组学里比较常用的选择。归一化方法没有绝对优劣关键是全流程统一并且在方法部分交代清楚。3. 多元统计与差异筛选PCA、PLS-DA 和火山图怎么配合用预处理做完下一步是看整体结构、找差异代谢物。PCA 用来看样本聚集和离群PLS-DA 用来放大组间差异差异筛选则结合统计检验和倍数变化。这三者不是替代关系是递进关系。3.1 PCA 与离群样本判断PCA 不需要分组信息先跑一遍能快速发现离群样本和批次效应。library(ggplot2) # 转置让行是样本、列是代谢物 pca_input - t(peak_table_scaled) # 跑 PCA pca_res - prcomp(pca_input, center TRUE, scale. FALSE) # 提取得分和前两个主成分的方差解释比例 scores - as.data.frame(pca_res$x[, 1:2]) scores$group - sample_meta$group var_exp - round(summary(pca_res)$importance[2, 1:2] * 100, 1) ggplot(scores, aes(PC1, PC2, color group)) geom_point(size 3) labs(x paste0(PC1 (, var_exp[1], %)), y paste0(PC2 (, var_exp[2], %))) theme_bw()scale. FALSE是因为前面已经做过帕累托缩放再标准化会重复处理。看 PCA 图时重点关注两件事组内样本是否聚在一起有没有样本跑到置信椭圆外面。如果有离群样本先查原始峰表里该样本的总峰面积和缺失率确认是技术问题还是生物学差异不要直接删。3.2 PLS-DA 与置换检验PLS-DA 是有监督方法能最大化组间分离但也容易过拟合。必须做置换检验来评估模型是否真的有效。library(mixOmics) # 构建 PLS-DA 模型ncomp 先设一个较大值 plsda_res - plsda(pca_input, sample_meta$group, ncomp 3) # 置换检验200 次 perm_res - perf(plsda_res, validation Mfold, folds 5, nrepeat 50, progressBar FALSE) # 查看分类错误率 plot(perm_res)ncomp是主成分数不是越多越好通常看分类错误率曲线选拐点。perf做交叉验证评估模型稳定性nrepeat越大结果越稳但耗时越长。如果置换检验的 Q2 截距明显小于零说明模型没有过拟合组间差异可信。这一步很多人跳过直接拿 PLS-DA 得分图说事审稿人一问模型验证就露馅。3.3 差异代谢物筛选VIP 加 Fold Change 加 P 值差异代谢物通常要同时满足三个条件VIP 大于 1、倍数变化超过阈值、统计检验显著。# 提取 VIP 值 vip_values - vip(plsda_res) # 对每个代谢物做 t 检验 p_values - apply(peak_table_log, 1, function(x) { t.test(x ~ sample_meta$group)$p.value }) # 计算倍数变化以两组为例 group_means - t(apply(peak_table_log, 1, function(x) { tapply(x, sample_meta$group, mean) })) fold_change - group_means[, 2] - group_means[, 1] # log2 尺度上的差值 # 汇总结果 diff_res - data.frame( metabolite rownames(peak_table_log), VIP vip_values, log2FC fold_change, p_value p_values ) # 筛选条件 diff_sig - subset(diff_res, VIP 1 abs(log2FC) 1 p_value 0.05)VIP 衡量每个代谢物对组间分离的贡献大于 1 是常用阈值。log2FC在 log2 尺度上直接相减就是倍数变化的对数值绝对值大于 1 对应原始倍数变化超过 2 倍。P 值这里没有做多重检验校正实际分析中建议加一列 FDR用p.adjust(p_values, method BH)筛选时用 FDR 小于 0.05 更稳妥。三个条件同时用能有效降低假阳性但也会漏掉一些变化幅度小却稳定的代谢物阈值要根据课题灵活调整。4. 通路富集与结果解释从代谢物列表到生物学故事差异代谢物列表本身没有生物学意义必须映射到通路才能讲出故事。富集分析的核心是看某些通路里的差异代谢物是否比随机期望更多。4.1 基于 KEGG 的富集分析流程常见做法是用代谢物名称或 KEGG ID 做超几何检验。R 里可以用clusterProfiler配合 KEGG 注释但代谢组学的 KEGG 富集和转录组略有不同需要自己整理代谢物到通路的映射。library(clusterProfiler) # 差异代谢物的 KEGG ID 列表 diff_kegg - c(C00031, C00186, C00041) # 示例 ID # 背景为所有检测到的代谢物的 KEGG ID background_kegg - c(C00031, C00186, C00041, C00064, C00065) # 富集分析 enrich_res - enrichKEGG(gene diff_kegg, universe background_kegg, organism hsa, pvalueCutoff 0.05) # 查看结果 head(as.data.frame(enrich_res))universe参数必须填所有检测到的代谢物对应的 KEGG ID不能省略否则背景不对富集结果会虚高。organism hsa是人类其他物种要换成对应的 KEGG 代码。富集结果里的Count是通路中命中的差异代谢物数量pvalue是超几何检验的原始 P 值建议同时看p.adjust。4.2 富集结果的可视化与解读边界富集气泡图是最常用的展示方式横轴是富集因子纵轴是通路名点的大小表示命中数量颜色表示显著性。library(ggplot2) enrich_df - as.data.frame(enrich_res) enrich_df$richFactor - enrich_df$Count / as.numeric(sub(/.*, , enrich_df$BgRatio)) ggplot(enrich_df, aes(x richFactor, y reorder(Description, richFactor))) geom_point(aes(size Count, color p.adjust)) scale_color_gradient(low red, high blue) labs(x Rich Factor, y NULL, color p.adjust, size Count) theme_bw()richFactor是命中数量除以背景中该通路的代谢物总数值越大说明该通路富集程度越高。解读时注意两点一是富集显著不代表因果关系代谢物变化可能是下游效应二是 KEGG 通路覆盖的代谢物有限很多脂质和未知代谢物没有注释富集结果只能反映有注释的那部分。我一般会把富集结果和差异代谢物的热图放在一起看确认通路里的代谢物变化方向是否一致。5. 避坑与排查代谢组学 R 分析里最容易翻车的五个地方这一章记录的是实际项目里反复出现的问题每一条都按现象、原因、解决来写。5.1 样本名对不上导致分组错位现象PCA 图里组间完全分开但置换检验 Q2 截距为正模型明显过拟合。原因峰表列名和元数据行名顺序不一致R 按位置匹配分组信息整体错位。解决读入数据后立刻用stopifnot(all(colnames(peak_table) rownames(sample_meta)))检查不一致就按元数据顺序重排不要靠肉眼核对。5.2 缺失值填补引入假信号现象差异代谢物列表里出现大量低丰度代谢物且这些代谢物在原始数据里缺失率很高。原因用最小值一半填补后缺失样本被赋予了一个固定的小值如果该代谢物在某一组里缺失集中填补后会人为制造组间差异。解决填补前先按缺失率过滤缺失率超过 20% 的代谢物直接剔除填补后做一次敏感性分析对比填补前后差异结果是否一致。5.3 归一化方法不统一导致批次效应残留现象PCA 图里样本按批次聚集而不是按分组聚集。原因只做了总峰面积归一化没有做批次校正或者不同批次用了不同的归一化参数。解决多批次数据建议用sva或limma的removeBatchEffect做校正校正后再跑 PCA 确认批次效应是否消除。归一化参数要全流程统一不要中途换方法。5.4 PLS-DA 过拟合没有做置换检验现象PLS-DA 得分图组间分离完美但用新样本预测时准确率很低。原因主成分数选得太多模型把噪声也拟合进去了。解决用perf做交叉验证选主成分数用permutest做置换检验Q2 截距小于零才认为模型可靠。主成分数一般不超过 3样本量小的时候更要保守。5.5 富集分析背景基因集设置错误现象富集结果里出现大量明显不相关的通路且 P 值极小。原因universe参数没有设置或设置成了全部 KEGG 代谢物而不是本次实验检测到的代谢物。解决背景必须用本次实验检测到的所有代谢物对应的 KEGG ID不能图省事用默认值。富集分析前先确认有多少代谢物有 KEGG 注释注释率太低的话富集结果参考价值有限。6. 把分析流程封装成可复现脚本的几个实用技巧走到这里单次分析已经能跑通了。但代谢组学项目往往不止一批数据下次换一组样本如果还要从头改代码效率太低。我的习惯是把整个流程拆成三个文件一个负责读入和预处理一个负责统计和差异筛选一个负责富集和可视化。每个文件开头用source调用上一个文件的输出中间结果存成 RDS 格式方便单独调试某一步而不用重跑全流程。# 保存中间结果方便后续单独加载 saveRDS(peak_table_log, file step1_normalized.rds) saveRDS(diff_sig, file step2_diff_metabolites.rds) # 下次调试直接从某一步开始 peak_table_log - readRDS(step1_normalized.rds)参数不要写死在代码里统一放到一个config.R文件比如缺失率阈值、VIP 阈值、P 值阈值、分组名称。换项目时只改配置文件主流程代码不动。这样做的另一个好处是方法部分写起来有据可查每个阈值都能追溯到配置文件里的具体数值。还有一个容易被忽略的点是随机种子。PLS-DA 的交叉验证、置换检验、KNN 填补都涉及随机过程不设种子的话每次跑结果略有差异别人复现时会对不上。在脚本开头加一行set.seed(123)所有随机过程就固定了。种子数字本身不重要重要的是全流程统一并且在方法里注明。最后说一个我踩过的坑不要等到分析全部做完再写注释和整理输出。代谢组学分析步骤多隔一周回头看自己的代码经常想不起来某个筛选条件为什么设成 0.05 而不是 0.01。每跑完一步就把参数含义、选择理由、输出文件说明写在代码块上方比事后补要省力得多。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/13 4:24:00 · 阅读 1.2k R语言代谢组学数据分析全流程：从数据清洗到通路富集的完整指南
简介这份资源面向具备一定R语言基础、希望进入代谢组学领域的科研人员与生物信息学初学者系统整理了从数据导入到结果报告的全流程分析思路。压缩包共61个文件以52个R脚本为核心辅以3个Rmd报告文档、1个README说明及若干辅助文件整体约277KB体积轻便却覆盖PCA、PLS-DA、OPLS-DA、通路分析、bootstrap与热图相似性等典型分析模块。内容围绕数据预处理、多元统计降维、差异代谢物筛选、共表达网络构建及KEGG、HMDB注释等关键环节展开并涉及MetaboAnalystR等专用包的调用方式。目前已有3066人学习下载适合作为代谢组学课题的脚本模板与流程参考帮助读者快速搭建可复用的分析框架理解各步骤的代码实现与结果解读逻辑。1. 代谢组学数据到手之后为什么 R 语言是绕不开的那条路拿到一份代谢组学数据很多人的第一反应是打开 Excel 做差异倍数和 T 检验几十个样本还能忍一旦样本量上百、代谢物几千个手工操作立刻变成灾难。R 语言代谢组学数据分析之所以成为这个领域的主流选择核心原因在于它把数据清洗、缺失值填补、归一化、多元统计、差异筛选、通路富集这一整条链路都变成了可复现的脚本。你不需要记住每一步点了哪个菜单只需要保存一份能重跑的代码。这套流程适合三类人刚拿到 LC-MS 或 GC-MS 峰表的研究生、需要批量处理多组学数据的分析人员、以及想把分析过程标准化沉淀下来的实验室成员。标题里的“R语言代谢组学数据分析”本质上是一套从原始峰表到生物学解释的完整工作流而不是某一个单独的包或函数。下面我按实际做项目的顺序把这条链路拆开讲清楚。2. 从峰表到分析就绪矩阵数据清洗与预处理的完整链路代谢组学分析翻车最多的地方不在统计模型而在数据进入模型之前。原始峰表通常包含大量缺失值、批次效应、量纲差异直接丢进 PCA 或 PLS-DA出来的图看着漂亮但结论站不住。这一章把预处理拆成可执行的步骤每一步都给出 R 代码和参数解释。2.1 读入峰表与样本元数据的对齐常见做法是把峰表存成 CSV行是代谢物含 m/z、保留时间、名称列是样本另一份元数据表记录样本分组、批次、性别等。两份表靠样本名对齐这一步出错后面全错。# 读入峰表第一列作为行名代谢物标识 peak_table - read.csv(peak_table.csv, row.names 1, check.names FALSE) # 读入样本元数据 sample_meta - read.csv(sample_meta.csv, row.names 1, check.names FALSE) # 检查样本名是否完全一致顺序也要对齐 stopifnot(all(colnames(peak_table) rownames(sample_meta))) # 如果顺序不一致按元数据顺序重排峰表列 peak_table - peak_table[, rownames(sample_meta), drop FALSE]row.names 1表示把第一列当作行名check.names FALSE防止 R 把代谢物名称里的特殊字符自动改写。stopifnot是一道保险样本名对不上直接报错比后面发现分组错位再回头查要省事得多。元数据里分组列建议用因子类型并显式指定水平顺序否则后续画图时对照组和实验组的排列可能和你预期相反。2.2 缺失值处理删、补还是留缺失值在代谢组学里分两类真缺失低于检出限和随机缺失仪器波动。处理策略取决于缺失比例。我一般按代谢物缺失率过滤缺失超过 20% 的代谢物直接剔除剩余缺失用最小值一半或 KNN 填补。# 计算每个代谢物的缺失比例 na_ratio - rowSums(is.na(peak_table)) / ncol(peak_table) # 剔除缺失率大于 20% 的代谢物 peak_table - peak_table[na_ratio 0.2, ] # 剩余缺失用每行最小值的一半填补 min_half - apply(peak_table, 1, function(x) min(x, na.rm TRUE) / 2) for (i in seq_len(nrow(peak_table))) { peak_table[i, is.na(peak_table[i, ])] - min_half[i] }20% 这个阈值不是铁律样本量小的时候可以放宽到 30%但要在方法里写清楚。用最小值一半填补是代谢组学里比较保守的做法避免引入过大的假信号。如果数据来自多批次填补前最好先做批次校正否则填补值会带着批次偏差进入下游。2.3 归一化与转换让样本之间可比峰强度受样本上样量、仪器灵敏度影响必须归一化。常见方法有总峰面积归一化、内标归一化、概率商归一化PQN。PQN 对整体差异的校正更稳健但计算稍复杂。# 总峰面积归一化 total_area - colSums(peak_table) peak_table_norm - sweep(peak_table, 2, total_area, /) * median(total_area) # 对数转换压缩量纲差异 peak_table_log - log2(peak_table_norm 1) # 帕累托缩放常用于 PCA 前处理 pareto_scale - function(x) { (x - mean(x)) / sqrt(sd(x)) } peak_table_scaled - t(apply(peak_table_log, 1, pareto_scale))sweep按列除以总峰面积再乘中位数让所有样本回到同一量级。log2(x1)里的加一是为了防止零值取对数出错。帕累托缩放除以标准差的平方根相比自动缩放除以标准差不会过度放大低丰度代谢物的噪声这是代谢组学里比较常用的选择。归一化方法没有绝对优劣关键是全流程统一并且在方法部分交代清楚。3. 多元统计与差异筛选PCA、PLS-DA 和火山图怎么配合用预处理做完下一步是看整体结构、找差异代谢物。PCA 用来看样本聚集和离群PLS-DA 用来放大组间差异差异筛选则结合统计检验和倍数变化。这三者不是替代关系是递进关系。3.1 PCA 与离群样本判断PCA 不需要分组信息先跑一遍能快速发现离群样本和批次效应。library(ggplot2) # 转置让行是样本、列是代谢物 pca_input - t(peak_table_scaled) # 跑 PCA pca_res - prcomp(pca_input, center TRUE, scale. FALSE) # 提取得分和前两个主成分的方差解释比例 scores - as.data.frame(pca_res$x[, 1:2]) scores$group - sample_meta$group var_exp - round(summary(pca_res)$importance[2, 1:2] * 100, 1) ggplot(scores, aes(PC1, PC2, color group)) geom_point(size 3) labs(x paste0(PC1 (, var_exp[1], %)), y paste0(PC2 (, var_exp[2], %))) theme_bw()scale. FALSE是因为前面已经做过帕累托缩放再标准化会重复处理。看 PCA 图时重点关注两件事组内样本是否聚在一起有没有样本跑到置信椭圆外面。如果有离群样本先查原始峰表里该样本的总峰面积和缺失率确认是技术问题还是生物学差异不要直接删。3.2 PLS-DA 与置换检验PLS-DA 是有监督方法能最大化组间分离但也容易过拟合。必须做置换检验来评估模型是否真的有效。library(mixOmics) # 构建 PLS-DA 模型ncomp 先设一个较大值 plsda_res - plsda(pca_input, sample_meta$group, ncomp 3) # 置换检验200 次 perm_res - perf(plsda_res, validation Mfold, folds 5, nrepeat 50, progressBar FALSE) # 查看分类错误率 plot(perm_res)ncomp是主成分数不是越多越好通常看分类错误率曲线选拐点。perf做交叉验证评估模型稳定性nrepeat越大结果越稳但耗时越长。如果置换检验的 Q2 截距明显小于零说明模型没有过拟合组间差异可信。这一步很多人跳过直接拿 PLS-DA 得分图说事审稿人一问模型验证就露馅。3.3 差异代谢物筛选VIP 加 Fold Change 加 P 值差异代谢物通常要同时满足三个条件VIP 大于 1、倍数变化超过阈值、统计检验显著。# 提取 VIP 值 vip_values - vip(plsda_res) # 对每个代谢物做 t 检验 p_values - apply(peak_table_log, 1, function(x) { t.test(x ~ sample_meta$group)$p.value }) # 计算倍数变化以两组为例 group_means - t(apply(peak_table_log, 1, function(x) { tapply(x, sample_meta$group, mean) })) fold_change - group_means[, 2] - group_means[, 1] # log2 尺度上的差值 # 汇总结果 diff_res - data.frame( metabolite rownames(peak_table_log), VIP vip_values, log2FC fold_change, p_value p_values ) # 筛选条件 diff_sig - subset(diff_res, VIP 1 abs(log2FC) 1 p_value 0.05)VIP 衡量每个代谢物对组间分离的贡献大于 1 是常用阈值。log2FC在 log2 尺度上直接相减就是倍数变化的对数值绝对值大于 1 对应原始倍数变化超过 2 倍。P 值这里没有做多重检验校正实际分析中建议加一列 FDR用p.adjust(p_values, method BH)筛选时用 FDR 小于 0.05 更稳妥。三个条件同时用能有效降低假阳性但也会漏掉一些变化幅度小却稳定的代谢物阈值要根据课题灵活调整。4. 通路富集与结果解释从代谢物列表到生物学故事差异代谢物列表本身没有生物学意义必须映射到通路才能讲出故事。富集分析的核心是看某些通路里的差异代谢物是否比随机期望更多。4.1 基于 KEGG 的富集分析流程常见做法是用代谢物名称或 KEGG ID 做超几何检验。R 里可以用clusterProfiler配合 KEGG 注释但代谢组学的 KEGG 富集和转录组略有不同需要自己整理代谢物到通路的映射。library(clusterProfiler) # 差异代谢物的 KEGG ID 列表 diff_kegg - c(C00031, C00186, C00041) # 示例 ID # 背景为所有检测到的代谢物的 KEGG ID background_kegg - c(C00031, C00186, C00041, C00064, C00065) # 富集分析 enrich_res - enrichKEGG(gene diff_kegg, universe background_kegg, organism hsa, pvalueCutoff 0.05) # 查看结果 head(as.data.frame(enrich_res))universe参数必须填所有检测到的代谢物对应的 KEGG ID不能省略否则背景不对富集结果会虚高。organism hsa是人类其他物种要换成对应的 KEGG 代码。富集结果里的Count是通路中命中的差异代谢物数量pvalue是超几何检验的原始 P 值建议同时看p.adjust。4.2 富集结果的可视化与解读边界富集气泡图是最常用的展示方式横轴是富集因子纵轴是通路名点的大小表示命中数量颜色表示显著性。library(ggplot2) enrich_df - as.data.frame(enrich_res) enrich_df$richFactor - enrich_df$Count / as.numeric(sub(/.*, , enrich_df$BgRatio)) ggplot(enrich_df, aes(x richFactor, y reorder(Description, richFactor))) geom_point(aes(size Count, color p.adjust)) scale_color_gradient(low red, high blue) labs(x Rich Factor, y NULL, color p.adjust, size Count) theme_bw()richFactor是命中数量除以背景中该通路的代谢物总数值越大说明该通路富集程度越高。解读时注意两点一是富集显著不代表因果关系代谢物变化可能是下游效应二是 KEGG 通路覆盖的代谢物有限很多脂质和未知代谢物没有注释富集结果只能反映有注释的那部分。我一般会把富集结果和差异代谢物的热图放在一起看确认通路里的代谢物变化方向是否一致。5. 避坑与排查代谢组学 R 分析里最容易翻车的五个地方这一章记录的是实际项目里反复出现的问题每一条都按现象、原因、解决来写。5.1 样本名对不上导致分组错位现象PCA 图里组间完全分开但置换检验 Q2 截距为正模型明显过拟合。原因峰表列名和元数据行名顺序不一致R 按位置匹配分组信息整体错位。解决读入数据后立刻用stopifnot(all(colnames(peak_table) rownames(sample_meta)))检查不一致就按元数据顺序重排不要靠肉眼核对。5.2 缺失值填补引入假信号现象差异代谢物列表里出现大量低丰度代谢物且这些代谢物在原始数据里缺失率很高。原因用最小值一半填补后缺失样本被赋予了一个固定的小值如果该代谢物在某一组里缺失集中填补后会人为制造组间差异。解决填补前先按缺失率过滤缺失率超过 20% 的代谢物直接剔除填补后做一次敏感性分析对比填补前后差异结果是否一致。5.3 归一化方法不统一导致批次效应残留现象PCA 图里样本按批次聚集而不是按分组聚集。原因只做了总峰面积归一化没有做批次校正或者不同批次用了不同的归一化参数。解决多批次数据建议用sva或limma的removeBatchEffect做校正校正后再跑 PCA 确认批次效应是否消除。归一化参数要全流程统一不要中途换方法。5.4 PLS-DA 过拟合没有做置换检验现象PLS-DA 得分图组间分离完美但用新样本预测时准确率很低。原因主成分数选得太多模型把噪声也拟合进去了。解决用perf做交叉验证选主成分数用permutest做置换检验Q2 截距小于零才认为模型可靠。主成分数一般不超过 3样本量小的时候更要保守。5.5 富集分析背景基因集设置错误现象富集结果里出现大量明显不相关的通路且 P 值极小。原因universe参数没有设置或设置成了全部 KEGG 代谢物而不是本次实验检测到的代谢物。解决背景必须用本次实验检测到的所有代谢物对应的 KEGG ID不能图省事用默认值。富集分析前先确认有多少代谢物有 KEGG 注释注释率太低的话富集结果参考价值有限。6. 把分析流程封装成可复现脚本的几个实用技巧走到这里单次分析已经能跑通了。但代谢组学项目往往不止一批数据下次换一组样本如果还要从头改代码效率太低。我的习惯是把整个流程拆成三个文件一个负责读入和预处理一个负责统计和差异筛选一个负责富集和可视化。每个文件开头用source调用上一个文件的输出中间结果存成 RDS 格式方便单独调试某一步而不用重跑全流程。# 保存中间结果方便后续单独加载 saveRDS(peak_table_log, file step1_normalized.rds) saveRDS(diff_sig, file step2_diff_metabolites.rds) # 下次调试直接从某一步开始 peak_table_log - readRDS(step1_normalized.rds)参数不要写死在代码里统一放到一个config.R文件比如缺失率阈值、VIP 阈值、P 值阈值、分组名称。换项目时只改配置文件主流程代码不动。这样做的另一个好处是方法部分写起来有据可查每个阈值都能追溯到配置文件里的具体数值。还有一个容易被忽略的点是随机种子。PLS-DA 的交叉验证、置换检验、KNN 填补都涉及随机过程不设种子的话每次跑结果略有差异别人复现时会对不上。在脚本开头加一行set.seed(123)所有随机过程就固定了。种子数字本身不重要重要的是全流程统一并且在方法里注明。最后说一个我踩过的坑不要等到分析全部做完再写注释和整理输出。代谢组学分析步骤多隔一周回头看自己的代码经常想不起来某个筛选条件为什么设成 0.05 而不是 0.01。每跑完一步就把参数含义、选择理由、输出文件说明写在代码块上方比事后补要省力得多。希望帮到你。本文还有配套的精品资源点击获取