--- title: "Projet" output: pdf_document: default html_document: default 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) ``` ```{r} set.seed(1234) 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) #centrer T head(T) summary(T) str(T) levels(T$ExpT1) ``` # A REFAIRE LES ANALYSES, JAVAIS PAS VU QUE T3 ETAIT UNE COMBINAISON DE T1 ET T2 # Contenu du jeu de données : - 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èeme heure lors du traitement $$T \in \{T1,T2,T3\}$$ - $$3*6 + 3*6 = 36$$ variables quantitatives continues représentant les effets des traitements sur l'expression des gênes T1 T2 et T3 à 1h,2h,3h,4h,5h,6h après l'administration pour les replicats R1 et 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 unidimentionnelle : ### LES 3 PROCHAINS PLOTS N'APPORTENT RIEN PAR RAPPORT AUX HISTOGRAMMES, JE LES AI DONC CACHES ### Expression des gênes lors du traitement T1 ```{r,echo=FALSE,results = FALSE,message=FALSE,warning=FALSE,fig.show='hide'} g1<-ggplot(T, aes(x=T$ExpT1))+ geom_bar()+ ylab("Effectifs")+ggtitle("Effectifs") g2<-ggplot(T, aes(x = T$ExpT1)) + geom_bar(aes(y = (..count..)/sum(..count..)))+ylab("")+ggtitle("Frequences") df <- data.frame(group = levels(T$ExpT1), value = as.vector(table(T$ExpT1))/nrow(T)) g3<-ggplot(df, aes(x="", y=value, fill=group))+ geom_bar(width = 1, stat = "identity")+ coord_polar("y", start=0)+ theme(legend.position="bottom") grid.arrange(g3,g1,g2,ncol=3,top=textGrob("Visualualisation des expressions relative des gènes lors du traitement T1")) ``` ### Expression des gênes lors du traitement T2 ```{r,echo=FALSE,results = FALSE,message=FALSE,warning=FALSE,fig.show='hide'} g1<-ggplot(T, aes(x=T$ExpT2))+ geom_bar()+ ylab("Effectifs")+ggtitle("Effectifs") g2<-ggplot(T, aes(x = T$ExpT2)) + geom_bar(aes(y = (..count..)/sum(..count..)))+ylab("")+ggtitle("Frequences") df <- data.frame(group = levels(T$ExpT2), value = as.vector(table(T$ExpT2))/nrow(T)) g3<-ggplot(df, aes(x="", y=value, fill=group))+ geom_bar(width = 1, stat = "identity")+ coord_polar("y", start=0)+ theme(legend.position="bottom") grid.arrange(g3,g1,g2,ncol=3,top=textGrob("Visualualisation des expressions relative des gènes lors du traitement T2")) ``` ### Expression des gênes lors du traitement T3 ```{r,echo=FALSE,results = FALSE,message=FALSE,warning=FALSE,fig.show='hide'} g1<-ggplot(T, aes(x=T$ExpT3))+ geom_bar()+ ylab("Effectifs")+ggtitle("Effectifs") g2<-ggplot(T, aes(x = T$ExpT3)) + geom_bar(aes(y = (..count..)/sum(..count..)))+ylab("")+ggtitle("Frequences") df <- data.frame(group = levels(T$ExpT3), value = as.vector(table(T$ExpT3))/nrow(T)) g3<-ggplot(df, aes(x="", y=value, fill=group))+ geom_bar(width = 1, stat = "identity")+ coord_polar("y", start=0)+ theme(legend.position="bottom") grid.arrange(g3,g1,g2,ncol=3,top=textGrob("Visualualisation des expressions relative des gènes lors du traitement T3")) ``` ### Analyse de de ces 3 variables 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 rapporrt 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 #### Traitement T1 ```{r, fig.height=5,fig.cap="Histogrammes de l'expression relative des gènes aux différents relevés lors du traitement T1"} #apply(T[-c(37:39)],2,function(col){ # which(T == col) #hist(col, main = paste("Histogram of", colnames(T)[which(T == col)[2]]), # xlab = "Values", col = "lightblue", border = "black") #}) 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) + labs(title = "Histogrammes de l'expression relative des gènes aux différents relevés lors du traitement T1", x = "Valeurs", y = "Effectifs") ``` ### Traitement T2 ```{r, fig.height=5,fig.cap="Histogrammes de l'expression relative des gènes aux différents relevés lors du traitement T3"} 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") ``` ### Traitement T3 ```{r, fig.height=5,fig.cap="Histogrammes de l'expression relative des gènes aux différents relevés lors du traitement T3"} 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") ``` Nous observons bien une concordance avec l'analyse des expressions des gênes figure . 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. ### boxplots pour faire joli ```{r} ggplot(melt(T[1:18]),aes(x=variable,y=value))+ geom_boxplot()+ theme(axis.text.x = element_text(angle=90, vjust = 0.5, hjust = 1)) ggplot(melt(T[19:36]),aes(x=variable,y=value))+ geom_boxplot() + theme(axis.text.x = element_text(angle=90, vjust = 0.5, hjust = 1)) ``` On voit que même sans réduire les données, chaque variable s'exprime environ avec la même intensité donc pas besoin de centrer réduire lors de l'ACP ## Boxplots par paire de variables (qualitative,quantitative) ```{r} # traitement 1 corrélation avec l'expression des genes du T1 T2 et T3 #ggplot(T,aes(y=T$T1_6H_R1,x=T$ExpT1))+ #geom_boxplot() #ggplot(T,aes(y=T$T1_6H_R2,x=T$ExpT1))+ #geom_boxplot() 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() 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 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() #ggplot(T,aes(y=T$T2_6H_R1,x=T$ExpT2))+ #geom_boxplot() #ggplot(T,aes(y=T$T2_6H_R2,x=T$ExpT2))+ #geom_boxplot() 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 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() 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() #ggplot(T,aes(y=T$T3_6H_R1,x=T$ExpT3))+ #geom_boxplot() #ggplot(T,aes(y=T$T3_6H_R2,x=T$ExpT3))+ #geom_boxplot() ``` ### Analyse des boxplots : - traitement 1 (réplicats 1 et 2) Les genes sur-exprimés au T1 sont non-exprimé durant le T2. Il est difficile d'observer une catégorie de genes de T1 qui se soient sous exprimés dans T2. De même pour la sur-expression dans T2. Ceux qui s'étaient sur-exprimés au T1 ont affiché aucun changement semblent ne pas avoir changé d'expression au T3 ( relativement à l'absence de traitement). Ceux qui n'avaient pas changés d'expression relative durant T1, se sont sous exprimé ou sur exprimé. Il faut bien se rappeler que T1 ne cible que très peu de genes donc ces observations sont cohérentes. - traitement 2 On observe une légère tendance des genes s'étant sous-exprimé avec T1 à se sous exprimer avec T2 mais le reste des boxplots ne permettent pas de relever quoi que ce soit par rapport à T2. En revanche, il est très clair que T2 et T3 ciblent les mêmes genes, toutes les expressions relevées par T2 concordent aux modalités qualitatives moyennes calculées sur T3. ## Analyse bi-dimensionnelle ### matrice de covariance des variables quantitatives ```{r, fig.height = 18,fig.cap="Visualisation de la matrice de covariance des variables quantitatives"} cr = cor(T[-c(37:39)]) corrplot(cr,method="ellipse", type="lower", bg = "lightgrey",title ="Visualisation de la matrice de covariance des variables quantitatives" ) ``` ### Regression linéaire des variables 2 à 2 ??? ### librairie biostatR ```{r} library(BioStatR) # Calculate eta² for Treatment 1 print("T1 vs T2") eta2(T$T1_6H_R1, T$ExpT2) eta2(T$T1_6H_R2, T$ExpT2) print("T1 vs T3") eta2(T$T1_6H_R1, T$ExpT3) eta2(T$T1_6H_R2, T$ExpT3) print("T2 vs T3") eta2(T$T2_6H_R1, T$ExpT3) eta2(T$T2_6H_R2, T$ExpT3) ``` Le calcul du rapport de correlation eta² 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. ### table de contingence pour les variables quali 2 à 2, mosaic plot ? ```{r} 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 genes. Plus finement, on peut confirmer l'observation faite sur les boxplots tendant à dire que le peu de genes 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. ## Menez une analyse en composantes principales où les Tt sH Rr sont les individus décrits par les gènes. 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 sur un jeu de données centrées réduites pour que chaque variable s'exprime avec la même force dans les résultats et qu'ils soient lisibles. ```{r, fig.width=10, fig.height=5, fig.cap="Participation des chaque valeur propre de la matrice de correlation à l'intertie totale des données"} donnees_transposees = t(T[-c(37:39)]) res_pca<-PCA(donnees_transposees,scale.unit=TRUE,graph=FALSE) res_pca$eig fviz_eig(res_pca,title="Participation des chaque valeur propre de la matrice de correlation à l'intertie totale des données") ``` Ce graphique représente les valeurs propres de la matrice de corrélation du jeu de données. L'inertie totale des données étant la somme des valeurs propres ( qui elles sont les inerties axiale associées à l'axe de vecteur directeur un vecteur propre associé ), chaque valeur propre est donc une fraction de l'inertie totale. On voit qu'on dépasse 80% de l'inertie totale des points rien qu'avec les deux premieres valeurs propres, on en prend donc les vecteurs propres associés respectivement à chacune de ces deux valeurs propres comme axes principaux de l'analyse. ```{r,fig.cap="Corrélations des variables avec les composantes principales"} fviz_pca_ind(res_pca,label="all",title="Projection des individus sur un plan factoriel") fviz_pca_var(res_pca,axes=c(1,2),label="none",title="Corrélations des variables avec les composantes principales") ``` Contexte : 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 ## Interprétation globale du couple de graphes On voit que les genes se polarisent principalement sur l'axe 1 dans un sens ou l'autre. Les 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. Il n'y a pas de tendance particulière sur la direction selon l'axe 2 des flèches : 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. 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 EST SENSE VOIR UN TRUC IMPORTANT D'APRES LA PROF MAIS JE VOIS RIEN ### Les axes générés par l'ACP sont purement virtuels, mais on peut y préter un sens plus ou moins défini selon les variables qui y sont corrélées axe 1: force du changement d'expression des gènes : si on sur-exprime/sous-exprime de beaucoup ( genre 0=>5) ou pas, polarisation ? axe 2 : vers le haut => ça va changer d'expression bientot, vers le bas => ça se stabilise ( donc les genes pointant vers le haut changent d'expression, et ceux vers le bas ont tendance à rester plutot pareil ) > On remarque que les traitements ne sont pas du tout décris de la même manière selon les horaires, mais il y a une régularité d'apparition dans chaque tranche horaire. ### hypothèses sur la signification : qu'est-ce qu'ils ont en commun ces gènes qui pourrait décrire ces relevés aux différentes heures et les différents traitements En regardant le graphe des individus (résultats aux heures de relevés), on a effectivement les heures groupées à des valeurs negatives de l'axe 1 correspondant aux relevés du traitement 1 qui, souvenons-nous toujours des histogrammes, ne change l'expression relative que de très peu de gènes. Pour l'interprétation du second axe, les gènes semblent y être positivement et negativement corrélés quel que soit leur correlation avec l'axe 1. En regardant les individus, on observe que plus l'heure est tardive, plus ils tendent vers des valeurs negatives sur l'axe 2. De plus, on observe que les points liés aux relevés du traitement 1 ne vont pas énormément vers les valeurs positives. Il semble donc que l'axe 2 soit indicateur des expressions des gènes sont susceptibles de changer. ## Clustering ### Clustering k-means ```{r,fig.cap="Visualisation de l'inertie intra-classe en fonction du nombre de classes demandées pour le clustering"} #centrage et réduction des données s = donnees_transposees 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") ``` nn pense dénoter un coude dans la courbe d'inertie intraclasse aux alentours de 4 clusters ```{r,fig.cap="Visualisation du critère de Silhouette et de Calinski-Harabasz en fonction du nombre de classes demandées pour le clustering"} 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) ggplot(df,aes(x=K,y=Silhouette))+ geom_point()+ geom_line()+theme(legend.position = "bottom") CH<-NULL Kmax<-20 for (k in 2:Kmax){ CH<-c(CH,index.G1(x=s , cl= kmeans(s,k)$cluster)) } daux<-data.frame(NbClust=2:Kmax,CH=CH) ggplot(daux,aes(x=NbClust,y=CH))+ geom_line()+ geom_point() ``` Silhouette et l'indice de Calinski-Harabasz ont un pic à 2 ### visualisation du clustering ```{r,fig.cap="Visualisation des clusters générés par la méthode kmeans dnas le plan factoriel 1,2"} #graphes des clusters res_kmeans = kmeans(s,2) fviz_cluster(res_kmeans,data=donnees_transposees, ellipse.type="norm",labelsize=8, geom=c("point"))+ggtitle("Visualisation des clusters générés par la méthode kmeans dnas le plan factoriel 1,2") # on choisis 2 clusters aux<-silhouette(reskmeanscl[,2-1], daisy(s)) fviz_silhouette(aux)+ theme(plot.title = element_text(size =9)) rm(df,Silhou,aux) ``` On remaque bien 2 clusters ### PAM ```{r, fig.cap="Visualisation des clusters générés par la méthode PAM dnas le plan factoriel 1,2"} res_pam = pam(s,2,metric="euclidean") fviz_cluster(res_pam,data=donnees_transposees, ellipse.type="norm",labelsize=8, geom=c("point"))+ggtitle("Visualisation des clusters générés par la méthode PAM dnas le plan factoriel 1,2") ``` ### clustering CAH ```{r, fig.width=9, fig.cap="Dendogramme du clustering de l'ACP des variables Tt en tant qu'individus, obtenu par méthode CAH"} dx<-dist(donnees_transposees) hward<-hclust(dx,method = "ward.D2") CH<-NULL Kmax<-20 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() ClustCH<-cutree(hward,k=2) fviz_dend(hward,k=2,show_labels=FALSE)+ggtitle("Dendogramme du clustering de l'ACP des variables Tt en tant qu'individus, obtenu par méthode CAH") fviz_pca_ind(res_pca, habillage = as.factor(ClustCH), geom = c("point")) ``` ### Comparaison des clusterings On voit bien qu'à 4 classes, les regroupements ne sont pas consistents entre chaque méthode de clustering. A 3 classes nous obtenons une classification qui ne change pas, ou presque, 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, celui composé des relevés de T2 et T3 aux heures des changements d'expression les plus brutaux, et finalement celui s'étalant sur la droite qui semble représenter la fin de l'évolution des traitements T2 et T3 où l'expression de gènes y est très polarisée. FAIRE LES GRAPHES RONDS CHELOUS POUR LES CLUSTERINGS 2 A 2 ```{r} # kmeans et PAM clust1f = paste("CLkm-",res_kmeans$cluster, sep="") clust2f = paste("CLCAH-", res_pam$clustering, sep="") chordDiagram(table(clust1f,clust2f)) # kmeans et CAH clust1f = paste("CLkm-",res_kmeans$cluster, sep="") clust2f = paste("CLCAH-", cutree(hward,2), sep="") chordDiagram(table(clust1f,clust2f)) # PAM ET CAH clust1f = paste("CLkm-",res_pam$clustering, sep="") clust2f = paste("CLCAH-", cutree(hward,2), sep="") chordDiagram(table(clust1f,clust2f)) ``` # ANALYSE DES GENES ## Generation de dataExpMoy ```{r} # 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)) colnames(DataExpMoy) <- paste(rep(traitements, each = length(heures)), heures, sep = "_") rownames(DataExpMoy) <- rownames(T) colnames(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) DataExpMoy$ExpT1 = T$ExpT1 DataExpMoy$ExpT2 = T$ExpT2 DataExpMoy$ExpT3 = T$ExpT3 head(DataExpMoy) ``` ## ACP de DataExpMoy ```{r, fig.width=13} # Effectuer l'ACP res_pca2 <- PCA(DataExpMoy, scale.unit = TRUE,quali.sup = c(19:21), graph = FALSE) # Création des graphiques plot1 <- fviz_eig(res_pca2, title = "Participation de chaque valeur propre à l'inertie totale des données") plot2 <- fviz_pca_ind(res_pca2, label = "point", title = "Projection des individus sur un plan factoriel") plot3 <- fviz_pca_var(res_pca2, axes = c(1, 2), label = "all", title = "Corrélations des variables avec les composantes principales") plot1 # Affichage côte à côte des graphiques grid.arrange(plot2,plot3, ncol = 2) ``` ### Ajout des variables qualitatives ```{r} fviz_pca_ind(res_pca2,habillage=19,geom=c("point")) fviz_pca_ind(res_pca2,habillage=20,geom=c("point")) fviz_pca_ind(res_pca2,habillage=21,geom=c("point")) ``` ## Clustering De L'ACP ```{r,fig.cap="Visualisation de l'inertie intra-classe en fonction du nombre de classes demandées pour le clustering"} #centrage et réduction des données s_2 = DataExpMoy[-c(19:21)] Kmax<-15 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) ggplot(df,aes(x=K,y=Iintra))+ geom_line()+ geom_point()+ xlab("Nombre de classes")+ ylab("Inertie intraclasse") ``` ```{r,fig.cap="Visualisation du critère de Silhouette en fonction du nombre de classes demandées pour le clustering"} 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) ggplot(df,aes(x=K,y=Silhouette))+ geom_point()+ geom_line()+theme(legend.position = "bottom") aux<-silhouette(reskmeanscl_2[,3-1], daisy(s_2)) fviz_silhouette(aux)+ theme(plot.title = element_text(size =9)) rm(df,Silhou,aux) ``` ```{r,fig.cap="Visualisation des clusters générés par la méthode kmeans dnas le plan factoriel 1,2"} #graphes des clusters res_kmeans_2 = kmeans(s_2,3) fviz_cluster(res_kmeans_2,data=s_2, ellipse.type="norm",labelsize=8, geom=c("point"))+ggtitle("Visualisation des clusters générés par la méthode kmeans dnas le plan factoriel 1,2") ``` ```{r, fig.cap="Visualisation des clusters générés par la méthode PAM dnas le plan factoriel 1,2"} res = pam(s_2,3,metric="euclidean") fviz_cluster(res,data=s_2, ellipse.type="norm",labelsize=8, geom=c("point"))+ggtitle("Visualisation des clusters générés par la méthode PAM dnas le plan factoriel 1,2") ``` ```{r, fig.width=9, fig.cap="Dendogramme du clustering de l'ACP, obtenu par méthode CAH"} dx=dist(s_2,method="euclidian") hward=hclust(dx,method = "ward.D2") fviz_dend(hward,k=3, show_labels = FALSE, rect=TRUE, rect_fill = TRUE, palette = "npg", rect_border = "npg", labels_track_height = 0.8)+ggtitle("Dendogramme du clustering de l'ACP, obtenu par méthode CAH") ``` ## Evolution de l'expression des genes en fonction de leur traitement et cluster : ```{r,fig.width=13} # Ajouter une colonne unique d'identification pour chaque gène DataExpMoy$GeneID = rownames(DataExpMoy) # Ajouter les identifiants des gènes DataExpMoy$Cluster = as.factor(res_kmeans_2$cluster) # Ajouter les clusters comme une colonne # Restructurer les données en format long data_long = DataExpMoy %>% pivot_longer( cols = -c(ExpT1, ExpT2, ExpT3, Cluster, GeneID), # Colonnes des expressions moyennes names_to = "Condition", values_to = "Expression" ) %>% separate(Condition, into = c("Traitement", "Temps"), sep = "_") %>% # Séparer Traitement et Temps mutate( Temps = as.numeric(gsub("H", "", Temps)) # Supprimer "H" et convertir en numérique ) #head(data_long) # Tracer les courbes ggplot(data_long, aes(x = Temps, y = Expression, group = GeneID, color = Cluster)) + geom_line(alpha = 0.6) + # Ajout des courbes avec transparence facet_wrap(~ Traitement, scales = "free_y") + # Un graphique par traitement labs(title = "Évolution de l'expression génique par traitement et cluster", x = "Temps (heures)", y = "Expression génique") + theme_minimal() + scale_color_manual(values = c("red", "blue", "green")) # Définir des couleurs fixes pour les clusters ``` ## Clustering des gènes à partir des variables ExpT1, ExpT2 et ExpT3. Comme ce sont des données qualitatives, on doit utiliser des méthodes alternatives de clustering comme les kmodes ou les kmeans mais sur les coordonnées de points projetés dans le plan d'une ACM ### K-modes ```{r,fig.cap=""} MatriceExp = T[,c(-1,-2,-3)] rownames(MatriceExp) = rownames(T) colnames(MatriceExp) = c("ExpT1","ExpT2","ExpT3") MatriceExp = as.data.frame(MatriceExp) MatriceExp$ExpT1 = T$ExpT1 MatriceExp$ExpT2 = T$ExpT2 MatriceExp$ExpT3 = T$ExpT3 clkmodes=kmodes(MatriceExp,3,iter.max=100,weight=FALSE) clusplot(MatriceExp, clkmodes$cluster, color=TRUE, shade=TRUE, labels=0, lines=0) ``` ### K-means sur ACM ```{r,fig.cap=""} head(MatriceExp) afcm=MCA(MatriceExp,graph = FALSE) coeff=afcm$ind$coord clkmeans=kmeans(coeff,3) # Visualisation des clusters fviz_cluster(clkmeans, data = coeff, geom = "point", ellipse.type = "norm", main = "Clustering des gènes basé sur l'ACM et k-means") ``` J'ai absolument rien compris à ce qu'il se passe là, il y a beaucoup moins de nuance dans les données donc les points sont empilés. On dénote peut-être des clusters vaguement ressemblant à ceux sur L'ACP des moyennes ? ExpT étant calculée sur les relevés à 6h, on s'éloigne tout de même des clustering précédents. # TODO : - refaire sans centrer-réduire - interpréter par rapport aux méta-variable - ne pas afficher le clustering dans les plans de l'acp pour les 2 dernières ACP - présenter chaque méthode/algo et pourquoi on l'utilise avant le code - analyser la seconde et troisième acp et les clustering qu'on en fait - refaire l'analyse de la première acp - choisir des graphes - voir si la troisième acp est bien ce qu'il faut psq wtf ??