###################################################
### chunk number 1: loadPacks
###################################################
library(cluster)
library(hopach)
library(ellipse)
library(lattice)
library(Biobase)
library(golubEsets)


###################################################
### chunk number 2: golub
###################################################
X<-exprs(golubTrain)

# Thresholding
X[X<100]<-100
X[X>16000]<-16000

# Filtering
M<-apply(X,1,max)
m<-apply(X,1,min)
which1<-(1:nrow(X))[(M/m)>5]
which2<-(1:nrow(X))[(M-m)>500]
geneSub<-intersect(which1,which2)
length(geneSub) ## Should get 3051
X <- X[geneSub,]

# Log transformation
X <- log2(X)

# New exprSet object
golub<-golubTrain[geneSub,]
gnames<-dimnames(X)[[1]]
slot(golub,"exprs")<-X
golub

# Class labels
Y <- golub$ALL.AML
Y<-paste(golubTrain$ALL.AML,golubTrain$T.B.cell)
Y<-sub("NA","",Y)
table(Y)


###################################################
### chunk number 3: top25F
###################################################
Fstat<-apply(exprs(golub), 1, function(z)  summary(lm(z ~ factor(Y)))$f[1])
notInf<-(1:length(Fstat))[Fstat<100]
top25F<-intersect(rev(order(Fstat)),notInf)[1:25]
affyID25<-geneNames(golub[top25F,])
affyID25


###################################################
### chunk number 4: cor
###################################################
corr<-cor(exprs(golub))
corr25<-cor(exprs(golub)[top25F,])
dimnames(corr)<-dimnames(corr25)<-list(as.vector(Y),as.vector(Y))


###################################################
### chunk number 5: corEllipseGolubAll
###################################################
plotcorr(corr, col="orange", 
main="Golub et al. (1999). Correlation matrix for n=38 learning cases \n G=3,051 genes")


###################################################
### chunk number 6: corEllipseGolubTop25
###################################################
plotcorr(corr25, col="purple", 
main="Golub et al. (1999). Correlation matrix for n=38 learning cases \n G=25 genes with largest F-statistics")


###################################################
### chunk number 7: corNumGolubAll
###################################################
plotcorr(corr, numbers=TRUE, 
main="Golub et al. (1999). Correlation matrix for n=38 learning cases \n G=3,051 genes")


###################################################
### chunk number 8: corNumGolubTop25
###################################################
plotcorr(corr25, numbers=TRUE, 
main="Golub et al. (1999). Correlation matrix for n=38 learning cases \n G=25 genes with largest F-statistics")


###################################################
### chunk number 9: corLevelGolubAll
###################################################
print(levelplot(corr, col.regions=rainbow(25), labels=as.vector(Y), 
main="Golub et al. (1999). Correlation matrix for n=38 learning cases \n G=3,051 genes"))


###################################################
### chunk number 10: corLevelGolubTop25
###################################################
print(levelplot(corr25, col.regions=heat.colors(25), labels=as.vector(Y), 
main="Golub et al. (1999). Correlation matrix for n=38 learning cases \n G=25 genes with largest F-statistics"))


###################################################
### chunk number 11: mds2and3
###################################################
Y<-paste(golub$ALL.AML,golub$T.B.cell)
Y<-sub("NA","",Y)
corr<-cor(exprs(golub))
d<-1-corr
mds2<- cmdscale(d, k=2, eig=TRUE)
mds3<- cmdscale(d, k=3, eig=TRUE)
names(mds2)
mds2$GOF
mds3$GOF


###################################################
### chunk number 12: mds2Golub
###################################################
plot(mds2$points, type="n", xlab="", ylab="", main="Golub et al. (1999). MDS, correlation matrix, G=3,051 genes, k=2")
text(mds2$points[,1],mds2$points[,2],Y, col=rank(unique(Y))[factor(Y)]+1)


###################################################
### chunk number 13: mds3Golub
###################################################
pairs(mds3$points,  main="Golub et al. (1999). MDS, correlation matrix, G=3,051 genes, k=3", pch=c("B","T","M")[rank(unique(Y))[factor(Y)]],col=rank(unique(Y))[factor(Y)]+1)


###################################################
### chunk number 14: pam2and3
###################################################
corr<-cor(exprs(golub))
d<-1-corr
dimnames(d)<-list(as.vector(Y),as.vector(Y))
set.seed(99)
pam2 <- pam(as.dist(d), k=2, diss=TRUE)
pam3 <- pam(as.dist(d), k=3, diss=TRUE)
table(pam2$clustering, Y)
table(pam3$clustering, Y)
summary(pam2)


###################################################
### chunk number 15: pamClusplot2Golub
###################################################
clusplot(d, pam2$clustering, diss=TRUE, labels=3, col.p=1, col.txt=rank(unique(Y))[factor(Y)]+1, main="Golub et al. (1999). Bivariate cluster plot for PAM \n Correlation matrix, K=2, G=3,051 genes")


###################################################
### chunk number 16: pamClusplot3Golub
###################################################
clusplot(d, pam3$clustering, diss=TRUE, labels=3, col.p=1, col.txt=rank(unique(Y))[factor(Y)]+1, main="Golub et al. (1999). Bivariate cluster plot for PAM \n Correlation matrix, K=3, G=3,051 genes")


###################################################
### chunk number 17: pamSil2Golub
###################################################
plot(pam2,which.plots=2,main="Golub et al. (1999). Silhouette plot for PAM \n Correlation matrix, K=2, G=3,051 genes") 


###################################################
### chunk number 18: pamSil3Golub
###################################################
plot(pam3,which.plots=2,main="Golub et al. (1999). Silhouette plot for PAM \n Correlation matrix, K=3, G=3,051 genes") 


###################################################
### chunk number 19: avgSilWidths
###################################################
K<-2:7
avgSil<-rep(NA, length(K))
for(k in K)
  avgSil[k-1]<-pam(as.dist(d), k=k, diss=TRUE)$silinfo$avg.width
pam4 <- pam(as.dist(d), k=4, diss=TRUE)
table(pam4$clustering, Y)


###################################################
### chunk number 20: avgSilBarplotGolub
###################################################
barplot(avgSil, names.arg=K, xlab="Number of clusters, K", ylab="Average silhouette width")


###################################################
### chunk number 21: avgSilPlotGolub
###################################################
plot(K, avgSil, pch=16, xlab="Number of clusters, K", ylab="Average silhouette width")


###################################################
### chunk number 22: hclust
###################################################
corr<-cor(exprs(golub))
d<-1-corr
dimnames(d)<-list(as.vector(Y),as.vector(Y))
hc <- hclust(as.dist(d), method="average")
hc
round(cor(cophenetic(hc),as.dist(d)),2)


###################################################
### chunk number 23: hclustPlotGolub
###################################################
plot(hc, labels=Y, main="Golub et al. (1999). Hierarchical clustering dendrogram",sub="Average linkage, correlation matrix, G=3,051 genes")
cutree(hc,k=5)
rect.hclust(hc,k=5,which=c(2,4,5),border=c("blue","green","red"))


###################################################
### chunk number 24: heatmapHeatGolub
###################################################
X<-exprs(golub)
cv<-apply(X,1, function(z) sd(z)/mean(z))
Xsub<-X[rev(order(cv))[1:50],]
dimnames(Xsub)[[2]]<-Y

args(heatmap)
heatmap(Xsub,ColSideColors=c("red","green","blue")[rank(unique(Y))[factor(Y)]])


###################################################
### chunk number 25: heatmapRainbowGolub
###################################################
heatmap(Xsub,col=rainbow(30),ColSideColors=c("red","green","blue")[rank(unique(Y))[factor(Y)]])


###################################################
### chunk number 26: cordist
###################################################
X<-exprs(golub)
dcor<-distancematrix(t(X), d="cor")        # sqrt of 1-cor
dimnames(dcor)<-list(as.vector(Y),as.vector(Y))


###################################################
### chunk number 27: hopach1
###################################################
hp1 <- hopach(t(X),dmat=dcor, K=1)
names(hp1)
names(hp1$clustering)
hp1$clustering$k # number of clusters
cl1<-hp1$clustering$labels # cluster labels
table(cl1,Y)


###################################################
### chunk number 28: dplot1Golub
###################################################
dplot(d, hp1, ord="final", showclusters=TRUE, main="Golub et al. (1999). Correlation distance matrix\n G=3,051 genes, one-level HOPACH ordering")
axis(2, labels=rev(Y[hp1$clustering$order]), at=1:38, cex.axis=0.75, col.axis=2, las=2)
axis(1, labels=Y[hp1$clustering$order], at=1:38, cex.axis=0.75, col.axis=2, las=2)


###################################################
### chunk number 29: boothopach1
###################################################
set.seed(99)
boot1<-boothopach(t(X), hp1, I=100)


###################################################
### chunk number 30: bootplot1Golub
###################################################
bootplot(boot1,hp1,ord="bootp", main="Golub et al. (1999). Barplot of bootstrap cluster membership probabilities \n Correlation distance, G=3,051 genes, one-level HOPACH ordering",showclusters=FALSE)


###################################################
### chunk number 31: hopach4
###################################################
hp4 <- hopach(t(X), d="cor")
cl4<-hp4$clustering$labels
table(cl4,Y)


###################################################
### chunk number 32: dplot4Golub
###################################################
dplot(d, hp4, col=topo.colors(15), ord="final", showclusters=TRUE, main="Golub et al. (1999). Correlation distance matrix\n G=3,051 genes, four-level HOPACH ordering")
axis(2, labels=rev(Y[hp4$clustering$order]), at=1:38, cex.axis=0.75, col.axis=2, las=2)
axis(1, labels=Y[hp4$clustering$order], at=1:38, cex.axis=0.75, col.axis=2, las=2)