【问题标题】：multiGSEA package in R not appropriately assigning metabolites to pathwaysR 中的 multiGSEA 包没有适当地将代谢物分配给途径
【发布时间】：2023-02-02 06:48:34
【问题描述】：

我一直在尝试使用 multiGSEA 包 Vignette for multiGSEA 为合并转录组学和代谢组学的途径生成组合 p 值。

即使在他们的小插图中，您也可以看到我遇到的问题——在我看来，代谢物映射没有适当地将代谢物分配给它们各自的途径。

下面，我使用 multiGSEA 小插图数据来演示我所看到的问题。有没有人对如何修复代谢物调用与实际途径的对齐有想法？我错过了一些明显的东西吗？

提前致谢！

library(multiGSEA)
library("org.Hs.eg.db")
library(magrittr)
library(AnnotationDbi)
library(AnnotationHub)

从小插图加载样本数据

data(transcriptome)
data(proteome)
data(metabolome)

下一节直接来自小插图，只是创建数据结构并用示例数据填充它

omics_data <- initOmicsDataStructure( layer = c("transcriptome", 
                                                "proteome",
                                                "metabolome"))
omics_data$transcriptome <- rankFeatures( transcriptome$logFC, 
                                          transcriptome$pValue)
names( omics_data$transcriptome) <- transcriptome$Symbol

omics_data$proteome <- rankFeatures(proteome$logFC, proteome$pValue)
names( omics_data$proteome) <- proteome$Symbol

omics_data$metabolome <- rankFeatures(metabolome$logFC, metabolome$pValue)
names( omics_data$metabolome) <- metabolome$HMDB
names( omics_data$metabolome) <- gsub( "HMDB", "HMDB00", 
                                       names( omics_data$metabolome))

下一节是自定义路径定义，我认为这是问题的根源

databases <- c( "kegg", "reactome")
layers <- names( omics_data)

pathways <- getMultiOmicsFeatures( dbs = databases, layer = layers,
                                   returnTranscriptome = "SYMBOL",
                                   returnProteome = "SYMBOL",
                                   returnMetabolome = "HMDB",
                                   useLocal = TRUE)

pathways_short <- lapply( names( pathways), function( name){
                          head( pathways[[name]], 2)
                        })
names( pathways_short) <- names( pathways)
pathways_short

在这里，您可以看到没有任何东西成功映射到代谢组通路——这是不正确的。我已经验证了许多 HMDB 值应该已经映射（其中超过 300 个与 KEGG 路径一致，特别是）。

接下来，我将运行浓缩分数，然后提取/更正 p 值。但是，由于代谢组的通路比对失败，我将在下面重点介绍我在进行扩充之前尝试过的一些故障排除。

我创建了一个注释中心文件来仔细查看我的代谢组学标识符并确保它们应该映射

## create a "data" file that shows a key for each HMDB to other identifiers, and merge with metabolome data
ah <- AnnotationHub()
datasets <- query( ah, "metaboliteIDmapping")
data <- ah[["AH83115"]]

metabolome$HMDB <- sub("HMDB","HMDB00",metabolome$HMDB)
merge(metabolome,data, by = "HMDB") -> test
## remove duplicated HMDB values from dataset
test[!duplicated(test$HMDB),] -> test

再试一次，但只使用代谢组和清理后的数据

omics_data <- initOmicsDataStructure( layer = c("metabolome"))
omics_data$metabolome <- rankFeatures(test$logFC, test$pValue)
names( omics_data$metabolome) <- test$HMDB

databases <- c( "kegg", "reactome")
layers <- names( omics_data)

pathways <- getMultiOmicsFeatures( dbs = databases, layer = layers,
                                   returnTranscriptome = "SYMBOL",
                                   returnProteome = "SYMBOL",
                                   returnMetabolome = "HMDB",
                                   useLocal = TRUE)

pathways_short <- lapply( names( pathways), function( name){
                          head( pathways[[name]], 2)
                        })
names( pathways_short) <- names( pathways)
pathways_short

我尝试了同样的事情，但将 returnMetabolome 输出更改为 KEGG 以查看它是否正确识别输入但随后未能输出它们

databases <- c( "kegg", "reactome")
layers <- names( omics_data)

pathways <- getMultiOmicsFeatures( dbs = databases, layer = layers,
                                   returnTranscriptome = "SYMBOL",
                                   returnProteome = "SYMBOL",
                                   returnMetabolome = "KEGG",
                                   useLocal = TRUE)

pathways_short <- lapply( names( pathways), function( name){
                          head( pathways[[name]], 2)
                        })
names( pathways_short) <- names( pathways)
pathways_short

现在，getMultiOmicsFeatures 至少将 KEGG 标识符分配给特定路径

因为我现在看到了路径值，所以我尝试运行扩充：

enrichment_scores <- multiGSEA( pathways, omics_data)

不幸的是，它没有注释我正确输入的任何 HMDB 值并将它们分配给任何 KEGG 或 recatome 通路

接下来，我尝试将输入重新映射到 KEGG 而不是 HMDB

omics_data <- initOmicsDataStructure( layer = c("metabolome"))

omics_data$metabolome <- rankFeatures(test$logFC, test$pValue)
names( omics_data$metabolome) <- test$KEGG

注意：映射的 KEGG ID 比 HMDB 少

我尝试了同样的事情，但将 returnMetabolome 输出更改为 KEGG 以查看它是否正确识别输入但随后未能输出它们

databases <- c( "kegg", "reactome")
layers <- names( omics_data)


pathways <- getMultiOmicsFeatures( dbs = databases, layer = layers,
                                   returnTranscriptome = "SYMBOL",
                                   returnProteome = "SYMBOL",
                                   returnMetabolome = "KEGG",
                                   useLocal = TRUE)`

pathways_short <- lapply( names( pathways), function( name){
                          head( pathways[[name]], 2)
                        })
names( pathways_short) <- names( pathways)
pathways_short

现在，getMultiOmicsFeatures 至少将 KEGG 标识符分配给特定路径

另一种丰富的尝试

enrichment_scores <- multiGSEA( pathways, omics_data)

看起来它有效，所以现在我将提取 pvalues 并更正

df <- extractPvalues( enrichmentScores = enrichment_scores,
                      pathwayNames = names( pathways[[1]]))

df$combined_pval <- combinePvalues( df)
df$combined_padj <- p.adjust( df$combined_pval, method = "BH")

df <- cbind( data.frame( pathway = names( pathways[[1]])), df)

它成功地将 KEGG 标识符链接到 KEGG 通路，但它在反应组完全失败（或者，如果我将数据库更改为“全部”，它几乎在除了 KEGG 之外的所有地方都失败）

我尝试将输入保持为 KEGG，但将 returnMetabolome 切换为 HMDB

databases <- c( "kegg", "reactome")
layers <- names( omics_data)

pathways \<- getMultiOmicsFeatures( dbs = databases, layer = layers,
returnTranscriptome = "SYMBOL",
returnProteome = "SYMBOL",
returnMetabolome = "HMDB",
useLocal = TRUE)



pathways_short <- lapply( names( pathways), function( name){
                          head( pathways[[name]], 2)
                        })
names( pathways_short) <- names( pathways)
pathways_short

但这也无法用 HMDB ID 注释任何内容

我尝试了将 HMDB 标识符链接到路径的不同方法。我尝试与 metaboliteIDmapping 合并并从 HMDB 切换到 KEGG，特别针对 KEGG 途径取得了一些成功，但对任何其他途径都没有。

【问题讨论】：

标签： r bioinformatics

【解决方案1】：

我联系了 multiGSEA 包开发人员，他回答说存在一个问题，即来自 metaboliteIDmapping 的 tibble 没有导入到当前版本的 multiGSEA 包的命名空间中。他已经在 Bioconductor 的开发和发布分支中推送了一个错误修复，所以如果其他人遇到类似的问题，请使用开发版本或等待未来的更新。

【讨论】：