---
title: "Exercices en R - Analyse des données de biopuces de Den Boer et al. (2009)"
author: "Jacques van Helden & Annie Broglio"
date: '`r Sys.Date()`'
output:
  pdf_document:
    fig_caption: no
    highlight: zenburn
    toc: yes
    toc_depth: 2
  md_document:
    variant: markdown_github
  html_document:
    code_folding: show
    fig_caption: no
    highlight: zenburn
    theme: cerulean
    toc: yes
    toc_depth: 3
    toc_float: yes
  word_document:
    toc: no
    toc_depth: 2
---

<!-- bibliography: ../../bibliography/ASG1_references.bib--->


```{r knitr setup, include=FALSE,  eval=TRUE, echo=FALSE, warning=FALSE}
# knitr::opts_chunk$set(echo=TRUE, eval=TRUE, cache=TRUE, message=FALSE, warning=FALSE, comment="")

library(knitr)
library(gplots)

options(width=300)
knitr::opts_chunk$set(
  fig.width = 7, 
  fig.height = 5, 
  fig.path='figures/',
  fig.align = "center", 
  size = "tiny", 
  echo = TRUE, eval=TRUE, 
  warning = FALSE, message = FALSE, 
  results = TRUE, comment = "")
# knitr::asis_output("\\footnotesize")

```

* * * * * * 


## Introduction: le jeu de Den Boer 2009


Comme cas d'étude, nous analyserons un jeu de données de Den Boer et collègues [@DenBoer:2009ik].

Les auteurs ont caractérisé le transcriptome de 190 patients souffrand de leucémie lymphoblastique aigue (*Acute Lymphoblastic Leukaemia*, *ALL*).

Cette leucémie est caractérisée par une prolifération clonale anormale, dans la moëlle osseuse, de progéniteurs de lymphocytes qui restent bloqués dans  une phase précise de leur différenciation. La prolifération  peut résulter de différentes mutations, dont certaines sont récurrentes chez différents patients.  Le jeu de données de Den Boer inclut l'information concernant les problèmes génétiques particuliers associés à chaque échantillon d'un patient atteint de leucémie. 

Dans leur article original, Den Boer et collègues ont défini une ***signature transcriptomique***, c'est-à-dire une sélection de gènes dont les profils de transcription sont spécifiquement affectés dans une ou plusieurs sous-classes des leucémies. Les profils d'expression de ces gènes-signature peuvent ensuite être utilisés pour classer de nouveaux échantillons de patients atteints de leucémie, et prédire le type de perturbation génétique associée à leur cas. 

Un avantage de ce jeu de données est qu'il inclut des groupes d'effectifs relativement élevés (4 groupes avec >30 patients), par rapport à d'autres publications dans le domaine de l'analyse transcriptomique. 

Le tableau de donnée sur lequel nous travaillerons comporte 190 colonnes (une par patient) et ~22.000 lignes (une par gène) et les valeurs indiquent le niveau d'expression du gène chez le patient considéré. 

On dispose également d'une table "pheno" qui fournit des informations sur chaque échantillon, et en particulier la classe de leucémie dont souffre ce patien (en fait cette classe est déterminée par son génotype plutôt que son phénotype, mais pour des raisons historiques les tableaux annotations des échantillons sont communément dénommés "pheno" dans le domaine des biopuces). 


## But du TP

Le but de ce TP sera de mettre en place une procédure qui combinera les tâches suivants: 

1. Télécharger les fichiers de données à partir d'un site Web.
2. Mesurer une série de paramètres statistiques sur les gènes (lignes) et échantillons (colonnes). 
3. Détecter les gènes différentiellement exprimés entre deux groupes d'échantillons (correspondant à deux classes de leucémie). 


Nous essaierons pour chaque étape 
- de modulariser le code en écrivant des fonctions ré-utilisables
- d'utiliser des techniques efficaces en évitant les boucles, et en utilisant plutôt la fonction `apply()`


## Exercices

### Créaction d'un dossier local pour les données et résultats

1. Ecrivez un bloc de code **R** qui définit (dans une variable dénommé `dir.output`) le chemin d'un dossier où vous sauvegarderez les données et les résultats de ce TP. Ce dossier sera défini comme un sous-dossier `STAT/TP/DenBoer2009` à partir de la racine de votre compte. Créez ce dossier (sauf s'il existe déjà) et vérifiez son contenu en fichiers.  

    Fonctions utiles: 
    
    - `Sys.getenv("HOME")` indiquera la racine de votre compte (Unix et Mac OS X)
    - `file.path()` permet de construire un chemin en concaténant des chaînes de caractères
    - `dir.create()` pour créer un dossier. Lisez l'aide pour choisir les options appropriées. 
    - `list.files()` pour obtenir la liste des fichiers dans un dossier donné

```{r denboer_directory, echo=TRUE, results=FALSE}
## Define the directory for the practical
dir.home <- Sys.getenv("HOME")
dir.output <- file.path(dir.home, "STAT", "TP", "DenBOer2009")
message(dir.output)

## Create the directory
dir.create(dir.output, recursive = TRUE, showWarnings = FALSE)

## Change directory to dir.output
setwd(dir.output)

## Print a message with output directory
message("Output directory for this practical: ", dir.output)

## Check it this directory already contains some files.
list.files(dir.output)
```

Si c'est la première fois que vous réalisez ce TP, la réponse devrait être

`character(0)`

Ceci indique que le dossier que vous venez de créer ne contient aucun fichier. Le résultat de `list.files()` est un vecteur contenant 0 chaînes de caractères.


### Téléchargement des données

Créez un sous-dossier `data` dans votre dossier TP, et téléchargez-y les deux fichiers de données ci-dessous, sauf si ces fichiers existent déjà à cet endroit. 

URL du dossier des données de Den Boer:

-  <http://pedagogix-tagc.univ-mrs.fr/courses/ASG1/data/marrays/>


Dans ce dossier nous utiliserons les deux tables suivantes. 

- Profils transcriptomiques: [GSE13425_Norm_Whole.tab](http://pedagogix-tagc.univ-mrs.fr/courses/ASG1/data/marrays/GSE13425_Norm_Whole.tab)

- Description de chaque échantillon: [phenoData_GSE13425.tab](http://pedagogix-tagc.univ-mrs.fr/courses/ASG1/data/marrays/phenoData_GSE13425.tab)


Vous pouvez éventuellement utiliser le fichier de  Description des sous-types de cancers (abréviation, couleurs d'affichage): 

- [GSE13425_group_descriptions.tab	1](http://pedagogix-tagc.univ-mrs.fr/courses/ASG1/data/marrays/GSE13425_group_descriptions.tab	1)


Le dossier contient également un fichier de qualité des mesures (absent / marginal / présent), que nous ignorerons pour ce TP-ci: 

- [GSE13425_AMP_Whole.tab](http://pedagogix-tagc.univ-mrs.fr/courses/ASG1/data/marrays/GSE13425_AMP_Whole.tab)

Vous trouverez ci-dessous le code pour télécharger une copie des fichiers utiles dans votre dossier personnel. 


Le bloc de code ci-dessous définit l'URL du dossier de données, à partir duquel nous allons télécharger la table de transcriptome + la description des échantillons (pheno). 

```{r data_url}
## Define the URL of the example data sets
url.course <- "http://pedagogix-tagc.univ-mrs.fr/courses/ASG1"
microarray.data.url <- file.path(url.course, "data", "marrays")

## Check the results
print(paste("Course URL: ", url.course))
print(paste("URL for the microarray data folder: ", microarray.data.url))

```

Nous allons maintenant charger la matrice d'expression, fournie dans le fichier `GSE13425_Norm_Whole.txt` dans le dossier Den Boer du serveu. Ce fichier contient un tableau avec 1 ligne par gène^[Plus précisément, chaque ligne correspond à un "probeset", qui mesure la concentration d'une série d'oligonucléotides représentatifs d'un gène. Certains gènes sont représentés par plusieurs probesets, mais par souci de simplicité nous utiliserons le mot "gène" pour les mesures ci-dessous. ] et 1 colonne par échantillon. Chaque échantillon provient d'un patient souffrant de leucémie lymphoblatique aiguë (ALL). 

We will once use the following functions: 

- *file.path()* to build the full URL to the expression table
- *download.file()* to download the file from the URL and obtain a local copy. This will avoid us to download the whole file each time we restart the analysis. 
- *read.table()* to read the table from this URL. Note that this functions allows you to load a table either from a local folder (on your computer), or from any Web site.

```{r load_expression_matrix}

## Get some help about the read.table fonction
# help(read.table)

## URL of the expression table on the course server
expr.url <- file.path(microarray.data.url, "GSE13425_Norm_Whole.txt") ## Build the full URL of the matrix

## Load expression values
dir.data <- file.path(dir.output, "data")
dir.create(dir.data, recursive = TRUE, showWarnings = FALSE)
expr.file <- file.path(dir.data, "DenBoer_2009_GSE13425_Norm_Whole.txt")

if (file.exists(expr.file)) {
  message("Expression file already downloaded: ", expr.file)
} else {
  message("Downloading matrix table from ", expr.url)
  download.file(url = expr.url, destfile = expr.file)
  message("Local expression file: ", expr.file)
}

## Load the expression table in memory
message("Loading expression table from ", expr.file)
expr.matrix <-  read.table(expr.file,sep="\t", header = T, row.names = 1)

```




#### Exercices

1. Adaptez le code ci-dessus pour télécharger une copie locale du fichier pheno. 

```{r load_pheno}

## Get some help about the read.table fonction
# help(read.table)

## URL of the pheno table on the course server
pheno.url <- file.path(microarray.data.url, "phenoData_GSE13425.tab") ## Build the full URL of the matrix

## Load pheno values
pheno.file <- file.path(dir.data, "phenoData_GSE13425.tab")

if (file.exists(pheno.file)) {
  message("pheno file already downloaded: ", pheno.file)
} else {
  message("Downloading matrix table from ", pheno.url)
  download.file(url = pheno.url, destfile = pheno.file)
  message("Local pheno file: ", pheno.file)
}

## Load the pheno table in memory
message("Loading pheno table from ", pheno.file)
pheno.matrix <-  read.table(pheno.file,sep="\t", header = T, row.names = 1)

# View(pheno.matrix)


## Check that the order ofthe pheno table rows is identifcal to the order of the columns in the expression table. For this, we count the number of differences
sum(colnames(expr.matrix) != rownames(pheno.matrix))

```


2. Chargez les fichiers d'expression et de description des échantillons dans des objets de type `data.frame()`, et récupérez dans un vecteur les sous-types de cancer associés à chaque échantillon. 

3. Comptez le nombre d'échantillons par sous-classe de cancer.


```{r samples_per_class}

# View(pheno.matrix)
sample.labels <- as.vector(pheno.matrix$sample.labels)

class.sizes <- sort(table(sample.labels), decreasing = TRUE)

```


4. Sélectionnez deux groupes comportant au moins 30 échantillons chacun, et créez une table d'expression et une table pheno limitées à ces deux groupes. A partir de ce moment vous travaillerez systématiquement avec ces deux tables. 

```{r group_selection}
##  Select the  labels of the two most populated classes
group1 <- names(class.sizes[1])
group2 <- names(class.sizes[2])

## Select the columns of the expression matrix for each group

group1.expr <- expr.matrix[, sample.labels == group1]
dim(group1.expr)
group2.expr <- expr.matrix[, sample.labels == group2]
dim(group2.expr)

selected.samples <- (sample.labels == group1) | (sample.labels == group2)
group.labels <- sample.labels[selected.samples]

```


5. En utilisant la fonction `apply()`, calculez pour chaque gène les paramètres suivants, et stockez-les dans une `data.frame()`.  

    - $m_1, m_2$ moyennes d'échantillons des deux groupes;
    - $s_2.est, s_2.est$ estimations des écarts-types des populations dont sont extraits les échantillons du 1er et 2d groupe, respectivement;
    - $d$ différence entre les moyennes des deux groupes.


```{r stats_per_gene}

gene.stats <- data.frame(
  m1 = apply(group1.expr, 1, mean),
  m2 = apply(group2.expr, 1, mean),
  s1.pop = apply(group1.expr, 1, sd),
  s2.pop = apply(group2.expr, 1, sd)
)
gene.stats$d <- gene.stats$m2 - gene.stats$m1


# View(gene.stats)

hist(gene.stats$d, breaks=100)

```

6. Sélectionnez un gène au hasard, récupérez dans un vecteur la ligne correspondante dans la table d'expression, et effectuez un test de comparaison de moyennes entre les deux groupes.

```{r t_test_one_gene}
## Number of genes
my.gene <- sample(x = row.names(expr.matrix), size = 1)

welch.result <- t.test(x = group1.expr[my.gene, ], y = group2.expr[my.gene,], alternative = "two.sided", var.equal = FALSE)

## Get all the attributes of the t.test() result
attributes(welch.result)

result <- c(welch.result$statistic, 
            p = welch.result$p.value, 
            welch.result$parameter,
            welch.result$estimate)
names(result) <- c("t", "p", "df", "m1", "m2")
result["d"]  <- result["m2"] - result["m1"] 

```


7. Ecrivez une fonction R dénommée `t.test.wrapper` qui prend en entrée 2 vecteurs (expression et classes de cancer, respectivement), effectue un test de Student, et retourene un vecteur avec les statistiques pertinentes. 

```{r t_test_vector}
TTestVector <- function(values, group.labels, group1 = group.labels[1], ...) {
  
  ## Select samples of the first and second group, respectively
  group1.values <- values[group.labels == group1]
  group2.values <- values[group.labels != group1]
  t.result <- t.test(x = group1.values,  y = group2.values, ...)
  
  ## Collect the relevant statistics in the t.test result
  result <- c(t.result$statistic, 
              p = t.result$p.value, 
              t.result$parameter,
              t.result$estimate)
  names(result) <- c("t", "p", "df", "m1", "m2")
  result["d"]  <- result["m2"] - result["m1"] 
  result["s1.est"] <- sd(group1.values)
  result["s2.est"] <- sd(group2.values)
  result["s.diff"] <- result["s2.est"] - result["s1.est"]
  result["s.diff.percent"] <- 100*result["s.diff"] /max(result["s2.est"], result["s1.est"])
  return(result)
}

## Run Welch test with my nice function
my.gene <- sample(x = row.names(expr.matrix), size = 1)
welch.result <- TTestVector(values = expr.matrix[my.gene,selected.samples], 
                            group.labels = sample.labels[selected.samples], 
                            alternative = "two.sided", var.equal = FALSE)


## Just for fun, run a Student test with the same gene
student.result <- TTestVector(values = expr.matrix[my.gene,selected.samples], 
                            group.labels = sample.labels[selected.samples], 
                            alternative = "two.sided", var.equal = TRUE)


## Compare the results obtained with Student and Welch t.test()
kable(data.frame(
  Welch= welch.result,
  Student = student.result), caption="Comparison between the results of a Welch and Student test on a randomly picked up gene. ")

```


8. Ecrivez une boucle `for`qui applique ce même test à chaque ligne de la table d'expression, et récupère dans un tableau les résultats du test, avec une ligne par gène, une colonne par statistique pertinente produite par t.test(). Mesurez le temps d'exécution des 22.000 tests de Student avec cette boucle. 

```{r welch_per_gene_loop}
## Prepare a data frame to host all the results, with one row per gene 
## and one column per field in the result of TTestVector 
## (we use the previous result to adapt the number of columns.
welch.loop <- data.frame(matrix(nrow = nrow(expr.matrix), ncol=length(welch.result)))
names(welch.loop) <- names(welch.result) ## Set column names from previous result
row.names(welch.loop) <- row.names(expr.matrix) ## Set row names as gene names
head(welch.loop) ## Check that the data frame has the right shape and contains NA values

## Run a
loop.time <- system.time(
  for (my.gene in rownames(group1.expr)) {
    welch.loop[my.gene,]  <- TTestVector(
      values=expr.matrix[my.gene, selected.samples], 
      group.labels = sample.labels[selected.samples], 
      alternative = "two.sided", var.equal = FALSE)
  }
)

print(loop.time)

```



9. Utilisez cette fonction et la fonction `apply` pour appliquer le test de Student à chaque gène, et mesurz le temps d'exécution. 



```{r welch_per_gene_apply}
## Run all the t tests with the apply function
apply.time <- system.time(
  welch.apply <- as.data.frame(t(apply(X = expr.matrix[, selected.samples], 1, TTestVector, group.labels = sample.labels[selected.samples], group1 = group1, alternative = "two.sided", var.equal = FALSE)))
)

print(apply.time)

kable(head(welch.apply), caption = "First rows of the table resulting from the application of Welch test to each gene. ")


## Nice formatting for the knit report (a bit tricky, because proc_time)
time.frame <- data.frame(
  loop = as.factor(loop.time),
  apply = as.factor(apply.time)
  )
time.frame$loop <- as.numeric(as.vector(as.factor(loop.time)))
time.frame$apply <- as.numeric(as.vector(as.factor(apply.time)))

kable(time.frame, caption = paste("Time spent for", nrow(group1.expr), "Welch tests with apply:", apply.time), digits=3)

```

**Exercice supplémentaire optionnel:** à titre de curiosité, nous pouvons également appliquer le test de Student au même jeu de données (en combinant `apply()` et `TTestVector()`) et comparer (sur un graphe aux axes logarithmiques) les p-valeurs obtenues par Welch et Student, respectivement. 



```{r student_per_gene_apply}
## Run all the t tests with the apply function
student.apply.time <- system.time(
  student.apply <- as.data.frame(t(apply(X = expr.matrix[, selected.samples], 1, TTestVector, group.labels = sample.labels[selected.samples], group1 = group1, alternative = "two.sided", var.equal = TRUE)))
)

kable(head(student.apply), caption = "First rows of the table resulting from the application of Student test to each gene. ")

print(apply.time)


plot(welch.apply[, "p"], student.apply[, "p"], log="xy")


```
10. Dessinez un volcano plot avec les résultats. 

    Le volcano plot indique en abscisse la "taille d'effet" (dans notre cas, la différence de moyennes), et l'ordonnée  $-log_{10}(p)$, où $p$ est la p-valeur.
    
```{r volcano_plot}
## Define the threshold on p-value
alpha <- 0.05 ## Adopté à l'unanimité par les étudiants

## Draw a XY plot with the difference between means (X) versus the p-value
plot(welch.apply$d, welch.apply$p, xlab="m2 - m1", ylab="p", col="grey")
abline(v=0)
abline(h=alpha, col="red")

## Draw a volcano plot with Welch results
plot(welch.apply$d, -log10(welch.apply$p), xlab="m2 - m1", ylab="-log10(p)", col="grey", main="Volcano plot of the nominal p-values")
abline(v=0)
abline(h=-log10(alpha), col="red")

```

    
11. Dessinez un histogramme des p-valeurs. 

```{r p_value_histogram}
## Histogram of p-values for Den Boer 2009 expression matrix
hist(welch.apply$p, breaks=seq(from=0, to=1, by=0.05), main="P-value histogram", xlab="P-value", ylab="Nb genes", las=1, col="#BBDDFF")
```


**Exercices supplémentaires**: 

    a. Générez une matrice de données aléatoires de la même taille que la matrice d'expression, où chaque ligne comporte des données normales sous hypothèse nulle ($H_0: \mu_1 = \mu_2$). Effectuez le test de Student et dessinez le volcano plot et l'histogramme des p-valeurs. A titre d'exemple, choisissons $\mu_1 = \mu_2 = 0$, $\sigma_1 = \sigma_2 = 1$. 
    b. Générez de la même façon une autre matrice de données sous hypothèse alternative, avec $\mu_1 = 0$, $\mu_2 = 2$, $\sigma_1 = \sigma_2 = 1$, effectuez un test de Student et dessinez le volcano plot et l'histogramme des p-valeurs. 
    c. Sélectionnez 15000 lignes de la première matrice (sous $H_0$), 5000 lignes de la seconde (sous $H_A$), effectuez les tests de Student et dessinez le volcano plot et l'histogramme des p-valeurs.
    
```{r random_data}
## Get the dimensions of the expresion matrix
g <- nrow(expr.matrix) ## number of genes
s <- ncol(expr.matrix) ## Number of samples

## Generate a matrix with random data under H0
rand.h0 <- matrix(nrow=g, ncol=s, data = rnorm(n = g * s))

student.under.h0 <- as.data.frame(t(apply(X = rand.h0[, selected.samples], 1, TTestVector, group.labels = sample.labels[selected.samples], group1 = group1, alternative = "two.sided", var.equal = TRUE)))

## Plot the histogram of p-values for random normal data under H0
hist(student.under.h0$p, breaks=seq(from=0, to=1, by=0.05), main="Random numbers under H0", xlab="P-value", ylab="Nb genes", las=1, col="#BBDDFF")


## Draw a volcano plot with random normal data under H0
plot(student.under.h0$d, -log10(student.under.h0$p), xlab="m2 - m1", ylab="-log10(p)", col="grey", main="Random normal under H0")
abline(v=0)
abline(h=-log10(alpha), col="red")

```



12. Dessinez un histogramme (en densités) avec les valeurs de la statistique du test ($t$). 

```{r t_statistic_histogram}
hist(welch.apply$t, breaks=100)
```

13. **Corrections de tests multiples**. Ajoutez au tableau de résultats les colonnes suivantes.

    - **e-valeur** ($E = p \cdot N$, où $N$ est le nombre de test)
    - **p-valeur corrigée par la méthode de Bonferonni**, calculé avec la fonction `p.adjust()` (tond vous lirez préalablement l'aide).
    - **taux de de fausse découverte**  (***False Discovery Rate***, ***FDR***), calculé avec la fonction `p.adjust()`.

```{r multiple_test_corrections}
## E-value
g <- nrow(expr.matrix)
welch.apply$E <- welch.apply$p * g
welch.apply$Bonferroni <- p.adjust(welch.apply$p, method = "bonferroni")
welch.apply$BH <- p.adjust(welch.apply$p, method = "BH")
welch.apply$fdr <- p.adjust(welch.apply$p, method = "fdr")
plot(welch.apply$p, welch.apply$Bonferroni, log="xy", col="red")
points(welch.apply$p, welch.apply$BH, col="darkgreen")
points(welch.apply$p, welch.apply$fdr, col="blue", pch=20, cex=0.8)
grid()

## Negative control: all data under H0
student.under.h0$E <- student.under.h0$p * g
student.under.h0$Bonferroni <- p.adjust(student.under.h0$p, method = "bonferroni")
student.under.h0$BH <- p.adjust(student.under.h0$p, method = "BH")
student.under.h0$fdr <- p.adjust(student.under.h0$p, method = "fdr")



## Draw a volcano plot with random normal data under H0
plot(student.under.h0$d, -log10(student.under.h0$E), xlab="m2 - m1", ylab="-log10(E)", col="grey", main="Random normal under H0")
abline(v=0)
abline(h=-log10(alpha), col="red")

## Histogram of E-values for Den Boer
hist(welch.apply$E, breaks=g *seq(from=0, to=1, by=0.05), main="E-value histogram", xlab="E-value", ylab="Nb genes", las=1, col="#BBDDFF")

## Plot the histogram of E-values for random normal data under H0
hist(student.under.h0$E, breaks=g*seq(from=0, to=1, by=0.05), main="Random numbers under H0", xlab="E-value", ylab="Nb genes", las=1, col="#BBDDFF")
par(mfrow=c(1,1))


```


14. Dessinez un graphe qui compare la p-valeur nominale (en abscisse) et les différentes corrections (en ordonnée). Utilisez les axes logarithmiques pour faire apparaître les valeurs pertinentes (très petites). 

```{r}

## Plot the volcano and p-value histograms of E-values
par(mfrow=c(2,2))
## Draw a volcano plot with Welch results
plot(welch.apply$d, -log10(welch.apply$E), xlab="m2 - m1", ylab="-log10(E)", col="grey", main="Den Boer, E-value volcano")
abline(v=0)
abline(h=-log10(alpha), col="red")## Histogram of p-values for Den Boer 2009 expression matrix

```

15. Sélectionez les gènes significatifs avec un seuil $\alpha=0.05$ **sur le taux de fausse découverte**. 


```{r}
selected.genes <- welch.apply$fdr < alpha
sum(selected.genes)
```

16. Sélectionnez une matrice restreinte aux profils d'expression des gènes différentiellement exprimés pour les deux groupes d'intérêt (les deux groupes dans la même matrice). 

```{r}
selected.expr.matrix <- expr.matrix[selected.genes, selected.samples]
dim(selected.expr.matrix)
```


***Ici s'arrête la matière à connaître pour l'examen: le reste du tuto est fourni pour le plaisir. ***


17. Calculez -- avec la fonction `cor()` -- une matrice de corrélation entre échantillons de cette matrice de gènes différentiellement exprimés. Générez un diagramme de cette matrice avec `heatmap()`. 

```{r sample_correlations}
## Compute Pearson correlation between samples
dim(selected.expr.matrix)

## Compute inter-sample correlation
sample.cor <- cor(selected.expr.matrix)
dim(sample.cor)
heatmap(sample.cor, scale = "none", col=bluered(100), symm=TRUE)

```




18. Utilisez la fonction `hclust()` pour effectuer un clustering des échantillons, sur base de cette matrice de corrélation. 

```{r gene_correlation_tree}

## Compute inter-gene correlation
gene.cor <- 1 - cor(t(selected.expr.matrix))
dim(gene.cor)

## Build a gene tree by running hierarchical clustering between genes
gene.tree <- hclust(d = as.dist(gene.cor), method = "complete", members = NULL)
plot(gene.tree)

```


19. Avec la fonction `heatmap()`, dessinez une carte thermique (***heat map***) des gènes différentiellement exprimés, en activant le clustering sur les colonnes (samples) et lignes (gènes). Raffinez ensuite le graphique avec la fonction `heatmap2()` de la librairie `gplots` (installez cette librairie si nécessaire).

```{r expression_heatmap}

## Gene-wise centering for the heatmap. The gene-wise (row-wise) mean is substracted fromeach expression value
selected.expr.matrix.gene.centered <- selected.expr.matrix - apply(selected.expr.matrix,1, mean)


# library(corrplot)
# corrplot(sample.cor[1:10, 1:10])
x <- 
selected.expr.matrix.gene.centered

## Add class label to each gene name
colnames(x) <- paste(colnames(x), sample.labels[selected.samples])

## Compute maximal absolute value of x for the range of the heatmap
absmax <- ceiling(max(abs(range(x))))

## Use 99th percentile as limit for the color scale
zlim <- quantile(x = abs(unlist(x)), probs = 0.99)

## Draw a heatmap with one row per gene and one column per sample
heatmap(as.matrix(x), cexCol = 0.5, col=redgreen(100), scale="none", labRow = FALSE, symm=FALSE, margins=c(5,1), zlim=c(-zlim, zlim), Rowv = as.dendrogram(gene.tree))

```


20. Au moyen de la fonction `strsplit()`, extrayez les noms des gènes significatifs à partir des étiquettes de lignes de la matrice d'expression. Ces étiquettes contiennent 2 identifiants séparés par un caractère $|$ (tube): nom de gène et identifiant de la sonde sur la biopuce. 

```{r}
gene.rownames <- rownames(expr.matrix)

gene.desc <- t(data.frame(strsplit(x = gene.rownames, split = '|', fixed=TRUE)))
rownames(gene.desc) <- rownames(expr.matrix)
colnames(gene.desc) <- c("gene.name", "affy.ID")
# dim(gene.desc)
kable(head(gene.desc))

## Selected gene names
selected.gene.names <- as.vector(gene.desc[selected.genes, "gene.name"])

selected.affy.IDs <- as.vector(gene.desc[selected.genes, "affy.ID"])
```


21. Analysez l'enrichissement fonctionnel de ces gènes avec les deux outils suivants: 

    - gProfiler: <http://biit.cs.ut.ee/gprofiler/>
    - DAVID: <https://david.ncifcrf.gov/>
    
22. ACP (PCA), sparse PCA

## References


