################################################### ### chunk: ReadAffy1 eval=FALSE ################################################### ## library("affy") ## myAB = ReadAffy() ################################################### ### chunk: ReadAffy2 eval=FALSE ################################################### ## myAB = ReadAffy(filenames=c("a1.cel", "a2.cel", "a3.cel")) ################################################### ### chunk: CLL ################################################### library("CLL") data("CLLbatch") CLLbatch ################################################### ### chunk: sns ################################################### sampleNames(CLLbatch) ################################################### ### chunk: disease ################################################### data("disease") head(disease) ################################################### ### chunk: rownamesdisease ################################################### rownames(disease) = disease$SampleID ################################################### ### chunk: removeCELsuffix ################################################### sampleNames(CLLbatch) = sub("\\.CEL$", "", sampleNames(CLLbatch)) ################################################### ### chunk: mt ################################################### mt = match(rownames(disease), sampleNames(CLLbatch)) ################################################### ### chunk: ADF ################################################### vmd = data.frame(labelDescription = c("Sample ID", "Disease status: progressive or stable disease")) ################################################### ### chunk: samplesCLLbatch ################################################### phenoData(CLLbatch) = new("AnnotatedDataFrame", data = disease[mt, ], varMetadata = vmd) ################################################### ### chunk: dropmissing ################################################### CLLbatch = CLLbatch[, !is.na(CLLbatch$Disease)] ################################################### ### chunk: qccalc ################################################### library("affyQCReport") saqc = qc(CLLbatch) ################################################### ### chunk: qcplot ################################################### plot(saqc) ################################################### ### chunk: dist2 ################################################### dd = dist2(log2(exprs(CLLbatch))) ################################################### ### chunk: 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: affyPLM ################################################### library("affyPLM") dataPLM = fitPLM(CLLbatch) ################################################### ### chunk: NUSE ################################################### boxplot(dataPLM, main="NUSE", ylim = c(0.95, 1.22), outline = FALSE, col="lightblue", las=3, whisklty=0, staplelty=0) ################################################### ### chunk: RLE ################################################### Mbox(dataPLM, main="RLE", ylim = c(-0.4, 0.4), outline = FALSE, col="mistyrose", las=3, whisklty=0, staplelty=0) ################################################### ### chunk: dropBadarrays ################################################### badArray = match("CLL1", sampleNames(CLLbatch)) CLLB = CLLbatch[, -badArray] ################################################### ### chunk: repRLENUSE ################################################### 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: rma ################################################### CLLrma = rma(CLLB) ################################################### ### chunk: exprsDemo ################################################### e = exprs(CLLrma) dim(e) dim(CLLrma) ################################################### ### chunk: howmanyprobesets ################################################### dim(e)[1] nrow(e) dim(exprs(CLLrma))[1] nrow(CLLrma) length(featureNames(CLLrma)) ################################################### ### chunk: pData1 ################################################### pData(CLLrma)[1:3,] ################################################### ### chunk: pData2 ################################################### table(CLLrma$Disease) ################################################### ### chunk: nsFilter ################################################### CLLf = nsFilter(CLLrma, remove.dupEntrez=FALSE, var.cutof =0.5)$eset ################################################### ### chunk: rowM ################################################### CLLtt = rowttests(CLLf, "Disease") names(CLLtt) ################################################### ### chunk: rowMeans ################################################### a = rowMeans(exprs(CLLf)) ################################################### ### chunk: 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: eBayesEx ################################################### library("limma") design = model.matrix(~CLLf$Disease) CLLlim = lmFit(CLLf, design) CLLeb = eBayes(CLLlim) ################################################### ### chunk: compTstats ################################################### plot(CLLtt$statistic, CLLeb$t[,2], pch=".") ################################################### ### chunk: volcano ################################################### lod = -log10(CLLtt$p.value) plot(CLLtt$dm, lod, pch=".", xlab="log-ratio", ylab=expression(-log[10]~p)) abline(h=2) ################################################### ### chunk: volcanoeb ################################################### plot(CLLtt$dm, -log10(CLLeb$p.value[,2]), pch=".", xlab="log-ratio", ylab=expression(log[10]~p)) abline(h=2) ################################################### ### chunk: volcano2 ################################################### 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") ################################################### ### chunk: table ################################################### sum(CLLtt$p.value<=0.01) sum(CLLeb$p.value[,2]<=0.01) ################################################### ### chunk: toptable ################################################### tab = topTable(CLLeb, coef=2, adjust.method="BH", n=10) genenames = as.character(tab$ID) ################################################### ### chunk: annotate1 ################################################### library("annotate") ################################################### ### chunk: annotate2 ################################################### annotation(CLLf) library("hgu95av2.db") ################################################### ### chunk: entrezGeneAndSymbol ################################################### ll = getEG(genenames, "hgu95av2") sym = getSYMBOL(genenames, "hgu95av2") ################################################### ### chunk: 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: browse ################################################### browseURL("GeneList1.html") ################################################### ### chunk: ################################################### library("KEGG.db") library("annaffy") atab = aafTableAnn(genenames, "hgu95av2.db", aaf.handler()) saveHTML(atab, file="GeneList2.html") ################################################### ### chunk: annaffy2 ################################################### atab = aafTableAnn(genenames, "hgu95av2.db", aaf.handler()[c(2,5,8,12)]) saveHTML(atab, file="GeneList3.html") ################################################### ### chunk: pmmm ################################################### pms = pm(CLLB) mms = mm(CLLB) ################################################### ### chunk: 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: PMminusMM ################################################### table(sign(pms-mms)) ################################################### ### chunk: mmhist ################################################### 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: rmabgcorrect ################################################### bgrma = bg.correct.rma(CLLB) exprs(bgrma) = log2(exprs(bgrma)) ################################################### ### chunk: vsn ################################################### library("vsn") bgvsn = justvsn(CLLB) ################################################### ### chunk: sel ################################################### sel = sample(unlist(indexProbes(CLLB, "pm")), 500) sel = sel[order(exprs(CLLB)[sel, 1])] ################################################### ### chunk: yoyryv ################################################### yo = exprs(CLLB)[sel,1] yr = exprs(bgrma)[sel,1] yv = exprs(bgvsn)[sel,1] ################################################### ### chunk: 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: vsnrma ################################################### CLLvsn = vsnrma(CLLB) ################################################### ### chunk: vsndiff ################################################### CLLvsnf = nsFilter(CLLvsn, remove.dupEntrez=FALSE, var.cutoff=0.5)$eset CLLvsntt = rowttests(CLLvsnf, "Disease") ################################################### ### chunk: inboth ################################################### inboth = intersect(featureNames(CLLvsnf), featureNames(CLLf)) ################################################### ### chunk: comparermavsn ################################################### plot(CLLtt[inboth, "statistic"], CLLvsntt[inboth, "statistic"], pch=".", xlab="RMA", ylab="VSN", asp=1) ################################################### ### chunk: summary1 ################################################### pns = probeNames(CLLB) indices = split(seq(along=pns), pns) ################################################### ### chunk: summary2 ################################################### length(indices) indices[["189_s_at"]] ################################################### ### chunk: 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: summary2 ################################################### newsummary = t(sapply(indices, function(j) rowMedians(t(pms[j,]-mms[j,])))) dim(newsummary) ################################################### ### chunk: ans ################################################### colMeans(newsummary<0)*100