###################################################
### 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()