#########################################################
## Clustering Demo                                     ## 
## Author: Han-Ming Wu (hmwu@stat.sinica.edu.tw)       ## 
## Institute of Statistical Science, Academia Sinica   ## 
## http://www.sinica.edu.tw/~hmwu/                     ## 
## 2005/07/27                                          ## 
#########################################################


## Read Data
setwd("C:\\Program Files\\R\\rw2001\\WorkingData")
library(stats)
cell.matrix  <- read.table("TradCellCycle103_alpha.txt", header=TRUE)
n <- dim(cell.matrix)[1]
p <- dim(cell.matrix)[2]-2
cell.data <- cell.matrix[,3:p+2]
gene.name <- cell.matrix[,1]
gene.phase <- cell.matrix[,2]
phase <- unique(gene.phase)
phase.name <- c("G1", "S", "S/G2", "G2/M", "M/G1")

## standardized data
cell.sdata <- (cell.data-apply(cell.data, 1, mean))/sqrt(apply(cell.data, 1, var)) 

no.group <- 5
no.iter <- 20
cell.kmeans <- kmeans(cell.sdata, no.group, no.iter)
cell.kmeans 

plot(cell.sdata[,1:4], col = cell.kmeans$cluster)

## PCA with k-means
plot(pca.dim1, pca.dim2,  
main="PCA for Cell Cycle Data with K-means Clustering", xlab="PCA-1", ylab="PCA-2", col=cell.kmeans$cluster)

## MDS with k-means
plot(mds.dim1, mds.dim2, xlab="MDS-1", ylab="MDS-2", 
main="MDS for Cell Cycle Data with K-means Clustering", col = cell.kmeans$cluster)

## SOM
library(som)
cell.som <- som(cell.sdata, xdim=5, ydim=4, topol="rect", neigh="gaussian")
plot(cell.som) 

## Hierarchical Clustering on genes
cell.gene.hc.ave <- hclust(dist(cell.sdata), method = "ave")
plot(cell.gene.hc.ave, hang = -1, cex=0.5, labels=gene.name)

## Hierarchical Clustering on experiments
cell.exp.hc.ave <- hclust(dist(t(cell.sdata)), method = "ave")
plot(cell.exp.hc.ave, cex=0.8)

## NOTE
#plot(cell.hc.ave)
#cell.hc.complete <- hclust(dist(cell.sdata), method = "complete")