inst/doc/npGSEA.R

In npGSEA: Permutation approximation methods for gene set enrichment analysis (non-permutation GSEA)

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

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: loadLibrariesAndData
###################################################
library(ALL)
library(hgu95av2.db)
library(genefilter)
library(limma)
library(GSEABase)
library(npGSEA)

data(ALL)

ALL <- ALL[, ALL$mol.biol %in% c('NEG','BCR/ABL') &
     !is.na(ALL$sex)]
ALL$mol.biol <- factor(ALL$mol.biol, 
     levels = c('NEG', 'BCR/ABL'))
ALL <- featureFilter(ALL)


###################################################
### code chunk number 3: adjustNames
###################################################
featureNames(ALL) <- select(hgu95av2.db, featureNames(ALL), 
                            "ENTREZID", "PROBEID")$ENTREZID


###################################################
### code chunk number 4: MakeSets
###################################################
xData <- exprs(ALL)
geneEids <- rownames(xData)
set.seed(12345)
set1 <- GeneSet(geneIds=sample(geneEids,15, replace=FALSE), 
             setName="set1", 
             shortDescription="This is set1")
set2 <- GeneSet(geneIds=sample(geneEids,50, replace=FALSE), 
             setName="set2",  
             shortDescription="This is set2")
set3 <- GeneSet(geneIds=sample(geneEids,100, replace=FALSE), 
             setName="set3",  
             shortDescription="This is set3")
set4 <- GeneSet(geneIds=sample(geneEids,500, replace=FALSE), 
            setName="set4",  
            shortDescription="This is set4")


###################################################
### code chunk number 5: RunLimma2NewSets
###################################################
model <- model.matrix(~mol.biol, ALL)
fit  <-  eBayes(lmFit(ALL, model))
tt <- topTable(fit, coef=2, n=200)
ttUp <- tt[which(tt$logFC >0), ]
ttDown <- tt[which(tt$logFC <0), ]

set5 <- GeneSet(geneIds=rownames(ttUp)[1:20], 
        setName="set5",  
        shortDescription="This is a true set of the top 20 DE 
        genes with a positive fold change")
set6 <- GeneSet(geneIds=rownames(ttDown)[1:20], 
        setName="set6",  
        shortDescription="This is a true set of the top 20 DE genes 
        with a negative fold change")
set7 <- GeneSet(geneIds=c(rownames(ttUp)[1:10], rownames(ttDown)[1:10]), 
        setName="set7",  
        shortDescription="This is a true set of the top 10 DE genes 
        with a positive and a negative fold change")


###################################################
### code chunk number 6: GSC
###################################################
gsc <- GeneSetCollection( c(set1, set2, set3, set4, set5, set6, set7) )
gsc


###################################################
### code chunk number 7: runIt
###################################################
yFactor <- ALL$mol.biol
res1 <- npGSEA(x = ALL, y = yFactor, set = set1)  ##with the eset
res1

res2_exprs <- npGSEA(xData, ALL$mol.biol, gsc[[2]])  ##with the expression data
res2_exprs 


###################################################
### code chunk number 8: accessorFunctions
###################################################
res3 <- npGSEA(ALL, yFactor, set3)  
res3
geneSetName(res3)
stat(res3)
sigmaSq(res3)
zStat(res3)
pTwoSided(res3)
pLeft(res3)
pValues(res3)
dim(xSet(res3))


###################################################
### code chunk number 9: plot3
###################################################
npGSEAPlot(res3)


###################################################
### code chunk number 10: runItNorm
###################################################
res5_norm <- npGSEA(ALL, yFactor, set5, approx= "norm")  
res5_norm
betaHats(res5_norm)
npGSEAPlot(res5_norm)


###################################################
### code chunk number 11: runItBeta
###################################################
res5_beta <- npGSEA(ALL, yFactor, set5, approx= "beta")
res5_beta
betaHats(res5_beta)
npGSEAPlot(res5_beta)


###################################################
### code chunk number 12: runItChiSq
###################################################
res5_chiSq <- npGSEA(ALL, yFactor, set5, approx= "chiSq")  
res5_chiSq
betaHats(res5_chiSq)
chiSqStat(res5_chiSq)
stat(res5_chiSq)
stat(res5_chiSq)/sigmaSq(res5_chiSq)
npGSEAPlot(res5_chiSq)


###################################################
### code chunk number 13: runItwts
###################################################
res7_nowts <- npGSEA(x = ALL, y= yFactor, set = set7)  
res7_nowts

wts <- apply(exprs(ALL)[match(geneIds(set7), featureNames(ALL)), ], 
						1, var)	
wts <- 1/wts				
res7_wts <- npGSEA(x = ALL, y = yFactor, set = set7,  w = wts, approx= "norm")  
res7_wts						


###################################################
### code chunk number 14: addCovariates
###################################################
res3_age <- npGSEA(x = ALL, y = yFactor, set = set3, covars = ALL$age) 
res3_age

res3_agesex <- npGSEA(x = ALL, y = yFactor, set = set3, covars = cbind(ALL$age, ALL$sex)) 
res3_agesex


###################################################
### code chunk number 15: runMultiplenpGSEA
###################################################
resgsc_norm <- npGSEA(x = ALL, y = yFactor, set = gsc)  
unlist( pLeft(resgsc_norm) )
unlist( stat (resgsc_norm) )
unlist( zStat (resgsc_norm) )


###################################################
### code chunk number 16: publishMultiplenpGSEA (eval = FALSE)
###################################################
## pvals <- p.adjust( unlist(pTwoSided(resgsc_norm)), method= "BH" )
## library(ReportingTools)
## npgseaReport <- HTMLReport (shortName = "npGSEA",
##         title = "npGSEA Results", reportDirectory = "./reports")
## publish(gsc, npgseaReport, annotation.db = "org.Hs.eg",
##         setStats = unlist(zStat (resgsc_norm)), setPValues = pvals)
## finish(npgseaReport)


###################################################
### code chunk number 17: sessInfo
###################################################
toLatex(sessionInfo())


###################################################
### code chunk number 18: resetOptions
###################################################
options(prompt="> ", continue="+ ")