tests/testReport.R

In ytdai/gfplots: Interactive visualization function of gene mutations and fusions

#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)