代谢组学数据分析¶

代谢组学实验类型：

靶向代谢组学（Targeted）：
  - 预先选定50-500个已知代谢物
  - 使用标准品建立定量方法
  - 定量准确（绝对定量）
  - 适合验证性研究

非靶向代谢组学（Untargeted）：
  - 不预设目标，检测所有可检测的代谢物
  - 通常检测到1000-10000个features
  - 半定量（相对定量）
  - 适合发现性研究

实验设计要点：
  1. 样本量：每组≥6个生物学重复（推荐≥10）
  2. QC样本：混合所有样本的等量混合物，每5-10个样本间插入一个QC
  3. 空白对照：溶剂空白，排除背景干扰
  4. 随机化：样本注射顺序随机化，避免批次效应

代谢组学数据的三种输入形式：

1. 原始数据（.mzML/.mzXML/.raw）
   → 需要XCMS/MZmine/MS-DIAL处理
   → 峰提取 → 对齐 → 缺失值填补 → 特征表

2. 特征表（Feature Table）
   → 行=features（代谢物特征），列=样本
   → 可直接用MetaboAnalyst分析

3. 浓度表（Concentration Table）
   → 靶向代谢组学的结果
   → 已知代谢物名称和绝对浓度
   → 直接统计分析

# 安装XCMS
BiocManager::install("xcms")  # 安装
library(xcms)                  # 加载

# 读取原始数据（mzML格式）
# 将.raw文件转换为mzML（用MSConvert/ProteoWizard）
raw_files <- list.files("mzml_files/",
                         pattern = "*.mzML",
                         full.names = TRUE)  # 列出所有mzML文件

# 读取数据
raw_data <- readMSData(
  files = raw_files,
  mode = "onDisk"              # 磁盘模式（节省内存）
)

# 查看数据信息
raw_data                        # 对象摘要
table(msLevel(raw_data))        # MS1和MS2扫描数

# 1. 峰提取（Peak Detection）——CentWave方法
cwp <- CentWaveParam(
  peakwidth = c(5, 30),        # 峰宽范围（秒）
  ppm = 15,                    # m/z容差（ppm）
  snthresh = 10,               # 信噪比阈值
  noise = 1000,                # 噪声水平
  prefilter = c(3, 1000)       # 至少3个扫描超过1000强度
)

xdata <- findChromPeaks(raw_data, param = cwp)  # 提取色谱峰
cat("检测到峰数量:", nrow(chromPeaks(xdata)), "\n")

# 2. 保留时间校正（RT Alignment）
# 校正不同样本间的保留时间漂移
ogp <- ObiwarpParam(
  binSize = 0.6                # m/z分箱大小
)
xdata <- adjustRtime(xdata, param = ogp)  # 保留时间校正

# 可视化校正效果
plotAdjustedRtime(xdata)       # 校正前后的保留时间偏移

# 3. 峰对应（Peak Correspondence）
# 把不同样本中的相同峰"对齐"到同一个feature
pdp <- PeakDensityParam(
  sampleGroups = xdata$group,  # 样本分组信息
  bw = 5,                      # 带宽（秒）
  minFraction = 0.5,           # 至少50%的样本中出现
  binSize = 0.025              # m/z分箱大小
)
xdata <- groupChromPeaks(xdata, param = pdp)  # 峰对应

# 4. 缺失值填补（Gap Filling）
xdata <- fillChromPeaks(xdata)  # 用原始数据回填缺失值

# 5. 提取特征表
feature_table <- featureValues(xdata, value = "into")  # 峰面积
feature_defs <- featureDefinitions(xdata)               # 特征定义（mz, rt）

# 导出为CSV
write.csv(feature_table, "feature_table.csv")  # 导出

# MetaboAnalystR 4.0是MetaboAnalyst的R包版本
# 统一了从原始数据处理到功能解读的完整流程

# 安装
# install.packages("MetaboAnalystR")
library(MetaboAnalystR)

# 方法1：从原始mzML数据开始
# MetaboAnalystR 4.0内置了auto-optimized峰检测
mSet <- PerformDataProcessing(
  rawfiles = list.files("mzml_files/", full.names = TRUE),  # 原始文件
  param_optimization = TRUE     # 自动优化参数（推荐）
)

# MetaboAnalystR 4.0 vs XCMS：
# - 真实峰数量增加109.1%
# - 总峰数量仅增加6.79%（噪声少）
# - 自动参数优化，不需要手动调参

# 方法2：从特征表开始
mSet <- InitDataObjects("conc", "stat", FALSE)  # 初始化
mSet <- Read.TextData(mSet, "feature_table.csv",  # 读取数据
                       "colu", "disc")             # 列格式，离散分组
mSet <- SanityCheckData(mSet)                      # 数据检查
mSet <- ReplaceMin(mSet)                           # 缺失值用最小值替代
mSet <- FilterVariable(mSet, "iqr", "F", 25)      # IQR过滤（去低变异特征）
mSet <- PreparePrenormData(mSet)                   # 准备归一化
mSet <- Normalization(mSet, "MedianNorm",          # 中位数归一化
                      "LogNorm", "AutoNorm")       # log转换+自动缩放

MetaboAnalyst 6.0（https://www.metaboanalyst.ca）

2024-2026年更新：
  - 支持MS2数据处理和注释（新增）
  - 因果分析模块（孟德尔随机化）
  - 剂量-响应分析模块
  - 扩展到~130个物种
  - GSEA和mummichog式富集（2026年5月新增）
  - 增强的PCA（2026年1月，支持大数据）
  - R 4.5.2版本升级（2025年12月）

支持的数据格式：
  1. 浓度表（行=样本，列=代谢物）
  2. 峰强度表（行=样本，列=m/z_RT）
  3. 光谱箱（NMR数据）
  4. MS峰列表

数据预处理选项：
  - 缺失值：最小值替代/KNN/中位数/0
  - 过滤：IQR/标准差/变异系数
  - 归一化：总量/中位数/参考样本
  - 转换：log/平方根/立方根
  - 缩放：自动缩放/Pareto/范围缩放

# 在R中模拟MetaboAnalyst的分析流程

# 1. PCA（主成分分析）——无监督
library(FactoMineR)
library(factoextra)

pca_result <- PCA(t(feature_table_norm), graph = FALSE)  # 转置后PCA
fviz_pca_ind(pca_result,
             col.ind = sample_groups,    # 按分组着色
             addEllipses = TRUE,         # 添加置信椭圆
             palette = c("red", "blue"), # 颜色
             legend.title = "Group")     # 图例标题

# 2. PLS-DA（偏最小二乘判别分析）——有监督
library(mixOmics)

plsda_result <- plsda(
  X = t(feature_table_norm),  # 特征矩阵
  Y = sample_groups,          # 分组标签
  ncomp = 3                   # 3个成分
)
plotIndiv(plsda_result, comp = c(1, 2),  # 画前2个成分
          group = sample_groups,          # 按组着色
          ind.names = FALSE,              # 不显示样本名
          legend = TRUE,                  # 显示图例
          title = "PLS-DA Score Plot")

# VIP分数（Variable Importance in Projection）
vip_scores <- vip(plsda_result)  # 计算VIP
vip_important <- vip_scores[vip_scores[, 1] > 1, ]  # VIP > 1的变量
cat("VIP > 1的特征数:", nrow(vip_important), "\n")

# 3. 差异代谢物分析
# t检验/Wilcoxon检验/ANOVA
library(tidyverse)

# 每个feature做t检验
pvalues <- apply(feature_table_norm, 1, function(x) {
  t.test(x[group == "Disease"], x[group == "Control"])$p.value  # t检验
})
fdr <- p.adjust(pvalues, method = "BH")  # FDR校正

# 计算Fold Change
fc <- apply(feature_table_norm, 1, function(x) {
  mean(x[group == "Disease"]) - mean(x[group == "Control"])  # log2差值
})

# 差异代谢物
diff_metabolites <- data.frame(
  feature = names(pvalues),
  pvalue = pvalues,
  FDR = fdr,
  log2FC = fc
)
sig_metabolites <- diff_metabolites %>%
  filter(FDR < 0.05 & abs(log2FC) > 1)  # 显著差异代谢物
cat("显著差异代谢物:", nrow(sig_metabolites), "\n")

# 4. 火山图
ggplot(diff_metabolites, aes(x = log2FC, y = -log10(FDR))) +
  geom_point(aes(color = FDR < 0.05 & abs(log2FC) > 1), size = 1) +  # 着色
  scale_color_manual(values = c("grey", "red")) +  # 灰色/红色
  geom_hline(yintercept = -log10(0.05), linetype = "dashed") +  # P值线
  geom_vline(xintercept = c(-1, 1), linetype = "dashed") +  # FC线
  labs(title = "Volcano Plot", x = "log2(Fold Change)", y = "-log10(FDR)") +
  theme_minimal()

# 非靶向代谢组学的关键步骤：把m/z值对应到具体代谢物

# 方法1：精确质量搜索
# 用m/z值在HMDB/KEGG/PubChem中搜索可能的代谢物
# MetaboAnalyst的MS Peaks to Paths模块

# 方法2：MS2谱图匹配（MetaboAnalyst 6.0新功能）
# 将实验MS2谱与参考数据库（GNPS/MassBank/NIST）比对

# 鉴定等级（Metabolomics Standards Initiative, MSI）：
# Level 1: 与标准品比对确认（m/z + RT + MS2）
# Level 2: 与数据库MS2谱匹配（无标准品）
# Level 3: 仅基于精确质量推测
# Level 4: 未知的特征

# R中的MS2匹配示例
library(CompoundDb)
# 读取MS2数据库
# db <- CompDb("MassBank_human.sqlite")
# 搜索匹配
# matches <- matchSpectra(query_spectra, db, mzTol = 10)

# 代谢通路分析——把差异代谢物映射到KEGG通路

# 方法1：MetaboAnalyst的Pathway Analysis
# 使用MetaboAnalystR
mSet <- InitDataObjects("conc", "pathinteg", FALSE)
mSet <- SetMetabolomeFilter(mSet, FALSE)
mSet <- SetCurrentMsetLib(mSet, "pathway", 2)  # 设置通路库

# 通路分析（ORA + 拓扑分析）
mSet <- Setup.MapData(mSet, sig_metabolite_names)  # 设置代谢物列表
mSet <- CrossReferencing(mSet, "name")              # 名称匹配
mSet <- CreateMappingResultTable(mSet)              # 创建映射结果
mSet <- SetKEGG.PathLib(mSet, "hsa")               # 人类KEGG通路
mSet <- SetMetabolomeFilter(mSet, FALSE)
mSet <- CalculateOraScore(mSet, "rbc", "hyperg")   # ORA分析

# 方法2：mummichog算法（不需要预先鉴定代谢物）
# 直接从m/z列表推断通路活性
# 2026年MetaboAnalyst新增了GSEA支持

# 代谢组学分析流程
原始数据(.mzML) → XCMS/MetaboAnalystR峰提取
→ 保留时间校正 → 峰对齐 → 缺失值填补
→ 归一化(中位数+log+自动缩放)
→ PCA/PLS-DA → 差异分析(t检验+FDR)
→ 代谢物鉴定(MS2匹配) → 通路分析(KEGG)

# XCMS关键参数（CentWave）
peakwidth: c(5,30)（窗口宽度，秒）
ppm: 10-25（质量精度）
snthresh: 5-10（信噪比）
prefilter: c(3,1000)（最少扫描和强度）
bw: 5-10（对齐带宽）

# MetaboAnalyst 6.0功能模块
Statistical Analysis: 单变量/多变量统计
Enrichment Analysis: 代谢物集富集分析
Pathway Analysis: KEGG通路分析（ORA+拓扑）
Biomarker Analysis: ROC曲线+交叉验证
Causal Analysis: 孟德尔随机化（2024新增）
Dose-Response: 剂量-响应分析（2024新增）
MS2 Annotation: MS2谱图注释（2024新增）

# 差异代谢物筛选标准
VIP > 1（PLS-DA重要性）
FDR < 0.05（t检验/Wilcoxon）
|log2FC| > 1（倍数变化）
推荐三者结合使用

# 代谢物鉴定等级（MSI）
Level 1: 标准品确认（最高）
Level 2: 数据库MS2匹配
Level 3: 精确质量推测
Level 4: 未知特征（最低）