This is a test version of Biostars. For the public version, visit https://www.biostars.org.
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.

rna-seq

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.

0 answers

No answers yet.

Log in to answer this question.