944 lines
40 KiB
Text
944 lines
40 KiB
Text
---
|
||
title: "Projet d'Analyse de données"
|
||
author: "JIBIDAR Blessing, MOUGNIBAS Théo"
|
||
output:
|
||
pdf_document :
|
||
toc : TRUE
|
||
toc_depth : 2
|
||
number_section : TRUE
|
||
header-includes:
|
||
- \usepackage{dsfont}
|
||
- \usepackage{color}
|
||
- \newcommand{\1}{\mathds{1}}
|
||
---
|
||
date: "2024-12-04"
|
||
|
||
---
|
||
|
||
```{r setup, include=FALSE}
|
||
knitr::opts_chunk$set(echo = FALSE)
|
||
library(ggplot2)
|
||
library(gridExtra)
|
||
library(grid)
|
||
library(reshape2)
|
||
library(corrplot)
|
||
library(FactoMineR)
|
||
library(factoextra)
|
||
library(cluster)
|
||
library(mclust)
|
||
library(dplyr)
|
||
library(tidyr)
|
||
library(dbscan)
|
||
library(klaR)
|
||
library(clusterSim)
|
||
library(circlize)
|
||
```
|
||
|
||
```{r,echo=FALSE,results = FALSE,message=FALSE,warning=FALSE,fig.show='hide'}
|
||
set.seed(1234)
|
||
# chargement de la table, /!\ séparateur point virgule car on dirait que les données sont dans un genre de csv
|
||
T = read.table("DataProjet3MIC-2425.txt",header=TRUE,sep=";")
|
||
T$ExpT1 = as.factor(T$ExpT1)
|
||
T$ExpT2 = as.factor(T$ExpT2)
|
||
T$ExpT3 = as.factor(T$ExpT3)
|
||
# petit regard sur les données
|
||
head(T)
|
||
summary(T)
|
||
|
||
str(T)
|
||
levels(T$ExpT1)
|
||
```
|
||
|
||
|
||
|
||
# Contenu du jeu de données :
|
||
Le jeu de données contient 3 variables qualitatives nominales représentant l'expression du gêne g dont les modalités sont {"sur","sous","non"}. Chaque variable correspond respectivement à la différence d'expression du gêne mesurée à la 6ème heure lors du traitement T appartenant à {T1,T2,T3} en moyenne, sur les réplicats {R1,R2}.
|
||
|
||
3x6 + 3x6 = 36 variables quantitatives continues représentant les effets des traitements T1 T2 et T3 sur l'expression des gênes à 1h,2h,3h,4h,5h,6h après l'administration pour les replicats {R1,R2}, par rapport à leur expression à T=0 ( sans traitement). Ce jeu de données contient des relevés sur 542 individus, ici des gênes.
|
||
|
||
|
||
# Analyse unidimensionnelle :
|
||
|
||
## Expression des gênes lors des traitements T1,T2 et T3
|
||
|
||
Nous avons réalisé 3 histogrammes pour visualiser les effectifs des gènes en fonction de leur expression relative moyenne à 6h suite aux traitements T1,T2 et T3, cette moyenne est représentée par les variables qualitatives nominales ExpT1,ExpT2 et ExpT3
|
||
|
||
```{r,fig.height=3,warning=FALSE,fig.cap="\\label{fig:histo1}Visualualisation des expressions relative des gènes lors du traitement T1,T2 et T3"}
|
||
g1<-ggplot(T, aes(x=T$ExpT1))+
|
||
geom_bar()+
|
||
ylab("Effectifs")
|
||
|
||
g2<-ggplot(T, aes(x=T$ExpT2))+
|
||
geom_bar()+
|
||
ylab("Effectifs")
|
||
|
||
g3<-ggplot(T, aes(x=T$ExpT3))+
|
||
geom_bar()+
|
||
ylab("Effectifs")
|
||
|
||
grid.arrange(g1,g2,g3,ncol=3)
|
||
```
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
|
||
### Analyse des résultats
|
||
Sur la Figure \ref{fig:histo1} on remarque que les traitements T2 et T3 semblent avoir un effet assez similaire sur l'expression des gênes relevée à la 6ème heure : Une polarisation entre la sous expression et la sur expression qui se partagent presque la totalité des relevés, avec un poids légèrement superieur à 55% pour la sur-expression au détriment de la sous-expression. Cela a peut être un rapport avec le fait que T3 est une combinaison des traitement T1 et T2.
|
||
|
||
T1 quant à lui se démarque grandement par une large majorité (Un peu plus de 80%), de gêne n'ayant pas changé d'expression après 6h de traitement.
|
||
|
||
## Expression relative des gênes mesurées à intervalle régulier
|
||
```{r, fig.height=5,fig.cap="Histogrammes de l'expression relative des gènes aux différents relevés lors du traitement T1",echo=FALSE,results = FALSE,message=FALSE,warning=FALSE,fig.show='hide'}
|
||
|
||
T_long = melt(T[c(1,2,3,4,5,6,19,20,21,22,23,24)])
|
||
ggplot(T_long, aes(x = value)) +
|
||
geom_histogram(binwidth = 1, fill = "blue", color = "black", alpha = 0.7) +
|
||
facet_wrap(~variable,scales = "free",ncol=6) + # va créer un sous graphique pour chaque variable, sans forcément la même échelle sur les 2 axes
|
||
labs(title = "Histogrammes de l'expression relative des gènes aux différents relevés lors du traitement T1", x = "Valeurs", y = "Effectifs")
|
||
```
|
||
```{r, fig.height=5,fig.cap="Histogrammes de l'expression relative des gènes aux différents relevés lors du traitement T2",echo=FALSE,results = FALSE,message=FALSE,warning=FALSE,fig.show='hide'}
|
||
T_long = melt(T[c(7,8,9,10,11,12,25,26,27,28,29,30)])
|
||
ggplot(T_long, aes(x = value)) +
|
||
geom_histogram(binwidth = 1, fill = "blue", color = "black", alpha = 0.7) +
|
||
facet_wrap(~variable,scales = "free",ncol=6) +
|
||
labs(title = "Histogrammes de l'expression relative des gènes aux différents relevés lors du traitement T2", x = "Valeurs", y = "Effectifs")
|
||
```
|
||
```{r, fig.height=5,fig.cap="Histogrammes de l'expression relative des gènes aux différents relevés lors du traitement T3",echo=FALSE,results = FALSE,message=FALSE,warning=FALSE,fig.show='hide'}
|
||
T_long = melt(T[c(13,14,15,16,17,18,31,32,33,34,35,36)])
|
||
ggplot(T_long, aes(x = value)) +
|
||
geom_histogram(binwidth = 1, fill = "blue", color = "black", alpha = 0.7) +
|
||
facet_wrap(~variable,scales = "free",ncol=6) +
|
||
labs(title = "Histogrammes de l'expression relative des gènes aux différents relevés lors du traitement T3", x = "Valeurs", y = "Effectifs")
|
||
```
|
||
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
Lors d'une analyse supplémentaire des fréquences d'expression des gènes à chaque heure et chaque traitement, nous observons bien une concordance avec l'analyse des expressions des gênes Figure \ref{fig:histo1}. En effet, les histogrammes en rapport avec le traitement 1 sont très nettement regroupés vers 0, soit une expression relative des gênes qui ne change peu. Les histogrammes pour les relevés des variables en lien avec T2 et T3 sont tout aussi similaires aux résultats précédents : La variance de l'expression relative des gênes est plus élevée et on observe bien une polarisation "sous-exprimé-"sur-exprimé" sur les relevés à 6h. Attention, ici on observe aussi que T2 et T3 n'ont pas leur effet caractéristique directement : à 2h, la distribution de l'expression des genes semble presque Gaussienne, et à 1h elle ne se distingue pas beaucoup du traitement 1 avec un regroupement sur 0.
|
||
|
||
c.f le document RMarkdown pour les histogrammes.
|
||
|
||
## Boxplots des différentes variables quantitatives du jeu de données
|
||
|
||
```{r,fig.height=4,message=FALSE,warning=FALSE,fig.cap="\\label{fig:boxplot1}Boxplots des différentes variables du jeu de données lors du réplicat R1"}
|
||
ggplot(melt(T[1:18]),aes(x=variable,y=value))+ # multiplie chaque gène par chaque variable, on va donc avoir plusieurs lignes avec le même gène mais une autre colonne va donner sa valeur sur une variable précise
|
||
geom_boxplot()+ theme(axis.text.x = element_text(angle=90, vjust = 0.5, hjust = 1))
|
||
```
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
Sur la figure \ref{fig:boxplot1} nous avons choisi de n'inclure que le graphe portant sur le réplicat R1 car il est porteur des mêmes informations que celui sur le réplicat R2.
|
||
L'analyse des boxplots montre que, même sans effectuer de réduction des données, chaque variable présente une dispersion similaire, ce qui suggère qu'il n'est pas nécessaire de procéder à un centrage et à une réduction avant de réaliser l'Analyse en Composantes Principales.
|
||
|
||
# Analyse bi-dimensionnelle
|
||
|
||
## Boxplots par paire de variables (qualitative,quantitative)
|
||
|
||
Nous avons ensuite générés des graphes contenants des boxplots d'expression des gènes moyennes d'un traitement par rapport à l'expression des gènes à 6h sur un autre traitement.
|
||
Cela nous a permis d'identifier les changements d'expressions de chaque groupes de gènes sur, sous ou non-exprimés entre un traitement et les deux autres.
|
||
Les deux réplicats ayant des résultats très similaires, pour ne pas surcharger le rapport nous avons décidé d'afficher seulement le réplicat R1.
|
||
|
||
```{r,fig.cap="\\label{fig:boxplot2}Boxplots par paire de variables (qualitative,quantitative)" ,warning=FALSE}
|
||
|
||
# traitement 1 corrélation avec l'expression des genes du T1 T2 et T3
|
||
|
||
p1 = ggplot(T,aes(y=T$T1_6H_R1,x=T$ExpT2))+
|
||
geom_boxplot()
|
||
#ggplot(T,aes(y=T$T1_6H_R2,x=T$ExpT2))+
|
||
# geom_boxplot()
|
||
p2 = ggplot(T,aes(y=T$T1_6H_R1,x=T$ExpT3))+
|
||
geom_boxplot()
|
||
#ggplot(T,aes(y=T$T1_6H_R2,x=T$ExpT3))+
|
||
# geom_boxplot()
|
||
|
||
# traitement 2 corrélation avec l'expression des genes du T1 T2 et T3
|
||
|
||
p3 = ggplot(T,aes(y=T$T2_6H_R1,x=T$ExpT1))+
|
||
geom_boxplot()
|
||
#ggplot(T,aes(y=T$T2_6H_R2,x=T$ExpT1))+
|
||
# geom_boxplot()
|
||
p4 = ggplot(T,aes(y=T$T2_6H_R1,x=T$ExpT3))+
|
||
geom_boxplot()
|
||
#ggplot(T,aes(y=T$T2_6H_R2,x=T$ExpT3))+
|
||
# geom_boxplot()
|
||
|
||
# traitement 3 corrélation avec l'expression des genes du T1 T2 et T3
|
||
p5 = ggplot(T,aes(y=T$T3_6H_R1,x=T$ExpT1))+
|
||
geom_boxplot()
|
||
#ggplot(T,aes(y=T$T3_6H_R2,x=T$ExpT1))+
|
||
# geom_boxplot()
|
||
p6 = ggplot(T,aes(y=T$T3_6H_R1,x=T$ExpT2))+
|
||
geom_boxplot()
|
||
#ggplot(T,aes(y=T$T3_6H_R2,x=T$ExpT2))+
|
||
# geom_boxplot()
|
||
# Utiliser grid.arrange pour afficher les graphiques dans une grille
|
||
grid.arrange(p1, p3, p6, p2, p4, p5, ncol = 3)
|
||
```
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
### Analyse des boxplots :
|
||
Nous allons analyser les boxplots obtenus à la Figure \ref{fig:boxplot2}
|
||
|
||
|
||
#### Traitement 1
|
||
Les genes sur-exprimés au T1 sont non-exprimé durant le T2 et T3 (Boxplots T1/ExpT2 et T1/ExpT3) .
|
||
Il est difficile d'observer une catégorie de genes de T1 qui se soient sous-exprimés ou sur-exprimés dans T2 et T3, ceux qui n'avaient pas changés d'expression relative durant T1, se sont soit sous-exprimé soit sur-exprimé (Boxplots T1/ExpT2 et T1/ExpT3).
|
||
|
||
#### Traitements 2 et 3
|
||
On observe une légère tendance des genes s'étant sous-exprimé avec T1 à se sous-exprimer avec T2 et T3 (Boxplots T2/ExpT1 et T3/expT1). En revanche, il est très clair que T2 et T3 ont les même effets sur les mêmes genes à 6h, toutes les expressions relevées par T2 concordent aux modalités qualitatives moyennes calculées sur T3 (Boxplots T2/ExpT3 et T3/ExpT2).
|
||
|
||
|
||
## matrice de covariance des variables quantitatives
|
||
Nous avons généré la matrice de covariance de nos données sans les variables quantitatives. Suite à cela nous avons remarqué qu'il était inutile d'afficher le réplicat R2 dans la matrice et que cela rendait le graphe moins lisible. Nous avons donc décidé de seulement afficher le graphe de la matrice de covariance avec le réplicat R1.
|
||
|
||
|
||
```{r, fig.height = 4,fig.cap="\\label{fig:covplot1}Visualisation de la matrice de covariance des variables quantitatives sur le replicat R1"}
|
||
# calcul de la matrice de covariance sur R1
|
||
cr = cor(T[c(1:18)])
|
||
# affichage
|
||
corrplot(cr,method="ellipse", type="lower", bg = "lightgrey")
|
||
```
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
|
||
### analyse de la matrice de covariance
|
||
Sur la Figure \ref{fig:covplot1} on observe clairement de grandes zones de de corrélation entre T2 et T3 délimitées par une correlation plus moyenne à la première heure des deux traitements sûrement dues au fait que T3 est influencé par T1 aux premières heures.
|
||
On remarque que T3 est plus fortement corrélé à T1 à 1H que T2 (qui ne l'est que légèrement), semblant indiquer que T1 a un effet avant T2; T3 étant la combinaison des deux traitements.
|
||
|
||
|
||
## Rapport de correlation Eta2
|
||
```{r,echo=FALSE,results = FALSE,message=FALSE,warning=FALSE,fig.show='hide'}
|
||
library(BioStatR)
|
||
|
||
print("T1 vs T2")
|
||
cat("Replicat 1 : ", eta2(T$T1_6H_R1, T$ExpT2), "\n")
|
||
cat("Replicat 2 : ", eta2(T$T1_6H_R2, T$ExpT2), "\n")
|
||
|
||
print("T1 vs T3")
|
||
cat("Replicat 1 : ", eta2(T$T1_6H_R1, T$ExpT3), "\n")
|
||
cat("Replicat 2 : ", eta2(T$T1_6H_R2, T$ExpT3), "\n")
|
||
|
||
print("T2 vs T3")
|
||
cat("Replicat 1 : ", eta2(T$T2_6H_R1, T$ExpT3), "\n")
|
||
cat("Replicat 2 : ", eta2(T$T2_6H_R2, T$ExpT3), "\n")
|
||
|
||
```
|
||
|
||
+ Pour la comparaison entre T1 et T2, le rapport de corrélation est de 0,187 pour le réplicat 1 et de 0,199 pour le réplicat 2.
|
||
+ Pour la comparaison entre T1 et T3, le rapport de corrélation est de 0,152 pour le réplicat 1 et de 0,142 pour le réplicat 2.
|
||
+ Pour la comparaison entre T2 et T3, le rapport de corrélation est de 0,902 pour le réplicat 1 et de 0,888 pour le réplicat 2.
|
||
Le calcul du rapport de correlation eta² confirme bien notre observation de la grande similarité d'expression des genes traités avec T2 et T3 et la dissimilarité des expression des genes lorsque la plante est traitée avec T1 comparée à T2 et T3, chose normale au vu du peu de genes affectés par T1.
|
||
|
||
|
||
## tables de contingence pour les variables qualitatives, 2 à 2
|
||
```{r,fig.cap="\\label{fig:table1}tables de contingence pour les variables qualitatives, 2 à 2"}
|
||
print("table de contingence entre T1 et T2")
|
||
table(T$ExpT1,T$ExpT2)
|
||
print("table de contingence entre T1 et T3")
|
||
table(T$ExpT1,T$ExpT3)
|
||
print("table de contingence entre T2 et T3")
|
||
table(T$ExpT2,T$ExpT3)
|
||
```
|
||
|
||
|
||
Nouvelle confirmation de nos résultats de manière encore plus précise, on observe que T1 ne change pas l'expression de la très grande majorité des gènes. Plus finement, on peut confirmer l'observation faite sur les boxplots tendant à dire que le peu de gènes s'étant sous exprimés avec T1 se sont aussi sous-exprimés avec T2 et T3.
|
||
|
||
La grande valeur des effectifs partiels sur la diagonale de la table de contingence entre T2 et T3 montre bien la similarité de l'effet de ces deux traitement sur l'expression des genes.
|
||
|
||
|
||
|
||
# Analyse en composantes principales où les Tt sH Rr sont les individus décrits par les gènes.
|
||
|
||
L'ACP permet d'obtenir une vision synthétique des données en réduisant leur complexité tout en préservant l'information essentielle. Elle facilite l'identification des tendances et profils d'individus, ainsi que la construction de méta-variables, plus simples à interpréter, tout en évitant les redondances entre les variables.
|
||
Pour faire cela, nous devons transposer la matrice de données originale qui elle décrivait les gènes (individus) en fonction des Tt sH Rr.
|
||
Nous décidons de faire directement une ACP car, comme mentionné plus tôt lors de l'analyse de la Figure \ref{fig:boxplot1}, nous n'avons pas besoin de centrer et réduire le jeu de données.
|
||
|
||
```{r, fig.width=10, fig.height=5, fig.cap="\\label{fig:vp1}Participation des chaque valeur propre de la matrice de correlation à l'intertie totale des données"}
|
||
donnees_transposees = t(T[-c(37:39)])
|
||
# ACP sans centrer ni réduire, on cache les graphes automatiques
|
||
res_pca<-PCA(donnees_transposees,scale.unit=FALSE,graph=FALSE)
|
||
# on montre d'abord le graphe des vp
|
||
fviz_eig(res_pca)
|
||
```
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
|
||
La Figure \ref{fig:vp1} représente la participation de chaque valeur propre de la matrice de corrélation du jeu de données dans l'inertie totale. L'inertie totale étant la somme des valeurs propres ( qui elles sont les inerties axiale associées à l'axe de vecteur directeur le vecteur propre associé ), chaque valeur propre est donc une fraction de l'inertie totale.
|
||
Nous voyons un dépassement des 80% de l'inertie totale rien qu'avec les deux premieres valeurs propres, nous en prenons donc les vecteurs propres associés comme axes principaux de l'analyse.
|
||
|
||
|
||
```{r,fig.width=8,fig.cap="\\label{fig:acp1}Projection des individus dans le plan factoriel et corrélations des variables avec les composantes principales"}
|
||
plot1=fviz_pca_ind(res_pca,label="all")
|
||
|
||
# on met un habillage en fonction des modalités des variables qualitatives, on voit beaucoup mieux
|
||
plot2=fviz_pca_var(res_pca,axes=c(1,2),label="none",habillage=T$ExpT1, geom=c("point"))
|
||
plot3=fviz_pca_var(res_pca,axes=c(1,2),label="none",habillage=T$ExpT2, geom=c("point"))
|
||
plot4=fviz_pca_var(res_pca,axes=c(1,2),label="none",habillage=T$ExpT3, geom=c("point"))
|
||
grid.arrange(plot1,plot2,plot3,plot4,ncol=2)
|
||
```
|
||
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
Contexte de la figure \ref{fig:acp1} : les relevés aux heures sont décrits par les gènes ( les gènes sont considérés comme les variables).
|
||
|
||
- les genes proches d'un axe sont très représentés par celui-ci
|
||
- les genes dont l'angle entre eux est petit sont corrélés entre eux
|
||
- plus le point représentant un gène est éloigné de l'origine du repère selon un axe, plus ce gène contribue à l'inertie de cet axe
|
||
|
||
## Interprétation globale du couple de graphes Figure \ref{fig:acp1}
|
||
|
||
On voit que les genes se polarisent principalement sur l'axe 1 dans un sens ou l'autre.
|
||
Certaines flèches sont d'une longueur presque du rayon du cercle, indiquant une participation très forte des genes dans la variance expliquée par ces dimensions.
|
||
Dans chaque "polarité" de fleches selon l'axe 1, il y a des fleches dont la direction est negative d'autres positive selon l'axe 2. Bien que l'on dénote une quantité plus grande de gènes corrélés positiviement à la dimension 2.
|
||
|
||
|
||
Gràce à l'habillage des graphes des variables (gènes) à l'aide des variables qualitatives, on remarque que les gènes se polarisent sur la dimension 1 du côté des valeurs negatives les gènes s'étant sous-exprimés à 6h en moyenne sur le traitement T2 et du côté des valeurs positives, ceux qui étaient sur-exprimés.
|
||
Pour la dimension 2, on observe les heures de relevés de traitement s'étaler selon cet axe, on suppose alors que les gènes corrélés à l'axe 2 très positivement s'exprime tôt, et inversement pour ceux qui s'expriment aux environs des 6h.
|
||
|
||
Le traitement 1 est entièrement groupé sur des valeurs très negatives de l'axe 1. On remarque dans ce groupement la présence des T2 et T3 à la première heure de relevés d'expression des gènes.
|
||
|
||
Pour le traitement 2 et 3, on les retrouves formant 2 groupements, 1 en haut à droite du graphe contenant les relevés à 2 et 3h puis un groupement s'étalant sur la droite du graphe du centre jusqu'en bas contenant les relevés à partir de 4h.
|
||
|
||
On en déduit alors que les gènes étant impliqués dans le traitement 1 s'expriment à vitesse moyenne, ni rapidement ni lentement mais qu'ils se sont sous-exprimés dans le traitement T2.
|
||
De plus, on observe que T2 et T3 visent durant la première heure les gènes qui vont se sous-exprimer et se déplacent petit-à-petit sur le graphe pour toucher les gènes qui vont se sur-exprimer.
|
||
|
||
## Clustering
|
||
|
||
L'objectif du clustering est de regrouper des individus (ici, les Tt_sH_Rr) en groupes homogènes selon leurs similarités, sans connaissance préalable des catégories. Cela permet d’identifier des profils ou des comportements communs, facilitant l’interprétation des données et la mise en évidence de structures sous-jacentes.
|
||
|
||
|
||
Nous avons commencé par effectuer un clustering avec la méthode K-means basée sur des estimations de nombre de clusters optimaux avec Silhouette et Calinski-Harabasz ainsi que l'inertie intra-classe, pour finir par effectuer un clustering avec la méthode CAH dont le nombre de classes est déterminé par l'indice de Calinski-Harabasz calculé sur des coupures du dendogramme.
|
||
|
||
Pour le dendogramme, nous avons utilisé la méthode d'agrégation de Ward car il y a trop de points pour utiliser la mesure euclidienne avec efficacité.
|
||
|
||
|
||
### Clustering k-means
|
||
|
||
On commence par afficher l'inertie intra classe en fonction du nombre de classe pour le clustering avec k-means, le but étant de la minimiser tout en gardant un nombre de classes raisonnable.
|
||
|
||
|
||
```{r,fig.width=5,fig.cap="\\label{fig:intra1}Visualisation de l'inertie intra-classe en fonction du nombre de classes demandées pour le clustering"}
|
||
s = donnees_transposees
|
||
|
||
set.seed(1234)
|
||
|
||
# code pour générer le graphe d'inertie intra issu du TP
|
||
Kmax<-15
|
||
reskmeanscl=matrix(0,nrow=nrow(s),ncol=Kmax-1)
|
||
Iintra<-NULL
|
||
for (k in 2:Kmax){
|
||
resaux<-kmeans(s,k)
|
||
reskmeanscl[,k-1]=resaux$cluster
|
||
Iintra<-c(Iintra,resaux$tot.withinss)
|
||
}
|
||
|
||
df<-data.frame(K=2:15,Iintra=Iintra)
|
||
ggplot(df,aes(x=K,y=Iintra))+
|
||
geom_line()+
|
||
geom_point()+
|
||
xlab("Nombre de classes")+
|
||
ylab("Inertie intraclasse")
|
||
```
|
||
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
Sur la Figure \ref{fig:intra1} on dénote un coude dans la courbe d'inertie intraclasse aux alentours de 4 clusters.
|
||
|
||
À présent on va afficher sur la Figure \ref{fig:sil1} l'évolution de l'indicateur Silouhette en fonction du nombre de classe ainsi que l'évolution de l'indice de Calinski-Harabasz en fonction du nombre de classe.
|
||
|
||
|
||
```{r,fig.cap="\\label{fig:sil1}Visualisation du critère de Silhouette et de Calinski-Harabasz en fonction du nombre de classes demandées pour le clustering",fig.width=5}
|
||
set.seed(1234)
|
||
|
||
# code pour générer le graphe Silhouette issu du TP
|
||
|
||
Silhou<-NULL
|
||
for (k in 2:Kmax){
|
||
aux<-silhouette(reskmeanscl[,k-1], daisy(s))
|
||
Silhou<-c(Silhou,mean(aux[,3]))
|
||
}
|
||
|
||
df<-data.frame(K=2:Kmax,Silhouette=Silhou)
|
||
p11=ggplot(df,aes(x=K,y=Silhouette))+
|
||
geom_point()+
|
||
geom_line()+theme(legend.position = "bottom")+
|
||
xlab("Nombre de classes")+
|
||
ylab("Silhouette")
|
||
|
||
CH<-NULL
|
||
for (k in 2:Kmax){
|
||
CH<-c(CH,index.G1(x=s , cl=reskmeanscl[,k-1]))
|
||
}
|
||
daux<-data.frame(NbClust=2:Kmax,CH=CH)
|
||
p12 =ggplot(daux,aes(x=NbClust,y=CH))+
|
||
geom_line()+
|
||
geom_point()+
|
||
xlab("Nombre de classes")+
|
||
ylab("CH")
|
||
|
||
grid.arrange(p11,p12,ncol=2)
|
||
|
||
```
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
Sur la Figure \ref{fig:sil1} Silhouette et l'indice de Calinski-Harabasz ont un pic à 2, mais cela est normal sachant que Silhouette et l'indice de Calinski-Harabasz ont tendance à sous-estimer le nombre de clusters.
|
||
|
||
### visualisation du clustering
|
||
```{r,fig.cap="\\label{fig:clust1}Visualisation des clusters générés par la méthode kmeans dans le plan factoriel 1,2",echo=FALSE,message=FALSE,results='hide'}
|
||
set.seed(1234)
|
||
#graphes des clusters
|
||
res_kmeans = kmeans(s,4)
|
||
p11 = fviz_cluster(res_kmeans,data=donnees_transposees,
|
||
ellipse.type="norm",labelsize=8,
|
||
geom=c("point"))
|
||
|
||
|
||
# on choisis 2 clusters
|
||
aux<-silhouette(reskmeanscl[,4-1], daisy(s))
|
||
p12 = fviz_silhouette(aux)
|
||
|
||
#graphes des clusters
|
||
res_kmeans = kmeans(s,2)
|
||
p21 = fviz_cluster(res_kmeans,data=donnees_transposees,
|
||
ellipse.type="norm",labelsize=8,
|
||
geom=c("point"))
|
||
|
||
|
||
# on choisis 2 clusters
|
||
aux<-silhouette(reskmeanscl[,2-1], daisy(s))
|
||
p22 = fviz_silhouette(aux)+
|
||
theme(plot.title = element_text(size =9))
|
||
|
||
grid.arrange(p11,p12,p21,p22,ncol=2)
|
||
```
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
Après comparaison sur la Figure \ref{fig:clust1}, on choisit 2 classes car cela nous semble plus optimal que 4.
|
||
|
||
### clustering CAH
|
||
|
||
À présent on va afficher l'évolution de l'indice de Calinski-Harabasz en fonction du nombre de classes utilisées pour découper le dendogramme.
|
||
|
||
|
||
```{r, fig.width=5, fig.cap="\\label{fig:cah1}Visualisation du critère de Calinski-Harabasz en fonction du nombre de classes demandées pour le clustering"}
|
||
set.seed(1234)
|
||
|
||
dx<-dist(donnees_transposees) # matrice de distance pour appliquer la méthode classification hierarchique dessus
|
||
hward<-hclust(dx,method = "ward.D2") # construis le dendogramme avec la méthode d'agregation de ward
|
||
|
||
# index.G1 est l'indice de Calinski-Harabasz, code repris du TP
|
||
CH<-NULL
|
||
for (k in 2:Kmax){
|
||
CH<-c(CH,index.G1(x=donnees_transposees, cl= cutree(hward,k)))
|
||
}
|
||
daux<-data.frame(NbClust=2:Kmax,CH=CH)
|
||
ggplot(daux,aes(x=NbClust,y=CH))+
|
||
geom_line()+
|
||
geom_point()+
|
||
xlab("Nombre de classes")+
|
||
ylab("CH")
|
||
```
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
|
||
En regardant la Figure \ref{fig:cah1}, on choisit alors de prendre 2 classes par observation du maximum du graphe.
|
||
|
||
```{r, fig.width=5, warning=FALSE, fig.cap="\\label{fig:dend1}Dendogramme du clustering de l'ACP des variables Tt en tant qu'individus, obtenu par méthode CAH"}
|
||
set.seed(1234)
|
||
|
||
# clustering par découpage du dendogramme
|
||
ClustCH<-cutree(hward,k=2)
|
||
|
||
# plot le dendogramme
|
||
plot1=fviz_dend(hward,k=2,show_labels=FALSE)
|
||
# plot les individus et affiche les clusters
|
||
plot2=fviz_pca_ind(res_pca, habillage = as.factor(ClustCH), geom = c("point"))
|
||
grid.arrange(plot1,plot2,ncol=2)
|
||
```
|
||
|
||
\vspace{1cm}
|
||
|
||
|
||
### Comparaison des clusterings
|
||
A 2 classes nous obtenons une classification qui ne change pas entre chaque méthode. Nous décidons donc qu'il s'agit donc d'un bon choix de nombre de classes.
|
||
|
||
La classification obtenue est en accord avec les observations faites lors de l'ACP, on y retrouve plus ou moins les mêmes groupements : celui majoritarement composé des relevés de T1 avec une majorité de gènes sans changement d'expression relative et celui composé des relevés de T2 et T3.
|
||
|
||
|
||
```{r,fig.width=5,fig.cap="\\label{fig:chord1}Diagramme mettant en valeur l'exacte correspondance des classes entre les deux clusterings"}
|
||
# comparaison kmeans et CAH, code issu de la correction d'un TP
|
||
clust1f = paste("CLkm-",res_kmeans$cluster, sep="")
|
||
clust2f = paste("CLCAH-", cutree(hward,2), sep="")
|
||
chordDiagram(table(clust1f,clust2f))
|
||
```
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
|
||
Sur la Figure \ref{fig:chord1} nous observons que les 2 clustering nous donnent les mêmes clusters peuplés des mêmes gènes.
|
||
|
||
|
||
# Analyse des gènes
|
||
|
||
## Generation de dataExpMoy
|
||
|
||
Nous construisons le jeu de données DataExpMoy contenant la moyenne des expressions sur les réplicats de chaque gène, pour chaque traitement et chaque heure. DataExpMoy est donc une matrice de taille 542×18.
|
||
|
||
```{r,echo=FALSE,results = FALSE,message=FALSE,warning=FALSE,fig.show='hide'}
|
||
# Liste des traitements, heures et réplicats
|
||
traitements <- c("T1", "T2", "T3")
|
||
heures <- c("1H", "2H", "3H", "4H", "5H", "6H")
|
||
replicats <- c("R1", "R2")
|
||
|
||
# Initialiser une matrice vide pour stocker les moyennes
|
||
G <- nrow(T) # Nombre de gènes
|
||
DataExpMoy <- matrix(NA, nrow = G, ncol = length(traitements) * length(heures))
|
||
|
||
#nommer les lignes et les colonnes comme celles des données originales
|
||
colnames(DataExpMoy) <- paste(rep(traitements, each = length(heures)), heures, sep = "_")
|
||
rownames(DataExpMoy) <- rownames(T)
|
||
|
||
|
||
# Calcul des moyennes
|
||
for (t in traitements) {
|
||
for (h in heures) {
|
||
# Colonnes correspondant au traitement et à l'heure pour R1 et R2
|
||
col_R1 <- paste(t, h, "R1", sep = "_")
|
||
col_R2 <- paste(t, h, "R2", sep = "_")
|
||
#print(col_R1)
|
||
#print(col_R2)
|
||
# Calculer la moyenne sur les deux réplicats
|
||
DataExpMoy[, paste(t, h, sep = "_")] <- rowMeans(T[, c(col_R1, col_R2)], na.rm = TRUE)
|
||
}
|
||
}
|
||
|
||
# Convertir en data.frame pour manipulation
|
||
DataExpMoy <- as.data.frame(DataExpMoy)
|
||
head(DataExpMoy)
|
||
```
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
|
||
## ACP de DataExpMoy
|
||
|
||
Pour l'ACP de DataExpMoy, nous avons centré et réduit les données car cela donnait des résultats sensiblement plus exploitables au niveau des indicateurs pour le clustering.
|
||
|
||
|
||
```{r, fig.width=5,fig.cap="\\label{fig:vp2}Participation de chaque valeur propre à l'inertie totale des données"}
|
||
|
||
# on fait l'acp
|
||
res_pca2 = PCA(DataExpMoy, scale.unit = TRUE, graph = FALSE)
|
||
|
||
plot1 = fviz_eig(res_pca2)
|
||
plot1
|
||
```
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
|
||
Sur la Figure \ref{fig:vp2} on voit qu'on dépasse 80% de l'inertie totale avec les deux premières valeurs propres, on prend donc les vecteurs propres associés à ces deux valeurs propres comme composantes principales de notre ACP.
|
||
|
||
|
||
|
||
```{r, fig.width=7,fig.cap="\\label{fig:acp2}Couple de graphes représentant la Projection des individus sur un plan factoriel (droite) et la Corrélations des variables avec les composantes principales (gauche)"}
|
||
plot2 = fviz_pca_ind(res_pca2, label = "point")
|
||
plot3 = fviz_pca_var(res_pca2, axes = c(1, 2), label = "all")
|
||
|
||
# Affichage côte à côte des graphiques
|
||
grid.arrange(plot2,plot3, ncol = 2)
|
||
|
||
```
|
||
|
||
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
|
||
|
||
En analysant la Figure \ref{fig:acp2}, il semble y apparaître deux groupes principaux de gènes (individus), ce qui pourrait refléter une différenciation claire entre les niveaux d'expression sous différents traitements.
|
||
Les points plus éloignés du centre vers le haut représentent des gènes ayant des comportements plus spécifiques ou atypiques.
|
||
|
||
Le premier axe (DIM 1) est fortement corrélé avec les traitements T2, T3 à des moments précis, en particulier pour les temps plus avancés (par exemple, T3 4H, T3 5H, T3 6H). Cela suggère que Dim1 reflète les différences d'expression liées à l'effet des traitements T2 au fil du temps, qui a la particularité de faire exprimer ses gènes sur les temps plus longs. Cela explique aussi la corrélation des variables de relevés de T3 à cette dimension, car T3 est un mélange de T1 et T2.
|
||
Le deuxième axe (DIM2) semble corrélé avec des effets à court terme, typiques de T1. Cela pourrait indiquer des gènes qui réagissent rapidement mais dont l'effet s'atténue à long terme.
|
||
|
||
En regardant la correlations des variables qui représentent les relevés sur T1, comme elles sont très positivement corrélées avec la dimension 2, on en déduis que les gènes vers les valeurs élevées de la dimension 2 sont ceux ciblés par le traitement 1, et 3.
|
||
|
||
Comme T3 est une combinaison des traitements T1 et T2, il semble structurer fortement l'axe principal (Dim1) pour des temps intermédiaires et avancés, mais est notablement corrélé positivement à la dimension 2 pour les premières heures de traitement, correspondant bien avec le fait que T1 commence à faire exprimer les gènes qu'il vise de manière plus rapide que les autres traitements. On retrouve bien le fait que T3 est la somme des deux autres traitements.
|
||
|
||
|
||
### Ajout des variables qualitatives
|
||
Nous allons maintenant afficher les points dans les dimensions de l'ACP, colorés en fonction des modalités des 3 variables qualitatives ExpT1, ExpT2 et ExpT3. Cela peut nous permettre d'affiner notre interprétation des résultats de l'ACP. Cela nous a donné les graphes Figure \ref{fig:proj1}
|
||
|
||
|
||
|
||
```{r, fig.height=7, fig.cap="\\label{fig:proj1}Données sous forme de points dans les dimensions de l'ACP, colorés en fonction des modalités des 3 variables qualitatives ExpT1, ExpT2 et ExpT3."}
|
||
|
||
# on change l'habillage des graphes des individus à chaque ligne pour montrer selon les différentes variables qualitatives quelles modalités ils prennent
|
||
p11 = fviz_pca_ind(res_pca2,habillage=T$ExpT1,geom=c("point"))
|
||
p12 = fviz_pca_ind(res_pca2,habillage=T$ExpT2,geom=c("point"))
|
||
p13 = fviz_pca_ind(res_pca2,habillage=T$ExpT3,geom=c("point"))
|
||
|
||
grid.arrange(p11,p12,p13)
|
||
```
|
||
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
|
||
### analyse
|
||
Sur la Figure \ref{fig:proj1},on observe donc 3 clusters de gènes :
|
||
- Les gènes étant poussés à la sur-expression par T1, affichés comme non-exprimés durant T2 et partiellement sur-exprimés avec T3
|
||
- Les gènes étant sous-exprimés (verts) durant T2 et T3, non-exprimés durant T1.
|
||
- Les gènes étant sur-exprimés (bleus) durant T2 et T3, non-exprimés durant T1.
|
||
|
||
On confirme donc bien notre analyse descriptive préliminaire, et nos suppositions de sens prétés aux dimensions 1 et 2 de l'ACP.
|
||
|
||
|
||
## Clustering
|
||
|
||
L'objectif de ce clustering est de regrouper des individus (ici, des gènes) en groupes homogènes selon leurs similarités.
|
||
|
||
```{r,fig.cap="\\label{fig:sil2}Visualisation de l'inertie intra-classe et Silouhette en fonction du nombre de classes demandées pour le clustering",fig.width=5}
|
||
set.seed(1234)
|
||
|
||
#centrage et réduction des données
|
||
s_2 = scale(DataExpMoy[-c(19:21)])
|
||
|
||
|
||
# code issu des TPs pour afficher le graphe d'inertie intra-classe en fonction du nombre de classe
|
||
reskmeanscl_2<-matrix(0,nrow=nrow(s_2),ncol=Kmax-1)
|
||
Iintra<-NULL
|
||
for (k in 2:Kmax){
|
||
resaux<-kmeans(s_2,k)
|
||
reskmeanscl_2[,k-1]<-resaux$cluster
|
||
Iintra<-c(Iintra,resaux$tot.withinss)
|
||
}
|
||
|
||
df<-data.frame(K=2:15,Iintra=Iintra)
|
||
p11=ggplot(df,aes(x=K,y=Iintra))+
|
||
geom_line()+
|
||
geom_point()+
|
||
xlab("Nombre de classes")+
|
||
ylab("Inertie intraclasse")
|
||
|
||
Silhou<-NULL
|
||
for (k in 2:Kmax){
|
||
aux<-silhouette(reskmeanscl_2[,k-1], daisy(s_2))
|
||
Silhou<-c(Silhou,mean(aux[,3]))
|
||
}
|
||
|
||
df<-data.frame(K=2:Kmax,Silhouette=Silhou)
|
||
p12=ggplot(df,aes(x=K,y=Silhouette))+
|
||
geom_point()+
|
||
geom_line()+theme(legend.position = "bottom")+
|
||
xlab("Nombre de classes")+
|
||
ylab("Silouhette")
|
||
|
||
grid.arrange(p11,p12,ncol=2)
|
||
|
||
```
|
||
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
En lisant le graphe de l'inertie intra-classe (Figure \ref{fig:sil2}), on observe le coude aux alentours de 3 classes. En observant l'indice de Silouhette, on remarque le pic recherché à 3 classes également. Les indices concordants, on en déduit que le choix optimal du nombre de classes pour le clustering est 3.
|
||
|
||
|
||
|
||
```{r,fig.cap="\\label{fig:clust2}Visualisation des clusters générés par la méthode kmeans dans le plan factoriel 1,2",fig.width=5}
|
||
set.seed(1234)
|
||
#kmeans definitif
|
||
res_kmeans_2 = kmeans(s_2,3)
|
||
#graphes des clusters
|
||
fviz_cluster(res_kmeans_2,data=s_2,
|
||
ellipse.type="norm",labelsize=8,
|
||
geom=c("point"))
|
||
|
||
```
|
||
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
|
||
|
||
```{r, fig.width=9,message=FALSE,warning=FALSE, fig.cap="\\label{fig:cah2}Dendogramme du clustering obtenu par méthode CAH",fig.width=5,echo=FALSE, }
|
||
# matrice de distances entre les points
|
||
dx=dist(s_2)
|
||
|
||
# calcule le clustering CAH avec la méthode d'agregation de ward
|
||
hward=hclust(dx,method = "ward.D2")
|
||
|
||
# affiche le dendogramme
|
||
fviz_dend(hward,k=3,
|
||
show_labels = FALSE,
|
||
palette = "npg",
|
||
rect_border = "npg",
|
||
labels_track_height = 0.8)
|
||
```
|
||
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
|
||
### analyse des clusterings
|
||
|
||
On remarque que le clustering de la figure \ref{fig:clust2} correspond aux groupes que nous avions discernés sur la Figure \ref{fig:proj1}. On en déduis donc la correlation des gènes de ces classes avec des expressions moyennes similaires dans chaque classe et distinctes entre les différentes classes.
|
||
|
||
|
||
|
||
### Evolution de l'expression des genes en fonction de leur traitement et cluster
|
||
|
||
Nous avons généré 3 graphiques (Figure \ref{fig:evol1}) , représentant l'évolution de l'expression moyenne des gènes par cluster et par traitement. Cela pourra nous permettre d'affiner nos observations sur les traitements.
|
||
|
||
|
||
|
||
```{r,fig.width=10, fig.cap="\\label{fig:evol1}Graphiques représentant l'évolution de l'expression moyenne des gènes par cluster et par traitement",warning=FALSE}
|
||
# on nomme les gènes dans la matrice
|
||
DataExpMoy$GeneID = rownames(DataExpMoy)
|
||
# stoquage plus utile du clustering
|
||
DataExpMoy$Cluster = as.factor(res_kmeans_2$cluster)
|
||
|
||
# on transforme les données de format large en format long,
|
||
# on renomme la colonne des variables => pratique pour ce qu'on fait après (mais pas indispensable car elle s'appelle "variable" par défaut)
|
||
data_long = melt(
|
||
DataExpMoy,
|
||
cols = -c(Cluster, GeneID),
|
||
variable.name = "Relevé",
|
||
value.name = "Expression"
|
||
)
|
||
|
||
# on sépare la colonne "Relevé" en deux nouvelles colonnes "Traitement" et "Temps", en utilisant "_" comme séparateur pour pouvoir faire un graphe en fonction du temps
|
||
# et 1 par traitement
|
||
data_long = separate(
|
||
data_long,
|
||
col = "Relevé",
|
||
into = c("Traitement", "Temps"),
|
||
sep = "_"
|
||
)
|
||
|
||
# convertit les valeurs de la colonne Temps en nombres après avoir supprimé le "H" du numéro des heures
|
||
# pour cela, on ustilise gsub qui permet de match des expressions regulières
|
||
data_long$Temps = as.numeric(gsub("H", "", data_long$Temps))
|
||
|
||
# groupe les données par Cluster, Traitement et Temps
|
||
data_grouped = group_by(data_long, Cluster, Traitement, Temps)
|
||
|
||
# calcule la moyenne de l'expression génique pour chaque combinaison de Cluster, Traitement et Temps. Chaque moyenne va donner un point du graphe
|
||
data_mean = summarise(
|
||
data_grouped,
|
||
Expression = mean(Expression),
|
||
.groups = "drop"
|
||
)
|
||
|
||
ggplot(data_mean, aes(x = Temps, y = Expression, group = Cluster, color = Cluster)) +
|
||
geom_line(size = 1.2) + # grosseur du trait
|
||
facet_wrap(~ Traitement, scales = "free_y") + # va créer un sous graphique pour chaque traitement, sans forcément la même échelle sur l'axe y
|
||
labs(
|
||
x = "Temps (heures)",
|
||
y = "Expression génique (moyenne)"
|
||
) +
|
||
scale_color_manual(values = c("red", "blue", "green"))
|
||
|
||
|
||
```
|
||
|
||
|
||
|
||
\vspace{1cm}
|
||
|
||
|
||
On lis sur les graphiques de la Figure \ref{fig:evol1} la liaison de toutes les observations que nous avons pu faire :
|
||
|
||
On remarque sur T1 le fait que le cluster 2 correspond au peu de gènes qu'il influence et qui vont se sur-exprimer plus tôt dans le traitement, pour qu'ensuite l'expression relative re-diminue.
|
||
Lors du traitement 2, les gènes de ce cluster ne vont pas/peu s'exprimer et lors du traitement 3, ils vont s'exprimer de la même façon que lors du traitement 1, ce qui est logique sachant que T3 est le mélange de T2 et T1.
|
||
|
||
Comme lors de l'analyse descriptive, on observe les deux autres clusters se sous-exprimer légèrement lors de T1 et se sur et sous exprimer pour T2 et T3, tout en conservant leur expression dans le temps.
|
||
Le cluster 1 correspond aux gènes sur-exprimés et le cluster 3 aux gènes sous-exprimés pour T2 et T3.
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
## Approfondissement de L'ACP avec un clustering de 7 classes
|
||
Sur la figure \ref{fig:sil2}, sur le graphe de l'inertie intraclasse on peut aussi remarquer un coude aux alentours de 7 classes, nous avons donc décidé de compléter notre analyse par un clustering à 7 classes dont on peut voir la projection dans le plan factoriel sur la Figure \ref{fig:clust3} et Figure \ref{fig:evol2} l'évolution de l'expression moyenne des gènes par cluster, en fonction du temps.
|
||
|
||
|
||
|
||
|
||
```{r,fig.cap="\\label{fig:clust3}Visualisation des clusters générés par la méthode kmeans dans le plan factoriel 1,2",fig.width=5}
|
||
set.seed(1234)
|
||
#graphes des clusters
|
||
res_kmeans_3 = kmeans(s_2,7)
|
||
fviz_cluster(res_kmeans_3,data=s_2,
|
||
ellipse.type="norm",labelsize=8,
|
||
geom=c("point"))
|
||
|
||
```
|
||
|
||
|
||
|
||
|
||
|
||
|
||
```{r,fig.width=10, fig.cap="\\label{fig:evol2}Graphiques représentant l'évolution de l'expression moyenne des gènes par cluster et par traitement",warning=FALSE}
|
||
set.seed(1234)
|
||
|
||
DataExpMoy$Cluster = as.factor(res_kmeans_3$cluster)
|
||
|
||
data_long = melt(
|
||
DataExpMoy,
|
||
cols = -c(Cluster, GeneID),
|
||
variable.name = "Relevé",
|
||
value.name = "Expression"
|
||
)
|
||
|
||
data_long = separate(
|
||
data_long,
|
||
col = "Relevé",
|
||
into = c("Traitement", "Temps"),
|
||
sep = "_"
|
||
)
|
||
|
||
data_long$Temps = as.numeric(gsub("H", "", data_long$Temps))
|
||
|
||
data_grouped = group_by(data_long, Cluster, Traitement, Temps)
|
||
|
||
data_mean = summarise(
|
||
data_grouped,
|
||
Expression = mean(Expression),
|
||
.groups = "drop"
|
||
)
|
||
|
||
ggplot(data_mean, aes(x = Temps, y = Expression, group = Cluster, color = Cluster)) +
|
||
geom_line(size = 1.2) +
|
||
facet_wrap(~ Traitement, scales = "free_y") +
|
||
labs(
|
||
x = "Temps (heures)",
|
||
y = "Expression génique (moyenne)"
|
||
)
|
||
|
||
```
|
||
|
||
|
||
|
||
|
||
|
||
### Analyse
|
||
|
||
Avec plus de clusters, on retrouve plus de précisions dans les comportements différents, avec différentes trajectoire d'expression (et donc comportement entre le premier et dernier relevé ) pour les gènes qui finissent sous-exprimés et ceux qui finissent sur exprimés dans T2 et T3 bien que le dernier relevé donne toujours les mêmes tendances globales. On observe par exemple, une classe de gènes qui se sur-exprime, se calme, puis se re-sur exprime. Chose que l'on ne distinguait pas en effectuant les moyennes des comportements sur uniquement 3 clusters : nous avions juste une courbe tendant vers la sur-expression à 6h. Il est tout de même intéressant de noter que la courbe d'expression moyenne pour les gènes sur-exprimés par T1 ne bouge pas par rapport à l'analyse sur 3 clusters, il semblerait donc que ce cluster n'ait sensiblement pas bougé. On peut aussi confirmer les tendances vues dans la première ACP, comme le fait que T1 vise des gènes qui s'expriment à un horizon temporel moyen.
|
||
|
||
|
||
|
||
|
||
## Question BONUS
|
||
|
||
Nous n'avons pas réussi à obtenir quelque-chose que nous étions capables d'interpréter dans la question bonus, mais elle est présente dans le Rmarkdown
|
||
```{r}
|
||
MatriceExp = T[,c(37,38,39)]
|
||
rownames(MatriceExp) = rownames(T)
|
||
colnames(MatriceExp) = c("ExpT1","ExpT2","ExpT3")
|
||
MatriceExp = as.data.frame(MatriceExp)
|
||
```
|
||
|
||
|
||
|
||
|
||
|
||
|
||
```{r, fig.cap="Visualisation des clusters générés par la méthode PAM, dans un plan de L'acp et dans les plans originaux des données",echo=FALSE,results = FALSE,message=FALSE,warning=FALSE,fig.show='hide'}
|
||
# Pour effectuer un clustering PAM avec des données qualitatives, nous avons généré en amont une matrice de dissimilarité grâce à la métrique de gower, pour ensuite lancer l'algorithme PAM sur cette matrice.
|
||
|
||
# La distance de Gower est adaptée aux données qualitatives
|
||
ds <- daisy(MatriceExp, metric = "gower")
|
||
|
||
|
||
resPAMcl<-matrix(0,nrow=nrow(MatriceExp),ncol=Kmax-1)
|
||
Silhou<-NULL
|
||
for (k in 2:Kmax){
|
||
resaux<-pam(ds,k,metric="euclidean")
|
||
resPAMcl[,k-1]<-resaux$clustering
|
||
aux<-silhouette(resPAMcl[,k-1], ds)
|
||
Silhou<-c(Silhou,mean(aux[,3]))
|
||
}
|
||
|
||
df<-data.frame(K=2:Kmax,Silhouette=Silhou)
|
||
ggplot(df,aes(x=K,y=Silhouette))+
|
||
geom_point()+
|
||
geom_line()+theme(legend.position = "bottom")
|
||
|
||
res_pam = pam(ds, k = 10)
|
||
# le champ data doit être remis dans les resultats de PAM
|
||
# car comme on a effectué l'algo en donnat directement la matrice de distance
|
||
# la variable de résultat ne contient pas les données originales
|
||
# Or elles sont necessaire car le champ data de fviz_cluster est utilisé
|
||
# que pour des résultats de kmeans et dbscan
|
||
res_pam$data = ds
|
||
|
||
fviz_cluster(
|
||
res_pam,
|
||
labelsize = 8
|
||
)
|
||
|
||
# visualisation des clusters dans le plan original
|
||
|
||
# Ajouter les clusters aux données d'origine
|
||
res = MatriceExp
|
||
res$cluster = as.factor(res_pam$clustering)
|
||
|
||
ggplot(res, aes(x = ExpT1, y = ExpT2, color = cluster)) +
|
||
geom_point(size = 3) +
|
||
labs(
|
||
x = "ExpT1",
|
||
y = "ExpT2",
|
||
color = "Cluster"
|
||
) +
|
||
theme_minimal()
|
||
ggplot(res, aes(x = ExpT1, y = ExpT3, color = cluster)) +
|
||
geom_point(size = 3) +
|
||
labs(
|
||
x = "ExpT1",
|
||
y = "ExpT2",
|
||
color = "Cluster"
|
||
) +
|
||
theme_minimal()
|
||
|
||
ggplot(res, aes(x = ExpT2, y = ExpT3, color = cluster)) +
|
||
geom_point(size = 3) +
|
||
labs(
|
||
x = "ExpT1",
|
||
y = "ExpT2",
|
||
color = "Cluster"
|
||
) +
|
||
theme_minimal()
|
||
```
|