Voici donc la suite de mon premier cours de R pour débutants pressés où nous avions vu les bases. Je vous invite à le lire si vous ne l’avez pas déjà fait.
Aujourd’hui le but est de s’approcher de la biologie et de traiter d’un cas concret. Nous avons des données de différents éléments fonctionnels du génome de D. melanogaster et leurs proportions de conservation et de contrainte évolutive. Nous souhaitons représenter la proportion dans laquelle chaque élément est contraint ou pas. Pour cela, nous allons faire un graphique en barre ou barplot, avec barres d'erreurs. Ce que nous allons voir aujourd'hui peut très bien s'appliquer à d'autres types de données : représenter la quantité de gènes suite à une qPCR ou encore le nombre de fois qu'une souris a choisi le bon chemin dans le labyrinthe suivant chaque condition choisie. L'important n'est donc pas ici le fond, mais bien la forme.
À la fin de ce tutoriel, vous devriez être capables de faire quelque chose comme ça :
R et les éditeurs de texte
Nous allons tout d'abord ouvrir R et un document (éditeur de texte intégré à R), qui va vous servir à sauvegarder le code de votre graphique, la console R n'étant pas forcément la plus pratique.
La console R :
Pour ouvrir un document, cliquez sur "Fichier/Nouveau Document" ou ctrl+N sur Windows et GNU/Linux, cmd+N sur mac. Une bonne habitude à prendre est d'écrire son code dans un éditeur de texte, nous permettant de l'éditer, de le sauvegarder, d'y revenir, sans à chaque fois avoir à remonter toutes les lignes. Je vous demanderai donc de suivre ce conseil. Surtout que si vous ne le faites pas, vous serez embêtés plus bas sur cette page !
Je voudrais que vous notiez avant de commencer que j'utiliserai la notation américaine du séparateur des décimales, donc le point et non la virgule.
Chargement du tableau de données
Nous allons commencer par créer une variable qui va contenir les données. Pour cela, nous allons charger dans R un tableau que je vous ai préparé : Table.txt
1 |
results = read.table ("Chemin_ou_vous_avez_sauve_Table/Table.txt", header = T); |
1 |
attach (results); |
Nous voyons ici plusieurs choses :
- comment lire un tableau ("read.table") ;
- comment sauver ces données dans une variable ("results =") ;
- comment dire à "read.table" que le tableau qu'on lui donne a ce qu'on appelle un header (un entête), c'est-à-dire que les colonnes ont un nom (la première ligne du fichier est composée des noms des colonnes).
[important]Quand on utilise la fonction "read.table", R s'attend à ce que le fichier qu'on lui donne ait les colonnes séparées par un espace ; il est possible de lui spécifier un autre séparateur avec l'argument "sep".[/important]
Enfin, "attach" sert à attacher les données contenues dans "results", dans le but de pouvoir appeler directement les colonnes de Table.txt, simplement en les nommant.
[notice]Essayez par exemple de taper "Conserved", il vous répondra les valeurs contenues dans la colonne "Conserved"[/notice]
1er barplot : La partie haute
Maintenant que nos valeurs sont bien accrochées, dessinons ! "barplot" est une fonction qui permet de faire... roulement de tambours... un barplot ! Soit en bon français, un graphique en barres. C'est exactement ce qui nous intéresse ici. Commençons donc par la ligne suivante :
1 |
barplot (Conserved, ylim = c(-1,1), main = "Bases", col = "yellow", axes = F); |
Ici, nous allons donc représenter les données de la colonne "Conserved" :
- en définissant les limites sur l'axe y (ylim) entre -1 et 1 ("c(-1,1)") ;
- en définissant un titre principal ("main") ;
- en colorant les barres ("col") en jaune ;
- le tout sans dessiner d'axes (axes = F) ;
Nous représentons ainsi la partie haute de notre premier graphe :
2ème barplot : La partie basse
Pour faire la partie basse, c'est aussi simple, il suffit de taper ceci :
1 |
barplot (-Unconserved, ylim = c(-1,1), col = "yellow4", add=T, axes = F); |
Si vous avez bien suivi, vous devez avoir :
- la partie du bas ("-Unconserved") ;
- dans les mêmes limites ;
- avec une couleur plus sombre ("yellow4") ;
- toujours sans les axes
[notice]Vous devez vous demander ce que c'est que ce "add=T". Eh bien, c'est un paramètre très utile que l'on peut trouver dans beaucoup de fonctions, qui permet très simplement d'ajouter ("add") un nouveau graphique à un premier, sans l'effacer ! C'est l'équivalent de la fonction points ou lines pour la fonction plot (je vous invite à taper "?lines", "?points" ou "?plot" si vous vous posez des questions). Cet "add" présente l'avantage de ne pas avoir à apprendre ou retenir le nom d'une autre fonction pour faire la même chose ![/notice]
Bon, trêve de bavardages, la suite !
Les axes
1 |
axis (2, label = c(1, 0.5, "0.0", 0.5, 1), at = c(-1, -.5, 0, .5, 1)) |
La fonction axis permet de dessiner un axe, que l'on n'avait pas dessiné avec barplot. Pourquoi faire ceci ? Essayez les lignes précédentes en mettant "axes = T", vous verrez que barplot vous dessine un axe y dont les valeurs vont de -1 à 1. L'intérêt étant ici non pas d'avoir un axe de -1 à 1, mais de 1 à 0, et de 0 à 1, car on représente 2 éléments qui sont inverses, en valeur absolue.
Le 2 est l'axe que l'on veut dessiner, en retenant le code suivant : 1 = bas, 2 = gauche, 3 = haut et 4 = droite.
[notice]Je vous invite à vous amuser un peu avec ces axes.[/notice]
"Label" est la liste de noms que l'on veut donner aux points affichés sur l'axe (une liste se définit en mettant des éléments séparés par une virgule dans "c()"), et "at" est la liste des positions respectives de ces points.
1 |
axis (1, label = c("FB CDS","Cel. CDS","FB 5'UTR","Cel. 5'UTR","FB 3'UTR","Cel. 3'UTR","FB Exons", |
1 |
"Cel. Exons","FB Introns","FB Pseudogenes","PolII binding sites","TF Sites","Histones", |
1 |
"CNV","Origins"), las = 3, at = 0.7+1.2*(0:14), pos = -1.1) |
Faisons de même pour l'axe des x, que l'on veut voir porter le nom des "Features" de Table.txt, mais de façon plus lisible :
- pour cela, on n'utilise pas directement Features, mais on redonne une liste (vous vous souvenez de "c()") ;
- on aurait aussi très bien pu définir d'abord le nom des Features et l'utiliser ;
- ensuite, on indique les positions de ces noms (tapez "0.7+1.2*(0:14)" dans R pour comprendre). Vous pouvez trouver les positions empiriquement en tâtonnant, mais sachez qu'en toute logique, il y a un espace avant les barres (ici le 0.7), puis que les barres sont également espacées (de 1.2), et qu'il y a 15 barres (on compte de 0 à 14 car la première est en position 0.7).
- enfin, "pos" détermine la position verticale de l'axe, que l'on veut ici en dessous du graphique du bas, donc placé à un peu plus de -1, par exemple -1.1.
Les barres d'erreurs
Plaçons maintenant les barres d'erreur, qui servent en général à représenter l'écart type, qui est en termes simples une indication de la dispersion des résultats. Vous pouvez cependant représenter ce que vous voulez, Variance ou autres. La fonction qui nous intéresse est "errbar", elle est contenue dans le package "Hmisc", qu'il vous faudra installer ("install.packages()") et appeler ("require()" ou "library()") :
1 2 3 4 |
install.packages("Hmisc") library(Hmisc) errbar (0.7+1.2*(0:14), Mean, yplus = Mean + 2*SD, yminus = Mean - 2*SD, add = T, pch = "") errbar (0.7+1.2*(0:14), -Mean2, yplus = -(Mean2 + 2*SD2), yminus = -(Mean2 - 2*SD2), add = T, pch = "") |
Nous donnons donc à errbar :
- en premier la position sur l'axe x où nous voulons nos barres d'erreurs ;
- ensuite les positions, qui sont les mêmes que celle des noms sur l'axe x ;
- puis la position du milieu de la barre ("Mean"), la position haute de la barre d'erreur ("yplus"), la position basse ("yminus") ;
- enfin on l'ajoute au graphique précédent ("add=T") ;
- et on utilise "pch" qui est le symbole utilisé pour le milieu de la barre (nous n'en voulons pas, donc on le met égal à un ensemble vide).
[error]Vous devez avoir un erreur qui apparaît quand vous tapez le code de la première barre d'erreur :
1 |
Erreur dans y[Type == 1] : objet de type 'closure' non indiçable |
Pour comprendre cette erreur, c'est très simple. Tapez dans votre console R :
1 |
Mean |
Il vous renvoie :
1 2 3 4 5 |
function (object, ...) UseMethod("Mean") <environment: namespace:Hmisc> |
Et pas la suite de moyennes attendu. Ceci est dû au fait que "Mean" est, en plus d'être notre liste de moyennes, une fonction interne à R et déjà utilisée. Donc R, quand on appelle "Mean", va tenter d'utiliser la fonction. Or on ne peut pas représenter graphiquement une fonction sans lui donner de données : on crée donc une erreur. Cette erreur nous permet de voir comment préciser que ce que l'on veut représenter appartient à un ensemble plus grand. Par exemple ici, notre ensemble est "results" duquel on veut représenter "Mean". Il faudra donc mettre dans notre commande non pas "Mean" mais "results$Mean", qui est en quelques sortes le chemin interne de notre liste.[/error]
[notice]Je vous invite à jouer avec "pch" en lui donnant des valeurs numériques, mais aussi des symboles ou des lettres ; pch = "+" est très utilisé.[/notice]
Les derniers détails
1 |
abline(h=0, lwd = 2) |
Ajoutons maintenant une ligne droite pour délimiter l'espace des deux graphiques, avec la fonction "abline" (qui signifie ligne entre le point A et le point B). Nous voulons qu'elle soit horizontale, au point 0 de l'axe des y et qu'elle ait une épaisseur ("lwd = line width", épaisseur de ligne) de 2. Enfin, ajoutons un peu de texte, avec la fonction "mtext" :
1 2 3 |
mtext("Fraction", side = 2, adj = 0.5, padj = -6) mtext("Conserved", side = 2, adj = 0.80, padj = -4) mtext("Unconserved", side = 2, adj = 0.2, padj = -4) |
On veut donc ajouter 3 mots, "Fraction", "Conserved" et "Unconserved", du côté gauche (2) du graphique :
- "adj" sert à placer le texte parallèlement à l'axe auquel le texte se rapporte, donc ici en hauteur
- "padj" place perpendiculairement à l'axe (d'où le p)
[important]Cette fois, cependant, on se place par rapport à l'extérieur du graphique et non plus par rapport aux valeurs de l'axe, car "mtext" place dans la marge (d'où le m). Alors là, vous allez me dire "Ouais mais moi ça sort de partout et ce n'est pas joli !"... et vous aurez raison. En effet, en écrivant des choses dans la marge, d'autant plus avec des valeurs inverses, sans avoir au préalable défini ces marges est dangereux, et vous fera bien souvent vous planter. [/important]
Nous définirons les marges ainsi :
1 |
par(oma = c(1,4,6,4), mar = c(9.1,4.1,1.1,4.1)); |
Ici, "oma" veut dire "outer margin", "mar" veut dire "margin", et les valeurs sont données toujours dans le même sens "bas, gauche, haut, droite". Maintenant, retapez votre code complet, ce qui devrait être facile si vous avez, au fur et à mesure de l'article, copié votre code dans votre éditeur de texte :
C'est-y pas beau tout ça ? Comment ? Vous voulez mettre ça dans un fichier .png ? Vous êtes bien exigeants... bon, eh bien c'est très simple, il suffit de taper ceci avant le graphe (personnellement, je le fais tout en haut) :
1 |
png ("chemin_ou_vous_voulez_le_sauver/nom.png", height = 1600, width = 1600, pointsize = 19); |
[important]Pour les valeurs de hauteur ("height") et largeur ("width"), par défaut elles sont données en pixels que je trouve être une bonne unité pour se rendre compte de la qualité et de la lisibilité d'un graphique. Mais vous pouvez bien sur choisir d'autres unités. Enfin, la taille des points du texte est à faire varier en fonction des gouts, en vérifiant bien à chaque changement la position des divers éléments texte, et en sachant qu'elle est calculée en fonction de la résolution de l'image ("?png" pour plus d'informations).[/important]
Alors, ça fonctionne ? Non ? Le fichier est vide ? Ah oui, il faut le fermer, sinon R est toujours entrain d'écrire dedans. Pour cela, on tape :
1 |
dev.off(); |
R doit alors vous indiquer qu'il est revenu dans le cadre de dessin précédent (quartz, x11 ou autres).
Enfin, vu qu'on a attaché "results" au début et qu'on veut prendre de bonnes habitudes, il faut détacher "results" :
1 |
detach(result); |
Si vous arrivez au même résultat, bravo ! Il ne vous reste plus qu'à vous attaquer à la dernière tâche ! Sinon, n'hésitez pas à poser des questions, je me ferai un plaisir d'y répondre !
À vous de jouer !
Votre défi, si vous l'acceptez, sera alors de faire comme sur la figure en haut de l'article, mettre 2 graphiques dans la même image, du texte en plus, et la tourner dans le bon sens. Pour cela 2 indices :
- "mfrow" est un paramètre très utile de "par " ;
- parfois, il vaut mieux se faciliter la vie et faire des choses simples avec un outil simple. Une rotation par exemple.
Le premier qui réussit aura gagné tout mon respect !
Allez, à la prochaine pour des choses encore plus folles !
Je tiens à remercier tous mes relecteurs, Hautbit, nallias, Guillaume Collet, Malicia, Yoann M. et nahoy !
Nasser
novembre 29, 2012 à 11:31
Vous êtes vraiment supers!!!! je suis sur qu'avec vous je serai un bioinformaticien avant d'avoir eu la licence! Je parlerai de vous sur mon blog!!!!!!!!!!
Merci!!!!!!!!!!!
Aurélien C.
novembre 29, 2012 à 11:44
Merci à toi Nasser, ça me fait plaisir de faire plaisir ! Et merci de parler de nous sur ton blog !
Nasser
novembre 30, 2012 à 4:52
Tous le plaisir est pour moi!
cortex2001
décembre 6, 2012 à 3:38
Super article 🙂
Si vous aviez pu être à la place de mon prof il y a 8 ans, j'aimerai p-e R au lieu de détester ce langage XD
raf
septembre 26, 2013 à 9:31
plus de tutorial depuis?
je suis biologiste et j'essai dans mes temps mort de me formé a R d'autre conseil?
site?
Yoann M.
septembre 26, 2013 à 9:38
Bonjour raf,
merci pour ton commentaire. Si tout va pour le mieux, le troisième épisode devrait arriver avant la fin de l'année. L'auteur habituel a eu un très gros planning ces derniers mois 🙂
Quant aux différents sites éventuels où tu pourrais trouver des informations, je lui laisse le soin de te répondre mieux que moi je ne le ferai puisque je ne manipule pas ce "langage" 🙂