################################################### ### chunk number 1: biocLiteEx eval=FALSE ################################################### ## source("http://bioconductor.org/biocLite.R") ## biocLite(c("graph", "xtable")) ################################################### ### chunk number 2: exSession ################################################### sessionInfo() ################################################### ### chunk number 3: ################################################### options(width = 40) ################################################### ### chunk number 4: loadAllLibsFirst ################################################### library("Biobase") library("hgu95av2") library("class") library("GO") library("RColorBrewer") library("BiocCaseStudies") library("CLL") library("matchprobes") library("hgu95av2probe") library("hgu95av2cdf") library("geneplotter") ################################################### ### chunk number 5: ex0 eval=FALSE ################################################### ## 2 * x + c(1,2) ################################################### ### chunk number 6: ex0sol1 ################################################### x <- c(0.1, 1.1, 2.5, 10) y <- 1:100 z <- y < 10 pets <- c(Rex="dog", Garfield="cat", Tweety="bird") ################################################### ### chunk number 7: ex0sol2 ################################################### 2 * x + c(1,2) ################################################### ### chunk number 8: ex0sol3 ################################################### ##logical y[z] ## integer y[1:4] y[-(1:95)] ## character pets["Garfield"] ################################################### ### chunk number 9: ex0sol4 ################################################### m <- matrix(1:12, ncol=4) m[1,3] ################################################### ### chunk number 10: ex0sol5 ################################################### l <- list(name="Paul", sex=factor("male"), age=35) l$name l[[3]] ################################################### ### chunk number 11: loadlib ################################################### library("Biobase") library("hgu95av2") ################################################### ### chunk number 12: showexSet ################################################### data(sample.ExpressionSet) exSet <- sample.ExpressionSet exSet ################################################### ### chunk number 13: exploringclass ################################################### class(exSet) slotNames(exSet) exSet$sex exSet[1,] exSet[,1] ################################################### ### chunk number 14: ex1Sol1 ################################################### is(exSet[1,], "ExpressionSet") ################################################### ### chunk number 15: ex1Sol2 ################################################### dim(pData(exSet[,1])) dim(exprs(exSet[,1])) ################################################### ### chunk number 16: ex1Sol3 ################################################### exSet[,exSet$sex=="Female"] ################################################### ### chunk number 17: loaddo ################################################### getwd() ls() dataPath = system.file("extdata", package="BiocCaseStudies") oldLoc = setwd(dataPath) dir() load("ALLmat.rda") ls() ################################################### ### chunk number 18: oldLoc ################################################### setwd(oldLoc) ################################################### ### chunk number 19: help1 eval=FALSE ################################################### ## help(load) ################################################### ### chunk number 20: classanddim ################################################### class(ALLmat) dim(ALLmat) ################################################### ### chunk number 21: setwd eval=FALSE ################################################### ## setwd(oldLoc) ################################################### ### chunk number 22: loadPhenoDataframe ################################################### samples <- read.table(file.path(dataPath, "ALL-sample-info.txt"), header=TRUE, check.names=FALSE) ################################################### ### chunk number 23: class ################################################### class(samples) ################################################### ### chunk number 24: colnames ################################################### names(samples) ################################################### ### chunk number 25: simpleSubsetting ################################################### samples[c(15, 30), c("sex", "age")] samples[samples$cod == "11005", c("sex", "age")] ################################################### ### chunk number 26: loadVarMeta ################################################### varInfo <- read.table(file.path(dataPath, "ALL-varMeta.txt"), header=TRUE, colClasses="character") varInfo[c("sex", "cod", "mol.biol"), , drop=FALSE] ################################################### ### chunk number 27: AnnotatedDataFrame ################################################### pd <- new("AnnotatedDataFrame", data=samples, varMetadata=varInfo) ################################################### ### chunk number 28: ExpressionSetFinally ################################################### ALLSet <- new("ExpressionSet", exprs=ALLmat, phenoData=pd, annotation="hgu95av2") ################################################### ### chunk number 29: usingDollar ################################################### ALLSet$sex[1:5] == "F" ALLSet$"t(9;22)"[1:5] ################################################### ### chunk number 30: featureNames ################################################### featureNames(ALLSet)[1:5] ################################################### ### chunk number 31: sampleNames ################################################### sampleNames(ALLSet)[1:5] varLabels(ALLSet) ################################################### ### chunk number 32: exprs ################################################### mat <- exprs(ALLSet) adf <- phenoData(ALLSet) ################################################### ### chunk number 33: simpleFun ################################################### sq1 <- function(x) return(x*x) sq2 <- function(x) x*x ################################################### ### chunk number 34: ppc ################################################### ppc <- function(x) paste("^", x, sep="") ################################################### ### chunk number 35: map ################################################### library("hgu95av2") hgu95av2MAP$"1001_at" ################################################### ### chunk number 36: eApplyex ################################################### myPos = eapply(hgu95av2MAP, function(x) grep("^17p", x, value=TRUE)) myPos = unlist(myPos) length(myPos) ################################################### ### chunk number 37: eApply2 ################################################### f17p = function(x) grep("^17p", x, value=TRUE) myPos2 = eapply(hgu95av2MAP, f17p) myPos2 = unlist(myPos2) length(myPos2) ################################################### ### chunk number 38: ans2 ################################################### myFindMap = function(mapEnv, which) { myg = ppc(which) a1 = eapply(mapEnv, function(x) grep(myg, x, value=TRUE)) unlist(a1) } ################################################### ### chunk number 39: envEx ################################################### e1 <- new.env(hash=TRUE) e1$a <- rnorm(10) e1$b <- runif(20) ls(e1) xx <- as.list(e1) names(xx) rm(a, envir=e1) ################################################### ### chunk number 40: ex2Sol1 ################################################### theEnv <- new.env(hash=TRUE) theEnv$exSet <- exSet ################################################### ### chunk number 41: ex2Sol2 ################################################### x <- 1:10 y <- 2 * x + rnorm(10, sd=0.25) model <- lm(y~x) theEnv$model <- model ################################################### ### chunk number 42: ex2Sol3 ################################################### myExtract <- function(env){ eset <- env$eset model <- env$model return(list(score=exSet$score, coeff=model$coeff)) } myExtract(theEnv) ################################################### ### chunk number 43: knn ################################################### library("class") apropos("knn") ################################################### ### chunk number 44: knnE1 ################################################### exprsExSet <- exprs(exSet) classExSet <- exSet$type esub <- exSet[, -1] pred1 <- knn(t(exprs(esub)), exprs(exSet)[,1], esub$type) classExSet[1] ################################################### ### chunk number 45: ex3Sol1 ################################################### predEset <- function(exSet, cov){ nc <- ncol(exSet) pred <- character(nc) for(i in 1:nc){ esub <- exSet[, -i] pred[i] <- as.character(knn(t(exprs(esub)), exprs(exSet)[, i], esub[[cov]])) } return(pred) } predEset(exSet, "type") ################################################### ### chunk number 46: ex3Sol2 ################################################### predEset <- function(exSet, cov, ...){ nc <- ncol(exSet) pred <- character(nc) pred <- character() for(i in 1:nc){ pred[i] <- as.character(knn(t(exprs(esub)), exprs(exSet)[, i], esub[[cov]], ...)) } return(pred) } ################################################### ### chunk number 47: GOexample ################################################### library("GO") library("hgu95av2") affyGO <- as.list(hgu95av2GO) #find the MF terms affyMF <- lapply(affyGO, function(x) { onts <- sapply(x, function(z) z$Ontology) if (is.null(unlist(onts)) || is.na(unlist(onts))) NA else unique(names(onts)[onts == "MF"]) }) ################################################### ### chunk number 48: ex4Sol1 ################################################### names(hgu95av2GO$"35889_at"[[1]]) ################################################### ### chunk number 49: ex4Sol2 ################################################### getGO <- function(ontology){ affyGO = as.list(hgu95av2GO) if(!ontology %in% c("MF", "BP", "CC")) stop("invalid ontology identifier") affyOnts = lapply(affyGO, function(x) { onts = sapply(x, function(z) z$Ontology) if( is.null(unlist(onts)) || is.na(unlist(onts)) ) NA else unique(names(onts)[onts==ontology]) }) return(affyOnts) } getGO("MF")[10:12] ################################################### ### chunk number 50: ex4Sol3 ################################################### getGO2 <- function(ontology, evidence){ affyGO = as.list(hgu95av2GO) if(!ontology %in% c("MF", "BP", "CC")) stop("invalid ontology identifier") affyOnts = lapply(affyGO, function(x) { onts = sapply(x, function(z) z$Ontology) evs = sapply(x, function(z) z$Evidence) if(is.null(unlist(onts)) || is.na(unlist(onts)) || is.null(unlist(evs)) || is.na(unlist(evs))) NA else unique(names(onts)[onts==ontology & evs==evidence]) }) return(affyOnts) } getGO2("MF", "IEA")[10:12] ################################################### ### chunk number 51: aprop ################################################### apropos("mean") find("mean") ################################################### ### chunk number 52: helpsearch eval=FALSE ################################################### ## help.search("mean") ################################################### ### chunk number 53: exHelp1print eval=FALSE ################################################### ## apropos("plot") ################################################### ### chunk number 54: exHelp1 ################################################### apropos("plot")[1:3] ################################################### ### chunk number 55: exHelp2 eval=FALSE ################################################### ## help.search("mann-whitney") ################################################### ### chunk number 56: dotplot ################################################### x <- exprs(exSet[,1]) y <- exprs(exSet[,3]) plot(x,y) ################################################### ### chunk number 57: ex5Sol1 ################################################### plot(x,y, xlab="gene expression sample #1", ylab="gene expression sample #3", main="scatterplot of expression intensities", pch=20) abline(a=0,b=1) ################################################### ### chunk number 58: pregcc ################################################### library("CLL") library("matchprobes") library("hgu95av2probe") library("hgu95av2cdf") library("RColorBrewer") data("CLLbatch") baseCt = basecontent(hgu95av2probe$sequence) gcCt = ordered(baseCt[, "C"]+baseCt[, "G"]) ################################################### ### chunk number 59: matchgcc ################################################### iab = get("xy2i", "package:hgu95av2cdf")(hgu95av2probe$x, hgu95av2probe$y) meanlogint = rowMeans(log2(exprs(CLLbatch)[iab, ])) ################################################### ### chunk number 60: boxplotgcc ################################################### mycol <- colorRampPalette(brewer.pal(9, "GnBu"))(nlevels(gcCt)) exlab <- expression(log[2]~intensity) boxplot(meanlogint ~ gcCt, col=mycol, outline=FALSE, xlab="Number of G and C", ylab=exlab) ################################################### ### chunk number 61: densplotgcc ################################################### library("geneplotter") tab = table(gcCt) ord = sort(order(tab, decreasing=TRUE)[1:10]) datsel = (as.character(gcCt) %in% names(tab)[ord]) gcCtsel = ordered(gcCt[datsel]) mycolsel = mycol[sort(ord)] ## why do I get strange errors for xlab and main arguments???? multidensity(meanlogint[datsel] ~ gcCtsel, xlim=c(6, 11), col=mycolsel, lwd=2) ################################################### ### chunk number 62: ex6Sol1 ################################################### multiecdf(meanlogint[datsel] ~ gcCtsel, xlim=c(6, 11), xlab=exlab, ylab="ECDF", main="", col=mycolsel, lwd=2) ################################################### ### chunk number 63: ################################################### sessionInfo()