---
title: "Distance metrics and clustering"
author: "Denis Puthier"
date: '`r Sys.Date()`'
output:
  html_document:
    fig_caption: yes
    highlight: zenburn
    theme: cerulean
    toc: yes
    toc_depth: 3
    toc_float: yes
  pdf_document:
    fig_caption: yes
    highlight: zenburn
    toc: no
    toc_depth: 3
  word_document: default
---


<script type="">
    $(function() {
        $(".hideshow").click(function() {
                $(this).parent().find(".exo").toggle();
        });
            
    })
</script>

<style>
.exo {
  display:none;
}

pre:not(.sourceCode) > code {
  background-color: white !important;
  color: grey !important;
} 

hr {
  color: grey !important;
}
</style>


## Introduction

Unsupervised classification methods aim at classifying objects without any prior knowledge of the class they belong to. This is a very common task that can be used to discover new structures or classes inside a given set of objects. Indeed, arranging objects based on a given criteria (that one needs to define) may lead to partitions that highlights the presence of objects of various classes. These unsupervized methods are widely used in genomics where objects (expression profiles, biological samples with various descriptor, sequences, motifs, genomic features,...)  need to be classified in order to mine large dataset in search for similar objects related to known or novel classes and displaying jointly particular properties. These approaches have been  applied in several reference papers in the context of biological sample classification leading to the discovery of novel tumor classes. 

The expected results of the classification is a partition. A partition can be defined as the division of a collection of objects into subsets or clusters. In clustering approaches, each cluster is expected to contain very similar object (low within-group variance) while object from different clusters are expected to differ (high inter-group variance).  These clusters can correspond to a simple division of the object collection but also to a  set of hierarchically nested subsets. In the most popular approaches the clusters produced are non overlapping, meaning that each object is assign to only one cluster. These approaches are called *hard* (or *crisp*) *clustering* methods in contrast to *fuzzy clustering* (or soft clustering) in which one object can fall into several clusters. Although rather complexe approaches have been proposed, the most popular methods for unsupervised clustering are based on rather simple algorithms that rely on simple mathematical basis. 
Whatever the method used, one important aspect is to choose a reasonable metric to assess the similarity between object. Although other parameters may strongly influence the partitioning result, the choice of the metric will have a major influence on clustering. They are lots of existing metrics that can be applied to a set of objects. In the next section, we will focus on some of them that are frequently encountered in the context of genomic data analysis. 

## Choosing a metric

### The euclidean distance.

One of the most classical metric is the well-known *euclidean distance*. Let's imagine two genes in a two dimensional space that could represent two biological samples. Thus these two genes could be represented as two points $a$ and $b$. The euclidean distance between these two genes can be represented in the sample space and corresponds to the physical distance between two points computed using pythagora's formula with $i=1,..,p$ corresponding to the samples:

$$d(a,b)=\sqrt{\sum_{i=1}^p (a_{i} - b_{i})^2}$$

We could also propose an alterntive representation in which the points would be the samples and the dimensions would correspond to genes. We could also propose an alternative representation in which the $x$ axis would represent the samples and the $y$ axis would represent the intensities. In this case we would represent the profiles of the two genes across the two samples $s1$ and $s2$. This three representations are depicted below.
In this particular case, the two points, $a$ and $b$ for which the distance is to be computed are in a two-dimensional space. However, this distance can be generalized to any space with p dimensions. Let's take an example with two points (genes) in a eight dimensional space. We will chose the third representation with samples names displayed on $x$ axis and intensities displayed on $y$ axis.The value for $a_{i}^2 - b_{i}^2$ are displayed with dashed lines

```{r, dpi=300}
# Preparation x window
col.vec <- c("black","red")
op <- par(no.readonly = TRUE)
par(mfrow=c(2,2), cex.main=0.7, mgp=c(2,1,0), mai=c(0.5,0.5,0.5,0.5))

# Genes as points and samples as dimensions
a <- c(1, 4)
b <- c(2, 3)
m <- rbind(a, b)
print(m)
colnames(m) <- c("s1", "s2")
plot(m , pch=16, xlim=c(0,3), 
     ylim=c(0,5), main="Genes as points and samples as dimensions",
     col=col.vec,
     panel.first=grid(lty=1))
suppressWarnings(arrows(a[1], a[2], b[1], b[2], angle=0))
suppressWarnings(arrows(a[1], a[2], b[1], a[2], angle=0, lty=3))
suppressWarnings(arrows(b[1], a[2], b[1], b[2], angle=0, lty=3))
text(m, lab=rownames(m), pos = 2, offset = 1)
text(7, 4,label="dist(a,b)")


# Sample as points and genes as dimensions
plot(t(m) , pch=16, 
     ylim=c(0,5), xlim=c(0,5), 
     main="Sample as points and genes as dimensions",
     panel.first=grid(lty=1))
text(t(m), lab=colnames(m), pos = 2, offset = 1)


# x axis correspond to samples and the y axis represent the intensities

matplot(t(m) , 
       ylim=c(0,5), xlim=c(0,3), 
       main="x axis for samples (n=2) and y axis for intensities", 
       ylab="Intensities",
       xlab="samples",
       xaxt = "n",
       type="n")
grid(lty=1)

axis(1, 0:4, c("", "s1", "s2", "", ""))
suppressWarnings(arrows(1, a[1], 2, a[2], angle=0))
suppressWarnings(arrows(1, b[1], 2, b[2], angle=0))

matpoints(t(m) , pch=16, 
       col=col.vec)


# 8 dimensions: x axis correspond to samples and the y axis represent the intensities
a <- c(7, 7, 7, 7, 6, 6, 10, 10)
b <- c(8, 2, 5,  6, 1, 6, 1, 4)
m <- rbind(a, b)
matplot(t(m),
        xlim=c(0,10), ylim=c(0,12),
        main="x axis for samples (n=8) and y axis for intensities",
        pch=16,
        col="black",
        ylab="Intensities",
        xlab="samples",
        lty=1,
        type="n")
grid(lty=1)

for(i in 1:length(a)){
  suppressWarnings(arrows(i, a[i] , i, b[i], angle=0, col=col.vec[i], lty=3))
}

matpoints(t(m),
        pch=16,
        type="b",
        lty=1)

points(a, type="p", col=col.vec[1],  pch=16)
points(b, type="p", col=col.vec[2],  pch=16)
```

**Exercise:** 

> - Compute the euclidean distance between a and b. 
> - Check that the same result is obtained with the **dist()** function used with default arguments.
> - Now assign the values of a to b.
>   - What is the computed distance.
> - Add successively 10, 20, 30, 40 to each element/dimension of $b$. 
>   - Compute the distance.
>   - What do you observe ?
> - Set $b$ values to $-a + mean(a)$ then a to $a - mean(a)$. 
>   - What is the distance obtained.
> - Now create a variable c with values $a + 3$
>   - What is the distance obtained with a.
>   - Plot the profile of a, b and c.
>   - Is the euclidean distance a good choice for comparing these profiles ? Discuss it. In which context could it be useful ?

<div class="hideshow"> << Hide | Show >> </div>
<div class="exo">
```{r}
par(op)

## The distance between a and b
(dist.euc  <- sqrt(sum((a-b)^2)))

# This can be also computed using the dist() function (requires) a matrix.
(dist.obj <- dist(rbind(a,b), method = "euclidean"))

# dist() return a dist object which can be converted to a matrix
(dist.mat <- as.matrix(dist.obj))

# Both lead to the same result
all.equal(dist.euc, dist.mat[1,2])

# Now assign the values of a to b.
b <- a
dist(rbind(a,b), method = "euclidean") # the distance is 0

# Add successively 10, 20, 30, 40 to each element/dimension of $b$.
range <- c(0, 10, 20, 30, 40)
b <- matrix(b, nc=8, nrow=5, byrow = TRUE)
b <- b + range
b <- rbind(b,a)
rownames(b) <- c("b+0","b+10", "b+20", "b+30", "b+40", "a")
plot(NA, xlim=c(1,8), 
     ylim=c(0,max(b)), type="n", 
     panel.first = grid(),
     ylab="Intensities",
     xlab="Samples")
matpoints(t(b), type="b", pch=16, lty=1, lwd=2)

# What we can observed is that although a and b have exactly the same profile 
# their distance increase as the absolute level of b is increasing.
# The euclidean distance has no upper limit and can thus increase indefinitely.
dist.a.b <- as.matrix(dist(b))[1:5,"a"]
text(2, range+2, labels=paste("d(a,", names(dist.a.b), ")=", round(dist.a.b,2), sep=""), cex=0.7, pos=3, offset=2)

# set $b$ values to -a + mean(a) then a to a - mean(a).
b <- -a + mean(a)
a <-  a - mean(a)
dist(rbind(a,b))

# Create vector c with values a + 6

c <- a + 3
dist(rbind(a,c))

matplot(t(rbind(a,b,c)), type="b", pch=1, lty=1, lwd=2, col=c("blue", "orange","red"))
grid()
legend(x="topleft", 
       legend = c("a","b","c"), 
       col=c("blue", "orange","red"), 
       lwd=2)

# As you can see the euclidean distance is unable to distinguish between profile displaying coordinate variations and profiles displaying clear anti-correlation.
```
</div>

### Pearson's correlation coefficient

The Pearson's correlation coefficient is a widely used statistical score to measure the co-variance of two variables. It is particularly used in the case of hierarchical clustering applied to microarray or RNA-Seq data. The Pearson's correlation coefficient can be computed by dividing the covariance of the two variables by the product of their standard deviations. 


$$r_{x,y}= \frac{1}{n}\sum\limits_{i=1}^n   \frac{(x_i - \bar{x})(y_i - \bar{y})}{s_xs_y} = \frac{cov(x,y)}{s_xs_y}$$


Let's go back to the example we provided above. We had three objects ($a$, $b$ and $c$) having values associated with a set of variables. These object could be viewed as a small set of genes having given expression values in a set of biological samples. We will construct a new diagram in which we will display for each gene (a, b, and c) through each sample (8 dimensions) the difference to the corresponding mean value ($x_i - \bar{x}$). The covariance is the sum of the products of the differences through the eight dimensions. If for in dimension the intensities of two gene goes jointly above their corresponding mean, this product will be positive. In the same way, if their intensities goes jointly above their corresponding mean, this product will be positive. It will be negative or zero in the other cases. The fact that the covariance is divided by the product of the standard deviation will ensure that the Pearson's $r$ coefficient is bound between -1 (the more negative corresponding to anti-correlation) and 1 (complete correlation). 

```{r}
col.vec <- c(col.vec, "green")
# The dataset
a <- c(7, 7, 7, 7, 6, 6, 10, 10)
a <- a - mean(a)
b <- -a + mean(a)
c <- a + 3
m <- rbind(a,b,c)

# The mean intensities
a.mean <- mean(a)
b.mean <- mean(b)
c.mean <- mean(c)

# The diagram with the three genes: a, b, c
matplot(t(m),
       main="x axis for samples (n=8) and y axis for intensities",
        pch=16,
        col="black",
        lty=1,
        type="n")
grid(lty=1)

mean.all <- c(a.mean, b.mean, c.mean)
abline(h=mean.all, lty=1, col="darkgrey")

gene.list <- list(a,b,c)

# The differences to the means
for(g in 1:length(gene.list)){
  for(i in 1:length(a)){
    arrows(i, mean.all[g] , i, gene.list[[g]][i], angle=0, lty=2)
  }
}

matpoints(t(m),
        pch=16,
        type="b",
        lty=1)

```

**Exercise:** 

> - Compute the Pearson's correlation coefficient between $a$ and $c$ using the formula.
>   - Is that results expected.
>   - Compute it using the cor() function. Is the same result obtained ?
> - What is the result obtained when comparing $a$ and $b$ ?
> - Is this metric a good choice for comparing two gene expression profile ?

<div class="hideshow"> << Hide | Show >> </div>
<div class="exo">

```{r}
n <- 8
# Compute sample standard deviation
sd.a <- sqrt(1/n * sum((a-a.mean)^2))
sd.c <- sqrt(1/n * sum((c-c.mean)^2))

# Compute Co-variance
cov.a.c <- 1/n*sum((a-a.mean) * (c-c.mean))
cor.a.c <- cov.a.c/(sd.a*sd.c)
cor.a.c 
all.equal(cor.a.c, cor(a,c))

# Correlation between a and b. 
# a and b are anti-correlated
cor(a,b)
```
</div>

This correlation coefficient can be simply transformed into a distance using the following formula:

$$d(x,y) = 1/2(1-r_{x,y})$$

For any gene pair with a Pearson's correlation coefficient of 1 the corresponding distance will be 0 while those displaying value of -1 will have an associated distance of 1. It should be noted that in contrast to euclidean distance, this distance does not follow the rule of  triangle inequality that states that for any triangle the sum of the length of two given sides  must be greater than the length of the third side. This type of distance can thus not be used to represent gene distances in a euclidean space. 

### Spearman's correlation coefficient $\rho$

Spearman's correlation coefficient can be also particularly useful to assess the similarity of two given expression profile. While this rank-based measure is more robust to outliers than Pearson's correlation coefficient, it is also less sensitive. While Pearson's correlation coefficient is constructed based on raw values, Spearman's correlation coefficient will assess the concordance between the ranks of these values. The formula of the Speaman 's correlation coefficient is provided below. In this formula $d$ is the difference between ranks.

$$\rho = {1- \frac {6 \sum d_i^2}{n(n^2 - 1)}}$$

As for the Pearson's correlation coefficient, $\rho$ can be transformed to a distance using:

$$d(x,y) = 1/2(1-\rho_{x,y})$$

### Comparing Spearman and Pearson correlation coefficient results

Also lot of other metric exist, Spearman and Pearson correlation coefficient (or their associated distance) are frequently a good choice when working with expression data. However, depending on the context, it is advisable to use the appropriate coefficient since their results may greatly differ. One can keep in mind that Pearson will capture linear relationships while Spearman will capture monotonic relationships. Let's take some example and compute the results of both coefficient.

```{r}
op <- par(no.readonly = TRUE)
par(mfrow = c(2, 2), mar = c(3, 4, 0, 0.2)) #it goes c(bottom, left, top, right) 
## An example where the result of both approaches are consistant
set.seed(1)
a <- rnorm(20)
b <- a + rnorm(20,sd=0.1)
plot(a,b, pch=16, las=1, panel.first=grid(lty=1))
cor.a.b.pear <- paste("r=", round(cor(a,b), 3))
cor.a.b.spear <- round(cor(a,b, method="spearman"), 3)
cor.a.b.spear <- paste("\u03C1=", cor.a.b.spear, sep="")
legend("topleft", leg=c(cor.a.b.spear, 
                           cor.a.b.pear ))

## Another example where the result of both approaches are consistant
set.seed(1)
a <- rnorm(20)
b <- -a + rnorm(20,sd=0.1)
plot(a,b, pch=16, las=1, panel.first=grid(lty=1))
cor.a.b.pear <- paste("r=", round(cor(a,b), 3))
cor.a.b.spear <- round(cor(a,b, method="spearman"), 3)
cor.a.b.spear <- paste("\u03C1=", cor.a.b.spear, sep="")
legend("topright", leg=c(cor.a.b.spear, 
                           cor.a.b.pear ))

## An example with non linear relationship
a <- 1:50
b <- exp(a)
plot(a,b, pch=16, las=1, panel.first=grid(lty=1))
cor.a.b.pear <- paste("r=", round(cor(a,b), 3))
cor.a.b.spear <- round(cor(a,b, method="spearman"), 3)
cor.a.b.spear <- paste("\u03C1=", cor.a.b.spear, sep="")
legend("topleft", leg=c(cor.a.b.spear, 
                           cor.a.b.pear ))

## An example with an outlier.
## Rank concordance is observed only in one dimension
set.seed(1)
a <- rnorm(10)
b <- rnorm(10)
rank(a) - rank(b)
a[5] <- a[5] + 10
b[5] <- b[5] + 10
plot(a,b, pch=16, las=1, panel.first=grid(lty=1))
cor.a.b.pear <- paste("r=", round(cor(a,b), 3))
cor.a.b.spear <- round(cor(a,b, method="spearman"), 3)
cor.a.b.spear <- paste("\u03C1=", cor.a.b.spear, sep="")
legend("topleft", leg=c(cor.a.b.spear, 
                           cor.a.b.pear ))

```


## Unsupervized clustering of microarray data



Here we will use the GSE13425 experiment which which was retrieved from the Gene Expression Omnibus ([GEO](http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE13425)) public database. In this experiment, the authors were interested in the molecular classification of acute lymphoblastic leukemia (ALL) that are characterized by the abnormal clonal proliferation, within the bone marrow, of lymphoid progenitors blocked at a precise stage of their differentiation.

Data were produced using Affymetrix geneChips (Affymetrix Human Genome U133A Array, HGU133A). Informations related to this platform are available on GEO website under identifier GPL96.

### Loading data into R

Start R. Have a look at the description of the *read.table* function. Load the expression matrix (GSE13425_Norm_Whole.txt), the A/P/M matrix (GSE13425_AMP_Whole.txt) and phenotypic data into R using the *read.table* function (assign the results to objects named data, amp, and pheno respectively).

| Name and URL |  Description |
|-----------------|-----------------------------------------------|
| [Normalized data](http://pedagogix-tagc.univ-mrs.fr/courses/ASG1/data/marrays/GSE13425_Norm_Whole.txt) | The file GSE13425_Norm_Whole.txt contains genes as rows and samples as columns. Data were previously normalized using rma algorithm (they are thus transformed in logarithm base 2). |
| [AMP calls](http://pedagogix-tagc.univ-mrs.fr/courses/ASG1/data/marrays/GSE13425_AMP_Whole.txt) | The GSE13425_APM_Whole.txt file contains information about A/P/M calls (genes as rows and samples as columns). |
| [Phenotypic data](http://pedagogix-tagc.univ-mrs.fr/courses/ASG1/data/marrays/phenoData_GSE13425.tab) | The phenoData_GSE13425.txt file contains phenotypic data about samples. |



<div class="hideshow"> << Hide | Show >> </div>
<div class="exo">
```{r}
#?read.table
data.url <- "http://pedagogix-tagc.univ-mrs.fr/courses/ASG1/data/marrays/GSE13425_Norm_Whole.txt"
amp.url <- "http://pedagogix-tagc.univ-mrs.fr/courses/ASG1/data/marrays/GSE13425_AMP_Whole.txt"
pheno.url <- "http://pedagogix-tagc.univ-mrs.fr/courses/ASG1/data/marrays/phenoData_GSE13425.tab"
data <-  read.table(data.url,sep="\t", head=T, row=1)
amp <- read.table(amp.url,sep="\t", head=T, row=1)
pheno <- read.table(pheno.url,sep="\t", head=T, row=1)
```
</div>



 

### Selecting a subset of genes

#### Selecting using the A/P/M criteria

First, we will select genes giving a significant signal in a given number of samples.

Select genes which are call present in at least 10% of the samples (n= 10724 genes).

<div class="hideshow"> << Hide | Show >> </div>
<div class="exo">
```{r}
isPresent <- amp == "P"
ind <- rowSums(isPresent) >= 19  
data <- data[ind, ]
```
</div>


#### Selecting using standard deviation. 

As the classification of the whole gene matrix is rather computer intensive we will select 30% of the genes based on standard deviation.

 - Select these genes.
 - Change column names so that the new matrix will contains information about sample types.
 - Write data onto disk (file GSE13425_sub_1.txt).



<div class="hideshow"> << Hide | Show >> </div>
<div class="exo">
```{r}
sd <- apply(data,1,sd)
summary(sd)
quantile(sd,0.7)
data <- data[sd > quantile(sd,0.7),]
colnames(data) <- paste(pheno$Sample_title, colnames(data), sep="| |")
dir.create("~/distance_and_clustering",showWarnings = FALSE)
setwd("~/distance_and_clustering")
write.table(data, "GSE13425_sub_1.txt", sep="\t", quote=F, col.names=NA)
```
</div>

### About phenotypic information

Here are the available phenotypic information. 

```{r}
colnames(pheno)
```

We will store some of these information in two vectors:

- One contain phenotype title and the GEO IDs as labels.
- One contain a set of colors related to the phenotypes and the GEO IDs as labels.

```{r}
pheno.vec <- as.character(pheno$Sample.title)
names(pheno.vec) <- as.character(rownames(pheno))
pheno.vec.col <- as.character(pheno$sample.colors)
names(pheno.vec.col) <- as.character(rownames(pheno))
```

### Hierarchical clustering with hclust

Euclidean distance is rarely used in the context of microarray analysis. A distance based on Pearson's correlation coefficient is most generally preferred (Spearman's rank correlation coefficient may also be used). Let's visualize the sample-sample correlation matrix using a heatmap.

<div class="hideshow"> << Hide | Show >> </div>
<div class="exo">
```{r}
pear <- cor(data, method="pearson")
palette <-colorRampPalette(c("yellow", "black","blueviolet"))
library(lattice)
levelplot(pear,col.regions=palette, scales=list(cex=0.2))

# we can also store the result as a high quality pdf file
pdf("coor.pdf"); levelplot(pear,col.regions=palette, scales=list(cex=0.2)); dev.off()

```
</div>



The Pearson's correlation coefficient is bounded between -1 and 1. We can transform it into a distance using the following command:


```{r}
pear.dist <- as.dist((1-pear)/2)
```

Using this distance matrix use the hclust function to perform hierarchical clustering of samples.

<div class="hideshow"> << Hide | Show >> </div>
<div class="exo">

```{r}
hp <- hclust(pear.dist, method="average")
plot(hp, hang=-1, labels=as.character(pheno$Sample.title), cex=0.2)

# we can also store the result as a high quality pdf file
pdf("hp.pdf")
plot(hp, hang=-1, labels=as.character(pheno$Sample.title), cex=0.2)
dev.off()
#system("evince hp.pdf&")
```

</div>

The objects resulting from *hclust* calls can also be used to construct a heatmap. This can be done for instance using the *heatmap.3* function. Again, take care to the arguments used to construct the tree.  

<div class="hideshow"> << Hide | Show >> </div>
<div class="exo">

```{r}

library("GMD")
library("RColorBrewer")

# Median center rows
data.med.center <- sweep(data, MARGIN=1, STATS=apply(data,1,median), FUN = "-",) 
# floor and ceil

data.med.center[ data.med.center > 5 ] <- 5
data.med.center[ data.med.center < -5 ] <- -5

# Compute hierarchical clustering trees for rows and columns
pear.cols <- cor(t(data.med.center), method="pearson")
pear.dist <- as.dist((1-pear.cols)/2)
hc.cols <- hclust(pear.dist, method="average")
pear.rows <- cor(data.med.center, method="pearson")
pear.dist <- as.dist((1-pear.rows)/2)
hc.rows <- hclust(pear.dist, method="average")

# Prepare a function that returns interpolated color sets
col.fun = colorRampPalette(c("green", "black", "red"))



library("GMD")
library("RColorBrewer")

x11(height=14, width=14)
heatmap.3(t(data.med.center), 
          scale="none", 
          RowIndividualColors = pheno.vec.col, 
          cluster.by.col=TRUE, 
          cluster.by.row=TRUE, 
          Rowv=TRUE,
          Colv=TRUE,
          hclust.col = hc.cols, 
          hclust.row = hc.rows, 
          labRow=pheno.vec,
          labCol=NULL,
          cexRow = 0.5,
          color.FUN = col.fun)

dev.off()
```

</div>

However, one can see the limits of using R in this context...

### Hierarchical clustering with the cluster and treeview software

R is not particularly well-suited to visualize classification results for very large datasets. We will thus install the cluster and treeview software that are very handy to browse the results of a hierarchical clustering.

<div class="hideshow"> << Hide | Show >> </div>
<div class="exo">
```{bash highlight=TRUE, eval=FALSE}
# These are shell/bash commands !
mkdir -p ~/bin
cd ~/bin
wget http://bonsai.hgc.jp/~mdehoon/software/cluster/cluster-1.50.tar.gz
tar xvfz cluster-1.50.tar.gz
cd cluster-1.50
./configure --without-x
make
echo -e "\nalias cluster=$PWD/src/cluster" >> ~/.bashrc
cd ..
wget http://sourceforge.net/projects/jtreeview/files/jtreeview/1.1.6r2/TreeView-1.1.6r2-bin.tar.gz
tar xvfz TreeView-1.1.6r2-bin.tar.gz
echo "alias javatreeview='java -jar  $PWD/TreeView-1.1.6r2-bin/TreeView.jar'" >> ~/.bashrc
source  ~/.bashrc
cd -
```
</div>

Using the cluster software, run hierarchical clustering on both genes and samples.


<div class="hideshow"> << Hide | Show >> </div>
<div class="exo">
```{r engine='bash', highlight=TRUE, eval=FALSE}
cluster -f GSE13425_sub_1.txt -g 2 -e 2 -m a  -cg m
javatreeview
```
</div>


Now you can open the .cdt file, that was produced by cluster using javatreeview. How are the samples classified ? What can you say about gene classification ?

### Kmeans clustering with R

Look at the tree produced by Treeview. How many classes would you guess from this tree. Use the *Kmeans()* function from the amap package to perform a kmeans with the number of classes you guessed from the hierachical clustering. Look carefully at the parameters. What is the default distance used here ? Is that the right distance to use in this context ? Use the *matplot()* function to draw the profiles of each classes. Compare you results with other students (that may have choosen a different number of classes).

<div class="hideshow"> << Hide | Show >> </div>
<div class="exo">

```{r}
## Reading the file (if required)
data <- read.table("GSE13425_sub_with_pheno.txt", sep="\t", head=T, row.names = 1)

## Performing Kmeans
library(amap)
k <- 9
km <- Kmeans(data, centers=k, method="pearson", iter.max = 100)

par(mfrow=c(3,3))
for(i in 1:k){

  cluster_k <- data[names(km$cluster[km$cluster==i]),]
  cluster_k <- sweep(cluster_k, 
                     MARGIN = 1, 
                     FUN = "-", 
                     STATS = apply(cluster_k, 1, median))
  matplot(t(cluster_k), 
          type="l", 
          col="gray", 
          ylab="Intensity",
          ylim=c(-8,8),
          main=paste("Cluster ", i, " #gene = ", nrow(cluster_k)))
  matpoints(colMeans(cluster_k), type="l", col="blue")
}

```
</div>



Perform the same clustering using the kmeans version implemented in cluster 3.0. Load the results in Java Treeview.

**NB**: if required, the filtered dataset can be obtained [here](GSE13425_sub_with_pheno.txt).

