################################################### ### chunk number 1: options ################################################### options(width=48) ################################################### ### chunk number 2: ReadAffy1 eval=FALSE ################################################### ## library("affy") ## myAB1 <- ReadAffy() ## myAB2 <- ReadAffy(filenames=c("a1.cel", "a2.cel", "a3.cel")) ################################################### ### chunk number 3: CLL ################################################### library("CLL") data("CLLbatch") CLLbatch ################################################### ### chunk number 4: sns ################################################### head(sampleNames(CLLbatch)) sampleNames(CLLbatch) <- sub("\\.CEL$", "", sampleNames(CLLbatch)) head(sampleNames(CLLbatch)) ################################################### ### chunk number 5: disease ################################################### data("disease") head(disease) table(disease$Disease, useNA = "always") rownames(disease) <- disease$SampleID ################################################### ### chunk number 6: ADF ################################################### mt <- match(rownames(disease), sampleNames(CLLbatch)) vmd <- data.frame(labelDescription = c("Sample ID", "Disease status: progressive or stable disease")) phenoData(CLLbatch) <- new("AnnotatedDataFrame", data = disease[mt, ], varMetadata = vmd) ################################################### ### chunk number 7: dropmissing ################################################### CLLbatch <- CLLbatch[, !is.na(CLLbatch$Disease)] ################################################### ### chunk number 8: qccalc ################################################### library("simpleaffy") saqc <- qc(CLLbatch) ################################################### ### chunk number 9: qcplot ################################################### plot(saqc) ################################################### ### chunk number 10: dist2 ################################################### library("genefilter") dd <- dist2(log2(exprs(CLLbatch))) ################################################### ### chunk number 11: setupdistplot ################################################### diag(dd) <- 0 dd.row <- as.dendrogram(hclust(as.dist(dd))) row.ord <- order.dendrogram(dd.row) library("latticeExtra") legend <- list(top=list(fun=dendrogramGrob, args=list(x=dd.row, side="top"))) lp <- levelplot(dd[row.ord, row.ord], scales=list(x=list(rot=90)), xlab="", ylab="", legend=legend) ################################################### ### chunk number 12: distplot ################################################### plot(lp) ################################################### ### chunk number 13: affyPLM ################################################### library("affyPLM") dataPLM <- fitPLM(CLLbatch) ################################################### ### chunk number 14: NUSE ################################################### boxplot(dataPLM, main="NUSE", ylim = c(0.95, 1.22), outline = FALSE, col="lightblue", las=3, whisklty=0, staplelty=0) ################################################### ### chunk number 15: RLE ################################################### Mbox(dataPLM, main="RLE", ylim = c(-0.4, 0.4), outline = FALSE, col="mistyrose", las=3, whisklty=0, staplelty=0) ################################################### ### chunk number 16: dropBadarrays ################################################### badArray <- match("CLL1", sampleNames(CLLbatch)) CLLB <- CLLbatch[, -badArray] ################################################### ### chunk number 17: repRLENUSE ################################################### dataPLMx <- fitPLM(CLLB) ################################################### ### chunk number 18: repNUSE ################################################### boxplot(dataPLMx, main="NUSE", ylim = c(0.95, 1.3), outline = FALSE, col="lightblue", las=3, whisklty=0, staplelty=0) ################################################### ### chunk number 19: repRLE ################################################### Mbox(dataPLMx, main="RLE", ylim = c(-0.4, 0.4), outline = FALSE, col="mistyrose", las=3, whisklty=0, staplelty=0) ################################################### ### chunk number 20: rma ################################################### CLLrma <- rma(CLLB) ################################################### ### chunk number 21: exprsDemo ################################################### exprs(CLLrma)[1:3,1:3] pData(CLLrma)[1:3,] ################################################### ### chunk number 22: howmanyprobesets ################################################### e <- exprs(CLLrma) dim(e)[1] nrow(e) dim(exprs(CLLrma))[1] nrow(CLLrma) length(featureNames(CLLrma)) ################################################### ### chunk number 23: nsFilter ################################################### CLLf <- nsFilter(CLLrma, remove.dupEntrez = FALSE, var.cutof = 0.5)$eset ################################################### ### chunk number 24: rowMeans ################################################### a <- rowMeans(exprs(CLLf)) ################################################### ### chunk number 25: rowttests ################################################### CLLtt <- rowttests(CLLf, "Disease") names(CLLtt) ################################################### ### chunk number 26: figdvsa ################################################### par(mfrow = c(1,2)) myPlot <- function(...){ plot(y = CLLtt$dm, pch = ".", ylim = c(-2,2), ylab = "log-ratio", ...) abline(h = 0, col = "blue") } myPlot(x = a, xlab="average intensity") myPlot(x = rank(a), xlab="rank of average intensity") ################################################### ### chunk number 27: eBayesEx ################################################### library("limma") design <- model.matrix(~CLLf$Disease) CLLlim <- lmFit(CLLf, design) CLLeb <- eBayes(CLLlim) ################################################### ### chunk number 28: compTstats ################################################### plot(CLLtt$statistic, - CLLeb$t[,2], xlim = c(-6, 6), ylim = c(-6, 6), xlab = "Classical Linear Model", ylab = "Empirical Bayes", main = "Comparison of t Statistics", pch=".") abline(a = 0, b = 1) ################################################### ### chunk number 29: volcano ################################################### plot(CLLtt$dm, -log10(CLLtt$p.value), pch = ".", xlab = "log-ratio", ylab = expression(-log[10]~p)) abline(h = 2) ################################################### ### chunk number 30: volcanoeb ################################################### plot(CLLtt$dm, -log10(CLLeb$p.value[,2]), pch = ".", xlab = "log-ratio", ylab = expression(log[10]~p)) abline(h = 2) ################################################### ### chunk number 31: volcano2 ################################################### plot(CLLtt$dm, -log10(CLLtt$p.value), pch = ".", xlab = "log-ratio", ylab = expression(log[10]~p)) o1 <- order(abs(CLLtt$dm), decreasing = TRUE)[1:25] points(CLLtt$dm[o1], -log10(CLLtt$p.value)[o1], pch = 18, col = "blue") ################################################### ### chunk number 32: table ################################################### sum(CLLtt$p.value <= 0.01) sum(CLLeb$p.value[,2] <= 0.01) ################################################### ### chunk number 33: toptable ################################################### tab <- topTable(CLLeb, coef = 2, adjust.method = "BH", n = 10) genenames <- as.character(tab$ID) genenames ################################################### ### chunk number 34: annotate ################################################### library("annotate") annotation(CLLf) suppressMessages(library("hgu95av2.db")) ################################################### ### chunk number 35: entrezGeneAndSymbol ################################################### ll <- getEG(genenames, "hgu95av2") sym <- getSYMBOL(genenames, "hgu95av2") ################################################### ### chunk number 36: output ################################################### tab <- data.frame(sym, signif(tab[,-1], 3)) htmlpage(list(ll), othernames = tab, filename = "GeneList1.html", title = "HTML report", table.center = TRUE, table.head = c("Entrez ID", colnames(tab))) ################################################### ### chunk number 37: browse1 eval=FALSE ################################################### ## browseURL("GeneList1.html") ################################################### ### chunk number 38: annaffy ################################################### library("annaffy") library("KEGG.db") atab <- aafTableAnn(genenames, "hgu95av2.db", aaf.handler()) saveHTML(atab, file = "GeneList2.html") atab <- aafTableAnn(genenames, "hgu95av2.db", aaf.handler()[c(2,5,8,12)]) saveHTML(atab, file = "GeneList3.html") ################################################### ### chunk number 39: browse2 eval=FALSE ################################################### ## browseURL("GeneList2.html") ## browseURL("GeneList3.html") ################################################### ### chunk number 40: pmmm ################################################### pms <- pm(CLLB) mms <- mm(CLLB) ################################################### ### chunk number 41: pmmmplot ################################################### smoothScatter(log2(mms[,1]), log2(pms[,1]), xlab = expression(log[2] * "MM values"), ylab = expression(log[2] * "PM values"), asp = 1) abline(a = 0, b = 1, col = "red") ################################################### ### chunk number 42: PMminusMM ################################################### table(sign(pms - mms)) ################################################### ### chunk number 43: mmhist ################################################### library("geneplotter") grouping <- cut(log2(pms)[,2], breaks = c(-Inf, log2(2000), Inf), labels = c("Low", "High")) multidensity(log2(mms)[,2] ~ grouping, main = "", xlab = "", col = c("red", "blue"), lwd = 2) legend("topright", levels(grouping), lty = 1, lwd = 2, col = c("red", "blue")) ################################################### ### chunk number 44: rmabgcorrect ################################################### bgrma <- bg.correct.rma(CLLB) exprs(bgrma) <- log2(exprs(bgrma)) ################################################### ### chunk number 45: vsn ################################################### library("vsn") bgvsn <- justvsn(CLLB) ################################################### ### chunk number 46: sel ################################################### sel <- sample(unlist(indexProbes(CLLB, "pm")), 500) sel <- sel[order(exprs(CLLB)[sel, 1])] ################################################### ### chunk number 47: yoyryv ################################################### yo <- exprs(CLLB)[sel,1] yr <- exprs(bgrma)[sel,1] yv <- exprs(bgvsn)[sel,1] ################################################### ### chunk number 48: bgrmavsn ################################################### par(mfrow = c(1,3)) plot(yo, yr, xlab = "Original", ylab = "RMA", log = "x", type = "l", asp = 1) plot(yo, yv, xlab = "Original", ylab = "VSN", log = "x", type = "l", asp = 1) plot(yr, yv, xlab = "RMA", ylab = "VSN", type = "l", asp = 1) ################################################### ### chunk number 49: vsnrma ################################################### vsnrma CLLvsn <- vsnrma(CLLB) ################################################### ### chunk number 50: vsndiff ################################################### CLLvsnf <- nsFilter(CLLvsn, remove.dupEntrez = FALSE, var.cutoff = 0.5)$eset CLLvsntt <- rowttests(CLLvsnf, "Disease") ################################################### ### chunk number 51: inboth ################################################### inboth <- intersect(featureNames(CLLvsnf), featureNames(CLLf)) ################################################### ### chunk number 52: comparermavsn ################################################### plot(CLLtt[inboth, "statistic"], CLLvsntt[inboth, "statistic"], pch = ".", xlab = "RMA", ylab = "VSN", asp = 1) ################################################### ### chunk number 53: summary1 ################################################### pns <- probeNames(CLLB) indices <- split(seq(along = pns), pns) ################################################### ### chunk number 54: summary2 ################################################### length(indices) indices[["189_s_at"]] ################################################### ### chunk number 55: matplot ################################################### colors <- brewer.pal(8, "Dark2") Index <- indices[["189_s_at"]][seq(along=colors)] matplot(t(pms[Index, 1:12]), pch = "P", log = "y", type = "b", lty = 1, main = "189_s_at", xlab = "samples", ylab = expression(log[2]~Intensity), ylim = c(50, 2000), col = colors) matplot(t(mms[Index, 1:12]), pch = "M", log = "y", type = "b", lty = 3, add = TRUE, col = colors) ################################################### ### chunk number 56: summary3 ################################################### newsummary <- t(sapply(indices, function(j) rowMedians(t(pms[j,] - mms[j,])))) dim(newsummary) ################################################### ### chunk number 57: ans ################################################### colMeans(newsummary < 0) * 100