################################################### ### chunk number 1: loadAllLibsFirst ################################################### library("Biobase") library("ALL") library("RColorBrewer") library("bioDist") library("genefilter") library("class") library("MLInterfaces") library("hgu95av2") ################################################### ### chunk number 2: 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 3: 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 4: distEx ################################################### dSub <- BNsub[1:60,] eucD <- dist(t(exprs(dSub))) eucD@Size eucM <- as.matrix(eucD) ################################################### ### chunk number 5: distHM ################################################### library("RColorBrewer") hmcol <- colorRampPalette(brewer.pal(10, "RdBu"))(256) heatmap(eucM, sym=TRUE, col=hmcol, distfun=function(x) as.dist(x)) ################################################### ### chunk number 6: distTau ################################################### tauD = tau.dist(t(exprs(dSub))) tauD@Size tauM <- as.matrix(tauD) ################################################### ### chunk number 7: tauFig ################################################### heatmap(tauM, sym=TRUE, col=hmcol, distfun=function(x) as.dist(x)) ################################################### ### chunk number 8: helperfuns ################################################### closestN = function(distM, label) { loc = match(label, row.names(distM)) names(which.min(distM[label,-loc])) } closestN(eucM, "03002") ################################################### ### chunk number 9: cM ################################################### library("bioDist") cD = MIdist(t(exprs(dSub))) cM = as.matrix(cD) closestN(cM, "03002") ################################################### ### chunk number 10: rowttests ################################################### library("genefilter") tt1 = rowttests(BNsub, "mol.biol") numToSel <- 50 ################################################### ### chunk number 11: topNumToSel ################################################### tt1ord = order(abs(tt1$statistic), decreasing=TRUE) top50 = tt1ord[1:numToSel] BNsub1 = BNsub[top50,] ################################################### ### chunk number 12: largest t-stat ################################################### largeT = which.max(tt1$statistic) tt1$statistic[largeT] tt1$p.value[largeT] featureNames(BNsub)[largeT] hgu95av2SYMBOL[[featureNames(BNsub)[largeT]]] ################################################### ### chunk number 13: rowIQRsFun ################################################### rowIQRs = function(eSet) { numSamp = ncol(eSet) lowQ = rowQ(eSet, floor(0.25 * numSamp)) upQ = rowQ(eSet, ceiling(0.75 * numSamp)) upQ - lowQ } ################################################### ### chunk number 14: useFun ################################################### byFun = rowIQRs(ALLBCRNEG) all(byFun == iqrs) ################################################### ### chunk number 15: standardize ################################################### standardize = function(x) (x - rowMedians(x)) / rowIQRs(x) exprs(BNsub1) = standardize(exprs(BNsub1)) ################################################### ### chunk number 16: 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 17: Samples ################################################### Negs = which(BNsub1$mol.biol == "NEG") Bcr = which(BNsub1$mol.biol == "BCR/ABL") S1 = sample(Negs, 20, replace=FALSE) S2 = sample(Bcr, 20, replace = FALSE) TrainInd = c(S1, S2) ################################################### ### chunk number 18: MLearn ################################################### kans = MLearn( mol.biol ~ ., BNsub1, "knn", TrainInd) confuMat(kans) dldans = MLearn( mol.biol ~ ., BNsub1, "dlda", TrainInd) confuMat(dldans) ldaans = MLearn( mol.biol ~ ., BNsub1, "dlda", TrainInd) confuMat(ldaans) ################################################### ### chunk number 19: xvalM ################################################### lk1 <- xvalML(mol.biol ~ ., BNsub1, "knn", xvalMethod="LOO") table(lk1, BNsub1$mol.biol) ################################################### ### chunk number 20: 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 21: findK ################################################### lk2 <- xvalML(mol.biol~., data=BNsub1, "knn", xvalMethod="LOO", k = 2) table(lk2, BNsub1$mol.biol) lk3 <- xvalML(mol.biol~., data=BNsub1, "knn", xvalMethod="LOO", k = 3) table(lk3, BNsub1$mol.biol) lk5 <- xvalML(mol.biol~., data=BNsub1, "knn", xvalMethod="LOO", k = 5) table(lk5, BNsub1$mol.biol) ################################################### ### chunk number 22: selectk ################################################### alist = list() 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 23: MLInterfaces ################################################### library("MLInterfaces") ################################################### ### chunk number 24: knn ################################################### knnResult <- MLearn(mol.biol ~ ., BNsub1, "knn", 1:50) knnResult ################################################### ### chunk number 25: knn-confusion ################################################### confuMat(knnResult) ################################################### ### chunk number 26: xval ################################################### knnXval <- xvalML(mol.biol~., data=BNsub1, "knn", xvalMethod="LOO") ################################################### ### chunk number 27: xval-table ################################################### table(given=BNsub1$mol.biol,predicted=knnXval) ################################################### ### chunk number 28: realXval ################################################### BNx = BNsub exprs(BNx) = standardize(exprs(BNx)) t.fun<-function(data, fac) { (abs(rowttests(data,data[[fac]], tstatOnly=FALSE)$statistic)) } lk3f <- xvalML(mol.biol~., data=BNx, "knn", xvalMethod="LOO", fsFun=t.fun, fsNum=50) table(given=BNx$mol.biol, predicted=lk3f$out) ################################################### ### chunk number 29: solXval ################################################### lk3f1 <- xvalML(mol.biol~., data=BNx, "knn", xvalMethod="LOO", fsFun=t.fun, fsNum=100) table(given=BNx$mol.biol, predicted=lk3f1$out) lk3f2 <- xvalML(mol.biol~., data=BNx, "knn", xvalMethod="LOO", 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(featureNames(BNsub)[varsused], hgu95av2SYMBOL)) ################################################### ### chunk number 30: loadlibs ################################################### library(randomForest) ################################################### ### chunk number 31: rf1 ################################################### set.seed(123) trainY = BNsub$mol.biol[TrainInd] Xm = t(exprs(BNsub)[,TrainInd]) rf1 <- randomForest(Xm, trainY, ntree=2000, mtry=55, importance=TRUE) rf1 rf2 <- randomForest(Xm, trainY, ntree=2000, mtry=35, importance=TRUE) rf2 vcrf1 = MLearn( mol.biol~., data=BNsub, "randomForest", TrainInd, ntree=2000, mtry=55, importance=TRUE) vcrf1 vcrf2 = MLearn( mol.biol~., data=BNsub, "randomForest", TrainInd, ntree=2000, mtry=35, importance=TRUE) vcrf2 ################################################### ### chunk number 32: rfpred ################################################### p1 <- predict(rf1, Xm, prox=TRUE) table(trainY, p1$pred) p2 <- predict(rf2, Xm, prox=TRUE) table(trainY, p2$pred) confuMat(vcrf1) confuMat(vcrf2) ################################################### ### chunk number 33: ################################################### varImpPlot(rf1, n.var=15) ################################################### ### chunk number 34: ################################################### varImpPlot(rf2, n.var=15) ################################################### ### chunk number 35: varimp ################################################### impvars <- function(x, which=2, k=10) { v1<-order(x$importance[,which]) l1 <- length(v1) x$importance[v1[(l1-k+1):l1], which] } iv.rf1 <- impvars(rf1, k=25) library("hgu95av2") library(annotate) isyms <- getSYMBOL(names(iv.rf1),data="hgu95av2") ################################################### ### chunk number 36: doipl ################################################### par(las=2) plot(getVarImp(vcrf1), resolveenv=hgu95av2SYMBOL) ################################################### ### chunk number 37: ################################################### sessionInfo()