简介这份资源面向具备一定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_modelsummaryDF) # 绘制得分图 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 或 2FDR 小于 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分组表里写的是sample01R 按字符串精确匹配差一个下划线就找不到。解决读入数据后先打印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检查路径和阈值再跑主控脚本省得跑到一半发现分组表列名又对不上。希望帮到你。本文还有配套的精品资源点击获取