Hi Dick, Glad to help! I have a newer version of the function, since the old one uses probe set names as unique identifiers. Dr. Gentleman pointed out to me that this is not appropriate since you end up double counting some. The function below takes care of that. Again, it was modified from GOHyperG and I've added a new "KEGG" option. The input 'x' is a set of Locus Link ID numbers. Please let me know if any problems. Thanks, ...Tao ##================================================================= GOKEGGHyperG <- function (x, lib = "hgu95av2", what = "MF") { getDataEnv <- function(name, lib) { get(paste(lib, name, sep = ""), mode = "environment") } require(lib, character.only = TRUE) || stop("need data package", lib) if (any(duplicated(x))) stop("input IDs must be unique") ################ MODIFIED ################# match.arg(what, c("MF", "BP", "CC", "KEGG")) ############################################ cLLs <- unlist(as.list(getDataEnv("LOCUSID", lib))) ourLLs <- cLLs[match(x, cLLs)] ################## MODIFIED ########################### if (what=="KEGG") { goV <- as.list(getDataEnv("PATH2PROBE", lib)) } else { goV <- as.list(getDataEnv("GO2ALLPROBES", lib)) } ###################################################### whWeHave <- sapply(goV, function(y) { if is.na(y) || length(y) == 0) return(FALSE) lls = unique(unlist(mget(y, getDataEnv("LOCUSID", lib), ifnotfound = NA))) any(x %in% lls) }) goV <- goV[whWeHave] ################## MODIFIED ######################### if (what!="KEGG") { goCat <- unlist(getGOOntology(names(goV))) goodGO <- goCat == what mm <- match(names(goCat), names(goV)) mm <- mm[goodGO] goV <- goV[mm] } ###################################################### goVcts = sapply(goV, function(x) { if (length(x) == 0 || is.na(x)) return(NA) lls <- unique(unlist(mget(x, getDataEnv("LOCUSID", lib)))) lls <- lls[!is.na(lls)] lls }) bad <- sapply(goVcts, function(x) (length(x) == 1 && is.na(x))) goVcts = goVcts[!bad] goV = goV[!bad] cLLs <- unique(unlist(goVcts)) nLL <- length(cLLs) goCounts <- sapply(goVcts, length) ourLLs <- unique(ourLLs[!is.na(ourLLs)]) ours <- ourLLs[!duplicated(ourLLs)] whGood <- ours[ours %in% cLLs] nInt = length(whGood) if (nInt == 0) warning("no interesting genes found") useCts <- sapply(goVcts, function(x) sum(whGood %in% x)) pvs <- phyper(useCts - 1, nInt, nLL - nInt, goCounts, lower.tail = FALSE) ord <- order(pvs) ################################ MODIFIED ################################################## if (what!="KEGG") { return(list(pvalues = pvs[ord], goCounts = goCounts[ord], chip = lib, go2Affy = goV, intCounts = useCts[ord], numLL = nLL, numInt = nInt, intLLs = x)) } else { return(list(pvalues = pvs[ord], keggCounts = goCounts[ord], chip = lib, kegg2Affy = goV, intCounts = useCts[ord], numLL = nLL, numInt = nInt, intLLs = x)) } ################################ MODIFIED ################################################## } #===================================================================== --- Dick Beyer <dbeyer at="" u.washington.edu=""> wrote: > Hi Tao, > > I tried out your KEGGHyperG function. Seems to work just great. Thanks very much, > Dick > > ******************************************************************** *********** > Richard P. Beyer, Ph.D. University of Washington > Tel.:(206) 616 7378 Env. & Occ. Health Sci. , Box 354695 > Fax: (206) 685 4696 4225 Roosevelt Way NE, # 100 > Seattle, WA 98105-6099 > http://depts.washington.edu/ceeh/ServiceCores/FC5/FC5.html > http://staff.washington.edu/~dbeyer > ******************************************************************** *********** > > --- "Shi, Tao" <shidaxia at="" yahoo.com=""> wrote: > > Message: 8 > Date: Tue, 6 Sep 2005 10:27:05 -0700 (PDT) > From: "Shi, Tao" <shidaxia at="" yahoo.com=""> > Subject: Re: [BioC] GOHyperG for KEGG > To: Gunnar Wrobel <bioc at="" gunnarwrobel.de=""> > Cc: bioconductor at stat.math.ethz.ch > Message-ID: <20050906172705.66173.qmail at web52703.mail.yahoo.com> > Content-Type: text/plain; charset=iso-8859-1 > > Thank you very much, Gunnar. I'll try that. > > At same time, I wrote a function by myself, which I totally stole from GOHyperG. > Just want to > share it with everybody. Please let me know if there are any bugs! > > ...Tao > > ==================================================================== ========= > > KEGGHyperG <- > function (geneIDs, lib = "hgu95av2") { > getDataEnv <- function(name, lib) { > get(paste(lib, name, sep = ""), mode = "environment") > } > require(lib, character.only = TRUE) || stop("need data package", lib) > if (any(duplicated(geneIDs))) stop("input IDs must be unique") > keggV <- as.list(getDataEnv("PATH2PROBE", lib)) > > whWeHave <- sapply(keggV, function(y) { > if is.na(y) || length(y) == 0) > return(FALSE) > ids = unique(unlist(y)) > any(geneIDs %in% ids) > }) > > keggV <- keggV[whWeHave] > keggV <- sapply(keggV, function(x) { > if(any(grep("AFFX",x))) { > return(x[-grep("AFFX",x)]) > } else { > return(x) > } > } ) ## get rid of control probes > > bad <- sapply(keggV, function(x) (length(x) == 1 && is.na(x))) > keggV <- keggV[!bad] > cIDs <- unique(unlist(keggV)) > nIDs <- length(cIDs) > keggCounts <- sapply(keggV, length) > ourIDs <- unique(geneIDs[!is.na(geneIDs)]) > ours <- ourIDs[!duplicated(ourIDs)] > whGood <- ours[ours %in% cIDs] > > nInt = length(whGood) > if (nInt == 0) { warning("no interesting genes found") } > useCts <- sapply(keggV, function(x) sum(whGood %in% x)) > > pvs <- phyper(useCts - 1, nInt, nIDs - nInt, keggCounts, lower.tail = FALSE) > ord <- order(pvs) > return(list(pvalues = pvs[ord], keggCounts = keggCounts[ord], > chip = lib, kegg2Affy = keggV, intCounts = useCts[ord], numIDs = nIDs, > numInt = nInt, intIDs = geneIDs)) > } > ==================================================================== ============ > == > > > > > --- Gunnar Wrobel <bioc at="" gunnarwrobel.de=""> wrote: > > > > Is there a similar function like GOHyperG that works on KEGG? It seems > there is no such thing > > > back in Feb. 05 > (https://stat.ethz.ch/pipermail/bioconductor/2005-February/007532.ht ml). Any > > > updates? > > Hi Tao, > > > > you might try to do this with goCluster. It does the same kind of > > calculation as GOHyperG but can use any kind of annotation. > > > > Cheers > > > > Gunnar > > > > > > > > >