domainenrichment.R
In domainsignatures: Geneset enrichment based on InterPro domain signatures

### R code from vignette source 'domainenrichment.Rnw'

###################################################
### code chunk number 1: KeggMapping
###################################################
library(domainsignatures)
keggMap <- mget(grep("^hsa", ls(KEGGPATHID2EXTID), value=TRUE), 
                KEGGPATHID2EXTID)
univGenes <-  unique(unlist(keggMap))
randGenes <- sample(univGenes, 30)
pwLength <- listLen(keggMap)
randPw <- "hsa04330"
set.seed(123)
pwGenes <- sample(keggMap[[randPw]], 10)


###################################################
### code chunk number 2: gseaprep
###################################################
geneset <- c(randGenes, pwGenes)
universe <- getKEGGdata(univGenes)
universe


###################################################
### code chunk number 3: gseashow (eval = FALSE)
###################################################
## res <- gseDomain(dataSource=universe, geneset=geneset, n=1000)


###################################################
### code chunk number 4: gseado
###################################################
## We only run this code for testing purpos. The actual object was pre-
## computed and we only load the serialized version here.
tmp <- gseDomain(universe, geneset, n=100, verbose=FALSE)
load(system.file("precomp/precomp.rda", package="domainsignatures"))


###################################################
### code chunk number 5: volcano
###################################################
volcData <- cbind(similarity=res$similarity, log_p=-log10(res$pvalue))
plot(volcData, main="volcano plot")
sel <- which.max(res$similarity)
points(volcData[sel,,drop=FALSE], pch=20, col="red")
text(volcData[sel,,drop=FALSE], names(sel), pos=2)


###################################################
### code chunk number 6: thesame
###################################################
names(sel) == randPw
getKEGGdescription(randPw)


###################################################
### code chunk number 7: biomart
###################################################
ensembl <- useMart("ensembl",dataset="hsapiens_gene_ensembl")
cMap <- getBM(attributes=c("entrezgene", "chromosome_name"), bmHeader=FALSE,
              filters="entrezgene", value=univGenes, mart=ensembl)
cMap <- cMap[cMap$chromosome_name %in% c(as.character(1:23), "X", "Y"),]
grouping <- split(cMap$chromosome_name, cMap$entrezgene)
universe2 <- dataSource(grouping)
geneset2 <- geneset[geneset %in% names(grouping)]


###################################################
### code chunk number 8: dsigshow (eval = FALSE)
###################################################
## ## Again, the object was precomputed for perfomance reasons
## res2 <- gseDomain(universe2, geneset2)


###################################################
### code chunk number 9: vplot2
###################################################
volcData2 <- cbind(similarity=res2$similarity, log_p=-log10(res2$pvalue))
plot(volcData2, main="volcano plot", 
     xlim=c(0,max(unlist(res$similarity))))