Code examples for Chapter 7: 'Differential Expression'

## Chunk 1
`library("ALL") data("ALL")`

## Chunk 2
`bcell = grep("^B", as.character(ALL$BT))`

## Chunk 3
`moltyp = which(as.character(ALL$mol.biol) %in% c("NEG", "BCR/ABL"))`

## Chunk 4
`ALL_bcrneg = ALL[, intersect(bcell, moltyp)] ALL_bcrneg$mol.biol = factor(ALL_bcrneg$mol.biol)`

## Chunk 5
`library("vsn") meanSdPlot(ALL_bcrneg)`

## Chunk 6
`sds = esApply(ALL, 1, sd) sel = (sds > quantile(sds, 0.8)) ALLset1 = ALL_bcrneg[sel,]`

## Chunk 7
`library("genefilter") tt = rowttests(ALLset1, "mol.biol")`

## Chunk 8
`names(tt)`

## Chunk 9
`plot(tt$dm, -log10(tt$p.value), pch=".", xlab = expression(mean~log[2]~fold~change), ylab = expression(-log[10](p)))`

(solution chunk) ## Chunk 10
`sum(tt$p.value<0.05 & abs(tt$dm)>0.5)`

## Chunk 11
`library("multtest") cl = as.numeric(ALLset1$mol.biol=="BCR/ABL") resT = mt.maxT(exprs(ALLset1), classlabel=cl, B=1000) ord = order(resT$index) ## the original gene order rawp = resT$rawp[ord] ## permutation p-values`

## Chunk 12
`hist(rawp, breaks=50, col="#B2DF8A")`

## Chunk 13
`sum(resT$adjp<0.05)`

## Chunk 14
`res = mt.rawp2adjp(rawp, proc = "BH") sum(res$adjp[,"BH"]<0.05)`

## Chunk 15
`library("limma") design = cbind(mean = 1, diff = cl)`

## Chunk 16
`fit = lmFit(exprs(ALLset1), design) fit = eBayes(fit)`

## Chunk 17
`library("hgu95av2.db") ALLset1Syms = unlist(mget(featureNames(ALLset1), env = hgu95av2SYMBOL)) topTable(fit, coef = "diff", adjust.method = "fdr", sort.by = "p", genelist = ALLset1Syms)`

## Chunk 18
`plot(-log10(tt$p.value), -log10(fit$p.value[, "diff"]), xlab = "-log10(p) from two-sample t-test", ylab = "-log10(p) from moderated t-test (limma)", pch=".") abline(c(0, 1), col = "red")`

## Chunk 19
`subs = c(35, 65, 75, 1, 69, 71) ALLset2 = ALL_bcrneg[, subs] table(ALLset2$mol.biol)`

## Chunk 20
`tt2 = rowttests(ALLset2, "mol.biol") fit2 = eBayes(lmFit(exprs(ALLset2), design=design[subs, ]))`

## Chunk 21
`plot(-log10(tt2$p.value), -log10(fit2$p.value[, "diff"]), xlab = "-log10(p) from two-sample t-test", ylab = "-log10(p) from moderated t-test (limma)", pch=".") abline(c(0, 1), col = "red")`

## Chunk 22
`g = which(tt2$p.value < 1e-4 & fit2$p.value[, "diff"] > 0.02)`

## Chunk 23
`sel = (ALLset2$mol.bio == "BCR/ABL")+1 col = c("black", "red")[sel] pch = c(1,16)[sel] plot(exprs(ALLset2)[g,], pch=pch, col=col, ylab="expression")`

## Chunk 24
`rocs = rowpAUCs(ALLset1, "mol.biol", p=0.2)`

## Chunk 25
`j = which.max(area(rocs)) plot(rocs[j], main = featureNames(ALLset1)[j])`

(solution chunk) ## Chunk 26
`mtyp = ALLset1$mol.biol sel = rep(1:2, each=rev(table(mtyp))) plot(exprs(ALLset1)[j, order(mtyp)], pch=c(1,15)[sel], col=c("black", "red")[sel], main=featureNames(ALLset1)[j], ylab=expression(log[2]~expression~level)) legend("bottomleft", col=c("black", "red"), pch=c(1,15), levels(mtyp), bty="n")`

## Chunk 27
`nrsel.ttest = function(x, pthresh=0.05) { pval = rowttests(x, "mol.biol")$p.value return(sum(pval < pthresh)) }`

## Chunk 28
`nrsel.pAUC = function(x, pAUCthresh=2.5e-2) { pAUC = area(rowpAUCs(x, fac="mol.biol", p=0.1)) return(sum(pAUC > pAUCthresh)) }`

## Chunk 29
`library(BiocCaseStudies) x = ALLset1[sample(nrow(ALLset1), 1000), ]`

## Chunk 30
`par(mfrow = c(1, 2)) resample(x, "nrsel.ttest") resample(x, "nrsel.pAUC")`

## Chunk 31
`## resample(x, "nrsel.ttest") ## resample(x, "nrsel.pAUC")`