###################################################
### chunk: loadAllLibsFirst
###################################################
library("Biobase")
library("genefilter")
library("cluster")
library("hgu95av2.db")
library("annotate")
library("MASS")
library("hopach")
library("kohonen")
library("RColorBrewer")


###################################################
### chunk: loadALLdata
###################################################
library("ALL")
data(ALL)
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)
ALLfilt_bcrneg = nsFilter(ALL_bcrneg, var.cutoff=0.75)$eset


###################################################
### chunk: restrictByGO
###################################################
GOTFfun = function(GOID) {
    x = hgu95av2GO2ALLPROBES[[GOID]]
    unique(x[ names(x) != "IEA"])
}

GOIDs = c("GO:0003700", "GO:0003702", "GO:0003709", 
    "GO:0016563", "GO:0016564")

TFs = unique(unlist(lapply(GOIDs, GOTFfun)))
inSel = match(TFs, featureNames(ALLfilt_bcrneg), nomatch=0)
es2 = ALLfilt_bcrneg[inSel,]


###################################################
### chunk: manh
###################################################
iqrs = esApply(es2, 1, IQR)
gvals = scale(t(exprs(es2)), rowMedians(es2), 
    iqrs[featureNames(es2)])
manDist = dist(gvals, method="manhattan")
hmcol = colorRampPalette(brewer.pal(10, "RdBu"))(256)
hmcol = rev(hmcol)
heatmap(as.matrix(manDist), sym=TRUE,  col=hmcol,
    distfun=function(x) as.dist(x))


###################################################
### chunk: mds
###################################################
cols = ifelse(es2$mol.biol == "BCR/ABL", "black", 
    "goldenrod")
sam1 = sammon(manDist, trace=FALSE)
plot(sam1$points, col=cols, xlab="Dimension 1", 
    ylab="Dimension 2")


###################################################
### chunk: samans
###################################################
sam2 = sammon(manDist, k=3,  trace=FALSE)


###################################################
### chunk: cmdscale
###################################################
cmd1 = cmdscale(manDist)


###################################################
### chunk: usett
###################################################
rtt = rowttests(ALLfilt_bcrneg, "mol.biol")
ordtt = order(rtt$p.value)
esTT = ALLfilt_bcrneg[ordtt[1:50],]
dTT = dist(t(exprs(esTT)), method="manhattan")
sTT = sammon(dTT, trace=FALSE)


###################################################
### chunk: silcheck
###################################################
mD = as.matrix(manDist)
silEst = silcheck(mD, diss=TRUE)
silEst
mD = as.hdist(mD)
mssEst = msscheck(mD)
mssEst

d2 = as.matrix(dist(t(gvals), method="man"))
silEstG = silcheck(d2, diss=TRUE)
silEstG
mssEstG = msscheck(as.hdist(d2))
mssEstG


###################################################
### chunk: silestSol
###################################################
dsol = as.matrix(dist(gvals), method="maximum")
silcheck(dsol, diss=TRUE)
msscheck(as.hdist(dsol))


###################################################
### chunk: hclust
###################################################
hc1 = hclust(manDist)
hc2 = hclust(manDist, method="single")
hc3 = hclust(manDist, method="ward")
hc4 = diana(manDist)


###################################################
### chunk: dendro
###################################################
par(mfrow=c(4,1))
plot(hc1, ann=FALSE) 
title(main="Complete Linkage", cex.main=2)
plot(hc2, ann=FALSE)
title(main="Single Linkage", cex.main=2)
plot(hc3, ann=FALSE)
title(main="Ward's Method", cex.main=2)
plot(hc4, ann=FALSE, which.plots=2)
title(main="Divisive Clustering", cex.main=2)
par(mfrow=c(1,1))


###################################################
### chunk: exCutree
###################################################
hc13 = cutree(hc1, k=3)
hc23 = cutree(hc2, k=3)
hc33 = cutree(hc3, k=3)
hc43 = cutree(hc4, k=3)


###################################################
### chunk: comp2dists
###################################################
table(hc13, hc33)


###################################################
### chunk: cophenetic1
###################################################
cph1 = cophenetic(hc1)
cor1 = cor(manDist, cph1)
cor1
plot(manDist, cph1, pch="|", col="blue")


###################################################
### chunk: cophensol
###################################################
cph2 = cophenetic(hc2)
cor2 = cor(manDist, cph2)
cor2

cph3 = cophenetic(hc3)
cor3 = cor(manDist, cph3)
cor3

cph4 = cophenetic(hc4)
cor4 = cor(manDist, cph4)
cor4


###################################################
### chunk: kmeans
###################################################
km2 = kmeans(gvals, centers=2, nstart=5)
kmx = kmeans(gvals, centers=2, nstart=25)


###################################################
### chunk: kmvals
###################################################
names(km2)
table(km2$cluster, kmx$cluster)


###################################################
### chunk: ALLfactors
###################################################

sapply(pData(es2), function(x) is.factor(x) || 
    is.logical(x) )



###################################################
### chunk: twowaytesting
###################################################

tt1 = table(es2$mdr, km2$cluster)
test1 = chisq.test(tt1)
test1$p.value


###################################################
### chunk: pamclustering
###################################################

pam2 = pam(manDist, k=2, diss=TRUE)
pam3 = pam(manDist, k=3, diss=TRUE)


###################################################
### chunk: comppamandkm
###################################################
all(names(km2$cluster) == names(pam2$clustering)) 

pam2km = table(km2$cluster, pam2$clustering)
pam2km


###################################################
### chunk: compKMPAM
###################################################
 km3 = kmeans(gvals, centers=3, nstart=25)


###################################################
### chunk: showcompKMPAM
###################################################
 table(km3$cluster, pam3$clustering)


###################################################
### chunk: tryouter
###################################################

 outSamekm3 = outer(km3$cluster, km3$cluster, 
                    FUN = function(x,y) x==y)

 outSamepam3 = outer(pam3$clustering, pam3$clustering, 
                     FUN = function(x,y) x==y)

 inSBoth = outSamekm3 & outSamepam3

 ##then we subtract 79, because an obs is in the same 
 ## cluster as itself this just means that the diagonal 
 ## is TRUE and divide by two
 sameBoth = (sum(inSBoth) - 79)/2

 ##not in the same one, in both are those FALSE entries
 notSBoth = (!outSamekm3) & (!outSamepam3)
 notSameBoth = sum(notSBoth)/2

 ##those that are different, are TRUE in one and FALSE 
 ## in the other or vice versa
 
 diffBoth = ((!outSamekm3) & outSamepam3) | 
    (outSamekm3 & (!outSamepam3))
 discordant = sum(diffBoth)/2



###################################################
### chunk: som
###################################################
set.seed(123)
s1 = som(gvals, grid=somgrid(4,4))
names(s1)
s2 = som(gvals, grid=somgrid(4,4), alpha=c(1,0.1), 
    rlen=1000)
s3 = som(gvals, grid=somgrid(4,4, topo="hexagonal"), 
    alpha=c(1,0.1), rlen=1000)
whGP = table(s3$unit.classif)
whGP


###################################################
### chunk: SOMsol1
###################################################
table(s1$unit.classif)
table(table(s2$unit.classif))
table(table(s3$unit.classif))


###################################################
### chunk: noname
###################################################
intOnes = s1$unit.classif == 13 | s1$unit.classif == 14
gvsub = gvals[intOnes,]
gvclasses = s1$unit.classif[intOnes]
sideC = ifelse(gvclasses==13, "yellow", "blue")
heatmap(t(gvsub), ColSideCol=sideC)



###################################################
### chunk: SOMBDR
###################################################
set.seed(777)
s4 = SOM(gvals, grid=somgrid(4,4, topo="hexagonal"))
SOMgp = knn1(s4$code, gvals, 1:nrow(s4$code))
table(SOMgp)


###################################################
### chunk: useMDS
###################################################
cD = dist(s4$code)
cD


###################################################
### chunk: dropCodes
###################################################
newCodes = s4$code[-c(2,5,6,10, 15, 9, 11, 13, 14),]

SOMgp2 = knn1(newCodes, gvals, 1:nrow(newCodes))
names(SOMgp2) = row.names(gvals)
table(SOMgp2)
cD2 = dist(newCodes)
cmdSOM = cmdscale(cD2)
#plot(cmdSOM)


###################################################
### chunk: km2sol
###################################################

km2sol = kmeans(gvals, centers=4, nstart=25)
table(km2sol$cluster, SOMgp2)



###################################################
### chunk: km3sol
###################################################
dropInds = SOMgp2 %in% c("1", "4", "5")
gvals2 = gvals[!dropInds,]
km3 = kmeans(gvals2, centers=4, nstart=50)
table(km3$cluster, SOMgp2[!dropInds])



###################################################
### chunk: hopachSamples
###################################################
  samp.hobj = hopach(gvals, dmat = manDist)
  samp.hobj$clust$k


###################################################
### chunk: hopachSamsizes
###################################################
samp.hobj$clust$sizes


###################################################
### chunk: hopachGenes
###################################################
gene.dist = distancematrix(t(gvals), d = "cosangle") 
gene.hobj = hopach(t(gvals), dmat = gene.dist) 
gene.hobj$clust$k


###################################################
### chunk: clustsize
###################################################
gene.hobj$clust$sizes


###################################################
### chunk: plotsil
###################################################
silpam2 = silhouette(pam2)


###################################################
### chunk: plotsilp2
###################################################
plot(silpam2, main="")


###################################################
### chunk: plotsil3
###################################################
silpam3 = silhouette(pam3)


###################################################
### chunk: plotsilp3
###################################################
plot(silpam3, main="")


###################################################
### chunk: findNegs
###################################################
silpam3[silpam3[,"sil_width"] < 0,]


###################################################
### chunk: negSol
###################################################
  ans = silpam2[silpam2[, "sil_width"] < 0, ]


###################################################
### chunk: make4grps
###################################################

 dcl4 = cutree(hc4, 4)
 table(dcl4)
 ## we presume the labels are in the order 
 ## given to the \indexTerm{clustering} algorithm

 sild4 =  silhouette(dcl4, manDist)


###################################################
### chunk: usettpc
###################################################
rtt = rowttests(ALLfilt_bcrneg, "mol.biol")
ordtt = order(rtt$p.value)
esTT = ALLfilt_bcrneg[ordtt[1:50],]


###################################################
### chunk: plotraw
###################################################
pairs(t(exprs(esTT)[1:5,]),
    col=ifelse(esTT$mol.biol=="NEG", "green", "blue"))


###################################################
### chunk: dopc
###################################################
pc = prcomp(t(exprs(esTT)))


###################################################
### chunk: plotpc
###################################################
pairs(pc$x[,1:5], col=ifelse(esTT$mol.biol=="NEG", 
    "green", "blue"))


###################################################
### chunk: plotbi
###################################################
biplot(pc)


###################################################
### chunk: dobt eval=FALSE
###################################################
## t.test(exprs(esTT)["33232_at",]~esTT$mol.biol)


###################################################
### chunk:  eval=FALSE
###################################################
## esTT.K = ALLfilt_bcrneg[ordtt[1:K],]