数据挖掘—使用MetaboAnalys进行代谢组分析及本地代码优化

原创

sheldor没耳朵

发布于 2026-08-21 13:48:29

1470

数据挖掘—使用MetaboAnalys进行代谢组分析及本地代码优化

Hello，大家好，好久不见。已经好久没写帖子记录了，不知道AI迅猛发展的今天还有没有必要记录一些代码什么的。毕竟和我目前使用的GPT-5.6 Sol相比，我两代码的能力是半斤八两。我是半斤废铁中的废铁，它是八两黄金中的黄金...

闲言少叙，最近在做代谢组学分析时，我用到了 MetaboAnalyst。这个网站对代谢组学分析非常友好，很多常见分析都可以直接在线完成，比如数据预处理、PCA、PLS-DA、差异代谢物筛选、富集分析、KEGG 通路分析等。甚至根据质谱图鉴定化合物都可以。对于不想自己从头写代码的人来说，MetaboAnalyst 的优势很明显：上传数据、选择参数、点击运行，就可以快速得到比较规范的结果和图。

我这里主要使用富集分析功能。我记得我最出做代谢组学富集分析的时候，还是自己手动搭了背景数据集，那时候，很多代谢物还并没有被KEGG收录，所以富集结果非常有限。用这个MetaboAnalyst。富集分析就极为方便了。此外，代谢组前期的数据预处理、PCA、PLS-DA、差异代谢物筛选等我还是基于我的代码进行了，并没有使用MetaboAnalyst，至于根据质谱图鉴定化合物，该功能我还并有测试。

1.MetaboAnalyst 的 Pathway Analysis 大致怎么用

这次我主要用的是它的 Pathway Analysis 模块，也就是根据一组差异代谢物，判断这些代谢物主要富集在哪些代谢通路中。

首页--> Pathway Analysis

在线分析时，大致流程是：

进入 MetaboAnalyst 的 Pathway Analysis；
上传或者粘贴差异代谢物；
指定代谢物 ID 类型，例如 KEGG ID；
选择物种；
设置背景代谢物；
选择富集分析方法和通路拓扑分析方法；
运行分析。

根据我的项目，我的参数是

富集方法：
Hypergeometric test
Topology：
Relative Betweenness Centrality
物种：拟南芥
Arabidopsis thaliana（ath）

注：另外，我没有直接使用数据库中的全部代谢物作为背景，而是使用本次实验中实际能够检测并注释到 KEGG ID 的全部代谢物作为 background。这个设计很重要。比如本次项目总共能够匹配到 431 个 KEGG 代谢物，那么通路富集分析的背景就应该尽量基于这 431 个代谢物，而不是整个 KEGG 数据库。因为真正有机会被筛选为“差异代谢物”的，本来就是实验中实际检测到的这些代谢物。

2.富集结果

上述点击proceed后，就能得到富集结果了，包括可以输出富集分析的气泡图，通路图等，可以进一步下载。此外，在match Status中，可以看到具体哪些代谢物被富集到了。遗憾的是这个结果只能

可以直接下载

局限性

3.解决方案及本地代码优化

在网页的右上端提供了R代码，直接点击即可。该代码主要使用了MetaboAnalystR包。我装这个包的时候比较难装。需要先手动装好依赖，根据提示一步步进行即可。

> remotes::install_github( 
+ "xia-lab/MetaboAnalystR", 
+ build = TRUE, 
+ build_vignettes = FALSE, 
+ dependencies = TRUE + )

#网页代码
# PID of current job: 944955
mSet<-InitDataObjects("conc", "pathora", FALSE, 150)
#差异分析代谢物列表
cmpd.vec<-c("C00655","C16417",...）#自己补全，需要向量格式
mSet<-Setup.MapData(mSet, cmpd.vec);
mSet<-CrossReferencing(mSet, "kegg");
mSet<-CreateMappingResultTable(mSet)
mSet<-SetKEGG.PathLib(mSet, "ath", "current")
mSet<-SetMetabolomeFilter(mSet, F);
mSet<-CalculateOraScore(mSet, "rbc", "hyperg")
mSet<-PlotPathSummary(mSet, F, "path_view_0_", "png", 150, width=NA, NA, NA )
PlotPathSummary(mSet, FALSE, "path_view_0_", "png", 150)
mSet<-SaveTransformedData(mSet)

安装好包，在本地实际跑上述过程中，mSet <- CalculateOraScore(mSet, "rbc", "hyperg")，这一步会因为线下API调用出错。网页版 MetaboAnalyst 很方便，但是在需要批量分析或者长期复现项目时，会有几个比较明显的问题。及其依赖服务器状态。因此我重新写了个本地运行的版本，且在最后保留了每天富集结果中对应的代谢物名称和ID，这样后续无论做什么样子的操作都更加方便了。

rm(list = ls())
options(stringsAsFactors = FALSE)
options(scipen = 20)
library(MetaboAnalystR)
library(qs2)
############################################################
# 1. 初始化 Pathway ORA
############################################################
mSet <- InitDataObjects(
  "conc",
  "pathora",
  FALSE,
  150
)

############################################################
# 2. 输入差异代谢物 KEGG ID
############################################################

cmpd.vec <- c(
  "C00655","C16417"...
)

cmpd.vec  <- as.vector(read.table("../result/step4/X50-XZ_vs_Y25-XZ_KEGG_ID.txt", 
                               header = FALSE,        # 无表头
                               stringsAsFactors = FALSE))

# 查看数据
head(data)
############################################################
# 3. KEGG ID mapping
############################################################

mSet <- Setup.MapData(
  mSet,
  cmpd.vec
)

mSet <- CrossReferencing(
  mSet,
  "kegg"
)

mSet <- CreateMappingResultTable(
  mSet
)


############################################################
# 4. 设置自定义背景代谢物
############################################################

mSet <- Setup.KEGGReferenceMetabolome(
  mSet,
  "../result/step4/All_metabolites_background_KEGG_ID.txt"
)

mSet <- SetMetabolomeFilter(
  mSet,
  TRUE
)


############################################################
# 5. 下载/读取本地 ath KEGG pathway library
############################################################

# MetaboAnalystR 官方服务器上的 KEGG pathway library
ath_url <- paste0(
  "https://www.metaboanalyst.ca/resources/libs/",
  "kegg/metpa/ath.qs"
)

ath_file <- "../result/step4/ath.qs"

# 第一次运行时下载，之后直接使用本地文件
if (!file.exists(ath_file)) {
  
  download.file(
    ath_url,
    destfile = ath_file,
    mode = "wb",
    method = "libcurl"
  )
}

current.kegglib <- qs2::qs_read(
  ath_file
)


############################################################
# 6. 检查 ath pathway library
############################################################

names(current.kegglib)

length(current.kegglib$mset.list)

length(current.kegglib$rbc)

head(current.kegglib$path.ids)


############################################################
# 7. 获取成功匹配的 KEGG ID
############################################################

nm.map <- GetFinalNameMap(
  mSet
)

# 与 MetaboAnalystR 官方逻辑一致：
# 去掉 NA 和重复 KEGG ID
valid.inx <- !(
  is.na(nm.map$kegg) |
    duplicated(nm.map$kegg)
)

ora.vec <- nm.map$kegg[
  valid.inx
]

ora.vec <- unique(
  ora.vec
)

cat(
  "Mapped query metabolites:",
  length(ora.vec),
  "\n"
)


############################################################
# 8. 获取 ath pathway compound sets
############################################################

current.mset <- current.kegglib$mset.list


############################################################
# 9. 使用自定义 background 过滤 pathway
############################################################

background.vec <- unique(
  mSet$dataSet$metabo.filter.kegg
)

cat(
  "Uploaded background metabolites:",
  length(background.vec),
  "\n"
)

# 每条 pathway 只保留本项目 background 中存在的代谢物
current.mset <- lapply(
  current.mset,
  function(x) {
    
    x[
      x %in% background.vec
    ]
  }
)

# 删除过滤后完全为空的 pathway
keep.path <- lengths(
  current.mset
) > 0

current.mset <- current.mset[
  keep.path
]

# 保存到 mSet，与官方对象结构保持一致
mSet$analSet$ora.filtered.mset <- current.mset


############################################################
# 10. 定义真正的 pathway universe
############################################################

# 注意：
# MetaboAnalystR 不是直接拿431作为 uniq.count，
# 而是取背景过滤后、至少属于一条 ath pathway 的代谢物并集

my.univ <- unique(
  unlist(
    current.mset,
    use.names = FALSE
  )
)

uniq.count <- length(
  my.univ
)

cat(
  "Effective pathway background:",
  uniq.count,
  "\n"
)


############################################################
# 11. query 也限制到有效 pathway universe
############################################################

ora.vec <- ora.vec[
  ora.vec %in% my.univ
]

q.size <- length(
  ora.vec
)

cat(
  "Query metabolites entering ORA:",
  q.size,
  "\n"
)

if (q.size < 3) {
  
  stop(
    "Less than 3 metabolites remain after pathway/background filtering."
  )
}


############################################################
# 12. 计算每条 pathway 的 hits
############################################################

hits <- lapply(
  current.mset,
  function(x) {
    
    x[
      x %in% ora.vec
    ]
  }
)

hit.num <- lengths(
  hits
)

set.num <- lengths(
  current.mset
)


############################################################
# 13. Hypergeometric ORA
############################################################

expected <- q.size * (
  set.num / uniq.count
)

raw.p <- phyper(
  q = hit.num - 1,
  m = set.num,
  n = uniq.count - set.num,
  k = q.size,
  lower.tail = FALSE
)


############################################################
# 14. 多重检验校正
############################################################

holm.p <- p.adjust(
  raw.p,
  method = "holm"
)

fdr.p <- p.adjust(
  raw.p,
  method = "fdr"
)


############################################################
# 15. Pathway topology：RBC impact
############################################################

# rbc = relative betweenness centrality
imp.list <- current.kegglib$rbc[
  names(current.mset)
]

impact <- mapply(
  function(imp, hit) {
    
    if (
      is.null(imp) ||
      length(hit) == 0
    ) {
      return(0)
    }
    
    # 官方 MetaboAnalystR：
    # sum(x[y])
    sum(
      imp[hit],
      na.rm = TRUE
    )
  },
  imp.list,
  hits
)


############################################################
# 16. 汇总 ORA 结果
############################################################

ora.mat <- data.frame(
  Total = set.num,
  Expected = expected,
  Hits = hit.num,
  `Raw p` = raw.p,
  `-log10(p)` = -log10(raw.p),
  `Holm adjust` = holm.p,
  FDR = fdr.p,
  Impact = impact,
  check.names = FALSE
)

rownames(
  ora.mat
) <- names(
  current.mset
)


############################################################
# 17. 只保留至少命中1个代谢物的 pathway
############################################################

ora.mat <- ora.mat[
  ora.mat$Hits > 0,
  ,
  drop = FALSE
]

ora.mat <- ora.mat[
  !is.na(ora.mat$Impact),
  ,
  drop = FALSE
]


############################################################
# 18. 排序
############################################################

if (nrow(ora.mat) > 1) {
  
  ora.mat <- ora.mat[
    order(
      ora.mat$`Raw p`,
      ora.mat$Impact
    ),
    ,
    drop = FALSE
  ]
}


############################################################
# 19. 保存到 mSet
############################################################

mSet$analSet$ora.mat <- signif(
  as.matrix(ora.mat),
  5
)

mSet$analSet$ora.hits <- hits

mSet$analSet$node.imp <- "rbc"

mSet$msgSet$topo.msg <-
  paste0(
    "Your selected node importance measure for ",
    "topological analysis is relative betweenness centrality."
  )

mSet$msgSet$rich.msg <-
  "The selected over-representation analysis method is Hypergeometric test."


############################################################
# 20. 将 pathway ID 转换成 pathway name
############################################################

pathway_id <- rownames(
  ora.mat
)

pathway_name <- names(
  current.kegglib$path.ids
)[
  match(
    pathway_id,
    current.kegglib$path.ids
  )
]


############################################################
# 21. 加入命中 KEGG ID
############################################################

hit_id <- sapply(
  pathway_id,
  function(x) {
    
    paste(
      hits[[x]],
      collapse = ";"
    )
  }
)


############################################################
# 22. 输出最终结果表
############################################################

pathway_result <- data.frame(
  Pathway = pathway_name,
  Pathway_ID = pathway_id,
  ora.mat,
  Hit_KEGG_ID = hit_id,
  row.names = NULL,
  check.names = FALSE
)

write.csv(
  pathway_result,
  "../result/step7/pathway_results_local.csv",
  row.names = FALSE
)

这样最终的结果中就保留代谢组名称和相应的ID了

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

生物基因

数据挖掘—使用MetaboAnalys进行代谢组分析及本地代码优化

数据挖掘—使用MetaboAnalys进行代谢组分析及本地代码优化

数据挖掘—使用MetaboAnalys进行代谢组分析及本地代码优化

1.MetaboAnalyst 的 Pathway Analysis 大致怎么用

2.富集结果

3.解决方案及本地代码优化

社区

活动

圈层

关于

腾讯云开发者

热门产品

热门推荐

更多推荐