Entering edit mode
Pierre-Yves Boëlle
▴
10
@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