download code
## Chunk 1
library("Biobase")
library("genefilter")
library("ALL")
data("ALL")
## Chunk 2
bcell = grep("^B", as.character(ALL$BT))
moltyp = which(as.character(ALL$mol.biol) 
    %in% c("NEG", "BCR/ABL"))
ALL_bcrneg = ALL[, intersect(bcell, moltyp)]
ALL_bcrneg$mol.biol = factor(ALL_bcrneg$mol.biol)
## Chunk 3
library("genefilter")
sds = rowSds(exprs(ALL_bcrneg))
sh = shorth(sds)
sh
## Chunk 4
hist(sds, breaks=50, col="mistyrose", xlab="standard deviation")
abline(v=sh, col="blue", lwd=3, lty=2)
## Chunk 5
ALLsfilt = ALL_bcrneg[sds>=sh, ]
dim(exprs(ALLsfilt))
## Chunk 6
table(ALLsfilt$mol.biol)
tt = rowttests(ALLsfilt, "mol.biol")
names(tt)
## Chunk 7
hist(tt$p.value, breaks=50, col="mistyrose", xlab="p-value", 
     main="Retained")
## Chunk 8
ALLsrest = ALL_bcrneg[sds ttrest = rowttests(ALLsrest, "mol.biol")
hist(ttrest$p.value, breaks=50, col="lightblue", xlab="p-value", 
     main="Removed")
## Chunk 9
library("multtest")
mt = mt.rawp2adjp(tt$p.value, proc="BH")
## Chunk 10
g = featureNames(ALLsfilt)[mt$index[1:10]]
## Chunk 11
library("hgu95av2.db")
links(hgu95av2SYMBOL[g])
## Chunk 12
mb  = ALLsfilt$mol.biol
y   = exprs(ALLsfilt)[g[1],]
ord = order(mb)
plot(y[ord], pch=c(1,16)[mb[ord]], 
    col=c("black", "red")[mb[ord]], 
    main=g[1], ylab=expression(log[2]~intensity), 
    xlab="samples")