简介：这份资源面向具备一定R语言基础、希望进入代谢组学领域的科研人员与生物信息学初学者，系统整理了从数据导入到结果报告的全流程分析思路。压缩包共61个文件，以52个R脚本为核心，辅以3个Rmd报告文件、1个README说明文档及若干辅助文件，整体约277KB，体积轻便却覆盖了PCA、PLS-DA、OPLS-DA、通路分析、bootstrap与热图相似性等典型分析模块。已有3066人学习下载，说明其在代谢组学入门与实战中具有较高参考价值。读者可借助脚本直接复现多元统计降维、差异代谢物筛选、共表达网络构建及ggplot2可视化等关键环节，并参考RMarkdown模板将代码、图表与解释整合为完整分析报告，从而快速搭建可复用的代谢组学分析流程，减少从零摸索的时间成本。
1. 从一份 R 语言代谢组学压缩包说起：它到底能帮你跑通哪几步
如果你手头有一批 LC-MS 或 GC-MS 的峰表数据，正卡在“用 Excel 处理太乱、用 Python 又得从头搭流程”的阶段，那这份 R 语言代谢组学数据分析压缩包值得先拆开看看。它本质上是一套围绕 R 语言构建的代谢组学分析脚本与配套示例数据的集合，核心解决的是从原始峰表到差异代谢物筛选、通路富集这一整条链路。适合谁？做代谢组学实验但编程基础一般的研究生、需要快速复现分析流程的检测机构技术员，以及想用 R 语言把代谢组学分析标准化的从业者。我见过太多人拿到峰表后第一反应是丢进 MetaboAnalyst 网页版，但一旦样本分组复杂、需要自定义统计模型，网页工具就捉襟见肘。这份资源的价值在于把数据预处理、缺失值填充、归一化、多元统计、单变量检验和通路映射串成了一条可修改的 R 脚本流水线，而不是一个黑匣子。你拿到手后，改改分组文件和峰表路径，就能跑出 PCA、PLS-DA 和火山图，省去大量查文档拼函数的时间。
2. 拆包先看结构：脚本、示例数据与依赖包怎么对应
2.1 压缩包目录的典型布局与文件职责
拿到压缩包后别急着 source 所有脚本，先解压看目录。常见做法是顶层放一个README.md或analysis_pipeline.R主控脚本，下面分data/、scripts/、output/三个文件夹。data/里通常包含示例峰表peak_table.csv、样本分组表sample_metadata.csv，有时还有一份feature_annotation.csv存代谢物注释信息。scripts/下按分析步骤拆成多个 R 文件，比如01_preprocessing.R、02_multivariate.R、03_univariate.R、04_pathway.R。output/是空目录，脚本跑完自动往里写图和数据表。我一般会先打开主控脚本看它 source 了哪些文件、用了哪些包，这样能快速判断这份资源的技术栈是偏 tidyverse 还是 base R，是用ropls做 PLS-DA 还是用mixOmics。
# 解压后先看目录树，确认文件层级 unzip R语言代谢组学数据分析.zip -d metabolomics_r cd metabolomics_r find . -maxdepth 2 -type f | sort这段命令只是把压缩包解开并列出两层内的文件，目的是让你在跑任何分析前先确认脚本和数据是否齐全。如果发现data/里只有峰表没有分组表，那说明这份资源需要你自己准备分组信息，常见格式是两列：样本名和分组标签，第一列列名通常要求为sample或Sample，具体看脚本里的read.csv参数。
2.2 依赖包清单与安装策略
代谢组学分析在 R 里依赖的包分两类：一类是 CRAN 上的通用包，比如tidyverse、ggplot2、ggrepel、pheatmap；另一类是 Bioconductor 上的专用包，比如ropls、mixOmics、MetaboAnalystR、pathview。压缩包里一般会附一个packages.R或直接在 README 里列出install.packages()和BiocManager::install()的命令。我习惯先检查 R 版本，代谢组学里不少包对 R 版本有要求，比如MetaboAnalystR在 R 4.2 以下容易出依赖冲突。
# 检查 R 版本和已安装包 R.version.string installed_packages <- rownames(installed.packages()) required_pkgs <- c("tidyverse", "ropls", "mixOmics", "pheatmap", "ggrepel") missing_pkgs <- setdiff(required_pkgs, installed_packages) if (length(missing_pkgs) > 0) { print(paste("缺失包：", paste(missing_pkgs, collapse = ", "))) } else { print("所有核心包已安装") }这段代码先打印 R 版本，再比对必需包列表和已安装包，输出缺失项。参数required_pkgs里的包名要根据你压缩包里实际用到的包来改，不要照抄。如果缺的是 Bioconductor 包，用BiocManager::install("ropls")安装，别直接用install.packages()，否则会报版本不匹配。常见坑是mixOmics和ropls同时加载时函数名冲突，比如plotIndiv在两个包里都有定义，脚本里最好用ropls::plotIndiv()这种命名空间前缀调用。
3. 数据预处理与归一化：从峰表到可分析矩阵的关键参数
3.1 缺失值填充与过滤阈值怎么定
代谢组学峰表里缺失值很常见，可能是低于检测限，也可能是峰提取失败。压缩包里的预处理脚本一般会先做缺失值过滤：保留在至少 80% 样本中都有值的特征。这个阈值不是固定的，如果样本分组多、每组样本量少，可以放宽到 70%；如果追求高置信度，可以收紧到 90%。填充方法常见的有最小值一半填充、KNN 填充和随机森林填充。我一般先用最小值一半填充跑一遍看分布，如果填充后 PCA 里某个组明显偏离，再换 KNN。
# 缺失值过滤与填充示例 library(tidyverse) peak_table <- read.csv("data/peak_table.csv", row.names = 1, check.names = FALSE) # 计算每个特征的非缺失比例 valid_ratio <- rowSums(!is.na(peak_table)) / ncol(peak_table) # 保留至少 80% 样本有值的特征 peak_filtered <- peak_table[valid_ratio >= 0.8, ] # 用最小值一半填充剩余缺失值 min_half <- min(peak_filtered, na.rm = TRUE) / 2 peak_filled <- peak_filtered %>% mutate(across(everything(), ~ replace_na(.x, min_half)))这段代码先读峰表，行是代谢物特征，列是样本。valid_ratio算每个特征的非缺失比例，valid_ratio >= 0.8是过滤阈值，你可以改成 0.7 或 0.9。min_half取全局最小值的一半，replace_na把所有 NA 替换成这个值。注意check.names = FALSE是为了防止 R 把样本名里的特殊字符自动改成点号，代谢组学样本名常带-或_，不改这个参数后面匹配分组表会出错。
3.2 归一化方法选择：总峰面积、内标还是 PQN
归一化是代谢组学里最容易被忽视但影响最大的步骤。压缩包通常提供三种归一化选项：总峰面积归一化、内标归一化和 PQN（概率商归一化）。总峰面积归一化假设所有样本的总代谢物量一致，适合尿液、血清这类组成相对稳定的样本；内标归一化需要你在峰表里指定内标特征列，适合有明确内标添加的实验；PQN 对样本间差异大的情况更稳健，但计算稍慢。我一般会先跑总峰面积归一化看 PCA 分组是否分开，如果分开不明显再换 PQN。
# 总峰面积归一化 peak_norm <- sweep(peak_filled, 2, colSums(peak_filled), FUN = "/") # 对数转换，减少量级差异 peak_log <- log2(peak_norm + 1) # 查看归一化前后样本总峰面积分布 boxplot(colSums(peak_filled), main = "归一化前", las = 2) boxplot(colSums(peak_log), main = "归一化后", las = 2)sweep函数把每个样本的峰强度除以该样本的总峰面积，实现总峰面积归一化。log2(x + 1)是常见的对数转换，加 1 是为了避免零值取对数报错。箱线图用来目视检查归一化效果，归一化后各样本的中位数应该大致对齐。如果某个样本的箱线图明显偏离，检查该样本的原始总峰面积是否异常低，可能是进样量不足或样品降解，这种样本可以考虑剔除，但要在报告里说明。
4. 多元统计与差异筛选：PCA、PLS-DA 和火山图的 R 实现
4.1 PCA 与 PLS-DA 的脚本参数与解读边界
PCA 是无监督的，用来观察样本自然聚集情况；PLS-DA 是有监督的，用来找组间差异最大的方向。压缩包里的多元统计脚本一般用ropls包实现，因为它同时支持 PCA、PLS-DA 和 OPLS-DA，而且自带置换检验。跑 PLS-DA 前要确认分组表里至少有两个组，每组样本数不少于 3，否则模型没有统计意义。ropls::opls()函数里predI = 1表示预测第一个分组变量，orthoI = NA表示不进行正交校正，如果要跑 OPLS-DA 就设orthoI = 1。
library(ropls) # 读入归一化后的矩阵和分组表 peak_log <- as.matrix(read.csv("output/peak_log.csv", row.names = 1)) metadata <- read.csv("data/sample_metadata.csv") # 确保样本顺序一致 peak_log <- peak_log[, metadata$sample] # 构建 PLS-DA 模型 plsda_model <- opls(x = t(peak_log), y = metadata$group, predI = 1, orthoI = NA, crossvalI = 7) # 查看模型 R2Y 和 Q2Y print(plsda_model@summaryDF) # 绘制得分图 plot(plsda_model, typeVc = "x-score")t(peak_log)是把特征×样本的矩阵转置成样本×特征，因为opls要求输入样本在行。crossvalI = 7是七折交叉验证，样本量少于 20 时建议用留一法，设crossvalI = 1。summaryDF里的 R2Y 表示模型解释分组变量的能力，Q2Y 表示预测能力，Q2Y 大于 0.5 算可接受，但不要盲目追求高 Q2Y，过拟合时 Q2Y 也会虚高。置换检验用permI = 200跑 200 次，看 Q2Y 的 p 值是否小于 0.05。
4.2 单变量检验与火山图：fold change 和 FDR 的阈值设定
多元统计告诉你组间有差异，单变量检验告诉你哪些代谢物具体差异显著。压缩包里的单变量脚本一般用limma或t.test加 BH 校正。limma更适合样本量小、方差齐性不确定的情况。火山图的横轴是 log2 fold change，纵轴是 -log10 FDR。阈值常见设法是 fold change 绝对值大于 1.5 或 2，FDR 小于 0.05。我一般先设 fold change > 1.5、FDR < 0.05 看差异代谢物数量，如果太多就收紧到 fold change > 2，太少就放宽 FDR 到 0.1，但要在方法里写清楚。
library(limma) # 构建设计矩阵 group <- factor(metadata$group) design <- model.matrix(~ 0 + group) colnames(design) <- levels(group) # 拟合线性模型 fit <- lmFit(peak_log, design) # 设置对比：实验组 vs 对照组 contrast_matrix <- makeContrasts(GroupA - GroupB, levels = design) fit2 <- contrasts.fit(fit, contrast_matrix) fit2 <- eBayes(fit2) # 提取差异结果 diff_results <- topTable(fit2, number = Inf, adjust.method = "BH") # 筛选显著差异代谢物 sig_metabolites <- diff_results[diff_results$adj.P.Val < 0.05 & abs(diff_results$logFC) > 1.5, ]model.matrix(~ 0 + group)生成不含截距的设计矩阵，makeContrasts里写你要比较的组名，注意组名要和分组表里的标签完全一致，大小写敏感。topTable的adjust.method = "BH"是 Benjamini-Hochberg 校正，代谢组学里常用。sig_metabolites就是最终差异代谢物列表，行名是特征 ID。火山图用ggplot2画，geom_point加geom_text_repel标注 top 10 代谢物，ggrepel包能避免标签重叠。
5. 避坑与排查：代谢组学 R 脚本跑不通的五个常见原因
5.1 样本名匹配翻车：分组表和峰表列名不一致
现象：脚本报错subscript out of bounds或undefined columns selected，或者跑出来的 PCA 图样本点全挤在一起。原因：峰表列名是Sample_01，分组表里写的是sample01，R 按字符串精确匹配，差一个下划线就找不到。解决：读入数据后先打印colnames(peak_table)和metadata$sample，用setdiff看差集，统一改成一致。我习惯在脚本开头加一句stopifnot(all(metadata$sample %in% colnames(peak_table)))，不满足直接停，比跑到一半报错好排查。
5.2 包版本冲突：ropls 和 mixOmics 同时加载后函数被覆盖
现象：调用plotIndiv时提示参数不对，或者opls函数找不到。原因：mixOmics和ropls都有plotIndiv，后加载的包会覆盖先加载的。解决：用命名空间前缀调用，比如ropls::plotIndiv()，或者在脚本里只加载一个包，另一个用requireNamespace按需调用。常见做法是在packages.R里按分析步骤分批加载，预处理阶段只加载tidyverse，多元统计阶段再加载ropls。
5.3 归一化后出现负值或零值导致对数转换报错
现象：log2报NaNs produced或-Inf。原因：PQN 归一化或某些内标归一化方法可能产生零值或负值，直接取对数就出问题。解决：先检查summary(peak_norm)，如果有零值，用log2(x + 1)或log2(x + min_half)；如果有负值，说明归一化方法选错了，换总峰面积归一化或检查内标列是否选错。我一般会在对数转换前加一句if (any(peak_norm <= 0)) stop("归一化矩阵存在非正值，请检查归一化方法")。
5.4 通路富集分析报错：KEGG 数据库连接失败或物种代码不对
现象：pathview或MetaboAnalystR跑通路富集时卡住或报cannot open URL。原因：KEGG 数据库需要在线访问，网络不稳定时会失败；另外物种代码写错也会返回空结果。解决：先确认物种代码，人类是hsa，小鼠是mmu，大鼠是rno，代谢组学里常见的是hsa和mmu。如果网络受限，可以下载 KEGG 的本地数据库文件，用pathview的kegg.native = FALSE参数走本地。常见做法是先把差异代谢物的 KEGG ID 整理成两列：代谢物名和 KEGG ID，再跑富集。
5.5 输出目录不存在导致脚本写入失败
现象：脚本跑完没有报错但output/里是空的。原因：脚本里用write.csv或ggsave写文件时，目标目录不存在，R 不会自动创建目录，静默失败。解决：在脚本开头加dir.create("output", showWarnings = FALSE)，或者用if (!dir.exists("output")) dir.create("output")。我习惯在每个写文件操作前都确认目录存在，尤其是把脚本拆成多个文件跑的时候，单独跑某个脚本容易忘了建目录。
6. 进阶技巧：把这份脚本改造成可复用的分析模板
跑通一遍之后，最有价值的动作是把这份压缩包里的脚本改造成你自己的分析模板。我一般会做三件事：第一，把分组表、峰表路径和输出目录提取成配置文件config.R，里面只放路径和阈值参数，分析脚本用source("config.R")读入，这样换一批数据只改一个文件。第二，把重复的绘图代码封装成函数，比如plot_pca()、plot_volcano()，参数里暴露point_size、label_top_n这些常用调整项。第三，加一个sessionInfo()输出到output/session_info.txt，记录每次跑分析时的 R 版本和包版本，投稿时审稿人问分析环境可以直接附上。
# config.R 示例 config <- list( peak_table_path = "data/peak_table.csv", metadata_path = "data/sample_metadata.csv", output_dir = "output", missing_threshold = 0.8, fc_cutoff = 1.5, fdr_cutoff = 0.05, species = "hsa" ) # 在分析脚本开头 source("config.R") dir.create(config$output_dir, showWarnings = FALSE)这个配置文件把路径和阈值集中管理，missing_threshold控制缺失值过滤比例，fc_cutoff和fdr_cutoff控制差异代谢物筛选阈值，species控制通路富集的物种代码。换项目时只改config.R，分析脚本一行不动。我还会在config.R里加一个set.seed(123)，保证随机森林填充和置换检验的结果可重复。从那以后我每次拿到新的峰表，都强制先跑一遍config.R检查路径和阈值，再跑主控脚本，省得跑到一半发现分组表列名又对不上。希望帮到你。
本文还有配套的精品资源，点击获取
