### R code from vignette source 'a4vignette.Rnw' ################################################### ### code chunk number 1: config ################################################### options(width = 50) options(continue=" ") options(prompt="R> ") set.seed(123) ################################################### ### code chunk number 2: loadPackage ################################################### library(a4) require(ALL) data(ALL, package = "ALL") ################################################### ### code chunk number 3: prepareSimulatedData ################################################### require(nlcv) esSim <- simulateData(nEffectRows=50, betweenClassDifference = 5, nNoEffectCols = 5, withinClassSd = 0.2) ################################################### ### code chunk number 4: prepareALL ################################################### ALL <- addGeneInfo(ALL) ################################################### ### code chunk number 5: show ALL (eval = FALSE) ################################################### ## # The phenotypic data ## head(pData(ALL)[, c(1:5, 13, 18, 21)]) ## # The gene expression data ## head(exprs(ALL)[, 1:5]) ## # The feature data ## fDat <- head(pData(featureData(ALL))) ## fDat[,"Description"] <- substr(fDat[,"Description"], 1, 30) ## fDat ################################################### ### code chunk number 6: Create bcrAblOrNeg ################################################### Bcell <- grep("^B", as.character(ALL$BT)) # create B-Cell subset for ALL subsetType <- "BCR/ABL" # other subsetType can be "ALL/AF4" bcrAblOrNegIdx <- which(as.character(ALL$mol) %in% c("NEG", subsetType)) bcrAblOrNeg <- ALL[, intersect(Bcell, bcrAblOrNegIdx)] bcrAblOrNeg$mol.biol <- factor(bcrAblOrNeg$mol.biol) ################################################### ### code chunk number 7: spectralMapALL ################################################### spectralMap(object = ALL, groups = "BT") # optional argument settings # plot.mpm.args=list(label.tol = 12, zoom = c(1,2), do.smoothScatter = TRUE), # probe2gene = TRUE) ################################################### ### code chunk number 8: spectralMapALLSubset ################################################### spectralMap(object = bcrAblOrNeg, groups = "mol.biol", probe2gene = TRUE) ################################################### ### code chunk number 9: filtering ################################################### selBcrAblOrNeg <- filterVarInt(object = bcrAblOrNeg) propSelGenes <- round((dim(selBcrAblOrNeg)[1]/dim(bcrAblOrNeg)[1])*100,1) ################################################### ### code chunk number 10: tTest ################################################### tTestResult <- tTest(selBcrAblOrNeg, "mol.biol") ################################################### ### code chunk number 11: tTestHist ################################################### histPvalue(tTestResult[,"p"], addLegend = TRUE) propDEgenesRes <- propDEgenes(tTestResult[,"p"]) ################################################### ### code chunk number 12: tabTTest ################################################### tabTTest <- topTable(tTestResult, n = 10) xtable(tabTTest, caption="The top 5 features selected by an ordinary t-test.", label ="tablassoClass") ################################################### ### code chunk number 13: tTestVolcanoPlot ################################################### volcanoPlot(tTestResult, topPValues = 5, topLogRatios = 5) ################################################### ### code chunk number 14: limmaTwoLevels ################################################### limmaResult <- limmaTwoLevels(selBcrAblOrNeg, "mol.biol") volcanoPlot(limmaResult) # histPvalue(limmaResult) # propDEgenes(limmaResult) ################################################### ### code chunk number 15: limma ################################################### tabLimma <- topTable(limmaResult, n = 10, coef = 2) # 1st is (Intercept) ################################################### ### code chunk number 16: annotationTableLimma ################################################### tabLimmaSel <- tabLimma[,c("SYMBOL", "logFC", "AveExpr","P.Value","adj.P.Val", "GENENAME" )] tabLimmaSel[, "GENENAME"] <- substr(tabLimmaSel[,"GENENAME"], 1, 38) dData <- data.frame(Gene = tabLimmaSel[, 1], tabLimmaSel[,-1], stringsAsFactors = FALSE, row.names = NULL) hData <- data.frame(Gene = generateEntrezIdLinks(tabLimma[,"ENTREZID"])) tabAnnot <- annotationTable(displayData = dData, hrefData = hData) xTabAnnot <- a4Reporting::xtable(tabAnnot, digits = 2, caption = "Top differentially expressed genes between disabled anf functional p53 cell lines.") print(xTabAnnot, include.rownames = FALSE, floating = FALSE) ################################################### ### code chunk number 17: limmaReg ################################################### ################################################### ### code chunk number 18: PAM ################################################### resultPam <- pamClass(selBcrAblOrNeg, "mol.biol") plot(resultPam) featResultPam <- topTable(resultPam, n = 15) xtable(head(featResultPam$listGenes), caption = "Top 5 features selected by PAM.") ################################################### ### code chunk number 19: randomForest ################################################### resultRF <- rfClass(selBcrAblOrNeg, "mol.biol") plot(resultRF, which = 2) featResultRF <- topTable(resultRF, n = 15) xtable(head(featResultRF$topList), caption = "Features selected by Random Forest variable importance.") ################################################### ### code chunk number 20: plotTop2_3genesRf ################################################### plotCombination2genes(probesetId1=rownames(featResultRF$topList)[1], probesetId2=rownames(featResultRF$topList)[2], object = selBcrAblOrNeg, groups = "mol.biol") ################################################### ### code chunk number 21: loadNlcvTT ################################################### # nlcvTT <- nlcv(selBcrAblOrNeg, classVar = 'mol.biol', # classdist = "unbalanced", # nRuns = 10, fsMethod = "t.test", verbose = TRUE) data(nlcvTT) ################################################### ### code chunk number 22: MCRPlot ################################################### mcrPlot_TT <- mcrPlot(nlcvTT, plot = TRUE, optimalDots = TRUE, layout = TRUE, main = "t-test selection") ################################################### ### code chunk number 23: tabmcrPlot ################################################### xtable(summary(mcrPlot_TT),rownames=TRUE, caption = "Optimal number of genes per classification method together with the respective misclassification error rate (mean and standard deviation across all CV loops).", label = "tabmcrPlot_TT") ################################################### ### code chunk number 24: ScoresPlot ################################################### scoresPlot(nlcvTT, tech = "svm", nfeat = 2) ################################################### ### code chunk number 25: lasso ################################################### resultLasso <- lassoClass(object = bcrAblOrNeg, groups = "mol.biol") plot(resultLasso, label = TRUE, main = "Lasso coefficients in relation to degree of penalization.") featResultLasso <- topTable(resultLasso, n = 15) ################################################### ### code chunk number 26: tabLasso ################################################### lassoTable <- xtable(featResultLasso, label = "tablassoClass", caption = "Features selected by Lasso, ranked from largest to smallest penalized coefficient.") print(lassoTable, include.rownames = FALSE) ################################################### ### code chunk number 27: plotTop2_3genesLasso ################################################### op <- par(mfrow=c(1,2)) plotCombination2genes(geneSymbol1 = featResultLasso$topList[1, 1], geneSymbol2 = featResultLasso$topList[2, 1], object = bcrAblOrNeg, groups = "mol.biol", main = "Combination of\nfirst and second gene", addLegend = TRUE, legendPos = "topright") plotCombination2genes(geneSymbol1 = featResultLasso$topList[1, 1], geneSymbol2 = featResultLasso$topList[3, 1], object = bcrAblOrNeg, groups = "mol.biol", main = "Combination of\nfirst and third gene", addLegend = FALSE) par(op) ################################################### ### code chunk number 28: LogisticRegression ################################################### logRegRes <- logReg(geneSymbol = "ABL1", object = bcrAblOrNeg, groups = "mol.biol") ################################################### ### code chunk number 29: LogisticRegressionPlot ################################################### probabilitiesPlot(proportions = logRegRes$fit, classVar = logRegRes$y, sampleNames = rownames(logRegRes), main = "Probability of being NEG") ################################################### ### code chunk number 30: LogisticRegressionPlotBars ################################################### probabilitiesPlot(proportions = logRegRes$fit, classVar = logRegRes$y, barPlot= TRUE, sampleNames = rownames(logRegRes), main = "Probability of being NEG") ################################################### ### code chunk number 31: ROC ################################################### ROCres <- ROCcurve(geneSymbol = "ABL1", object = bcrAblOrNeg, groups = "mol.biol") ################################################### ### code chunk number 32: plotProfile ################################################### plot1gene(probesetId = rownames(tTestResult)[1], object = selBcrAblOrNeg, groups = "mol.biol", legendPos = "topright") ################################################### ### code chunk number 33: otherSampleIDsInPlot1gene ################################################### plot1gene(probesetId = rownames(tTestResult)[1], object = selBcrAblOrNeg, groups = "mol.biol", sampleIDs = "mol.biol", legendPos = "topright") ################################################### ### code chunk number 34: plot1gene2vars ################################################### plot1gene(probesetId = rownames(tTestResult)[1], object = selBcrAblOrNeg, groups = "mol.biol", colgroups = 'BT', legendPos = "topright") ################################################### ### code chunk number 35: boxPlot ################################################### boxPlot(probesetId = rownames(tTestResult)[1], object = selBcrAblOrNeg, boxwex = 0.3, groups = "mol.biol", colgroups = "BT", legendPos = "topright") ################################################### ### code chunk number 36: plotTop2_3genesLasso ################################################### plotCombination2genes(geneSymbol1 = featResultLasso$topList[1, 1], geneSymbol2 = featResultLasso$topList[2, 1], object = bcrAblOrNeg, groups = "mol.biol", main = "Combination of\nfirst and second gene", addLegend = TRUE, legendPos = "topright") ################################################### ### code chunk number 37: profilesPlot ################################################### myGeneSymbol <- "LCK" probesetPos <- which(myGeneSymbol == featureData(ALL)$SYMBOL) myProbesetIds <- featureNames(ALL)[probesetPos] profilesPlot(object = ALL, probesetIds = myProbesetIds, orderGroups = "BT", sampleIDs = "BT") ################################################### ### code chunk number 38: plotComb2Samples ################################################### plotComb2Samples(ALL, "11002", "01003", xlab = "a T-cell", ylab = "another T-cell") ################################################### ### code chunk number 39: plotComb2Samples ################################################### require(Cairo) CairoPNG(filename="plotComb2Samples.png",width=500,height=500) plotComb2Samples(ALL, "11002", "01003", xlab = "a T-cell", ylab = "another T-cell") dev.off() ################################################### ### code chunk number 40: plotComb2SamplesWithAnnotation ################################################### plotComb2Samples(ALL, "84004", "01003", trsholdX = c(10,12), trsholdY = c(4,6), xlab = "a B-cell", ylab = "a T-cell") ################################################### ### code chunk number 41: plotComb2SamplesWithAnnotation2 ################################################### CairoPNG(filename="plotComb2SamplesWithAnnotation.png",width=500,height=500) plotComb2Samples(ALL,"84004", "01003", trsholdX = c(10,12), trsholdY = c(4,6), xlab = "a B-cell", ylab = "a T-cell") dev.off() ################################################### ### code chunk number 42: plotCombMultipleSamples ################################################### plotCombMultSamples(exprs(ALL)[,c("84004", "11002", "01003")]) # text.panel= function(x){x, labels = c("a B-cell", "a T-cell", "another T-cell")}) ################################################### ### code chunk number 43: plotCombMultipleSamples ################################################### CairoPNG(filename="plotCombMultipleSamples.png", width=500, height=500) plotCombMultSamples(exprs(ALL)[, c("84004", "11002", "01003")]) dev.off() ################################################### ### code chunk number 44: GeneLRlist ################################################### ALL$BTtype <- as.factor(substr(ALL$BT,0,1)) ALL2 <- ALL[,ALL$BT != 'T1'] # omit subtype T1 as it only contains one sample ALL2$BTtype <- as.factor(substr(ALL2$BT,0,1)) # create a vector with only T and B # Test for differential expression between B and T cells tTestResult <- tTest(ALL, "BTtype", probe2gene = FALSE) topGenes <- rownames(tTestResult)[1:20] # plot the log ratios versus subtype B of the top genes LogRatioALL <- computeLogRatio(ALL2, reference = list(var="BT", level="B")) a <- plotLogRatio(e = LogRatioALL[topGenes,], openFile = FALSE, tooltipvalues = FALSE, device = "pdf", filename = "GeneLRlist", colorsColumnsBy = "BTtype", main = 'Top 20 genes most differentially between T- and B-cells', orderBy = list(rows = "hclust"), probe2gene = TRUE) ################################################### ### code chunk number 45: plotLogRatioComplex ################################################### load(system.file("extdata", "esetExampleTimeCourse.rda", package = "a4")) logRatioEset <- computeLogRatio(esetExampleTimeCourse, within = "hours", reference = list(var = "compound", level = "DMSO")) # re-order idx <- order(pData(logRatioEset)$compound, pData(logRatioEset)$hours) logRatioEset <- logRatioEset[,idx] # plot LogRatioEset across all cl <- "TEST" compound <- "COMPOUND" shortvarnames <- unique(interaction(pData(logRatioEset)$compound, pData(logRatioEset)$hours)) shortvarnames <- shortvarnames[-grep("DMSO", shortvarnames), drop=TRUE] plotLogRatio(e = logRatioEset, mx = 1, filename = "logRatioOverallTimeCourse.pdf", gene.fontsize = 8, orderBy = list(rows = "hclust", cols = NULL), colorsColumnsBy = c('compound'), within = "hours", shortvarnames = shortvarnames, exp.width = 1, main = paste("Differential Expression (trend at early time points) in", cl, "upon treatment with", compound), reference = list(var = "compound", level = "DMSO"), device = 'pdf') ################################################### ### code chunk number 46: MLP ################################################### require(MLP) # create groups labels <- as.factor(ifelse(regexpr("^B", as.character(pData(ALL)$BT))==1, "B", "T")) pData(ALL)$BT2 <- labels # generate p-values limmaResult <- limmaTwoLevels(object = ALL, group = "BT2") pValues <- limmaResult@MArrayLM$p.value pValueNames <- fData(ALL)[rownames(pValues), 'ENTREZID'] pValues <- pValues[,2] names(pValues) <- pValueNames pValues <- pValues[!is.na(pValueNames)] ################################################### ### code chunk number 47: geneSet ################################################### geneSet <- getGeneSets(species = "Human", geneSetSource = "GOBP", entrezIdentifiers = names(pValues) ) tail(geneSet, 3) ################################################### ### code chunk number 48: MLP ################################################### mlpOut <- MLP( geneSet = geneSet, geneStatistic = pValues, minGenes = 5, maxGenes = 100, rowPermutations = TRUE, nPermutations = 50, smoothPValues = TRUE, probabilityVector = c(0.5, 0.9, 0.95, 0.99, 0.999, 0.9999, 0.99999), df = 9) ################################################### ### code chunk number 49: GOgraph ################################################### pdf(file = "GOgraph.pdf") plot(mlpOut, type = "GOgraph", nRow = 25) dev.off() ################################################### ### code chunk number 50: sessionInfo ################################################### toLatex(sessionInfo())