Hi, I am trying to use Gostats to calculate enrichment or depletion of
GO terms. The below is the code adapted from vignette. Assume we start
from a mouse 4302 expressionset object "exprSet1000". It always returns
empty table no matter what data I try.

Can somebody help me out?

Thanks.

## General Filtering, remove probe sets without Entrez ID
entrezIds <- mget(featureNames(exprSet1000), envir = mouse4302ENTREZID)
haveEntrezId <- names(entrezIds)[sapply(entrezIds, function(x)
!is.na(x))]
numNoEntrezId <- length(featureNames(exprSet1000)) -
length(haveEntrezId)
bcrAblOrNeg <- exprSet1000[haveEntrezId, ]

## Remove probe sets for which we have no GO annotation.
haveGo <- sapply(mget(featureNames(bcrAblOrNeg), mouse4302GO),
function(x) {
if (length(x) == 1 && is.na(x))
FALSE
else TRUE
})
numNoGO <- sum(!haveGo)
bcrAblOrNeg <- bcrAblOrNeg[haveGo, ]

## remove redundant probe set ids
bcrAblOrNegIqr <- apply(exprs(bcrAblOrNeg), 1, IQR)
numNsWithDups <- length(featureNames(bcrAblOrNeg))
uniqGenes <- findLargest(featureNames(bcrAblOrNeg), bcrAblOrNegIqr,
"mouse4302")
bcrAblOrNeg <- bcrAblOrNeg[uniqGenes, ]

entrezUniverse <- unlist(mget(featureNames(bcrAblOrNeg),
mouse4302ENTREZID))
if (any(duplicated(entrezUniverse))) stop("error in gene universe: can't
have duplicate Entrez")  
selectedEntrezIds <- unlist(mget(featureNames(bcrAblOrNeg),
mouse4302ENTREZID))

hgCutoff <- 1
params <- new("GOHyperGParams", geneIds = selectedEntrezIds, 
universeGeneIds = entrezUniverse, annotation = "mouse4302",
ontology = "BP", pvalueCutoff = hgCutoff, conditional = FALSE,
testDirection = "over")

paramsCond <- params
conditional(paramsCond) <- TRUE
hgOver <- hyperGTest(params)
hgCondOver <- hyperGTest(paramsCond)



-------------------------------------------------
Dongxiao Zhu
Stowers Institute for Medical Research
Phone: 816-926-4456
FAX: 816-926-4674


	[[alternative HTML version deleted]]

