简介：这份资源面向具备一定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(x+1)里的加一是为了防止零值取对数出错。帕累托缩放除以标准差的平方根，相比自动缩放（除以标准差）不会过度放大低丰度代谢物的噪声，这是代谢组学里比较常用的选择。归一化方法没有绝对优劣，关键是全流程统一，并且在方法部分交代清楚。
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。每跑完一步就把参数含义、选择理由、输出文件说明写在代码块上方，比事后补要省力得多。希望帮到你。
本文还有配套的精品资源，点击获取
