###################################################
### chunk number 1: load library
###################################################
#line 97 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
library(lumi)


###################################################
### chunk number 2: load example data
###################################################
#line 130 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
## load example data (a LumiBatch object)
data(example.lumi)
## summary of the example data
example.lumi


###################################################
### chunk number 3: summary of example data
###################################################
#line 145 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
## summary of the quality control
summary(example.lumi, 'QC')


###################################################
### chunk number 4: densityPlot
###################################################
#line 164 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
plot(example.lumi, what='density')	## plot the density


###################################################
### chunk number 5: densityPlot
###################################################
#line 183 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
plot(example.lumi, what='density')	## plot the density


###################################################
### chunk number 6: CDF
###################################################
#line 200 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
 plotCDF(example.lumi)


###################################################
### chunk number 7: pairPlot
###################################################
#line 220 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
plot(example.lumi, what='pair')		## pairwise plots


###################################################
### chunk number 8: 
###################################################
#line 229 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
plot(example.lumi, what='pair', smoothScatter=T)


###################################################
### chunk number 9: MAplot
###################################################
#line 250 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
## pairwise MAplot
plot(example.lumi, what='MAplot')	


###################################################
### chunk number 10: 
###################################################
#line 261 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
## pairwise MAplot
plot(example.lumi, what='MAplot', smoothScatter=T)	


###################################################
### chunk number 11: 
###################################################
#line 275 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
## density plot of coefficient of varience
plot(example.lumi, what='cv')	


###################################################
### chunk number 12: 
###################################################
#line 288 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
plot(example.lumi, what='sampleRelation')


###################################################
### chunk number 13: 
###################################################
#line 306 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
plot(example.lumi, what='sampleRelation', method='mds', color=c("01", "02", "01", "02"))


###################################################
### chunk number 14: VST transform
###################################################
#line 327 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
## Do default VST variance stabilizing transform
lumi.T <- lumiT(example.lumi)


###################################################
### chunk number 15: plotVST
###################################################
#line 343 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
trans <- plotVST(lumi.T)


###################################################
### chunk number 16: 
###################################################
#line 352 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
matplot(log2(trans$untransformed), trans$transformed, type='l', main='compare VST and log2 transform', xlab='log2 transformed', ylab='vST transformed')
abline(a=0, b=1, col=2)


###################################################
### chunk number 17: Default normalization
###################################################
#line 375 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
## Do quantile between microarray normaliazation
lumi.N <- lumiN(lumi.T)


###################################################
### chunk number 18: Customized normalization eval=FALSE
###################################################
## #line 381 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
## ## Do RSN between microarray normaliazation
## lumi.N <- lumiN(lumi.T, method='rsn')


###################################################
### chunk number 19: QC after normalization
###################################################
#line 390 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
## Do quality control estimation after normalization
lumi.N.Q <- lumiQ(lumi.N)

## summary of the quality control
summary(lumi.N.Q, 'QC')		## summary of QC


###################################################
### chunk number 20: 
###################################################
#line 401 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
plot(lumi.N.Q, what='density')		## plot the density


###################################################
### chunk number 21: 
###################################################
#line 411 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
plot(lumi.N.Q, what='boxplot')		## box plot
# boxplot(lumi.N.Q)


###################################################
### chunk number 22: 
###################################################
#line 422 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
plot(lumi.N.Q, what='pair')			## pairwise plots


###################################################
### chunk number 23: 
###################################################
#line 432 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
plot(lumi.N.Q, what='MAplot')		## plot the pairwise MAplot


###################################################
### chunk number 24: 
###################################################
#line 442 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
## plot the sampleRelation using hierarchical clustering
plot(lumi.N.Q, what='sampleRelation')


###################################################
### chunk number 25: 
###################################################
#line 453 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
## plot the sampleRelation using MDS
plot(lumi.N.Q, what='sampleRelation', method='mds', color=c("01", "02", "01", "02"))


###################################################
### chunk number 26: Encapsulate the processing steps
###################################################
#line 465 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
## Do all the default preprocessing in one step
lumi.N.Q <- lumiExpresso(example.lumi)


###################################################
### chunk number 27: 
###################################################
#line 472 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
## Do all the preprocessing with customized settings
# lumi.N.Q <- lumiExpresso(example.lumi, normalize.param=list(method='rsn'))


###################################################
### chunk number 28: inverse VST
###################################################
#line 481 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
## Parameters of VST transformed LumiBatch object
names(attributes(lumi.T))
## VST parameters: "vstParameter"  and  "transformFun"
attr(lumi.T, 'vstParameter')
attr(lumi.T, 'transformFun')
## Parameters of VST transformed and RSN normalized LumiBatch object
names(attributes(lumi.N.Q))
## VSN "targetArray" , VST parameters: "vstParameter"  and  "transformFun"
attr(lumi.N.Q, 'vstParameter')
attr(lumi.N.Q, 'transformFun')
## After doing statistical analysis of the data, users can recover to the raw scale for the fold-change estimation.
## Inverse VST to the raw scale
lumi.N.raw <- inverseVST(lumi.N.Q)


###################################################
### chunk number 29: retrieve normalized data
###################################################
#line 609 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
dataMatrix <- exprs(lumi.N.Q)


###################################################
### chunk number 30: filtering
###################################################
#line 614 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
presentCount <- detectionCall(example.lumi)
selDataMatrix <- dataMatrix[presentCount > 0,]
probeList <- rownames(selDataMatrix)


###################################################
### chunk number 31: Identify differentially expressed genes
###################################################
#line 621 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
## Specify the sample type
sampleType <- c('100:0', '95:5', '100:0', '95:5')
if (require(limma)) {
	##  compare '95:5' and '100:0'
	design <- model.matrix(~ factor(sampleType))
	colnames(design) <- c('100:0', '95:5-100:0')
	fit <- lmFit(selDataMatrix, design)
	fit <- eBayes(fit)
	## Add gene symbols to gene properties
	if (require(lumiHumanAll.db) & require(annotate)) {
               geneSymbol <- getSYMBOL(probeList, 'lumiHumanAll.db')
               geneName <- sapply(lookUp(probeList, 'lumiHumanAll.db', 'GENENAME'), function(x) x[1])
               fit$genes <- data.frame(ID= probeList, geneSymbol=geneSymbol, geneName=geneName, stringsAsFactors=FALSE)
          }
	## print the top 10 genes
	print(topTable(fit, coef='95:5-100:0', adjust='fdr', number=10))
	
	## get significant gene list with FDR adjusted p.values less than 0.01
	p.adj <- p.adjust(fit$p.value[,2])		
	sigGene.adj <- probeList[ p.adj < 0.01]
	## without FDR adjustment
	sigGene <- probeList[ fit$p.value[,2] < 0.01]
} 


###################################################
### chunk number 32: GO analysis eval=FALSE
###################################################
## #line 660 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
## if (require(GOstats) & require(lumiHumanAll.db)) {
## 
## 	## Convert the probe Ids as  Entrez Ids and make them unique
## 	sigLL <- unique(unlist(lookUp(sigGene,'lumiHumanAll.db','ENTREZID')))
## 	sigLL <- as.character(sigLL[!is.na(sigLL)])
## 	params <- new("GOHyperGParams",
##               geneIds= sigLL,
##               annotation="lumiHumanAll.db",
##               ontology="BP",
##               pvalueCutoff= 0.01,
##               conditional=FALSE,
##               testDirection="over")
##          
##          hgOver <- hyperGTest(params)
##          	
## 	## Get the p-values of the test
## 	gGhyp.pv <- pvalues(hgOver)
## 	
## 	## Adjust p-values for multiple test (FDR)
## 	gGhyp.fdr <- p.adjust(gGhyp.pv, 'fdr')
## 
## 	## select the Go terms with adjusted p-value less than 0.01
## 	sigGO.ID <- names(gGhyp.fdr[gGhyp.fdr < 0.01])
## 	
## 	## Here only show the significant GO terms of BP (Molecular Function)
## 	## 	For other categories, just follow the same procedure.
## 	sigGO.Term <- getGOTerm(sigGO.ID)[["BP"]]
## }


###################################################
### chunk number 33: Create selected GO table eval=FALSE
###################################################
## #line 693 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
## if (require(GOstats) & require(lumiHumanAll.db) & require(xtable)) {
## 
## 	##get gene counts at each GO category
## 	gg.counts <- geneCounts(hgOver)[sigGO.ID]
## 	total.counts <- universeCounts(hgOver)[sigGO.ID]
## 
## 	ggt <- unlist(sigGO.Term)
## 	numCh <- nchar(ggt)
## 	ggt2 <- substr(ggt, 1, 17)
## 	ggt3 <- paste(ggt2, ifelse(numCh > 17, "...", ""), sep="")
## 	
## 	## output the significant GO categories as a table
## 	ggMat <- matrix(c(names(sigGO.Term), ggt3, signif(gGhyp.pv[sigGO.ID],5), gg.counts, total.counts),
##     		byrow=FALSE, nc=5, dimnames=list(1:length(sigGO.Term), c("GO ID",
##    		"Term", "p-value","Significant Genes No.", "Total Genes No.")))
## 	xtable.matrix(ggMat,
##   		caption="GO terms, p-values and counts.", label="ta:GOggterms")
## }


###################################################
### chunk number 34: sessionInfo
###################################################
#line 756 "E:/biocbld/bbs-2.7-bioc/tmpdir/RtmpFUupme/R.INSTALL15e72ac7/lumi/inst/doc/lumi.Rnw"
toLatex(sessionInfo())


