################################################### ### chunk number 1: load library ################################################### library(lumi) ################################################### ### chunk number 2: load example data ################################################### ## load example data (a LumiBatch object) data(example.lumi) ## summary of the example data example.lumi ################################################### ### chunk number 3: summary of example data ################################################### ## summary of the quality control summary(example.lumi, 'QC') ################################################### ### chunk number 4: densityPlot ################################################### plot(example.lumi, what='density') ## plot the density ################################################### ### chunk number 5: boxplot ################################################### plot(example.lumi, what='boxplot') ## box plot ################################################### ### chunk number 6: pairPlot ################################################### plot(example.lumi, what='pair') ## pairwise plots ################################################### ### chunk number 7: ################################################### plot(example.lumi, what='pair', smoothScatter=T) ################################################### ### chunk number 8: MAplot ################################################### ## pairwise MAplot plot(example.lumi, what='MAplot') ################################################### ### chunk number 9: ################################################### ## pairwise MAplot plot(example.lumi, what='MAplot', smoothScatter=T) ################################################### ### chunk number 10: ################################################### ## density plot of coefficient of varience plot(example.lumi, what='cv') ################################################### ### chunk number 11: ################################################### plot(example.lumi, what='sampleRelation') ################################################### ### chunk number 12: ################################################### plot(example.lumi, what='sampleRelation', method='mds', color=c("01", "02", "01", "02")) ################################################### ### chunk number 13: VST transform ################################################### ## Do default VST variance stabilizing transform lumi.T <- lumiT(example.lumi) ################################################### ### chunk number 14: plotVST ################################################### trans <- plotVST(lumi.T) ################################################### ### chunk number 15: ################################################### 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 16: Default normalization ################################################### ## Do quantile between microarray normaliazation lumi.N <- lumiN(lumi.T) ################################################### ### chunk number 17: Customized normalization ################################################### ## Do RSN between microarray normaliazation # lumi.N <- lumiN(lumi.T, method='rsn') ################################################### ### chunk number 18: 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 ################################################### ### chunk number 19: ################################################### plot(lumi.N.Q, what='density') ## plot the density ################################################### ### chunk number 20: ################################################### plot(lumi.N.Q, what='boxplot') ## box plot # boxplot(lumi.N.Q) ################################################### ### chunk number 21: ################################################### plot(lumi.N.Q, what='pair') ## pairwise plots ################################################### ### chunk number 22: ################################################### plot(lumi.N.Q, what='MAplot') ## plot the pairwise MAplot ################################################### ### chunk number 23: ################################################### ## plot the sampleRelation using hierarchical clustering plot(lumi.N.Q, what='sampleRelation') ################################################### ### chunk number 24: ################################################### ## plot the sampleRelation using MDS plot(lumi.N.Q, what='sampleRelation', method='mds', color=c("01", "02", "01", "02")) ################################################### ### chunk number 25: Encapsulate the processing steps ################################################### ## Do all the default preprocessing in one step lumi.N.Q <- lumiExpresso(example.lumi) ################################################### ### chunk number 26: ################################################### ## Do all the preprocessing with customized settings # lumi.N.Q <- lumiExpresso(example.lumi, normalize.param=list(method='rsn')) ################################################### ### chunk number 27: 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) ################################################### ### chunk number 28: retrieve normalized data ################################################### dataMatrix <- exprs(lumi.N.Q) ################################################### ### chunk number 29: filtering ################################################### presentCount <- detectionCall(example.lumi) selDataMatrix <- dataMatrix[presentCount > 0,] probeList <- rownames(selDataMatrix) ################################################### ### chunk number 30: 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] } ################################################### ### chunk number 31: 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"]] ## } ################################################### ### chunk number 32: 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") ## } ################################################### ### chunk number 33: sessionInfo ################################################### toLatex(sessionInfo())