Astuce :
Représenter rapidement une ACP avec R et ggplot2

Je ne sais pas pour vous, mais moi, à chaque fois que j'assiste à une réunion de labo, il y a quasi systématiquement un graphique d'ACP pour montrer les données. Et à chaque fois, il s'agit d'un graphique de base, généré avec R, avec la fonction plot(), des couleurs qui piquent les yeux et des axes et légendes illisibles. La critique est facile me direz-vous, j'avoue avoir moi aussi présenté ce genre de graphique assez souvent. Mais au bout d'un moment, vu que les ACP ça devient très très routinier dans le travail d'un "data analyste", la nécessité de produire rapidement et facilement un graphique ACP se fait sentir.

À cette occasion, j'ai retroussé mes manches et j'ai écrit un petit bout de code maison qui génère tout seul les graphiques pour un nombre donné de composantes principales, et ce en une seule fonction. J'en ai profité pour utiliser ggplot2, dont je vous invite à lire la présentation, pour rendre ces graphiques un peu plus esthétiques. Je partage donc avec vous mon secret pour faire des graphiques sexy. Ce n'est pas très compliqué en soit, c'est juste des graphiques, mais je vous donne mon code si jamais ça peut aussi vous servir et vous épargner du temps.

Aperçu du type de graphique que produit le code qui va suivre

Première étape : chargement de la librairie

Comme je l'ai dit plus haut, nous allons utiliser la librairie ggplot2. Cette fonction permet d'installer automatiquement le paquet s'il n'est pas déjà installé.

Deuxième étape : le screeplot

Lorsque l'on fait une ACP, on aime bien savoir quelles composantes principales (PC) accumulent le plus de variance (pour les détails sur ce qu'est une ACP, c'est par là). Pour aller vite, on peut simplement lancer la commande screeplot(), mais le rendu n'est pas hyper sexy. En plus il n'y a pas de légende sur l'axe des abscisses, ce qui est fâcheux quand on a plus de 20 PC. D'ailleurs, quand on lance une ACP sur un set de données avec plus de 100 échantillons (si si, ça arrive), on se retrouve avec un barplot tout serré où en réalité l'information qui nous intéresse est compressée sur la gauche du graphique.

Pour résoudre cela, la fonction suivante prend en entrée le résultat de l'ACP (calculée avec prcomp), et le nombre de PC à afficher. C'est pas le plus "hot" des graphiques mais ça peut être pratique.

Troisième étape : l'ACP

Passons au nerf de la guerre : l'ACP en elle même. Pour apprécier la complexité d'un set de données, il est souvent nécessaire de regarder un peu plus loin que juste les 2 premières composantes. Ce qui me fatigue le plus (il ne me faut pas grand chose), c'est de relancer la commande de graphique pour chaque combinaison. Pour satisfaire ma fainéantise, il me faut donc une fonction avec le nombre de PC voulues et hop, ça me sort toutes les combinaisons possibles.

La fonction suivante prend en entrée le résultat d'une ACP (calculée avec prcomp), le nombre de PC à regarder, les conditions des échantillons (une liste qui fait la même taille que le nombre d'échantillons), et une palette de couleur en hexadécimale (par exemple: "#fb7072"), avec autant de couleurs que de conditions différentes (e.g. deux conditions, cas et contrôle, donc 2 couleurs).

Et enfin : la démo

Et voilà le résultat :

Screeplot et ACP

Résultat des commandes précédentes. À gauche le screeplot, à droite l'ACP avec les deux premières composantes et le pourcentage de variance entre parenthèse, le tout coloré en fonction des conditions (1 à 3).

 

Juste une dernière remarque : si vous spécifiez plus de 2 PC à afficher, ça va générer les graphiques chacun leur tour, donc je vous conseille d'appeler la fonction plot_pca() dans un pdf pour les sauvegarder dans un seul fichier. Peut-être un jour j'essayerai le paquet gridExtra pour afficher plusieurs graphiques sur une même page... Un jour...

 

Merci à Jnsll, Nico M., et Mat Blum pour leurs commentaires et relectures.

  • À propos de
  • Mi-bio, mi-bioinfo, et re-mi-bio derrière. Je suis actuellement en postdoc en épigénomique développementale de la Drosophile à l'IGFL (Lyon). J'ai suivi une licence de biologie cellulaire et génétique et un master de bioinformatique à l'université de Rennes, puis j'ai travaillé comme ingé d'étude en développement web pendant 1 an et demi. Enfin j'ai effectué un doctorat en bioinformatique à l'université de Genève (single-cell RNA-seq de la gonade de souris en développement).

Catégorie: Astuce, Didacticiel, Suivez l'guide | Tags: , , , ,

7 commentaires sur “Représenter rapidement une ACP avec R et ggplot2

  1. Super article, merci 🙂

    Pour les figures multi-panneaux, j'aime bien cowplot, que je trouve plus facile d'utilisation que juste gridExtra:
    https://cran.r-project.org/web/packages/cowplot/vignettes/introduction.html

    Sinon, le package ggfortify permet de ploter les objet issus de prcomp() directement:
    https://cran.r-project.org/web/packages/ggfortify/vignettes/plot_pca.html
    Il y a quelques features en plus, comme ploter la projection des axes d'origines sur le nouveau plan, et des features en moins, comme l'absence du pourcentage de variances sur les axes.

  2. Bonjour et merci pour ton article.

    Pour que le code fonctionne chez moi j'ai été obligé de modifier les deux lignes suivantes:

    # pour le barplot:
    geom_bar(stat="identity") au lieu de geom_col()

    # pour la fonction apply():
    scale_colour_manual(values = col, name = "") au lieu de scale_fill_manual(values=colours, name="")

    • Merci pour ce commentaire.
      J'ai en effet oublié de préciser que j'utilise ggplot2 version 2.2.1. Il peut donc y avoir de petites subtilités comme celle-ci si vous utilisez une version différente du paquet.

  3. Merci ! Pour compléter : l'excellent package explor (https://cran.r-project.org/web/packages/explor/vignettes/introduction_fr.html) permet d'"explorer" de façon interactive (interface shiny) les resultats d'acps (ou acms) provenant de differents packages R. Indispensable !

  4. Bonjour ! J'ai utilisé votre script pour mes données qui consistent en un dataframe de 9 variables numériques qui peuvent prendre des valeurs nulles. Voici ce que je remplace au script pour composer mes groupes :

    library(FactoMineR)
    hcpc_db<-HCPC(db)
    hcpc_db$data.clust
    db_clust<-as.data.frame(hcpc_db$data.clust[,"clust"])
    db.clust<-cbind(db,db_clust)
    group<-as.factor(db.clust$clust)

    et lors de l'application de la fonction du graphique ACP "plot_pca" voici ce que R me retourne :
    Error: Aesthetics must be either length 1 or the same as the data (61): fill

    J'ai épluché un peu les forums pour modifier le script en conséquence mais je n'ai pas réussi à me défaire de cette erreur. Je pense que je dois caler un length=length(db) quelque part...

    Sauriez vous m'aider s'il vous plait ?

    • Bonjour Rachel,

      ce script est fait pour prendre le résultat d'une ACP calculée avec prcomp. HCPC a un format de sortie différent de prcomp, il va donc falloir adapter le code pour remplacer toutes les allusions au format de sorti de prcomp (tous les "pca$...") par leurs équivalents HCPC.

      Si tu n'y arrives pas, je m'y collerai de mon côté pour t'aider.

      • Ah d'accord ! Merci beaucoup, je vais tenter et dans ce cas si j'ai encore un souci je reviendrais vers vous !

Laisser un commentaire