inst/doc/HTSanalyzeR-Vignette.R

In HTSanalyzeR: Gene set over-representation, enrichment and network analyses for high-throughput screens

### R code from vignette source 'HTSanalyzeR-Vignette.Rnw'

###################################################
### code chunk number 1: Ropts
###################################################
options(width=70)


###################################################
### code chunk number 2: setup
###################################################
library(HTSanalyzeR)
library(GSEABase)
library(cellHTS2)
library(org.Dm.eg.db)
library(GO.db)
library(KEGG.db)


###################################################
### code chunk number 3: cellHTSRead
###################################################
experimentName <- "KcViab"
dataPath <- system.file(experimentName, package = "cellHTS2")
x <- readPlateList("Platelist.txt", name = experimentName, 
path = dataPath,verbose=TRUE)


###################################################
### code chunk number 4: cellHTSConfig
###################################################
x <- configure(x, descripFile = "Description.txt", confFile = 
"Plateconf.txt", logFile = "Screenlog.txt", path = dataPath)


###################################################
### code chunk number 5: cellHTSNorm
###################################################
xn <- normalizePlates(x, scale = "multiplicative", log = FALSE, 
method = "median", varianceAdjust = "none")


###################################################
### code chunk number 6: cellHTSNorm
###################################################
xn <- annotate(xn, geneIDFile = "GeneIDs_Dm_HFA_1.1.txt", 
path = dataPath)


###################################################
### code chunk number 7: cellHTSScore
###################################################
xsc <- scoreReplicates(xn, sign = "-", method = "zscore")


###################################################
### code chunk number 8: cellHTSSum
###################################################
xsc <- summarizeReplicates(xsc, summary = "mean")	


###################################################
### code chunk number 9: showXSC
###################################################
xsc


###################################################
### code chunk number 10: getDataFromCellHTS
###################################################
data4enrich <- as.vector(Data(xsc))
names(data4enrich) <- fData(xsc)[, "GeneID"]
data4enrich <- data4enrich[which(!is.na(names(data4enrich)))]


###################################################
### code chunk number 11: selectHits
###################################################
hits <- names(data4enrich)[which(abs(data4enrich) > 2)]


###################################################
### code chunk number 12: gscList
###################################################
GO_MF <- GOGeneSets(species="Dm", ontologies=c("MF"))
GO_BP <- GOGeneSets(species="Dm", ontologies=c("BP"))
PW_KEGG <- KeggGeneSets(species="Dm")
ListGSC <- list(GO_MF=GO_MF, GO_BP=GO_BP, PW_KEGG=PW_KEGG)


###################################################
### code chunk number 13: gscaInit
###################################################
gsca <- new("GSCA", listOfGeneSetCollections=ListGSC, 
geneList=data4enrich, hits=hits)
gsca <- preprocess(gsca, species="Dm", initialIDs="FlybaseCG", 
keepMultipleMappings=TRUE, duplicateRemoverMethod="max", 
orderAbsValue=FALSE)


###################################################
### code chunk number 14: gscaAnalyses
###################################################
gsca<-analyze(gsca, para=list(pValueCutoff=0.05, pAdjustMethod
="BH", nPermutations=100, minGeneSetSize=180, exponent=1))


###################################################
### code chunk number 15: HTSanalyzeR-Vignette.Rnw:214-216 (eval = FALSE)
###################################################
## library(snow)
## options(cluster=makeCluster(4, "SOCK"))


###################################################
### code chunk number 16: HTSanalyzeR-Vignette.Rnw:221-225 (eval = FALSE)
###################################################
## if(is(getOption("cluster"), "cluster")) {
## 	stopCluster(getOption("cluster"))
## 	options(cluster=NULL)
## }


###################################################
### code chunk number 17: gscaSum
###################################################
summarize(gsca)


###################################################
### code chunk number 18: selectSigGS
###################################################
topGS_GO_MF <- getTopGeneSets(gsca, "GSEA.results", c("GO_MF", 
"PW_KEGG"), allSig=TRUE)


###################################################
### code chunk number 19: printTopGSs
###################################################
topGS_GO_MF


###################################################
### code chunk number 20: viewGSEARandWalk
###################################################
viewGSEA(gsca, "GO_MF", topGS_GO_MF[["GO_MF"]][1])


###################################################
### code chunk number 21: HTSanalyzeR-Vignette.Rnw:263-265 (eval = FALSE)
###################################################
## plotGSEA(gsca, gscs=c("GO_BP","GO_MF","PW_KEGG"), 
## ntop=1, filepath=".")


###################################################
### code chunk number 22: viewEnrichMap1_a
###################################################
data("KcViab_GSCA")
viewEnrichMap(KcViab_GSCA, resultName="HyperGeo.results", 
gscs=c("PW_KEGG"), allSig=FALSE, ntop=30, gsNameType="id", 
displayEdgeLabel=FALSE, layout="layout.fruchterman.reingold")


###################################################
### code chunk number 23: viewEnrichMap1_b
###################################################
data("KcViab_GSCA")
viewEnrichMap(KcViab_GSCA, resultName="GSEA.results", 
gscs=c("PW_KEGG"), allSig=FALSE, ntop=30, gsNameType="id", 
displayEdgeLabel=FALSE, layout="layout.fruchterman.reingold")


###################################################
### code chunk number 24: viewEnrichMap2_a
###################################################
KcViab_GSCA<-appendGSTerms(KcViab_GSCA, goGSCs=c("GO_BP",
"GO_MF","GO_CC"), keggGSCs=c("PW_KEGG"))
viewEnrichMap(KcViab_GSCA, resultName="HyperGeo.results",
gscs=c("PW_KEGG"), allSig=FALSE, ntop=30, gsNameType="term", 
displayEdgeLabel=FALSE, layout="layout.fruchterman.reingold")


###################################################
### code chunk number 25: viewEnrichMap2_b
###################################################
KcViab_GSCA<-appendGSTerms(KcViab_GSCA, goGSCs=c("GO_BP",
"GO_MF","GO_CC"), keggGSCs=c("PW_KEGG"))
viewEnrichMap(KcViab_GSCA, resultName="GSEA.results", 
gscs=c("PW_KEGG"), allSig=FALSE, ntop=30, gsNameType="term", 
displayEdgeLabel=FALSE, layout="layout.fruchterman.reingold")


###################################################
### code chunk number 26: HTSanalyzeR-Vignette.Rnw:335-339 (eval = FALSE)
###################################################
## plotEnrichMap(KcViab_GSCA, gscs=c("PW_KEGG"), allSig=TRUE, 
## ntop=NULL, gsNameType="id", displayEdgeLabel=FALSE,
## layout="layout.fruchterman.reingold", filepath=".", 
## filename="PW_KEGG.map.pdf",output="pdf", width=8, height=8)


###################################################
### code chunk number 27: gscaReport (eval = FALSE)
###################################################
## report(object=gsca, experimentName=experimentName, species="Dm", 
## allSig=TRUE, keggGSCs="PW_KEGG", goGSCs=c("GO_BP", "GO_MF"), 
## reportDir="HTSanalyzerGSCAReport")


###################################################
### code chunk number 28: gscaReportStruct (eval = FALSE)
###################################################
## print(dir("HTSanalyzerGSCAReport",recursive=TRUE))


###################################################
### code chunk number 29: gscaSaveAndLoad (eval = FALSE)
###################################################
## save(gsca, file="./gsca.RData")
## load(file="./gsca.RData")


###################################################
### code chunk number 30: nwaGetPval
###################################################
test.stats <- cellHTS2OutputStatTests(cellHTSobject=xn, 
annotationColumn="GeneID", alternative="two.sided", 
tests=c("T-test"))
library(BioNet)
pvalues <- aggrPvals(test.stats, order=2, plot=FALSE)


###################################################
### code chunk number 31: nwaInit2 (eval = FALSE)
###################################################
## data("Biogrid_DM_Interactome")
## nwa <- new("NWA", pvalues=pvalues, interactome=
## Biogrid_DM_Interactome, phenotypes=data4enrich)


###################################################
### code chunk number 32: nwaInit
###################################################
nwa <- new("NWA", pvalues=pvalues, phenotypes=data4enrich)
nwa <- preprocess(nwa, species="Dm", initialIDs="FlybaseCG", 
keepMultipleMappings=TRUE, duplicateRemoverMethod="max")


###################################################
### code chunk number 33: getInteractome2 (eval = FALSE)
###################################################
## nwa<-interactome(nwa, species="Dm", reportDir="HTSanalyzerReport",
## genetic=FALSE)


###################################################
### code chunk number 34: getInteractome
###################################################
data("Biogrid_DM_Mat")
nwa<-interactome(nwa, interactionMatrix=Biogrid_DM_Mat,
genetic=FALSE)


###################################################
### code chunk number 35: showNWAGraph
###################################################
nwa@interactome


###################################################
### code chunk number 36: fitBUMplot
###################################################
nwa<-analyze(nwa, fdr=0.001, species="Dm")


###################################################
### code chunk number 37: nwaSum
###################################################
summarize(nwa)


###################################################
### code chunk number 38: viewSubNet
###################################################
viewSubNet(nwa)


###################################################
### code chunk number 39: HTSanalyzeR-Vignette.Rnw:470-471 (eval = FALSE)
###################################################
## plotSubNet(nwa, filepath=".", filename="subnetwork.png")


###################################################
### code chunk number 40: HTSanalyzeR-Vignette.Rnw:478-481 (eval = FALSE)
###################################################
## report(object=nwa, experimentName=experimentName, species="Dm", 
## allSig=TRUE, keggGSCs="PW_KEGG", goGSCs=c("GO_BP", "GO_MF"), 
## reportDir="HTSanalyzerNWReport")


###################################################
### code chunk number 41: HTSanalyzeR-Vignette.Rnw:486-489 (eval = FALSE)
###################################################
## reportAll(gsca=gsca, nwa=nwa, experimentName=experimentName, 
## species="Dm", allSig=TRUE, keggGSCs="PW_KEGG", goGSCs=
## c("GO_BP", "GO_MF"), reportDir="HTSanalyzerReport")


###################################################
### code chunk number 42: HTSanalyzeR-Vignette.Rnw:494-496 (eval = FALSE)
###################################################
## save(nwa, file="./nwa.RData")
## load("./nwa.RData")


###################################################
### code chunk number 43: HTSanalyzeR-Vignette.Rnw:505-509 (eval = FALSE)
###################################################
## data("KcViab_Norm")
## GO_CC<-GOGeneSets(species="Dm",ontologies=c("CC"))
## PW_KEGG<-KeggGeneSets(species="Dm")
## ListGSC<-list(GO_CC=GO_CC,PW_KEGG=PW_KEGG)


###################################################
### code chunk number 44: HTSanalyzeR-Vignette.Rnw:514-526 (eval = FALSE)
###################################################
## HTSanalyzeR4cellHTS2(
## 	normCellHTSobject=KcViab_Norm,
## 	annotationColumn="GeneID",
## 	species="Dm",
## 	initialIDs="FlybaseCG",
## 	listOfGeneSetCollections=ListGSC,
## 	cutoffHitsEnrichment=2,
## 	minGeneSetSize=200,
## 	keggGSCs=c("PW_KEGG"),
## 	goGSCs=c("GO_CC"),
## 	reportDir="HTSanalyzerReport"
## ) 


###################################################
### code chunk number 45: HTSanalyzeR-Vignette.Rnw:535-538 (eval = FALSE)
###################################################
## c2<-getGmt(con="c2.all.v2.5.symbols.gmt.txt",geneIdType=
## SymbolIdentifier(), collectionType=
## BroadCollection(category="c2"))


###################################################
### code chunk number 46: HTSanalyzeR-Vignette.Rnw:542-543 (eval = FALSE)
###################################################
## c2entrez<-mapIdentifiers(c2, EntrezIdentifier('org.Hs.eg.db'))


###################################################
### code chunk number 47: HTSanalyzeR-Vignette.Rnw:548-550 (eval = FALSE)
###################################################
## collectionOfGeneSets<-geneIds(c2entrez)
## names(collectionOfGeneSets)<-names(c2entrez)


###################################################
### code chunk number 48: sessionInfo
###################################################
toLatex(sessionInfo())