################################################### ### chunk number 1: Ropts ################################################### options(width=70) ################################################### ### chunk number 2: setup ################################################### library(HTSanalyzeR) library(GSEABase) library(cellHTS2) library(org.Dm.eg.db) library(GO.db) library(KEGG.db) ################################################### ### chunk number 3: cellHTSRead ################################################### experimentName <- "KcViab" dataPath <- system.file(experimentName, package = "cellHTS2") x <- readPlateList("Platelist.txt", name = experimentName, path = dataPath,verbose=TRUE) ################################################### ### chunk number 4: cellHTSConfig ################################################### x <- configure(x, descripFile = "Description.txt", confFile = "Plateconf.txt", logFile = "Screenlog.txt", path = dataPath) ################################################### ### chunk number 5: cellHTSNorm ################################################### xn <- normalizePlates(x, scale = "multiplicative", log = FALSE, method = "median", varianceAdjust = "none") ################################################### ### chunk number 6: cellHTSNorm ################################################### xn <- annotate(xn, geneIDFile = "GeneIDs_Dm_HFA_1.1.txt", path = dataPath) ################################################### ### chunk number 7: cellHTSScore ################################################### xsc <- scoreReplicates(xn, sign = "-", method = "zscore") ################################################### ### chunk number 8: cellHTSSum ################################################### xsc <- summarizeReplicates(xsc, summary = "mean") ################################################### ### chunk number 9: showXSC ################################################### xsc ################################################### ### chunk number 10: getDataFromCellHTS ################################################### data4enrich <- as.vector(Data(xsc)) names(data4enrich) <- fData(xsc)[, "GeneID"] data4enrich <- data4enrich[which(!is.na(names(data4enrich)))] ################################################### ### chunk number 11: selectHits ################################################### hits <- names(data4enrich)[which(abs(data4enrich) > 2)] ################################################### ### 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) ################################################### ### 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) ################################################### ### chunk number 14: gscaAnalyses ################################################### gsca<-analyze(gsca, para=list(pValueCutoff=0.05, pAdjustMethod ="BH", nPermutations=100, minGeneSetSize=180, exponent=1)) ################################################### ### chunk number 15: eval=FALSE ################################################### ## library(snow) ## options(cluster=makeCluster(4, "SOCK")) ################################################### ### chunk number 16: eval=FALSE ################################################### ## if(is(getOption("cluster"), "cluster")) { ## stopCluster(getOption("cluster")) ## options(cluster=NULL) ## } ################################################### ### chunk number 17: gscaSum ################################################### summarize(gsca) ################################################### ### chunk number 18: selectSigGS ################################################### topGS_GO_MF <- getTopGeneSets(gsca, "GSEA.results", c("GO_MF", "PW_KEGG"), allSig=TRUE) ################################################### ### chunk number 19: printTopGSs ################################################### topGS_GO_MF ################################################### ### chunk number 20: viewGSEARandWalk ################################################### viewGSEA(gsca, "GO_MF", topGS_GO_MF[["GO_MF"]][1]) ################################################### ### chunk number 21: eval=FALSE ################################################### ## plotGSEA(gsca, gscs=c("GO_BP","GO_MF","PW_KEGG"), ## ntop=1, filepath=".") ################################################### ### 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") ################################################### ### 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") ################################################### ### 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") ################################################### ### 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") ################################################### ### chunk number 26: 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) ################################################### ### 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") ################################################### ### chunk number 28: gscaReportStruct eval=FALSE ################################################### ## print(dir("HTSanalyzerGSCAReport",recursive=TRUE)) ################################################### ### chunk number 29: gscaSaveAndLoad eval=FALSE ################################################### ## save(gsca, file="./gsca.RData") ## load(file="./gsca.RData") ################################################### ### 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) ################################################### ### chunk number 31: nwaInit2 eval=FALSE ################################################### ## data("Biogrid_DM_Interactome") ## nwa <- new("NWA", pvalues=pvalues, interactome= ## Biogrid_DM_Interactome, phenotypes=data4enrich) ################################################### ### chunk number 32: nwaInit ################################################### nwa <- new("NWA", pvalues=pvalues, phenotypes=data4enrich) nwa <- preprocess(nwa, species="Dm", initialIDs="FlybaseCG", keepMultipleMappings=TRUE, duplicateRemoverMethod="max") ################################################### ### chunk number 33: getInteractome2 eval=FALSE ################################################### ## nwa<-interactome(nwa, species="Dm", reportDir="HTSanalyzerReport", ## genetic=FALSE) ################################################### ### chunk number 34: getInteractome ################################################### data("Biogrid_DM_Mat") nwa<-interactome(nwa, interactionMatrix=Biogrid_DM_Mat, genetic=FALSE) ################################################### ### chunk number 35: showNWAGraph ################################################### nwa@interactome ################################################### ### chunk number 36: fitBUMplot ################################################### nwa<-analyze(nwa, fdr=0.001, species="Dm") ################################################### ### chunk number 37: nwaSum ################################################### summarize(nwa) ################################################### ### chunk number 38: viewSubNet ################################################### viewSubNet(nwa) ################################################### ### chunk number 39: eval=FALSE ################################################### ## plotSubNet(nwa, filepath=".", filename="subnetwork.png") ################################################### ### chunk number 40: eval=FALSE ################################################### ## report(object=nwa, experimentName=experimentName, species="Dm", ## allSig=TRUE, keggGSCs="PW_KEGG", goGSCs=c("GO_BP", "GO_MF"), ## reportDir="HTSanalyzerNWReport") ################################################### ### chunk number 41: eval=FALSE ################################################### ## reportAll(gsca=gsca, nwa=nwa, experimentName=experimentName, ## species="Dm", allSig=TRUE, keggGSCs="PW_KEGG", goGSCs= ## c("GO_BP", "GO_MF"), reportDir="HTSanalyzerReport") ################################################### ### chunk number 42: eval=FALSE ################################################### ## save(nwa, file="./nwa.RData") ## load("./nwa.RData") ################################################### ### chunk number 43: 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) ################################################### ### chunk number 44: 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" ## ) ################################################### ### chunk number 45: eval=FALSE ################################################### ## c2<-getGmt(con="c2.all.v2.5.symbols.gmt.txt",geneIdType= ## SymbolIdentifier(), collectionType= ## BroadCollection(category="c2")) ################################################### ### chunk number 46: eval=FALSE ################################################### ## c2entrez<-mapIdentifiers(c2, EntrezIdentifier('org.Hs.eg.db')) ################################################### ### chunk number 47: eval=FALSE ################################################### ## collectionOfGeneSets<-geneIds(c2entrez) ## names(collectionOfGeneSets)<-names(c2entrez) ################################################### ### chunk number 48: sessionInfo ################################################### toLatex(sessionInfo())