################################################### ### chunk number 1: loadlib ################################################### library("Biobase") library("ALL")

data(ALL, package="ALL")

ALLBs = ALL[,grep("^B", as.character(ALL$BT))]

ALLBCRNEG = ALLBs[, ALLBs$mol == "BCR/ABL" | ALLBs$mol =="NEG"]

ALLBCRNEG$mol.biol = factor(ALLBCRNEG$mol.biol)

numBN = length(ALLBCRNEG$mol.biol)

ALLBCRALL1 = ALLBs[, ALLBs$mol == "BCR/ABL" | ALLBs$mol == "ALL1/AF4"]

ALLBCRALL1$mol.biol = factor(ALLBCRALL1$mol.biol) numBA = length(ALLBCRALL1$mol.biol)

################################################### ### chunk number 2: IQRfilter ################################################### lowQ = rowQ(ALLBCRNEG, floor(0.25 numBN)) upQ = rowQ(ALLBCRNEG, ceiling(0.75 numBN)) iqrs = upQ - lowQ ##you can also do: ##iqrs = esApply(ALLBCRNEG, 1, IQR)

giqr = iqrs > quantile(iqrs, probs=.75) sum(giqr)

BNsub = ALLBCRNEG[giqr,]

################################################### ### chunk number 3: distEx ################################################### dSub <- BNsub[1:60,] eucD <- dist(t(exprs(dSub))) eucD@Size eucM <- as.matrix(eucD)

################################################### ### chunk number 4: distHM ################################################### library("RColorBrewer") hmcol <- colorRampPalette(brewer.pal(10, "RdBu"))(256) heatmap(eucM, sym=TRUE, col=hmcol, distfun=function(x) as.dist(x))

################################################### ### chunk number 5: helperfuns ################################################### closestN = function(distM, label) { loc = match(label, row.names(distM)) names(which.min(distM[label,-loc])) }

closestN(eucM, "03002")

################################################### ### chunk number 6: cM ################################################### library("bioDist") cD = MIdist(t(exprs(dSub))) cM = as.matrix(cD) closestN(cM, "03002")

################################################### ### chunk number 7: rowttests ################################################### library("genefilter") tt1 = rowttests(BNsub, "mol.biol") numToSel <- 50

################################################### ### chunk number 8: topNumToSel ################################################### tt1ord = order(abs(tt1$statistic), decreasing=TRUE)

top50 = tt1ord[1:numToSel]

BNsub1 = BNsub[top50,]

################################################### ### chunk number 9: rowIQRsFun ################################################### rowIQRs = function(eSet) { if (is(eSet, "exprSet")) numSamp = length(sampleNames(eSet)) else numSamp = ncol(eSet) lowQ = rowQ(eSet, floor(0.25 numSamp)) upQ = rowQ(eSet, ceiling(0.75 numSamp)) upQ - lowQ }

################################################### ### chunk number 10: standardize ################################################### standardize = function(x) (x - rowMedians(x)) / rowIQRs(x)

exprs(BNsub1) = standardize(exprs(BNsub1))

################################################### ### chunk number 11: heatmap ################################################### library("RColorBrewer") hmcol <- colorRampPalette(brewer.pal(10, "RdBu"))(256) spcol <- ifelse(BNsub1$mol.biol=="BCR/ABL", "goldenrod", "skyblue") heatmap(exprs(BNsub1), col=hmcol, ColSideColors=spcol)

################################################### ### chunk number 12: knn1 ###################################################

library("class")

a1 = knn.cv(t(exprs(BNsub1)), BNsub1$mol.biol)

ctab1 = table(a1, BNsub1$mol.biol)

errrate = (ctab1["BCR/ABL", "NEG"] + ctab1["NEG", "BCR/ABL"])/sum(ctab1)

################################################### ### chunk number 13: selectk ################################################### alist = NULL for(i in 1:4) alist[[i]] = knn.cv(t(exprs(BNsub1)), BNsub1$mol.biol, k=i)

sapply(alist, function(x) { ct1 = table(x, BNsub1$mol.biol) (ct1["BCR/ABL", "NEG"] + ct1["NEG", "BCR/ABL"]) / sum(ct1) })

################################################### ### chunk number 14: MLInterfaces ################################################### library("MLInterfaces")

################################################### ### chunk number 15: knn ################################################### knnResult <- MLearn("mol.biol", BNsub1, "knn", 1:50) knnResult

################################################### ### chunk number 16: knn-confusion ################################################### confuMat(knnResult)

################################################### ### chunk number 17: xval ################################################### knnXval <- xval(BNsub1, "mol.biol", knnB, xvalMethod="LOO")

################################################### ### chunk number 18: xval-table ################################################### table(given=BNsub1$mol.biol,predicted=knnXval)

################################################### ### chunk number 19: xval-feature-selection-function ################################################### t.fun<-function(data, fac) { (abs(rowttests(data,data[[fac]], tstatOnly=FALSE)$statistic)) }

################################################### ### chunk number 20: realXval ################################################### BNx = BNsub exprs(BNx) = standardize(exprs(BNx))

lk3f <- xval(BNx, "mol.biol", knnB, xvalMethod="LOO", 0:0, fsFun=t.fun, fsNum=50) table(given=BNx$mol.biol, predicted=lk3f$out)

################################################### ### chunk number 21: solXval ################################################### lk3f1 <- xval(BNx, "mol.biol", knnB, xvalMethod="LOO", 0:0, fsFun=t.fun, fsNum=100) table(given=BNx$mol.biol, predicted=lk3f1$out)

lk3f2 <- xval(BNx, "mol.biol", knnB, xvalMethod="LOO", 0:0, fsFun=t.fun, fsNum=5) table(given=BNx$mol.biol, predicted=lk3f2$out)

table(lk3f2$fs.memory)

##for the small example we can do a bit of investigating varsused = sort(unique(lk3f2$fs.memory)) heatmap(exprs(BNsub)[varsused,], ColSide = ifelse(BNx$mol=="NEG", "skyblue", "salmon"), col = topo.colors(255), keep.dendro=TRUE)

library("hgu95av2") unlist(mget(geneNames(BNsub)[varsused], hgu95av2SYMBOL))

################################################### ### chunk number 22: ################################################### sessionInfo()