vignettes/lumi.R

In lumi: BeadArray Specific Methods for Illumina Methylation and Expression Microarrays

### R code from vignette source 'lumi.Rnw'

###################################################
### code chunk number 1: load library
###################################################
library(lumi)


###################################################
### code chunk number 2: load example data
###################################################
## load example data (a LumiBatch object)
data(example.lumi)
## summary of the example data
example.lumi


###################################################
### code chunk number 3: summary of example data
###################################################
## summary of the quality control
summary(example.lumi, 'QC')


###################################################
### code chunk number 4: densityPlot
###################################################
plot(example.lumi, what='density')	## plot the density


###################################################
### code chunk number 5: densityPlot
###################################################
plot(example.lumi, what='density')	## plot the density


###################################################
### code chunk number 6: CDF
###################################################
 plotCDF(example.lumi)


###################################################
### code chunk number 7: pairPlot
###################################################
plot(example.lumi, what='pair')		## pairwise plots


###################################################
### code chunk number 8: lumi.Rnw:233-234
###################################################
plot(example.lumi, what='pair', smoothScatter=T)


###################################################
### code chunk number 9: MAplot
###################################################
## pairwise MAplot
plot(example.lumi, what='MAplot')	


###################################################
### code chunk number 10: lumi.Rnw:265-267
###################################################
## pairwise MAplot
plot(example.lumi, what='MAplot', smoothScatter=T)	


###################################################
### code chunk number 11: lumi.Rnw:279-281
###################################################
## density plot of coefficient of varience
plot(example.lumi, what='cv')	


###################################################
### code chunk number 12: lumi.Rnw:292-293
###################################################
plot(example.lumi, what='sampleRelation')


###################################################
### code chunk number 13: lumi.Rnw:310-311
###################################################
plot(example.lumi, what='sampleRelation', method='mds', color=c("01", "02", "01", "02"))


###################################################
### code chunk number 14: VST transform
###################################################
## Do default VST variance stabilizing transform
lumi.T <- lumiT(example.lumi)


###################################################
### code chunk number 15: plotVST
###################################################
trans <- plotVST(lumi.T)


###################################################
### code chunk number 16: lumi.Rnw:356-358
###################################################
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)


###################################################
### code chunk number 17: Default normalization
###################################################
## Do quantile between microarray normaliazation
lumi.N <- lumiN(lumi.T)


###################################################
### code chunk number 18: Customized normalization (eval = FALSE)
###################################################
## ## Do RSN between microarray normaliazation
## lumi.N <- lumiN(lumi.T, method='rsn')


###################################################
### code chunk number 19: QC after normalization
###################################################
## 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


###################################################
### code chunk number 20: lumi.Rnw:405-406
###################################################
plot(lumi.N.Q, what='density')		## plot the density


###################################################
### code chunk number 21: lumi.Rnw:415-417
###################################################
plot(lumi.N.Q, what='boxplot')		## box plot
# boxplot(lumi.N.Q)


###################################################
### code chunk number 22: lumi.Rnw:426-427
###################################################
plot(lumi.N.Q, what='pair')			## pairwise plots


###################################################
### code chunk number 23: lumi.Rnw:436-437
###################################################
plot(lumi.N.Q, what='MAplot')		## plot the pairwise MAplot


###################################################
### code chunk number 24: lumi.Rnw:446-448
###################################################
## plot the sampleRelation using hierarchical clustering
plot(lumi.N.Q, what='sampleRelation')


###################################################
### code chunk number 25: lumi.Rnw:457-459
###################################################
## plot the sampleRelation using MDS
plot(lumi.N.Q, what='sampleRelation', method='mds', color=c("01", "02", "01", "02"))


###################################################
### code chunk number 26: Encapsulate the processing steps
###################################################
## Do all the default preprocessing in one step
lumi.N.Q <- lumiExpresso(example.lumi)


###################################################
### code chunk number 27: lumi.Rnw:476-478
###################################################
## Do all the preprocessing with customized settings
# lumi.N.Q <- lumiExpresso(example.lumi, normalize.param=list(method='rsn'))


###################################################
### code chunk number 28: inverse VST
###################################################
## 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)


###################################################
### code chunk number 29: retrieve normalized data
###################################################
dataMatrix <- exprs(lumi.N.Q)


###################################################
### code chunk number 30: filtering
###################################################
presentCount <- detectionCall(example.lumi)
selDataMatrix <- dataMatrix[presentCount > 0,]
if (require(lumiHumanAll.db) & require(annotate)) {
  selDataMatrix <- selDataMatrix[!is.na(getSYMBOL(rownames(selDataMatrix), 'lumiHumanAll.db')),]
}
probeList <- rownames(selDataMatrix)


###################################################
### code chunk number 31: Identify differentially expressed genes
###################################################
## 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]
} 


###################################################
### code chunk number 32: GO analysis (eval = FALSE)
###################################################
## 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"]]
## }


###################################################
### code chunk number 33: Create selected GO table (eval = FALSE)
###################################################
## 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")
## }


###################################################
### code chunk number 34: sessionInfo
###################################################
toLatex(sessionInfo())