Annotation with KEGGrest
The first few lines of my file look like this
Gene adj.P.Val
ENSDARG00000087345 1.39E-128
ENSDARG00000058371 5.68E-50
ENSDARG00000026403 1.77E-48
ENSDARG00000006427 2.36E-48
ENSDARG00000041381 8.10E-39
ENSDARG00000071601 1.81E-33
ENSDARG00000054191 2.73E-32
ENSDARG00000098488 5.67E-31
ENSDARG00000076129 2.31E-29.
I wanted to perform the KEGG enrichment test for selected genes in ZebraFIsh.
Code snippets used are
# Pull all pathways for AT
pathways.list <- keggList("pathway", "dre")
head(pathways.list)
# Pull all genes for each pathway
pathway.codes <- sub("path:", "", names(pathways.list))
genes.by.pathway <- sapply(pathway.codes,
function(pwid){
pw <- keggGet(pwid)
if (is.null(pw[[1]]$GENE)) return(NA)
pw2 <- pw[[1]]$GENE[c(TRUE,FALSE)] # may need to modify this to c(FALSE, TRUE) for other organisms
pw2 <- unlist(lapply(strsplit(pw2, split = ";", fixed = T), function(x)x[1]))
return(pw2)
}
)
head(genes.by.pathway)
My dataframe is read as DE.table <- read.delim(df, stringsAsFactors = F)
geneList <- DE.table$P.Value
names(geneList) <- DE.table$Gene
head(geneList)
# Wilcoxon test for each pathway
pVals.by.pathway <- t(sapply(names(genes.by.pathway),
function(pathway) {
pathway.genes <- genes.by.pathway[[pathway]]
list.genes.in.pathway <- intersect(names(geneList), pathway.genes)
list.genes.not.in.pathway <- setdiff(names(geneList), list.genes.in.pathway)
scores.in.pathway <- geneList[list.genes.in.pathway]
scores.not.in.pathway <- geneList[list.genes.not.in.pathway]
if (length(scores.in.pathway) > 0){
p.value <- wilcox.test(scores.in.pathway, scores.not.in.pathway, alternative = "less")$p.value
} else{
p.value <- NA
}
return(c(p.value = p.value, Annotated = length(list.genes.in.pathway)))
}
))
# Assemble output table
outdat <- data.frame(pathway.code = rownames(pVals.by.pathway))
outdat$pathway.name <- pathways.list[outdat$pathway.code]
outdat$p.value <- pVals.by.pathway[,"p.value"]
outdat$Annotated <- pVals.by.pathway[,"Annotated"]
outdat <- outdat[order(outdat$p.value),]
head(outdat)
The output is like this pathway.code pathway.name p.value Annotated
1 dre00010 <NA> NA 0
2 dre00020 <NA> NA 0
3 dre00030 <NA> NA 0
4 dre00040 <NA> NA 0
and what is expected is
pathway.code pathway.name p.value Annotated
109 ath03010 <NA> 5.378589e-35 301
96 ath00966 <NA> 1.517953e-06 24
131 ath04141 <NA> 5.032367e-06 189
10 ath00071 <NA> 8.804424e-06 45
110 ath03013 <NA> 2.771177e-04 155
56 ath00592 <NA> 1.097775e-03 3
8
the expected data is built with respect to the Arabidopsis database and my database of interest is Zebrafish. Preferably pathway name inclusion will also be helpful.
• 1,136 views
•
link
0 answers
No answers yet.
Log in to answer this question.
Hello aranyak111!
Despite being asked to format your post multiple times in the past, this post is still not formatted. Please format it properly and your post will be reopened.
WIll you please provide some outlines about how I can format my posts? I do the best I can.
Your edits look good, so I'll open the post.