Code examples for Chapter 3: 'Processing Affymetrix Expression Data'

## Chunk 1
`## library("affy") ## myAB = ReadAffy()`

## Chunk 2
`## myAB = ReadAffy(filenames=c("a1.cel", "a2.cel", "a3.cel"))`

## Chunk 3
`library("CLL") data("CLLbatch") CLLbatch`

## Chunk 4
`sampleNames(CLLbatch)`

## Chunk 5
`data("disease") head(disease)`

## Chunk 6
`rownames(disease) = disease$SampleID`

## Chunk 7
`sampleNames(CLLbatch) = sub("\\.CEL$", "", sampleNames(CLLbatch))`

## Chunk 8
`mt = match(rownames(disease), sampleNames(CLLbatch))`

## Chunk 9
`vmd = data.frame(labelDescription = c("Sample ID", "Disease status: progressive or stable disease"))`

## Chunk 10
`phenoData(CLLbatch) = new("AnnotatedDataFrame", data = disease[mt, ], varMetadata = vmd)`

## Chunk 11
`CLLbatch = CLLbatch[, !is.na(CLLbatch$Disease)]`

## Chunk 12
`library("affyQCReport") saqc = qc(CLLbatch)`

## Chunk 13
`plot(saqc)`

## Chunk 14
`dd = dist2(log2(exprs(CLLbatch)))`

## Chunk 15
`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 16
`library("affyPLM") dataPLM = fitPLM(CLLbatch)`

## Chunk 17
`boxplot(dataPLM, main="NUSE", ylim = c(0.95, 1.22), outline = FALSE, col="lightblue", las=3, whisklty=0, staplelty=0)`

## Chunk 18
`Mbox(dataPLM, main="RLE", ylim = c(-0.4, 0.4), outline = FALSE, col="mistyrose", las=3, whisklty=0, staplelty=0)`

## Chunk 19
`badArray = match("CLL1", sampleNames(CLLbatch)) CLLB = CLLbatch[, -badArray]`

(solution chunk) ## Chunk 20
`dataPLMx = fitPLM(CLLB) boxplot(dataPLM, main="NUSE", ylim = c(0.95, 1.3), outline = FALSE, col="lightblue", las=3, whisklty=0, staplelty=0) Mbox(dataPLM, main="RLE", ylim = c(-0.4, 0.4), outline = FALSE, col="mistyrose", las=3, whisklty=0, staplelty=0)`

## Chunk 21
`CLLrma = rma(CLLB)`

## Chunk 22
`e = exprs(CLLrma) dim(e) dim(CLLrma)`

(solution chunk) ## Chunk 23
`dim(e)[1] nrow(e) dim(exprs(CLLrma))[1] nrow(CLLrma) length(featureNames(CLLrma))`

## Chunk 24
`pData(CLLrma)[1:3,]`

## Chunk 25
`table(CLLrma$Disease)`

## Chunk 26
`CLLf = nsFilter(CLLrma, remove.dupEntrez=FALSE, var.cutof =0.5)$eset`

## Chunk 27
`CLLtt = rowttests(CLLf, "Disease") names(CLLtt)`

## Chunk 28
`a = rowMeans(exprs(CLLf))`

(solution chunk) ## Chunk 29
`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 30
`library("limma") design = model.matrix(~CLLf$Disease) CLLlim = lmFit(CLLf, design) CLLeb = eBayes(CLLlim)`

(solution chunk) ## Chunk 31
`plot(CLLtt$statistic, CLLeb$t[,2], pch=".")`

## Chunk 32
`lod = -log10(CLLtt$p.value) plot(CLLtt$dm, lod, pch=".", xlab="log-ratio", ylab=expression(-log[10]~p)) abline(h=2)`

(solution chunk) ## Chunk 33
`plot(CLLtt$dm, -log10(CLLeb$p.value[,2]), pch=".", xlab="log-ratio", ylab=expression(log[10]~p)) abline(h=2)`

(solution chunk) ## Chunk 34
`plot(CLLtt$dm, lod, pch=".", xlab="log-ratio", ylab=expression(log[10]~p)) o1 = order(abs(CLLtt$dm), decreasing=TRUE)[1:25] points(CLLtt$dm[o1], lod[o1], pch=18, col="blue")`

(solution chunk) ## Chunk 35
`sum(CLLtt$p.value<=0.01) sum(CLLeb$p.value[,2]<=0.01)`

## Chunk 36
`tab = topTable(CLLeb, coef=2, adjust.method="BH", n=10) genenames = as.character(tab$ID)`

## Chunk 37
`library("annotate")`

## Chunk 38
`annotation(CLLf) library("hgu95av2.db")`

## Chunk 39
`ll = getEG(genenames, "hgu95av2") sym = getSYMBOL(genenames, "hgu95av2")`

## Chunk 40
`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 41
`browseURL("GeneList1.html")`

## Chunk 42
`library("KEGG.db") library("annaffy") atab = aafTableAnn(genenames, "hgu95av2.db", aaf.handler()) saveHTML(atab, file="GeneList2.html")`

## Chunk 43
`atab = aafTableAnn(genenames, "hgu95av2.db", aaf.handler()[c(2,5,8,12)]) saveHTML(atab, file="GeneList3.html")`

## Chunk 44
`pms = pm(CLLB) mms = mm(CLLB)`

(solution chunk) ## Chunk 45
`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")`

(solution chunk) ## Chunk 46
`table(sign(pms-mms))`

(solution chunk) ## Chunk 47
`grouping = cut(log2(pms)[,1], breaks=c(-Inf, log2(2000), Inf), labels=c("Low", "High")) library(geneplotter) multidensity(log2(mms)[,1] ~ grouping, main="", xlab="", col=c("red", "blue"), lwd=2) legend("topright", levels(grouping), lty=1, lwd=2, col=c("red", "blue"))`

## Chunk 48
`bgrma = bg.correct.rma(CLLB) exprs(bgrma) = log2(exprs(bgrma))`

## Chunk 49
`library("vsn") bgvsn = justvsn(CLLB)`

(solution chunk) ## Chunk 50
`sel = sample(unlist(indexProbes(CLLB, "pm")), 500) sel = sel[order(exprs(CLLB)[sel, 1])]`

(solution chunk) ## Chunk 51
`yo = exprs(CLLB)[sel,1] yr = exprs(bgrma)[sel,1] yv = exprs(bgvsn)[sel,1]`

(solution chunk) ## Chunk 52
`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 53
`CLLvsn = vsnrma(CLLB)`

## Chunk 54
`CLLvsnf = nsFilter(CLLvsn, remove.dupEntrez=FALSE, var.cutoff=0.5)$eset CLLvsntt = rowttests(CLLvsnf, "Disease")`

(solution chunk) ## Chunk 55
`inboth = intersect(featureNames(CLLvsnf), featureNames(CLLf))`

(solution chunk) ## Chunk 56
`plot(CLLtt[inboth, "statistic"], CLLvsntt[inboth, "statistic"], pch=".", xlab="RMA", ylab="VSN", asp=1)`

## Chunk 57
`pns = probeNames(CLLB) indices = split(seq(along=pns), pns)`

## Chunk 58
`length(indices) indices[["189_s_at"]]`

(solution chunk) ## Chunk 59
`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 60
`newsummary = t(sapply(indices, function(j) rowMedians(t(pms[j,]-mms[j,])))) dim(newsummary)`

(solution chunk) ## Chunk 61
`colMeans(newsummary<0)*100`