###################################################
### chunk number 1: 
###################################################
options(width=100) 
 


###################################################
### chunk number 2: 
###################################################
library(ALL) 
data(ALL) 
help(ALL) 
cl <- substr(as.character(ALL$BT), 1, 1)  
table(cl)
dat <- exprs(ALL)  
 


###################################################
### chunk number 3: 
###################################################
d <- dist(t(dat)) 
image(as.matrix(d))  
hc <- hclust(d, method="complete")  
plot(hc, labels=cl)  
 


###################################################
### chunk number 4: 
###################################################
plot(hc, labels = cl)  


###################################################
### chunk number 5: 
###################################################
groups <- cutree(hc, k = 2)  
table(groups,cl)  
fisher.test(groups,cl)$p.value  
 


###################################################
### chunk number 6: 
###################################################
genes.var <- apply(dat,1,var)  
genes.var.select <- order(genes.var, decreasing=T)[1:100]  
dat.s <- dat[genes.var.select, ]  
d.s <- dist(t(dat.s)) 
hc.s <- hclust(d.s,method="complete")  
plot(hc.s,labels=cl)  
 


###################################################
### chunk number 7: 
###################################################
plot(hc.s,labels=cl)  


###################################################
### chunk number 8: 
###################################################
groups.s <- cutree(hc.s,k=2)  
table(groups.s,cl)  
fisher.test(groups.s,cl)$p.value 
 


###################################################
### chunk number 9: 
###################################################
genes.random.select <- sample(nrow(dat),100)  
dat.r <- dat[genes.random.select,]  
d.r <- dist(t(dat.r)) 
image(as.matrix(d.r))  
hc.r <- hclust(d.r,method="complete")  
plot(hc.r,labels=cl)  
groups.r <- cutree(hc.r,k=2)  
table(groups.r,cl)  
fisher.test(groups.r,cl)$p.value  
 


###################################################
### chunk number 10: 
###################################################
par(mfrow=c(1,3))  
image(as.matrix(d)) 
image(as.matrix(d.r)) 
image(as.matrix(d.s)) 


###################################################
### chunk number 11: 
###################################################
library(cluster) 
par(mfrow=c(2,1))  
clusplot(t(dat.r),groups.r,main="100 random genes",color=T)  
clusplot(t(dat.s),groups.s,main="100 high-variance genes",color=T)  


###################################################
### chunk number 12: 
###################################################
library(cluster) 
par(mfrow=c(2,1))  
clusplot(t(dat.r),groups.r,main="100 random genes",color=T)  
clusplot(t(dat.s),groups.s,main="100 high-variance genes",color=T)  


###################################################
### chunk number 13: 
###################################################
k <- 2  
withinss <- Inf  
for (i in 1:10) {  
  kmeans.run <- kmeans(t(dat),k) 
  print(sum(kmeans.run$withinss)) 
  print(table(kmeans.run$cluster,cl)) 
  cat("----\n") 
  if (sum(kmeans.run$withinss) < withinss) {  
    result <- kmeans.run  
    withinss <- sum(result$withinss) 
        } 
}  
table(result$cluster,cl)  
 


###################################################
### chunk number 14: 
###################################################
kmeans.s <- kmeans(t(dat.s),k)  
table(kmeans.s$cluster,cl)  


###################################################
### chunk number 15: 
###################################################
result <- pam(t(dat),k)  
table(result$clustering,cl)  
 


###################################################
### chunk number 16: 
###################################################
ngenes <- 2:50  
o <- order(genes.var,decreasing=T)  
miscl <- NULL  
for (k in ngenes) {  
  dat.s2 <- dat[o[1:k],]  
  pam.result <- pam(t(dat.s2),k=2)  
  ct <- table(pam.result$clustering,cl)  
  miscl[k] <- min(ct[1,2]+ct[2,1],ct[1,1]+ct[2,2]) }
xl = "# genes"
yl = "# misclassification"
plot(ngenes,miscl[ngenes],type="l",xlab=xl,ylab=yl)  
 


###################################################
### chunk number 17: 
###################################################
totalsum <- sum(diag((ncol(dat.s)-1)*cov(t(dat.s))))  
withinss <- numeric() 
withinss[1] <- totalsum 
for (k in 2:20) {  
  withinss[k] <- sum(kmeans(t(dat.s),k)$withinss) }  
plot(1:20,withinss,xlab="# clusters",ylab="objective function",type="b")  
 


###################################################
### chunk number 18: 
###################################################
asw <- numeric()  
for (k in 2:20) {  
  asw[k] <- pam(t(dat.s),k)$silinfo$avg.width }  
plot(1:20,asw,xlab="# clusters",ylab="average silhouette width",type="b")  
plot(silhouette(pam(t(dat.s),2)))  
 


###################################################
### chunk number 19: 
###################################################
plot(silhouette(pam(t(dat.s),2)))  


###################################################
### chunk number 20: 
###################################################
plot(silhouette(pam(t(dat.s),3)))  


###################################################
### chunk number 21: 
###################################################
d.s <- dist(t(dat.s))  
cl.true <- as.numeric(cl == "B") + 1  
asw <- mean(silhouette(cl.true,d.s)[,3])  


###################################################
### chunk number 22: 
###################################################
asw.random <- rep(0,20)  
for (sim in 1:20) {  
  cl.random = sample(cl.true)  
  asw.random[sim] = mean(silhouette(cl.random,d.s)[,3]) }  
symbols(1,asw,circles=0.01,ylim=c(-0.1,0.4),inches=F,bg="red")  
text(1.2,asw,"Average silhouette value\n for true labels")  
boxplot(data.frame(asw.random),col="blue",add=T)  
 


###################################################
### chunk number 23: 
###################################################
library(hgu95av2) 
gene.names <- sapply(featureNames(ALL),get,hgu95av2SYMBOL) 
gene.names <- gene.names[genes.var.select] 
d.g = dist(dat.s)  
hc = hclust(d.g,method="complete")  
plot(hc,labels=gene.names)  
 


###################################################
### chunk number 24: 
###################################################
plot(hc,labels=gene.names)  
 


###################################################
### chunk number 25: 
###################################################
asw = numeric()  
for (k in 2:20) {  
  asw[k] = pam(dat.s,k)$silinfo$avg.width }  
plot(1:20,asw,xlab="# clusters",ylab="average silhouette width",type="b")  
plot(silhouette(pam(dat.s,2)))  
 


###################################################
### chunk number 26: 
###################################################
plot(1:20,asw,xlab="# clusters",ylab="average silhouette width",type="b")  


###################################################
### chunk number 27: 
###################################################
plot(silhouette(pam(dat.s,2)))  


###################################################
### chunk number 28: 
###################################################
help(heatmap) 
heatmap(dat.s) 
 


###################################################
### chunk number 29: 
###################################################
heatmap(dat.s) 


###################################################
### chunk number 30: 
###################################################
select = c(which(cl=="B")[1:5],which(cl=="T")[1:5])  
dat.small = dat[,select]  
#  
cl.random = sample(cl[select])  
library(multtest)  
tstat = mt.teststat(dat.small,classlabel=(cl.random=="B"),test="t")  
genes = order(abs(tstat),decreasing=T)[1:100]  
dat.split = dat.small[genes,]  
#  
d = dist(t(dat.split))  
hc = hclust(d,method="complete")  
par(mfrow=c(1,2))  
plot(hc,labels=cl.random,main="Labels used for gene selection",xlab="")  
plot(hc,labels=cl[select],main="True labels",xlab="")  


###################################################
### chunk number 31: 
###################################################
par(mfrow=c(1,2))  
plot(hc,labels=cl.random,main="Labels used for gene selection",xlab="")  
plot(hc,labels=cl[select],main="True labels",xlab="")  


###################################################
### chunk number 32: 
###################################################
list.distances = as.matrix(d.g)[5,]
o = order(list.distances)
list.distances[o][1:10]
gene.names[o][1:10]



###################################################
### chunk number 33: 
###################################################
hc = hclust(d.g, method = "complete")
cutree(hc,3)[o][1:10]



###################################################
### chunk number 34: 
###################################################
library(xtable)
used.version = sapply(sub("package:","", grep("package:", search(), value=TRUE)), package.version)
xtable(data.frame(used.version))