UE5-ACP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analyse en Composantes Principales sur des transcriptomes de patients atteints de diabète de type 1\n",
    "------------------------------------------------------------------------------------------------------\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dans ce TP vous allez réaliser une Analyse en Composantes Principales (ACP). Cette méthode permet d'évaluer la variabilité de vos données. Elle repose sur une transformation des variables (ici: niveaux de transcrits) partiellement corrélées entre elles en de nouvelles variables dé-corrélées les unes des autres. L'ACP vous permet ainsi d'optimiser la représentation de données massives pour en tirer un maximum d'information.\n",
    "\n",
    "Vous allez analyser des données de transcriptomes générées sur puces d'expression, les mêmes que celles étudiées au cours des TP 6 et 7 de l'UE1.\n",
    "\n",
    "__Les ARN ont été extraits à partir de cultures de lignées lymphoblastoïdes (lymphocytes B immortalisés par le virus EBV) de patients et apparentés (germains partageant les mêmes allèles HLA), stimulés 6h ou 24h avec du Phorbol Myristate Acetate (PMA), ou non stimulés.__\n",
    "\n",
    "Vous allez :\n",
    "- Réaliser l'ACP en utilisant le package FactoMineR de R.\n",
    "- Analyser les graphiques des individus et repérer si les nouveaux axes séparent des groupes d'individus pertinents sur le plan clinique\n",
    "- Analyser les cercles des corrélations et identifier les transcrits les mieux représentés par les axes de l'ACP.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "---\n",
    "\n",
    "## 1. Préparation de l'environnement R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Choix du répertoire de travail et chargement du paquet FactoMineR (déjà installé)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "getwd()\n",
    "list.files() # on verifie que l'objet for_factominer.RData est bien présent sinon on décommente la ligne suivante:\n",
    "#setwd('~/m1_meg-acp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Importation du package FactoMineR:\n",
    "options(width=160)\n",
    "installed.packages()[\"FactoMineR\",c(1,2,3)] # pour vérifier que le paquet est bien installé\n",
    "\n",
    "if (!require(\"FactoMineR\", quietly = T)) {\n",
    "    install.packages(\"FactoMineR\")}  # commande pour installer factominer au cas où -> peut-être utile si vous travailler sur un autre ordinateur/serveur\n",
    "\n",
    "library(FactoMineR)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat(\"Voici mon environnement de travail avec les paquets de R chargés:\\n\")\n",
    "sessionInfo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Importation des données de transcriptome:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "load(\"for_factominer.RData\")\n",
    "ls()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'objet for_factominer contient les variables en colonnes, ainsi que des variables supplémentaires, qualitatives et quantitatives."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Combien ce dataframe contient-elle de lignes et de colonnes ? A quoi correspondent-elles ?\n",
    "Faites attention notamment aux 6 dernières colonnes du dataframe.\n",
    "\n",
    "*__Fonctions à utiliser__*: `str()` ou `dim()` ou `nrow()` et `ncol()`; `head()` et `tail()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ligne de code\n",
    "str(for_factominer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "N'affichez la structure que des 10 dernières colonnes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "str(for_factominer[,(ncol(for_factominer)-9):ncol(for_factominer)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Vos commentaires :"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Autre façon de visualiser le dataframe:\n",
    "-> ici en le transposant\n",
    "-> en utilisant la fonction datatable de la librairie \"DT\" qui permet de naviguer avec les ascenceurs dans la table d'un notebook et de filtrer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if (!require(\"DT\", quietly = T)) {\n",
    "    install.packages(\"DT\")}  \n",
    "options(warn = -1)\n",
    "transposed_data <- as.data.frame(t(for_factominer), stringsAsFactors=F)\n",
    "str(transposed_data)\n",
    "transposed_data$Gene_name <- row.names(transposed_data)\n",
    "transposed_data <- transposed_data[,c(ncol(transposed_data), 1: ncol(transposed_data)-1)]\n",
    "datatable(transposed_data, rownames=F, caption=\"Dataset for ACP analysis\", filter=\"top\", options = list(lengthMenu=c(10,20, 30, 40, 50)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## 2. Analyse en Composantes Principales par FactoMineR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1. Création des résultats de l'ACP:\n",
    "\n",
    "La fonction `PCA()` de FactoMineR crée un objet de type `list` contenant tous les résultats de l'ACP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Utilisation de la fonction PCA() de FactoMineR\n",
    "#?PCA\n",
    "PCAres <- PCA(for_factominer, quali.sup = 1921:1925, quanti.sup = 1926, graph = FALSE)\n",
    "#str(PCAres) # longue structure!\n",
    "print(class(PCAres))\n",
    "names(PCAres)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Analyse visuelle des deux premières composantes de l'ACP\n",
    "Les deux premières composantes rendent comptent au maximum de l'inertie du nuage de points d'origine. Visualisons les résultats grâce à un graphique affichant la position des individus en fonction des deux premières composantes principes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Création de graphiques indiquant la position des échantillons par rapport aux deux premières composantes et enregistrement dans un fichier PDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#pdf(\"PCA_individus.pdf\")\n",
    "\n",
    "plot(PCAres, choix = \"ind\", autoLab = \"yes\", invisible = \"quali\", cex=0.3)\n",
    "plot(PCAres, choix = \"ind\", autoLab = \"yes\", invisible = \"quali\", habillage = \"Stim\", cex=0.3)\n",
    "plot(PCAres, choix = \"ind\", autoLab = \"yes\", invisible = \"quali\", cex=0.3, habillage = \"PedID\")\n",
    "plot(PCAres, choix = \"ind\", autoLab = \"yes\", invisible = \"quali\", cex=0.3, habillage = \"ID\")\n",
    "plot(PCAres, choix = \"ind\", autoLab = \"yes\", invisible = \"quali\", cex=0.3, habillage = \"Status\")\n",
    "plot(PCAres, choix = \"ind\", autoLab = \"yes\", invisible = \"quali\", cex=0.3, habillage = \"Sex\")\n",
    "\n",
    "#dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Pensez à ajouter un titre et une légende!\n",
    "\n",
    "Figure 1. Titre\n",
    "Légende"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Commentez les éléments communs à ces différents graphiques, notamment  :__\n",
    "- à quoi correspondent les axes de ces graphiques ?\n",
    "- que signifient les pourcentages associés à chacun de ces graphiques ?\n",
    "- que représentent les points de ce graphique et les coordonnées de ces points ?\n",
    "- les points son-ils étalés sur ces deux axes ? Constate-t-on des regroupements de certains points ?\n",
    "\n",
    "__Commentez ensuite pour chacun de ce ces graphiques :__\n",
    "- Constatez-vous un lien entre les annotations indiquées en couleurs de ces points et leur position sur chacun des axes créés par l'ACP ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vos commentaires :"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Eléments communs à tous les graphiques :"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Eléments spécifiques à chaque graphique :"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3. Qualité de l'ACP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plus les axes de l'ACP rendent comptent de l'inertie du nuage de points d'origine, meilleure est la qualité de l'ACP (cela signifie que l'on a maintenu au mieux la forme du nuage de points d'origine).\n",
    "L'inertie exprimée par chacune des composantes est enregistrée dans `PCAres$eig`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Inertie traduite par les 5 premières composantes\n",
    "PCAres$eig[1:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Commentez ces valeurs :__\n",
    "- Quel pourcentage de l'inertie du nuage de point d'origine est traduit par les composantes 1 et 2 ? \n",
    "- Quel pourcentage de l'intertie est traduit par les cinq premières composantes ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Vos commentaires :"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il est possible d'afficher sous forme graphique le pourcentage d'inertie traduit par chaque dimension :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Graphique Inertia and dimensions\n",
    "\n",
    "#pdf(\"PCA_inertia.pdf\")\n",
    "eig.val <- PCAres$eig\n",
    "barplot(eig.val[, 2], \n",
    "        names.arg = 1:nrow(eig.val), \n",
    "        main = \"Variances Explained by PCs (%)\",\n",
    "        xlab = \"Principal Components\",\n",
    "        ylab = \"Percentage of variances\",\n",
    "        ylim = c(0,15),\n",
    "        col =\"steelblue\")\n",
    "\n",
    "barplot(eig.val[, 3], \n",
    "        names.arg = 1:nrow(eig.val), \n",
    "        main = \"Variances Explained by PCs (%)\",\n",
    "        xlab = \"Principal Components\",\n",
    "        ylab = \"Cumulative Percentage of variance\",\n",
    "        ylim = c(0,100),\n",
    "        col =\"steelblue\")\n",
    "\n",
    "#dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Comment évolue l'inertie avec les différentes composantes ?__\n",
    "- Combien de composantes sont finalement nécessaires pour rendre compte de 50% de l'inertie du nuage d'origine ?\n",
    "- Comment évolue le pourcentage d'inertie traduit par les composantes ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Vos commentaires :"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## 3. analyse de corrélations entre les variables du dataframe et/ou de l'ACP\n",
    "\n",
    "### 3.1. Corrélation entre les différentes composantes et les variables supplémentaires (traitement, sexe, âge, famille, individu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plutôt que d'étudier de façon visuelle chacune des composantes, il est possible de tester la corrélation entre les différentes composantes et chacune des variables. On peut ainsi mettre en évidence que les échantillons situés sur les valeurs basses de la composante 1 auront par exemple été traitées avec le PMA, tandis que les échantillons situés dans les valeurs hautes n'auront pas été stimulés."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test de corrélation entre les variables supplémentaires qualitatives et les composantes\n",
    "\n",
    "round(PCAres$quali.sup$eta2,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test de corrélation entre les variables supplémentaires quantitatives et les composantes\n",
    "\n",
    "round(PCAres$quanti.sup$cor,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Commentez les résultats de ces tests de corrélation__\n",
    "- Quelle(s) composante(s) vous semble corrélée(s) avec quelle(s) variable(s) supplémentaire(s) ?\n",
    "- Cela va-t-il dans le même sens que vos observations sur les deux premières composantes ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Vos commentaires :"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*__Vous pouvez observer les éventuelles corrélations entre variables qualitatives et composantes via des graphiques !__*\n",
    "\n",
    "En particulier, utilisez l'argument `axes` pour spécifier les dimensions d'intérêt comme dans l'exemple ci-dessous."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#pdf(\"nom_fichier.pdf\")\n",
    "# les composantes par défaut sont les 1 et 2\n",
    "# vous pouvez les modifier via l'argument \"axes\" (cf. exemple ci-dessous pour les composantes3 et 4)\n",
    "\n",
    "plot(PCAres, choix = \"ind\", autoLab = \"yes\", invisible = \"quali\", cex=0.3, habillage = \"Status\", axes = c(3,4))\n",
    "\n",
    "#dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Vos commentaires :"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Corrélation entre les différentes composantes et les variables utilisées par l'ACP (quantité des ARNm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un gène ayant fortement contribué à une composante sera corrélé avec elle ! \n",
    "\n",
    "Cela signifie que :\n",
    "- les individus situés dans les valeurs hautes de la composante expriment fortement ce gène (corrélation, 0 < coeff. de corrélation < 1)\n",
    "\n",
    "- OU les individus situés dans les valeurs hautes de la composante expriment faiblement ce gène (corrélation inverse, -1 < coeff. de corrélation < 0)\n",
    " \n",
    " Le __cercle de corrélation__ révèle la corrélation entre les différentes variables (ici les niveaux de transcrits) et les composantes. Chaque transcrit est représenté par une flèche. La position de l'extrémité de la flèche indique le coefficient de corrélation entre le transcrit et les deux composantes affichées."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cercle de corrélation\n",
    "# L'argument \"contrib 30\" indique que vous choisisser de n'afficher que les 30 transcrits les plus corrélés avec les composantes 1 et 2\n",
    "\n",
    "#pdf(file = \"PCA_diab_contrib30.pdf\")\n",
    "plot(PCAres, choix = \"var\", cex = 0.6, select = \"contrib 30\", unselect = 1)\n",
    "#dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Identifiez quelques transcrits dont les niveaux vous semblent corrélés avec les composantes 1 et 2__\n",
    "- Identifiez les transcrits représenrés par des flèches proches de l'un des axes\n",
    "- S'agit-il d'une corrélation ou d'une corrélation inverse ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Vos commentaires :"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vous pouvez calculer le niveau de corrélation de ces transcrits avec la position de ces échantillons sur les axes qui vous intéressent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lignes de codes pour calculer les coefficients de corrélation entre le niveau des transcrits et la position des échantillons sur les axes\n",
    "cor.test(for_factominer$SLAMF1,PCAres$ind$coord[,1])\n",
    "cor.test(for_factominer$SLAMF1,PCAres$ind$coord[,2])\n",
    "cor.test(for_factominer$RGS16,PCAres$ind$coord[,1])\n",
    "cor.test(for_factominer$RGS16,PCAres$ind$coord[,2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Vos commentaires : (lien entre la flèche du cercle de corrélation et les valeurs que vous venez de calculer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Cette analyse a-t-elle répondu à votre question scientifique ? Quelles analyses pourriez vous imaginer afin d'aller plus loin ?__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "---\n",
    "\n",
    "# Consignes pour votre compte-rendu de TP\n",
    "\n",
    "Vous allez réaliser un compte-rendu de TP __sous la forme d'un nouveau Notebook__ récapitulant les informations importantes ce TP.\n",
    "Les consignes sont les suivantes :\n",
    "\n",
    "- compte-rendu à remettre sur Moodle pour le 14 décembre 2020\n",
    "- un compte-rendu par binôme\n",
    "- les plagiats sont interdits. Des comptes-rendus considérés comme excessivement proches par le correcteur pourront entraîner des procédures disciplinaires\n",
    "- Nombre limite de caractères : 4000 en dehors des légendes des figures\n",
    "- Les figures doivent être légendées. La légende inclut le titre de la figure et doit être suffisamment explicite, en quelques lignes, pour que votre lecteur soit capable de la comprendre (et de se forger sa propre opinion) sans avoir à lire votre texte.\n",
    "- Le compte-rendu inclura : __un titre__, __une introduction__ situant  votre  analyse  dans  les  recherches sur le diabète. Elle amène à la question posée par cette étude(0.5 pt avec le titre) ; __les résultats__ qui décriventles résultats obtenus et concluent de façon factuelle ce que vous avez découvert via l'ACP. Ils seront illustrés par __une/des figure(s)__. Vous devez donc citer dans le  texte de façon pertinente les  figures ou parties de figures illustrant les résultats que vous êtes en train de décrire. (2,5 pts texte + 2 pts figure(s)) ; __une discussion/conclusion__ dans laquelle vous pourrez donner 1) une interprétation des résultats allant au delà d'une pure description, 2) des idées d'expériences possibles pour aller au-delà et 3) une phrase d'ouverture. Dans certains articles la discussion peut contenir un schéma de modèle (mais cela ne vous est pas demandé pour ce TP).(0,5 pt)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}