microarrayAnalysis.R
In ReportingTools: Tools for making reports in various formats

### R code from vignette source '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")
## publish(fit, deReport, eSet=ALL, factor=ALL$mol.biol, coef=2, n=100)
## finish(deReport)


###################################################
### code chunk number 5: making the DE report with new images (eval = FALSE)
###################################################
## library(hwriter)
## makeNewImages <- function(df,...){
## 	imagename <- c()
## 	for (i in 1:nrow(df)){
## 		probeId <- df$ProbeId[i]
## 		y_at <- pretty(exprs(ALL)[probeId,])
## 		y_labels <- formatC(y_at, digits = 1, format = 'f')
## 		imagename[i] <- paste0("plot", probeId, ".png")
## 		png(filename = paste0("./reports/figuresde_analysis/", 
##                       imagename[i]))
## 		print(stripplot(exprs(ALL)[probeId,]~ALL$mol.biol|ALL$sex))
## 		dev.off()
## 	}
## 	df$Image <- hwriteImage(paste0("figuresde_analysis/", imagename), 
##                                 link=paste0("figuresde_analysis/", imagename), 
##                                 table=FALSE, width=100)
## 	return(df)
## }
## deReport2 <- HTMLReport(shortName='de_analysis2',
##                         title = 'Analysis of BCR/ABL translocation differential expression with new plots',
##                         reportDirectory = "./reports")
## publish(fit, deReport2, eSet = ALL, factor = ALL$mol.biol, coef=2, 
##         n=100, 
##        ##.modifyDF=list(modifyReportDF, makeNewImages) ) ##to add new images to default RT output
##        .modifyDF=makeNewImages)
## finish(deReport2)


###################################################
### code chunk number 6: Do GO analysis (eval = FALSE)
###################################################
## library(GOstats)
## tt <- topTable(fit, coef = 2, n = 100)
## selectedIDs <- unlist(mget(rownames(tt), 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 7: make the GO report (eval = FALSE)
###################################################
## goReport <- HTMLReport(shortName = 'go_analysis',
##                        title = 'GO analysis of BCR/ABL translocation',
##                        reportDirectory = "./reports")
## publish(goResults, goReport)
## finish(goReport)


###################################################
### code chunk number 8: 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 9: make the PFAM report (eval = FALSE)
###################################################
## PFAMReport <- HTMLReport(shortName = 'pfam_analysis',
##                          title = 'PFAM analysis of BCR/ABL translocation',
##                          reportDirectory = "./reports")
## publish(PFAMResults, PFAMReport, categorySize = 3)
## finish(PFAMReport)


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


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


###################################################
### code chunk number 15: make the GSEA report with new columns (eval = FALSE)
###################################################
## runGSEA <- function(df,...){
##   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)	
##   lm <- lmPerGene(ALL, ~mol.biol+sex, na.rm=TRUE)
##   GSNorm <- GSNormalize(lm$tstat[2,], mat)
##   pVals <- gsealmPerm(ALL,~mol.biol+sex, mat, nperm = 100)  
##   bestPval <- apply(pVals,1, min)
##   df <- cbind(df, GSNorm, bestPval)
##   return(df)
## }
## 
## gseaReport2 <- HTMLReport(shortName = "gsea_analysis2",
##                           title = "GSEA analysis", 
##                           reportDirectory = "./reports")
## publish(geneSets, gseaReport2, annotation.db = "org.Hs.eg", 
##         .modifyDF = runGSEA)
## finish(gseaReport2)


###################################################
### code chunk number 16: make the index page (eval = FALSE)
###################################################
## indexPage <- HTMLReport(shortName = "index",
##                         title = "Analysis of ALL Gene Expression",
##                         reportDirectory = "./reports")
## publish(Link(list(deReport, goReport), report = indexPage), indexPage)
## publish(Link(PFAMReport, report = indexPage), indexPage)
## publish(Link("GSEA report has a new title", "gsea_analysis.html"), indexPage)
## finish(indexPage)