R中的multiGSEA软件包未将代谢物适当分配至通路

oknrviil  于 2023-02-06  发布在  其他
关注(0)|答案(1)|浏览(237)

我一直在尝试使用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))

下一部分是定制路径定义,我认为这就是问题的根源

一个三个三个一个

在这里,您可以看到没有成功Map到代谢组途径-这是不正确的。我已经验证了许多HMDB值应该已经Map(其中〉300个与KEGG途径一致)。

接下来,我将运行富集分数,然后提取/校正p值。但是,由于代谢物组的途径比对失败,我将在继续富集之前强调我在下面尝试的一些故障排除。

我创建了一个注解中心文件,以更仔细地查看我的代谢组学标识符,并确保它们应该Map

## 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路径

接下来,我尝试将输入重新Map到KEGG而不是HMDB

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

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

注意:Map的KEGG ID比HMDB少

我尝试了同样的操作,但将returnMetabolome输出更改为KEGG,以查看它是否正确识别输入,但随后无法输出它们

一个12b1x一个13b1x

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

另一次致富的尝试

enrichment_scores <- multiGSEA( pathways, omics_data)

看起来它起作用了,所以现在我将提取p值并更正

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途径,但在reactome中完全失败(或者,如果我将数据库更改为"所有",则在除KEGG外的几乎所有位置失败)

我尝试将输入保留为KEGG,但将returnMetabolome切换为HMDB

一个15b1x一个16b1x一个17b1x

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

我尝试了不同的方法将HMDB标识符与途径联系起来。我尝试了与代谢物IDMap合并,并从HMDB切换到KEGG,在KEGG途径中取得了一些成功,但在其他途径中没有成功。

yi0zb3m4

yi0zb3m41#

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

相关问题