GOstats v2.2.3 . bug in inducedTermGraph
1
0
Entering edit mode
@pierre-yves-boelle-2185
Last seen 10.3 years ago
Using R : 2.5.0 / BioC 2.0 / GOStats 2.2.3 / GO : 1.16 on Windows XP. the following code reproduces the bug. (quite lengthy... I took the example from the vignette) ---------------------------------- library("ALL") library("hgu95av2") library("annotate") library("genefilter") library("GOstats") library("RColorBrewer") library("xtable") library("Rgraphviz") data(ALL, package = "ALL") Bcell <- grep("^B", as.character(ALL$BT)) bcrAblOrNegIdx <- which(as.character(ALL$mol) %in% c("NEG", "ALL1/AF4")) bcrAblOrNeg <- ALL[, intersect(Bcell, bcrAblOrNegIdx)] bcrAblOrNeg$mol.biol = factor(bcrAblOrNeg$mol.biol) entrezIds <- mget(featureNames(bcrAblOrNeg), envir = hgu95av2ENTREZID) haveEntrezId <- names(entrezIds)[sapply(entrezIds, function(x) !is.na(x))] numNoEntrezId <- length(featureNames(bcrAblOrNeg)) - length(haveEntrezId) bcrAblOrNeg <- bcrAblOrNeg[haveEntrezId, ] haveGo <- sapply(mget(featureNames(bcrAblOrNeg), hgu95av2GO), function(x) { if (length(x) == 1 && is.na(x)) FALSE else TRUE}) numNoGO <- sum(!haveGo) bcrAblOrNeg <- bcrAblOrNeg[haveGo, ] iqrCutoff <- 0.5 bcrAblOrNegIqr <- apply(exprs(bcrAblOrNeg), 1, IQR) selected <- bcrAblOrNegIqr > iqrCutoff chrN <- mget(featureNames(bcrAblOrNeg), envir = hgu95av2CHR) onY <- sapply(chrN, function(x) any(x == "Y")) onY[is.na(onY)] <- FALSE selected <- selected & !onY nsFiltered <- bcrAblOrNeg[selected, ] nsFilteredIqr <- bcrAblOrNegIqr[selected] uniqGenes <- findLargest(featureNames(nsFiltered), nsFilteredIqr, "hgu95av2") nsFiltered <- nsFiltered[uniqGenes, ] numSelected <- length(featureNames(nsFiltered)) BCRcols = ifelse(nsFiltered$mol == "ALL1/AF4", "goldenrod", "skyblue") cols = brewer.pal(10, "RdBu") affyUniverse <- featureNames(nsFiltered) entrezUniverse <- unlist(mget(affyUniverse, hgu95av2ENTREZID)) if (any(duplicated(entrezUniverse))) stop("error in gene universe: can't have duplicate EntrezIds") chipAffyUniverse <- featureNames(bcrAblOrNeg) chipEntrezUniverse <- mget(chipAffyUniverse, hgu95av2ENTREZID) chipEntrezUniverse <- unique(unlist(chipEntrezUniverse)) ttestCutoff <- 0.05 ttests = rowttests(nsFiltered, "mol.biol") smPV = ttests$p.value < ttestCutoff pvalFiltered <- nsFiltered[smPV, ] selectedEntrezIds <- unlist(mget(featureNames(pvalFiltered), hgu95av2ENTREZID)) hgCutoff <- 0.001 params <- new("GOHyperGParams", geneIds = selectedEntrezIds, universeGeneIds = entrezUniverse, annotation = "hgu95av2", ontology = "BP", pvalueCutoff = hgCutoff, conditional = FALSE, testDirection = "over") hgOver <- hyperGTest(params) inducedTermGraph(hgOver,names(pvalues(hgOver)[95]),children=FALSE) ---------------------------- fails with : Erreur dans checkValidNodeName(node) : invalid node names: missing value NA not allowed looking at the code (in GOHyperGResults-accessors.R) , it appears that parents are looked for in the "KidsEnv" instead of "ParentsEnv" Also, the graphs is reversed? (ancestor at the bottom, children on top) -- Pierre-Yves Bo?lle Universite Pierre et Marie Curie-Paris6, UMR S 707, Paris, F-75012 INSERM, UMR S 707, Paris, F-75012 AP-HP, hopital St Antoine, Paris, F-75012 Tel : +33 (0)1 49 28 32 26 Fax : +33 (0)1 49 28 32 33
Annotation GO GOstats Annotation GO GOstats • 1.1k views
ADD COMMENT
0
Entering edit mode
Seth Falcon ★ 7.4k
@seth-falcon-992
Last seen 10.3 years ago
Hi, Thanks for the nice reproducible bug report. Pierre-Yves Bo?lle <boelle at="" u707.jussieu.fr=""> writes: > Using R : 2.5.0 / BioC 2.0 / GOStats 2.2.3 / GO : 1.16 on Windows XP. > > inducedTermGraph(hgOver,names(pvalues(hgOver)[95]),children=FALSE) > > ---------------------------- > fails with : > Erreur dans checkValidNodeName(node) : invalid node names: missing value > NA not allowed > > looking at the code (in GOHyperGResults-accessors.R) , it appears that > parents are looked for in the "KidsEnv" instead of "ParentsEnv" Yes, that was part of the bug. Also, the code needed to handle NA's explicitly since this could occur for children=TRUE. Below is a new version of inducedTermGraph which you can test and which I plan to include in a patch release for GOstats. > Also, the graphs is reversed? (ancestor at the bottom, children on > top) The edges of the returned graph connect children to their parents. This matches the behavior of the GOGraph function as well as matches a convention in UML of edges going from subclasses to parent classes. As for top/bottom, that is a display issue and I believe there are ways of asking graphviz to flip things if that is what you want. Since these graphs are tree-like, one could argue that rendering a tree with root at the bottom and leaves at the top has a certain appeal in that it matches the way most trees choose to orient themselves ;-) Best, + seth Here is a new version of the function to test out: inducedTermGraph <- function(r, id, children=TRUE, parents=TRUE) { if (!children && !parents) stop("children and parents can't both be FALSE") ## XXX: should use more structure here goName <- paste(testName(r), collapse="") goKidsEnv <- get(paste(goName, "CHILDREN", sep="")) goParentsEnv <- get(paste(goName, "PARENTS", sep="")) goIds <- character(0) wantedNodes <- id ## children if (children) { wantedNodes <- c(wantedNodes, unlist(edges(goDag(r))[id], use.names=FALSE)) } ## parents g <- reverseEdgeDirections(goDag(r)) if (parents) { wantedNodes <- c(wantedNodes, unlist(edges(g)[id], use.names=FALSE)) } wantedNodes <- unique(wantedNodes) g <- subGraph(wantedNodes, g) ## expand; add children and/or parents that are not present in g, ## but are definedin the GO data. if (children) { for (goid in id) { kids <- unique(goKidsEnv[[goid]]) for (k in kids) { if is.na(k)) next if (!(k %in% nodes(g))) { g <- addNode(k, g) g <- addEdge(k, goid, g) } } } } if (parents) { for (goid in id) { elders <- unique(goParentsEnv[[goid]]) for (p in elders) { if is.na(p)) next if (!(p %in% nodes(g))) { g <- addNode(p, g) g <- addEdge(goid, p, g) } } } } g } -- Seth Falcon | Computational Biology | Fred Hutchinson Cancer Research Center http://bioconductor.org
ADD COMMENT

Login before adding your answer.

Traffic: 432 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6