### R code from vignette source 'vignettes/ReportingTools/inst/doc/microarrayAnalysis.Rnw'

###################################################
### code chunk number 1: load ALL data (eval = FALSE)
###################################################
## library(ReportingTools)
## library(ALL)
## library(hgu95av2.db)
## library(genefilter)
## 
## data(ALL)


###################################################
### code chunk number 2: filter ALL data (eval = FALSE)
###################################################
## 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: limma linear model (eval = FALSE)
###################################################
## library(limma)
## 
## model <- model.matrix(~mol.biol+sex, ALL)
## fit <- eBayes(lmFit(ALL, model))


###################################################
### code chunk number 4: making the DE report (eval = FALSE)
###################################################
## library(lattice)
## rep.theme <- reporting.theme()
## lattice.options(default.theme = rep.theme)
## 
## deReport <- HTMLReport(shortName = 'de_analysis',
##     title = 'Analysis of BCR/ABL translocation differential expression',
##     reportDirectory = "./reports/",
##     baseUrl = "")
## publish(fit, deReport, eSet=ALL, factor=ALL$mol.biol, coef=2, n=100)
## finish(deReport)


###################################################
### code chunk number 5: Do GO analysis (eval = FALSE)
###################################################
## library(GOstats)
## tt <- topTable(fit, coef = 2, n = 100)
## selectedIDs <- unlist(mget(tt$ID, hgu95av2ENTREZID))
## universeIDs <- unlist(mget(featureNames(ALL), hgu95av2ENTREZID))
## goParams <- new("GOHyperGParams", 
##     geneIds = selectedIDs, 
##     universeGeneIds = universeIDs, 
##     annotation = annotation(ALL), 
##     ontology = "BP", 
##     pvalueCutoff = 0.01,
##     conditional = TRUE, 
##     testDirection = "over")
## goResults <- hyperGTest(goParams)


###################################################
### code chunk number 6: make the GO report (eval = FALSE)
###################################################
## goReport <- HTMLReport(shortName = 'go_analysis',
##     title = 'GO analysis of BCR/ABL translocation',
##     reportDirectory = "./reports/",
##     baseUrl = "")
## publish(goResults, goReport, selectedIDs=selectedIDs, annotation.db="org.Hs.eg")
## finish(goReport)


###################################################
### code chunk number 7: Do PFAM analysis (eval = FALSE)
###################################################
## library(Category)
## pfamParams <- new("PFAMHyperGParams", 
##     geneIds = selectedIDs, 
##     universeGeneIds = universeIDs, 
##     annotation = annotation(ALL),  
##     pvalueCutoff = 0.01,
##     testDirection = "over")
## PFAMResults <- hyperGTest(pfamParams)


###################################################
### code chunk number 8: make the PFAM report (eval = FALSE)
###################################################
## PFAMReport <- HTMLReport(shortName = 'pfam_analysis',
##     title = 'PFAM analysis of BCR/ABL translocation',
##     reportDirectory = "./reports/",
##     baseUrl = "")
## publish(PFAMResults, PFAMReport, selectedIDs=selectedIDs, annotation.db="org.Hs.eg",categorySize=3)
## finish(PFAMReport)


###################################################
### code chunk number 9: Make Gene Sets (eval = FALSE)
###################################################
## library(GSEAlm)
## library(GSEABase)
## mapped_genes <- mappedkeys(org.Hs.egSYMBOL)
## eidsAndSymbols <- as.list(org.Hs.egSYMBOL[mapped_genes])
## geneEids<-names(eidsAndSymbols)
## set.seed(123)
## set1<-GeneSet(geneIds=sample(geneEids,100, replace=FALSE), setName="set1", 
## 	shortDescription="This is set1")
## set2<-GeneSet(geneIds=sample(geneEids,10, replace=FALSE), setName="set2",  
## 	shortDescription="This is set2")
## set3<-GeneSet(geneIds=sample(geneEids,37, replace=FALSE), setName="set3",  
## 	shortDescription="This is set3")
## set4<-GeneSet(geneIds=sample(geneEids,300, replace=FALSE), setName="set4", 
## 	shortDescription="This is set4")
## geneSets<-GeneSetCollection(c(set1,set2,set3,set4))


###################################################
### code chunk number 10: Make Simple Gene Set Table (eval = FALSE)
###################################################
## geneSetsReport <- HTMLReport(shortName = "gene_sets",
## 	title = "Gene Sets", 
## 	reportDirectory = "./reports",
## 	baseUrl = "")
## publish(geneSets, geneSetsReport, annotation.db="org.Hs.eg")
## finish(geneSetsReport)


###################################################
### code chunk number 11: Get incidence matrix (eval = FALSE)
###################################################
## mat <- matrix(data=0, ncol=length(universeIDs),nrow=length(geneSets))
## for(i in 1:length(geneSets)){
## 	geneIdEntrez<-unlist(geneIds(geneSets[[i]]))
## 	mat[i,match(geneIdEntrez, universeIDs)] <- 1
## }
## colnames(mat) <- universeIDs
## rownames(mat) <- sapply(geneSets, function(x) x@setName)


###################################################
### code chunk number 12: Run GSEA (eval = FALSE)
###################################################
## lm<-lmPerGene(ALL, ~mol.biol+sex, na.rm=TRUE)
## GSNorm<-GSNormalize(lm$tstat[2,], mat)
## #one-sided p-values
## pVals<-gsealmPerm(ALL,~mol.biol+sex, mat, nperm=100)  
## bestPval<-apply(pVals,1, min)


###################################################
### code chunk number 13: make the GSEA report (eval = FALSE)
###################################################
## gseaReport <- HTMLReport(shortName = "gsea_analysis",
## 	title = "GSEA analysis", 
## 	reportDirectory = "./reports",
## 	baseUrl = "")
## publish(geneSets, gseaReport, annotation.db="org.Hs.eg", setStats=GSNorm, setPValues=2*bestPval)
## finish(gseaReport)


###################################################
### code chunk number 14: make the index page (eval = FALSE)
###################################################
## indexPage <- HTMLReport(shortName = "index",
##     title = "Analysis of ALL Gene Expression",
##     reportDirectory = "./reports",
##     baseUrl = "")
## publish(c(deReport, goReport, PFAMReport, gseaReport), indexPage)
## finish(indexPage)