### R code from vignette source 'piano-vignette.Rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: piano-vignette.Rnw:55-56 ################################################### library(piano) ################################################### ### code chunk number 2: piano-vignette.Rnw:129-130 ################################################### library(piano) ################################################### ### code chunk number 3: piano-vignette.Rnw:179-192 ################################################### # An example setup: oxygen <- c("aerobic","anaerobic") limitation <- c("Clim","Nlim") mySetup <- cbind(oxygen[c(1,1,1,2,2,2,1,1,1,2,2,2)], limitation[c(1,1,1,1,1,1,2,2,2,2,2,2)]) # The rownames correspond to the CEL-file names (CAE1.CEL etc): rownames(mySetup) <- c("CAE1","CAE2","CAE3","CAN1","CAN2","CAN3", "NAE1","NAE2","NAE3","NAN1","NAN2","NAN3") colnames(mySetup) <- c("oxygen","limitation") # The final setup object can look like this: mySetup ################################################### ### code chunk number 4: piano-vignette.Rnw:219-226 ################################################### # Get path to example data and setup files: dataPath <- system.file("extdata", package="piano") # Load pre-normalized data: myArrayData <- loadMAdata(datadir=dataPath, dataNorm="norm_data.txt.gz", platform="yeast2") ################################################### ### code chunk number 5: piano-vignette.Rnw:231-232 ################################################### myArrayData ################################################### ### code chunk number 6: piano-vignette.Rnw:237-241 ################################################### # Check the setup: myArrayData$setup # Check the annotation (top 10 rows): myArrayData$annotation[1:10,] ################################################### ### code chunk number 7: piano-vignette.Rnw:250-251 ################################################### runQC(myArrayData) ################################################### ### code chunk number 8: piano-vignette.Rnw:256-262 ################################################### # To only run the PCA: runQC(myArrayData, rnaDeg=FALSE, nuseRle=FALSE, hist=FALSE, boxplot=FALSE, pca=TRUE) # Additionally, for the PCA you can specify other colors: runQC(myArrayData, rnaDeg=FALSE, nuseRle=FALSE, hist=FALSE, boxplot=FALSE, pca=TRUE, colors=c("cyan","orange")) ################################################### ### code chunk number 9: piano-vignette.Rnw:279-280 ################################################### extractFactors(myArrayData) ################################################### ### code chunk number 10: piano-vignette.Rnw:285-287 ################################################### pfc <- diffExp(myArrayData, contrasts=c("aerobic_Clim - anaerobic_Clim", "aerobic_Nlim - anaerobic_Nlim")) ################################################### ### code chunk number 11: piano-vignette.Rnw:297-303 ################################################### # Sort genes in "aerobic_Clim - anaerobic_Clim" according to adjusted p-value: ii <- sort(pfc$resTable[[1]]$adj.P.Val, index.return=TRUE)$ix pfc$resTable[[1]][ii[1:5],] # Sort genes in "aerobic_Nlim - anaerobic_Nlim" according to adjusted p-value: ii <- sort(pfc$resTable[[2]]$adj.P.Val, index.return=TRUE)$ix pfc$resTable[[2]][ii[1:5],] ################################################### ### code chunk number 12: piano-vignette.Rnw:332-336 ################################################### # Get p-values from the aerobic_Clim vs anaerobic_Clim comparison: myPval <- pfc$pValues["aerobic_Clim - anaerobic_Clim"] # Display the first values and gene IDs: head(myPval) ################################################### ### code chunk number 13: piano-vignette.Rnw:343-353 ################################################### # Custom gene to gene set mapping: genes2genesets <- cbind(paste("gene",c("A","A","A","B","B","C","C","C","D"),sep=""), paste("set",c(1,2,3,1,3,2,3,4,4),sep="")) genes2genesets # Load into correct format: myGsc <- loadGSC(genes2genesets) # View summary: myGsc # View all gene sets: myGsc$gsc ################################################### ### code chunk number 14: piano-vignette.Rnw:358-361 ################################################### myStats <- c(-1.5,-0.5,1,2) names(myStats) <- paste("gene",c("A","B","C","D"),sep="") myStats ################################################### ### code chunk number 15: piano-vignette.Rnw:395-396 ################################################### gsaRes <- runGSA(myStats, gsc=myGsc) ################################################### ### code chunk number 16: piano-vignette.Rnw:402-403 ################################################### gsaRes ################################################### ### code chunk number 17: piano-vignette.Rnw:406-407 ################################################### names(gsaRes) ################################################### ### code chunk number 18: piano-vignette.Rnw:414-415 ################################################### c(-4, -3, 2, 3.5) ################################################### ### code chunk number 19: piano-vignette.Rnw:418-419 ################################################### mean(c(-4, -3, 2.5, 4.5)) ################################################### ### code chunk number 20: piano-vignette.Rnw:422-423 ################################################### mean(abs(c(-4, -3, 2.5, 4.5))) ################################################### ### code chunk number 21: piano-vignette.Rnw:426-428 ################################################### mean(abs(c(-4, -3))) mean(abs(c(2.5, 4.5))) ################################################### ### code chunk number 22: piano-vignette.Rnw:433-436 ################################################### tmp<-c(0.0001,0.0001,0.9999,0.0001,0.005) names(tmp) <- c("p (non.dir)","p (dist.dir.up)","p (dist.dir.dn)","p (mix.dir.up)","p (mix.dir.dn)") tmp ################################################### ### code chunk number 23: piano-vignette.Rnw:439-442 ################################################### tmp<-c(100,95,5) names(tmp)<-c("Genes (tot)","Genes (up)","Genes (down)") tmp ################################################### ### code chunk number 24: piano-vignette.Rnw:473-479 ################################################### # Get the p-values for the contrast aerobic_Clim - anaerobic_Clim myPval <- pfc$pValues["aerobic_Clim - anaerobic_Clim"] head(myPval) # Get the fold-changes for the contrast aerobic_Clim - anaerobic_Clim myFC <- pfc$foldChanges["aerobic_Clim - anaerobic_Clim"] head(myFC) ################################################### ### code chunk number 25: piano-vignette.Rnw:482-494 (eval = FALSE) ################################################### ## library("biomaRt") ## # Select ensembl database and S. cerevisiae dataset: ## ensembl <- useMart("ENSEMBL_MART_ENSEMBL", dataset="scerevisiae_gene_ensembl" ,host="www.ensembl.org") ## # Map Yeast 2.0 microarray probeset IDs to GO: ## mapGO <- getBM(attributes=c('affy_yeast_2', 'name_1006'), ## #filters = 'affy_yeast_2', ## #values = rownames(myPval), ## mart = ensembl) ## # Remove blanks ("") ## mapGO <- mapGO[mapGO[,2]!="",] ## # Check the 10 first rows to see what we got: ## mapGO[1:10,] ################################################### ### code chunk number 26: piano-vignette.Rnw:497-498 (eval = FALSE) ################################################### ## myGsc <- loadGSC(mapGO) ################################################### ### code chunk number 27: piano-vignette.Rnw:535-536 (eval = FALSE) ################################################### ## GSAsummaryTable(gsaRes, save=TRUE, file="gsaResTab.xls") ################################################### ### code chunk number 28: piano-vignette.Rnw:569-570 (eval = FALSE) ################################################### ## nw <- networkPlot(gsaRes,class="non") ################################################### ### code chunk number 29: piano-vignette.Rnw:579-580 (eval = FALSE) ################################################### ## nw$geneSets ################################################### ### code chunk number 30: piano-vignette.Rnw:586-588 (eval = FALSE) ################################################### ## # An example usage: ## myGsc <- loadGSC("myModel.xml") ################################################### ### code chunk number 31: piano-vignette.Rnw:591-600 ################################################### # To load iTO977: metMap <- system.file("extdata", "probe2metabolites_iTO977.txt.gz", package="piano") # To load iIN800: metMap <- system.file("extdata", "probe2metabolites_iIN800.txt.gz", package="piano") # Convert into piano format: metMap <- read.delim(metMap) myGsc <- loadGSC(metMap) ################################################### ### code chunk number 32: piano-vignette.Rnw:603-607 ################################################### gsaRes <- runGSA(myPval,myFC,gsc=myGsc, geneSetStat="reporter", signifMethod="nullDist", nPerm=1000, gsSizeLim=c(5,100)) ################################################### ### code chunk number 33: piano-vignette.Rnw:609-610 ################################################### gsaRes ################################################### ### code chunk number 34: piano-vignette.Rnw:615-617 ################################################### nw <- networkPlot(gsaRes, class="distinct", direction="both", significance=0.005, label="numbers") ################################################### ### code chunk number 35: piano-vignette.Rnw:626-627 ################################################### nw$geneSets ################################################### ### code chunk number 36: piano-vignette.Rnw:631-635 ################################################### gsaResTab <- GSAsummaryTable(gsaRes) # Which columns contain p-values: grep("p \\(",colnames(gsaResTab),value=TRUE) grep("p \\(",colnames(gsaResTab)) ################################################### ### code chunk number 37: piano-vignette.Rnw:638-640 ################################################### ii <- which(gsaResTab[,10]<0.0001) gsaResTab$Name[ii] ################################################### ### code chunk number 38: piano-vignette.Rnw:643-650 ################################################### # Get minimum p-value for each gene set: minPval <- apply(gsaResTab[,c(4,7,10,14,18)],1,min,na.rm=TRUE) # Select significant gene sets: ii <- which(minPval<0.0001) gsaResTabSign <- gsaResTab[ii,c(1,4,7,10,14,18)] # Look at the first 10 gene sets: gsaResTabSign[1:10,] ################################################### ### code chunk number 39: piano-vignette.Rnw:666-668 (eval = FALSE) ################################################### ## myTval <- pfc$resTable[["aerobic_Clim - anaerobic_Clim"]]$t ## names(myTval) <- pfc$resTable[["aerobic_Clim - anaerobic_Clim"]]$ProbesetID ################################################### ### code chunk number 40: piano-vignette.Rnw:671-686 (eval = FALSE) ################################################### ## myGsc <- loadGSC(mapGO) ## gsaRes1 <- runGSA(myTval,geneSetStat="mean",gsc=myGsc, ## nPerm=1000,gsSizeLim=c(10,800)) ## gsaRes2 <- runGSA(myTval,geneSetStat="median",gsc=myGsc, ## nPerm=1000,gsSizeLim=c(10,800)) ## gsaRes3 <- runGSA(myTval,geneSetStat="sum",gsc=myGsc, ## nPerm=1000,gsSizeLim=c(10,800)) ## gsaRes4 <- runGSA(myTval,geneSetStat="maxmean",gsc=myGsc, ## nPerm=1000,gsSizeLim=c(10,800)) ## gsaRes5 <- runGSA(myPval,myFC,geneSetStat="fisher",gsc=myGsc, ## nPerm=1000,gsSizeLim=c(10,800)) ## gsaRes6 <- runGSA(myPval,myFC,geneSetStat="stouffer",gsc=myGsc, ## nPerm=1000,gsSizeLim=c(10,800)) ## gsaRes7 <- runGSA(myPval,myFC,geneSetStat="tailStrength",gsc=myGsc, ## nPerm=1000,gsSizeLim=c(10,800)) ################################################### ### code chunk number 41: piano-vignette.Rnw:689-693 (eval = FALSE) ################################################### ## resList <- list(gsaRes1,gsaRes2,gsaRes3,gsaRes4,gsaRes5,gsaRes6,gsaRes7) ## names(resList) <- c("mean","median","sum","maxmean","fisher", ## "stouffer","tailStrength") ## ch <- consensusHeatmap(resList,cutoff=30,method="mean") ################################################### ### code chunk number 42: piano-vignette.Rnw:696-697 (eval = FALSE) ################################################### ## ch$pMat