我一直在尝试使用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途径中取得了一些成功,但在其他途径中没有成功。
1条答案
按热度按时间yi0zb3m41#
我联系了multiGSEA包的开发人员,他回答说,在当前版本中,来自metaboliteIDmapping的tibble没有导入到multiGSEA包的命名空间中。他已经在Bioconductor的devel和release分支中推出了一个bug修复,所以如果其他人遇到类似的问题,请使用devel版本或等待未来的更新。