Hello all. I am using bioconductor pathview package, and have a question. My goal is to input 'up-regulated genes' and get output of enriched pathway maps. I generated one pathview map as an example compared it to my Differentially Expressed Genes statistics, just to be sure that the map is drawn from my input DEGs.
Then I faced some doubts : - Some genes are represented in the pathway map, but not mapped (all.mapped = NA). How can this be possible? - Mapped genes are colored (ex - mol.data = 6.85822838471096, mol.col = #FF0000). Red colors seem to be based on foldchanges, but 1 gene data does not match between pathway map datat and DEG statistics. (ex - mol.data = 48.21641641, log2FoldChange = 5.477477252)
Is this kind of a bug?
Here is my script :
## Let's annotate entrez gene id, for KEGG pathway is based on entrez id
library("AnnotationDbi")
library("org.Hs.eg.db")
columns(org.Hs.eg.db)
res$symbol = mapIds(org.Hs.eg.db,
keys=row.names(res),
column="SYMBOL",
keytype="ENSEMBL",
multiVals="first")
res$entrez = mapIds(org.Hs.eg.db,
keys=row.names(res),
column="ENTREZID",
keytype="ENSEMBL",
multiVals="first")
res$name = mapIds(org.Hs.eg.db,
keys=row.names(res),
column="GENENAME",
keytype="ENSEMBL",
multiVals="first")
head(res, 10)
## setup the KEGG data-sets
library(pathview)
library(gage)
library(gageData)
data(kegg.sets.hs)
data(sigmet.idx.hs)
kegg.sets.hs = kegg.sets.hs[sigmet.idx.hs]
head(kegg.sets.hs, 3)
##### Input up DEG pathview #####################################################
## get significant res : logFoldchange >=2, p.val <= 0.05
## subset removes NAs.
res_sig = subset(res, (log2FoldChange >= 2) & pvalue <= 0.05)
foldchanges = res_sig$log2FoldChange
names(foldchanges) = res_sig$entrez
head(foldchanges)
##### Run GAGE ################################################################
keggres = gage(foldchanges, gsets=kegg.sets.hs, same.dir=TRUE)
##### Draw an example pathway map ################################################
## 1) Get data that are used to draw a pathway map
pv.out <- pathview(gene.data = foldchanges, pathway.id = '00040', speices = 'hsa', same.layer = F, limit = c(-10,10))
list.files(pattern = 'hsa04110', full.names = T)
write.csv(pv.out$plot.data.cpd, 'hsa00040_cpd_data.csv')
write.csv(pv.out$plot.data.gene, 'hsa00040_gene_data.csv')
## 2) Draw example pathway map
pathview(gene.data = foldchanges, pathway.id = '00040', speices = 'hsa', same.layer = F, limit = c(-10,10))