silico.biotoul.fr
 

M1 Traitement de Donnees Biologiques - TP 2 R

From silico.biotoul.fr

(Difference between revisions)
Jump to: navigation, search
m (Lecture et exploration du tableau de données)
 
(60 intermediate revisions not shown)
Line 6: Line 6:
* Estimation d'un intervalle de confiance d'un paramètre (moyenne, proportion)
* Estimation d'un intervalle de confiance d'un paramètre (moyenne, proportion)
 +
<span style='color: #990000'>
 +
Créer un répertoire de travail sur le bureau (par exemple TDB-TP2) et commencez par télécharger le fichier source que vous allez utiliser et compléter pour générer le compte rendu de TP : [[Media:M1.TDB.TP_reglin_proba_ic_R.Rmd|M1.TDB.TP_reglin_proba_ic_R.Rmd]] (click droit de la souris -- enregistrer la cible sous...).
 +
</span>
 +
Ouvrez le logiciel RStudio et chargez ce fichier puis lancez sa compilation pour voir le compte rendu. Pour cela cliquez sur le bouton Knit HTML ou bien utilisez la combinaison de touches Ctrl + shift + K.
= Régression linéaire =
= Régression linéaire =
-
<span style='color: #990000'>On cherche à détecter une corrélation entre deux variables quantitatives mesurées sur les mêmes individus.</span>
+
On cherche à détecter une corrélation entre deux variables quantitatives mesurées sur les mêmes individus.
Notre exemple: pour différentes espèces de bactéries, nous connaissons (i) la taille du génome grâce au séquençage, (ii) le nombre de séquences codantes prédites par un algorithme bioinformatique.
Notre exemple: pour différentes espèces de bactéries, nous connaissons (i) la taille du génome grâce au séquençage, (ii) le nombre de séquences codantes prédites par un algorithme bioinformatique.
-
Quelle question peut-on se poser ?
+
Téléchargez et sauvegardez le fichier [[Media:bacterial_genomes.txt|bacterial_genomes.txt]] dans le répertoire crée pour ce TP (click droit de la souris -- enregistrer la cible sous...).
 +
 
 +
<span style='color: #990000'>
 +
Quelle(s) question(s) peut-on se poser ?
 +
</span>
== Lecture et exploration du tableau de données ==
== Lecture et exploration du tableau de données ==
Line 19: Line 27:
</source>
</source>
-
Accéder directement aux variables simplement en donnant leurs noms:
+
Accéder directement aux colonnes du data.frame <tt>genomes</tt> simplement à partir de leurs noms (fonction <tt>attach</tt>):
 +
<source lang='rsplus'>
attach(genomes)
attach(genomes)
names(genomes)
names(genomes)
Genome_size
Genome_size
ORF_number
ORF_number
 +
</source>
 +
 +
<span style='color: #990000'>
 +
Quelle est l'unité de mesure pour <tt>Genome_size</tt> ?
 +
</span>
Faire un graphique:
Faire un graphique:
 +
<source lang='rsplus'>
plot(Genome_size,ORF_number,pch=16)
plot(Genome_size,ORF_number,pch=16)
 +
</source>
Quantifier la relation entre ces 2 variables:
Quantifier la relation entre ces 2 variables:
* covariance
* covariance
-
cov(Genome_size,ORF_number)
+
<source lang='rsplus'>
-
* coefficient de corrélation r de Pearson
+
cov(Genome_size,ORF_number)
-
cor(Genome_size,ORF_number)
+
</source>
 +
* coefficient de corrélation ''r'' de Pearson
 +
<source lang='rsplus'>
 +
cor(Genome_size,ORF_number)
 +
</source>
-
== Régression linéaire ==
+
== Régression ==
-
Calculer une fonction linéaire qui relie les 2 variables, avec la commande <span style='color: #990000;'>lm()</span>:
+
<span style='color: #990000;'>
 +
Calculer une fonction linéaire qui relie les 2 variables, avec la commande <tt>lm()</tt>
 +
</span>:
 +
<source lang='rsplus'>
lm(ORF_number ~ Genome_size) # remarque: c'est une régression de "y" sur "x", d'où lm(y~x)
lm(ORF_number ~ Genome_size) # remarque: c'est une régression de "y" sur "x", d'où lm(y~x)
 +
</source>
Gardons en mémoire le résultat de la régression:
Gardons en mémoire le résultat de la régression:
 +
<source lang='rsplus'>
reglin=lm(ORF_number ~ Genome_size)
reglin=lm(ORF_number ~ Genome_size)
 +
</source>
On peut vérifier la significativité des coefficients de la droite avec:
On peut vérifier la significativité des coefficients de la droite avec:
 +
<source lang='rsplus'>
summary(reglin)
summary(reglin)
 +
</source>
 +
<span style='color: #990000;'>
Quelle est l'équation de la droite de régression ?
Quelle est l'équation de la droite de régression ?
 +
</span>
Calculez le coefficient de détermination R2 (% de variance expliquée par le modèle linéaire ==> bien si > 70%):
Calculez le coefficient de détermination R2 (% de variance expliquée par le modèle linéaire ==> bien si > 70%):
-
cor(Genome_size,ORF_number)^2 (stocké aussi dans summary(reglin))
+
<source lang='rsplus'>
-
+
cor(Genome_size,ORF_number)^2 #(stocké aussi dans summary(reglin))
 +
</source>
 +
 
Représenter le nuage de points avec la droite de régression:
Représenter le nuage de points avec la droite de régression:
 +
<source lang='rsplus'>
plot(Genome_size,ORF_number,pch=16)
plot(Genome_size,ORF_number,pch=16)
abline(reglin,col="red",lwd=2)
abline(reglin,col="red",lwd=2)
 +
</source>
<span style='color: #990000;'>Ajoutez ces parties à votre compte rendu.</span>
<span style='color: #990000;'>Ajoutez ces parties à votre compte rendu.</span>
 +
= Calculs de probabilité =
 +
== Calculs basés sur des données ==
 +
L'exemple porte sur la taille des arbres dans une forêt de séquoias.
 +
Téléchargez et sauvegardez le fichier [[Media:sequoia.txt|sequoia.txt]] (click droit de la souris -- enregistrer la cible sous...).
 +
<span style='color: #990000;'>
 +
Chargez le fichier avec la commande <tt>read.table</tt> dans une variable <tt>sequoia</tt>, puis utilisez la commande <tt>attach</tt> pour accéder directement aux noms des colonnes. Utilisez également la commande <tt>names</tt> pour connaître le nom des colonnes dans le fichier.
 +
</span>
-
 
+
<!--
-
 
+
Lecture du tableau de données:
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
[[Image:M1.TDB.RStudio-screenshot.jpg]]
+
-
 
+
-
 
+
-
 
+
-
<span style='color: #990000'>Créer un répertoire de travail sur le bureau (par exemple <tt>TDB-TP1_introduction</tt>) et commencez par télécharger le fichier source que vous allez utiliser et compléter pour générer le compte rendu de TP : [[Media:M1.TDB.TP_introduction_R.Rmd|M1.TDB.TP_introduction_R.Rmd]] (click droit de la souris -- enregistrer la cible sous...).</span>
+
-
 
+
-
<span style='color: #990000'>Ouvrez le logiciel RStudio et chargez ce fichier puis lancez sa compilation pour voir le rendu comme dans la capture d'écran ci-dessus. Pour cela cliquez sur le bouton '''Knit HTML''' ou bien utilisez la combinaison de touches <tt>Ctrl + shift + K</tt>.</span>
+
-
 
+
-
Vous verrez que si la compilation est réussie, un fichier <tt>M1.TDB.TP_introcution_R.html</tt> va être généré dans le même répertoire que le fichier <tt>M1.TDB.TP_introduction_R.Rmd</tt> que vous avez téléchargé.
+
-
 
+
-
= Utilisation et calculs avec du code R =
+
-
 
+
-
== Utilisation de variables ==
+
-
 
+
-
Les variables (ou objets) permettent de stocker des données qui peuvent être :
+
-
* une valeur simple de type numérique (<tt>numeric</tt>), logique (<tt>logical</tt>), chaîne de caractères (<tt>character</tt>) ou qualitative (<tt>factor</tt>).
+
-
* une liste (appelée <tt>vector</tt>).
+
-
* un tableau à 2 dimensions, les colonnes pouvant avoir des types différents (<tt>data frame</tt>). Ce sont les plus utilisés en statistiques.
+
-
* un tableau à 2 dimensions, toutes les cases ayant le même type (<tt>matrix</tt>).
+
-
* un tableau à plus de 2 dimension (<tt>array</tt>).
+
-
* une combinaison des précédents (<tt>list</tt>).
+
-
 
+
-
Une variable a un nom (défini par l'utilisateur) qui permet d'accéder à son contenu.
+
-
 
+
-
<span style='color: #990000;'>Créer une variable appelée <tt>x</tt> qui stocke le résultat du calcul précédent <tt>1 + 2 + 3</tt>.</span> L'affectation d'une valeur à une variable se fait avec <tt><-</tt> ou le signe <tt>=</tt> comme suit dans la console :
+
<source lang='rsplus'>
<source lang='rsplus'>
-
x = 1 + 2 + 3
+
sequoia=read.table("sequoia.txt", sep="\t", header=TRUE)
 +
attach(sequoia);names(sequoia)
</source>
</source>
 +
-->
-
Vous devriez voir s'afficher <tt>x</tt> dans l'environnement en haut à droite comme sur la capture d'écran précédente.
+
Représentation graphique des données:  
-
 
+
-
Pour afficher le contenu d'une variable, il suffit de taper son nom (dans la console en bas à gauche) :
+
<source lang='rsplus'>
<source lang='rsplus'>
-
x
+
hist(taille_arbre)
</source>
</source>
-
<span style='color: #990000;'>Rajoutez ces commandes à votre compte rendu de TP (script en haut à gauche) et assurez-vous que cela fonctionne (en recompilant le fichier).</span>
+
<span style='color: #990000;'>
 +
Quelle est la probabilité qu'un arbre mesure 80m ?
 +
</span>
-
<span style='color: #990000;'>Essayez avec les autres types de variables simples. Par exemple :</span>
 
<source lang='rsplus'>
<source lang='rsplus'>
-
texte = "Phrase qu'il faut mettre entre guillemets."
+
length(taille_arbre[taille_arbre==80])/length(taille_arbre)
-
logique = TRUE
+
-
faux = FALSE
+
-
# ce qui est après un # correspond à du commentaire et n'est pas interprété comme du code
+
-
x = x * 2 # on multiplie x de toute à l'heure par 2
+
-
y = x**2  # y prend la valeur de x au carré (puissance 2)  
+
</source>
</source>
-
 
+
ou mieux :  
-
Les valeurs logiques peuvent être le résultats de tests :
+
<source lang='rsplus'>
<source lang='rsplus'>
-
x < 10  # strictement inférieur à 10 ?
+
length( which(taille_arbre==80) ) / length(taille_arbre)
-
x <= 11 # inférieur ou égal
+
-
x == 12 # test d'égalité
+
-
x >= 13
+
-
x > 0
+
-
x != 0 # test si x est différent de 0
+
-
! FALSE  # inverse une valeur logique
+
</source>
</source>
-
<span style='color: #990000;'>Ajoutez ces parties à votre compte rendu.</span>
 
-
== Vecteurs ==
+
<span style='color: #990000;'>
-
<span style='color: #990000;'>Créez un vecteur <tt>x</tt> contenant les valeurs 2,3,5,8,4,6 rassemblées avec la fonction <tt>c()</tt> :</span>
+
Quelle est la probabilité qu'un arbre mesure plus de 100 ?
-
<source lang='rsplus'>
+
</span>
-
x = c(2,3,5,8,4,6)
+
-
x
+
-
</source>
+
-
<span style='color: #990000;'>Calculez sa longueur avec la fonction <tt>length()</tt> :</span>
 
<source lang='rsplus'>
<source lang='rsplus'>
-
length(x)
+
length( which(taille_arbre>100) ) / length(taille_arbre)
</source>
</source>
-
Accès à certains éléments :
+
== Calculs basés sur des lois (ou distributions) de probabilité  ==
-
<span style='color: #990000;'>Affichez la deuxième valeur de x.</span> Pour cela, on utilise les crochets :
+
-
<source lang='rsplus'>
+
-
x[2]
+
-
</source>
+
-
Pour accéder à plusieurs éléments, on utilise un vecteur contenant les positions :
+
Différentes lois sont programmées (binomiale, Poisson, normale, chi2 ou <math>\chi^2</math>, ...)
-
<source lang='rsplus'>
+
Il y a plusieurs fonctions pour chaque loi. Par exemple, pour la loi normale:
-
x[ c(2,4,1) ]
+
* <tt>dnorm()</tt> : fonction densité ('''d'''ensity)
-
</source>
+
* <tt>pnorm()</tt> : fonction de répartition ('''p'''robability)
 +
* <tt>qnorm()</tt> : fonction quantile ('''q'''uantile)
 +
* <tt>rnorm()</tt> : générateur aléatoire ('''r'''andom)
-
Si les positions se suivent, on peut utiliser les <tt>:</tt>
+
=== Calculs basés sur une loi normale centrée réduite <math>N(0,1)</math> ===
-
<source lang='rsplus'>
+
-
2:4
+
-
x[2:4]
+
-
</source>
+
-
ou dans l'ordre inverse :
+
'''Rappels :'''
-
<source lang='rsplus'>
+
* Une '''loi normale''' est définie par '''2 paramètres''' : la '''moyenne''' et l''''écart-type'''.
-
4:2
+
* La '''loi normale centrée réduite''' est la loi normale dont la moyenne est '''0''' et l'écart-type '''1'''. Elle est notée généralement <math>N(0,1)</math> (0 pour la moyenne et 1 pour l'écart-type de la variable X).
-
x[4:2]
+
-
</source>
+
-
On peut aussi accéder à tous les éléments sauf certains :
+
[[Image:loi_normale.jpeg|thumb|400px|link=|Fonction de densité de la loi normale centrée réduite <math>N(0,1)</math>]]
-
<source lang='rsplus'>
+
-
x[ -2 ] # tous, sauf le 2ème
+
-
x[ -2:-3 ] # tous, sauf le 2ème et le 3ème
+
-
</source>
+
-
Exploration suivant un critère: vrai/faux, indices, ou valeurs
+
Probabilité de X = 0, notée <math>P(X = 0)</math>
<source lang='rsplus'>
<source lang='rsplus'>
-
x>4
+
dnorm(0,0,1)
-
which( x>4 )
+
-
x[ x>4 ]
+
</source>
</source>
-
Effectuer une opération sur chaque élément :
+
<math>P(X \le 0)</math> :
<source lang='rsplus'>
<source lang='rsplus'>
-
x
+
pnorm(0,0,1)
-
1 / x
+
-
# en arrondissant 2 chiffres après la virgule)
+
-
round(1/x, 2)
+
-
x + 1
+
-
x * 2
+
</source>
</source>
-
<span style='color: #990000;'>Ajoutez ces parties à votre compte rendu.</span>
+
<math>P(X > 0)</math> :
-
 
+
-
== Fonctions ==
+
-
 
+
-
R charge certaines librairies au démarrage selon comment il est installé et configuré. Les librairies mettent à disposition des fonctions qui peuvent être très simples, jusqu'à très très compliquées.
+
-
 
+
-
<span style='color: #990000;'>Essayez les suivantes sur votre vecteur <tt>x</tt>. Et notez dans votre compte rendu ce qu'elles permettent de faire.</span>
+
-
* <tt>summary(x)</tt>
+
-
* <tt>min</tt>
+
-
* <tt>max</tt>
+
-
* <tt>median</tt>
+
-
* <tt>mean</tt>
+
-
* <tt>var</tt>
+
-
* <tt>sd</tt>
+
-
 
+
-
Afin de mieux comprendre ce qu'elles font et quels sont les paramètres que l'on peut leur passer, <span style='color: #990000;'>utilisez l'aide :</span>
+
<source lang='rsplus'>
<source lang='rsplus'>
-
help(median)
+
1-pnorm(0,0,1) # ou pnorm(0,0,1, lower.tail=F)
</source>
</source>
-
ou plus simplement :
+
Histogramme pour 100 valeurs échantillonnées aléatoirement suivant une <math>N(0,1)</math> :
<source lang='rsplus'>
<source lang='rsplus'>
-
?mean
+
hist(rnorm(100,0,1))
</source>
</source>
-
Quand on ne connaît pas exactement le nom, on peut faire une recherche par mot :
 
-
<source lang='rsplus'>
 
-
??mean
 
-
</source>
 
-
== Lecture et écriture de fichiers contenant des données ==
 
-
Les fonctions permettant de lire et d'écrire des données dans des fichiers sont :
+
'''Application à la forêt de séquoias'''
-
* '''read.table'''
+
-
* read.csv et read.csv2
+
-
* read.delim et read.delim2
+
-
* '''write.table'''
+
-
* write.csv et write.csv2
+
-
* write.delim et write.delim2
+
-
Ce sont à peu près les mêmes fonctions mais elles prennent des paramètres par défaut légèrement différents. <span style='color: #990000;'>Commencez par consulter leur documentation.</span>
+
<span style='color: #990000;'>
 +
Utilisez la moyenne et l'écart-type des tailles des séquoias chargées précédemment pour calculer la probabilité qu'un arbre mesure 80m selon la loi normale paramétrée ainsi.
 +
</span>
-
<span style='color: #990000;'>Téléchargez et sauvegardez le fichier suivant dans le même répertoire que toute à l'heure.</span> (click droit puis enregistrer la cible sous...) : [[Media:croissance_plantes.txt|croissance_plantes.txt]]
+
<span style='color: #990000;'>
 +
Quelle est la probabilité qu'un séquoia mesure plus de 110m ?
 +
</span>
-
Afin de visualiser son contenu, vous pouvez l'ouvrir d'abord avec un éditeur de texte (Textpad, Notepad++ ou geany), ou un tableur LibreOffice Calc (ou M office excel).
+
<!--span style='color: #990000;'>
 +
Quelle est la probabilité qu'un arbre pris au hasard dans cette forêt et qui mesure 70m soit un séquioa  ?
 +
</span-->
-
'''Remarque :''' il faut en général éviter de mettre des espaces dans les noms de variables, ainsi que des accents
+
=== Calculs basés sur une loi Binomiale <math>B(n,p)</math> ===
-
Ensuite, avant de charger les données, il faut '''dire à R de se placer dans votre répertoire de travail'''. Sinon, il ne trouvera pas le fichier et vous aurez une erreur à la suite de la commande <tt>read.table</tt>. Pour cela, <span style='color: #990000;'>cliquez sur le menu <tt>Session</tt> puis <tt>Set Working Directory > To Source File Location</tt>.</span> Il y aussi la possibilité de sélectionner le répertoire avec l'option <tt>Choose Directory...</tt> ou la combinaison de touches <tt>Ctrl + shift + H</tt>.
+
Exemple: pour un lot de n=100 graines d'Arabidopsis ayant chacune une probabilité de germination p=0.8 donc selon <math>B(100,0.8)</math>
-
Après cela, <span style='color: #990000;'>essayez la fonction <tt>read.table</tt> sur le fichier <tt>croissance_plantes.txt</tt> que vous avez récupéré.</span> Il vous faudra utiliser les paramètres <tt>sep</tt> et <tt>header</tt> pour désigner le séparateur de colonnes (tabulation notée <tt>\t</tt>) ainsi que le fait que la première ligne correspond aux noms des colonnes (<tt>header=TRUE</tt>):
+
La probabilité que k=80 graines germent (k=nombre de succès)
<source lang='rsplus'>
<source lang='rsplus'>
-
read.table("croissance_plantes.txt", sep="\t",header=TRUE)  
+
dbinom(80,100,prob=0.8)
</source>
</source>
-
 
+
-
Avec la commande ci-dessus, si le fichier est dans le répertoire de travail, il va être chargé. Comme il n'est pas affecté à une variable, il va être affiché. Pour le garder en mémoire dans une variable, il faut le spécifier :
+
La probabilité que au maximum 80 graines germent
<source lang='rsplus'>
<source lang='rsplus'>
-
croissance = read.table("croissance_plantes.txt", sep="\t",header=TRUE)  
+
pbinom(80,100,prob=0.8)
</source>
</source>
-
Après cela, la variable <tt>croissance</tt> doit apparaître dans l'environnement. Si vous cliquez dessus, son contenu doit être ouvert dans un nouvel onglet. Vous verrez aussi dans la console la commande qui a été exécutée pour l'ouvrir : <tt>View(croissance)</tt>.
+
La probabilité que plus de 80 graines germent
-
 
+
-
Ce type de tableau est appelé <tt>data.frame</tt> en R. Chaque colonne a un type. C'est ce qui est le plus utilisé en statistiques avec les individus en lignes et les variables en colonnes. <span style='color: #990000;'>Essayez la fonction <tt>summary</tt> ainsi que la fonction <tt>dim</tt>.</span>
+
-
 
+
-
Comme pour un vecteur, on peut afficher le contenu entier, ou accéder à certaines lignes et/ou colonnes :
+
<source lang='rsplus'>
<source lang='rsplus'>
-
croissance
+
pbinom(80,100,prob=0.8, lower.tail=F)
-
croissance[1 , ] # ligne 1
+
-
croissance[  , 1] # colonne 1
+
-
croissance[1, 1] # première case
+
-
croissance$poids # colonne qui s'appelle poids (la première)
+
-
croissance[ , 'poids'] # encore une autre manière d'y accéder
+
</source>
</source>
-
Pour obtenir le nom des colonnes :
+
<span style='color: #990000'>
 +
En moyenne, combien de graines devraient germer ?
 +
</span>
 +
 
 +
La distribution Binomiale <math>B(100,0.8)</math>:
<source lang='rsplus'>
<source lang='rsplus'>
-
names(croissance)
+
k=seq(0,100,1)    # k=nombre de graines qui germent
 +
plot(dbinom(k,100,prob=0.8),xlab="k",ylab="Prob(X=k)")
</source>
</source>
-
Il est aussi possible d'utiliser <tt>colnames</tt> et <tt>rownames</tt>.
 
-
Pour y accéder sans avoir à chaque fois à mettre le nom <tt>croissance</tt>, on utilise <tt>attach</tt> :
 
-
<source lang='rsplus'>
 
-
attach(croissance)
 
-
poids
 
-
</source>
 
-
Pour récupérer une partie des données, on peut faire comme précédemment avec un vecteur. Par exemple pour ne récupérer que les lignes correspondant à des données d'origine des pyrénées :
+
<span style='color: #990000;'>
-
<source lang='rsplus'>
+
Ajoutez ces parties à votre compte rendu.
-
origine_geo == 'pyr'    # test pour obtenir VRAI/FAUX
+
</span>
-
which(origine_geo == 'pyr') # test pour obtenir les numéros des lignes
+
-
# et ensuite
+
-
croissance[ origine_geo == 'pyr' , 1:2 ]
+
-
# ou bien
+
-
croissance[ which(origine_geo == 'pyr') , 1:2 ]
+
-
# pour stocker le résultats dans une autre variable
+
-
pyr = croissance[ origine_geo == 'pyr' , 1:2 ]
+
-
</source>
+
-
Ensuite, <span style='color: #990000;'>sauvegardez cette partie des données avec la fonction <tt>write.table</tt> et ouvrez le résultat avec Textpad ou LibreOffice :</span>
+
= Estimer un intervalle de confiance (IC) d'un paramètre populationnel (µ, p,...) à partir d'un échantillon =
-
<source lang='rsplus'>
+
-
write.table(pyr, "croissance_plantes_pyr.txt", quote=F, col.names=T, row.names=F, sep="\t")  
+
-
</source>
+
-
Vous remarquez au passage que l'on peut écrire TRUE avec seulement T, (et FALSE avec F).
+
-
== Fonctions graphiques ==
+
== Intervalle de confiance de la moyenne théorique µ de la population, à partir de la moyenne m calculée sur un échantillon ==
-
Pour commencer, un '''camembert''' avec le nombre d'individus par origine ('''variable qualitative''') :
+
Calculer l'intervalle de confiance (IC) du poids moyen théorique d'une récolte de tomates cerises, à partir d'un échantillon.
 +
<span style='color: #990000;'>
 +
Téléchargez et sauvegardez le fichier le fichier [[Media:tomates_cerises.txt|tomates_cerises.txt]] (click droit de la souris -- enregistrer la cible sous...).
 +
</span>
 +
 
 +
<span style='color: #990000;'>
 +
Chargez le fichier avec la commande <tt>read.table</tt> dans une variable <tt>tomates</tt>, puis utilisez la commande <tt>attach</tt> pour accéder directement aux noms des colonnes. Utilisez également la commande <tt>names</tt> pour connaître le nom des colonnes dans le fichier. Affichez le contenu de la colonne <tt>poids_tomate</tt>.
 +
</span>
 +
<!--
 +
Lecture du tableau de données:
<source lang='rsplus'>
<source lang='rsplus'>
-
summary(origine_geo)                                                   # les effectifs
+
tomates=read.table("tomates_cerises.txt", sep="\t", header=TRUE)
-
pie( summary(origine_geo) )                                            # le camembert
+
attach(tomates);names(tomates)
-
pie (summary(origine_geo), main="origine géographique des plantes")  # avec un titre
+
poids_tomate
-
</source>
+
</source> -->
-
Un '''histogramme''' ('''variable quantitative'''):
+
Script pour calculer l'IC:
<source lang='rsplus'>
<source lang='rsplus'>
-
hist(taille, xlim=c(40,90), xlab="taille (cm)", ylab="effectif", freq=T, main="histogramme de la taille des plantes", col="orange")
+
m = mean(poids_tomate)  # moyenne calculée sur l'échantillon
 +
s = sd(poids_tomate)   # écart-type calculé sur l'échantillon
 +
n = length(poids_tomate)# taille de l'échantillon
</source>
</source>
-
Une '''boite à moustaches''' ('''variable quantitative'''):
+
Si n&le;100, on utilise la distribution de Student pour calculer la marge d'erreur. Ici, nous allons calculer l'intervalle ayant 95% de chance de contenir la vraie moyenne de la population (et pas celle de l'échantillon).
<source lang='rsplus'>
<source lang='rsplus'>
-
boxplot(taille, main="boxplot de la taille des plantes", ylab="taille")  
+
error <- qt(0.975,df=n-1)*s/sqrt(n) # qt(0.975) = quantile 97,5% de la distribution de Student ; s/sqrt(n) = erreur standard de la moyenne
 +
left  <- m-error
 +
right <- m+error
</source>
</source>
-
 
+
L'intervalle ayant 95% de chance de contenir la vraie moyenne µ de la population est:
-
Plusieurs boites à moustaches (une variable quantitative en fonction des modalités d'une variable qualitative):
+
<source lang='rsplus'>
<source lang='rsplus'>
-
plot(origine_geo, taille, las=3, main="boxplot de la taille des plantes en fonction de l'origine géographique")
+
left;right
-
stripchart(taille~origine_geo,las=1)
+
</source>
</source>
-
Un graphique à 2 dimensions ou '''nuage de points''' (''' 2 variables quantitatives mesurées sur les mêmes individus''') :
+
Si n>100 on utilise la distribution normale pour calculer l'erreur.
<source lang='rsplus'>
<source lang='rsplus'>
-
plot( taille, poids)
+
error <- qnorm(0.975)*s/sqrt(n) # qnorm(0.975) = quantile 97,5% de la distribution normale
 +
left  <- m-error
 +
right <- m+error
</source>
</source>
-
 
+
L'intervalle ayant 95% de chance de contenir la vraie moyenne µ de la population est alors:
-
Affichage de plusieurs graphiques dans la même fenêtre :
+
<source lang='rsplus'>
<source lang='rsplus'>
-
par(mfrow=c(2,2))  # 2 en lignes et 2 en colonnes
+
left;right
-
hist(taille,xlim=c(40,90),xlab="taille (cm)",ylab="effectif",freq=T,main="histogramme de la taille des plantes",col="orange")
+
-
boxplot(taille,main="boxplot de la taille des plantes",ylab="taille")
+
-
plot(taille~origine_geo,las=3)
+
-
stripchart(taille~origine_geo,las=1)
+
</source>
</source>
-
Ouverture d'une nouvelle fenêtre graphique :
+
'''REMARQUE IMPORTANTE :''' si l'on veut un IC à 99% ==> remplacer 0.975 par 0.995.
 +
 
 +
== IC de la proportion théorique p d'un caractère dans la population, à partir de la fréquence f calculée sur un échantillon ==
 +
 
 +
Calcul de l'IC de la proportion théorique de pommes rouges dans une récolte comportant des pommes rouges et vertes, à partir d'un échantillon de n=125 pommes contenant f=0.4 de pommes rouges (40%).
 +
On utilise une approximation normale de la distribution d'une proportion :
<source lang='rsplus'>
<source lang='rsplus'>
-
x11()
+
f <- 0.4 # fréquence de pommes rouges dans l'échantillon
 +
n <- 125 # taille de l'échantillon
 +
error <- qnorm(0.975)*sqrt(f*(1-f)/n)
 +
left  <- f-error
 +
right <- f+error
 +
</source>
 +
L'intervalle ayant 95% de chance de contenir la vraie proportion p de la population est:
 +
<source lang='rsplus'>
 +
left;right
</source>
</source>
-
Sauvegarde et/ou exportation des graphiques. Dans l'onglet <tt>Plots</tt>, au choix dans le menu <tt>Export</tt> :
+
<span style='color: #990000;'>
-
* <tt>Save as Image...</tt> : différents format disponibles (PNG, JPEG, TIFF, SVG, ...)
+
Faites varier la taille de l'échantillon ainsi que la confiance.
-
* <tt>Save as PDF...</tt>
+
</span>
-
* <tt>Copy to Clipboard...</tt> : copie en mémoire (presse papier). Essayez de le coller par exemple dans MsWord, Powerpoint ou LibreOffice.
+
-
<span style='color: #990000;'>Ajoutez tout cela au compte rendu de TP avant de l'envoyer à votre enseignant par mail ([mailto:bonhomme@lrsv.ups-tlse.fr bonhomme@lrsv.ups-tlse.fr] ou [mailto:barriot@biotoul.fr barriot@biotoul.fr]). Le compte rendu est à envoyer '''avant de commencer le TP2'''. Envoyez les 2 fichiers (.Rmd et .html). Envoyez-vous aussi le mail en copie pour pouvoir vérifier que tout est bien passé. Mettez un titre tel que "Compte rendu TP1 TDB de -et votre Nom et Prénom-".</span>
+
 
 +
<span style='color: #990000;'>Ajoutez tout cela au compte rendu de TP avant de l'envoyer à votre enseignant par mail ([mailto:bonhomme@lrsv.ups-tlse.fr bonhomme@lrsv.ups-tlse.fr] ou [mailto:roland.barriot@univ-tlse3.fr roland.barriot@univ-tlse3.fr]). Le compte rendu est à envoyer '''avant de commencer le TP3'''. Envoyez les 2 fichiers (.Rmd et .html). Envoyez-vous aussi le mail en copie pour pouvoir vérifier que tout est bien passé. Mettez un titre tel que "Compte rendu TP2 TDB de -et votre Nom et Prénom-".</span>
= Liens =
= Liens =
Line 359: Line 302:
* Site de R : http://www.r-project.org et sites miroirs (dont ceux en France) pour télécharger le logiciel et les librairies : https://cran.r-project.org/mirrors.html
* Site de R : http://www.r-project.org et sites miroirs (dont ceux en France) pour télécharger le logiciel et les librairies : https://cran.r-project.org/mirrors.html
* RStudio : https://www.rstudio.com
* RStudio : https://www.rstudio.com
 +
* Utilisation de R depuis un navigateur : http://www.r-fiddle.org
 +
Chargement des données avec l'adresse des fichiers
 +
<source lang='rsplus'>
 +
genomes=read.table("http://silico.biotoul.fr/site/images/d/de/Bacterial_genomes.txt", sep="\t", header=TRUE)
 +
sequoia=read.table("http://silico.biotoul.fr/site/images/e/e1/Sequoia.txt", sep="\t", header=TRUE)
 +
tomates=read.table("http://silico.biotoul.fr/site/images/3/3d/Tomates_cerises.txt", sep="\t", header=TRUE)
 +
</source>

Current revision as of 08:10, 21 September 2021

Contents

Régression linéaire, probabilités et intervalles de confiance

Ce TP comporte trois parties:

  • Régression linéaire
  • Calculs simples de probabilité
  • Estimation d'un intervalle de confiance d'un paramètre (moyenne, proportion)

Créer un répertoire de travail sur le bureau (par exemple TDB-TP2) et commencez par télécharger le fichier source que vous allez utiliser et compléter pour générer le compte rendu de TP : M1.TDB.TP_reglin_proba_ic_R.Rmd (click droit de la souris -- enregistrer la cible sous...). Ouvrez le logiciel RStudio et chargez ce fichier puis lancez sa compilation pour voir le compte rendu. Pour cela cliquez sur le bouton Knit HTML ou bien utilisez la combinaison de touches Ctrl + shift + K.

Régression linéaire

On cherche à détecter une corrélation entre deux variables quantitatives mesurées sur les mêmes individus. Notre exemple: pour différentes espèces de bactéries, nous connaissons (i) la taille du génome grâce au séquençage, (ii) le nombre de séquences codantes prédites par un algorithme bioinformatique. Téléchargez et sauvegardez le fichier bacterial_genomes.txt dans le répertoire crée pour ce TP (click droit de la souris -- enregistrer la cible sous...).

Quelle(s) question(s) peut-on se poser ?

Lecture et exploration du tableau de données

Lecture:

genomes=read.table("bacterial_genomes.txt", sep="\t", header=TRUE)

Accéder directement aux colonnes du data.frame genomes simplement à partir de leurs noms (fonction attach):

attach(genomes)
names(genomes)
Genome_size
ORF_number

Quelle est l'unité de mesure pour Genome_size ?

Faire un graphique:

plot(Genome_size,ORF_number,pch=16)

Quantifier la relation entre ces 2 variables:

  • covariance
cov(Genome_size,ORF_number)
  • coefficient de corrélation r de Pearson
cor(Genome_size,ORF_number)

Régression

Calculer une fonction linéaire qui relie les 2 variables, avec la commande lm() :

lm(ORF_number ~ Genome_size) # remarque: c'est une régression de "y" sur "x", d'où lm(y~x)

Gardons en mémoire le résultat de la régression:

reglin=lm(ORF_number ~ Genome_size)

On peut vérifier la significativité des coefficients de la droite avec:

summary(reglin)

Quelle est l'équation de la droite de régression ?

Calculez le coefficient de détermination R2 (% de variance expliquée par le modèle linéaire ==> bien si > 70%):

cor(Genome_size,ORF_number)^2 #(stocké aussi dans summary(reglin))

Représenter le nuage de points avec la droite de régression:

plot(Genome_size,ORF_number,pch=16)			
abline(reglin,col="red",lwd=2)

Ajoutez ces parties à votre compte rendu.


Calculs de probabilité

Calculs basés sur des données

L'exemple porte sur la taille des arbres dans une forêt de séquoias. Téléchargez et sauvegardez le fichier sequoia.txt (click droit de la souris -- enregistrer la cible sous...).

Chargez le fichier avec la commande read.table dans une variable sequoia, puis utilisez la commande attach pour accéder directement aux noms des colonnes. Utilisez également la commande names pour connaître le nom des colonnes dans le fichier.


Représentation graphique des données:

hist(taille_arbre)

Quelle est la probabilité qu'un arbre mesure 80m ?

length(taille_arbre[taille_arbre==80])/length(taille_arbre)

ou mieux :

length( which(taille_arbre==80) ) / length(taille_arbre)


Quelle est la probabilité qu'un arbre mesure plus de 100 ?

length( which(taille_arbre>100) ) / length(taille_arbre)

Calculs basés sur des lois (ou distributions) de probabilité

Différentes lois sont programmées (binomiale, Poisson, normale, chi2 ou χ2, ...) Il y a plusieurs fonctions pour chaque loi. Par exemple, pour la loi normale:

  • dnorm() : fonction densité (density)
  • pnorm() : fonction de répartition (probability)
  • qnorm() : fonction quantile (quantile)
  • rnorm() : générateur aléatoire (random)

Calculs basés sur une loi normale centrée réduite N(0,1)

Rappels :

  • Une loi normale est définie par 2 paramètres : la moyenne et l'écart-type.
  • La loi normale centrée réduite est la loi normale dont la moyenne est 0 et l'écart-type 1. Elle est notée généralement N(0,1) (0 pour la moyenne et 1 pour l'écart-type de la variable X).
Fonction de densité de la loi normale centrée réduite N(0,1)

Probabilité de X = 0, notée P(X = 0)

dnorm(0,0,1)

P(X \le 0) :

pnorm(0,0,1)

P(X > 0) :

1-pnorm(0,0,1) # ou pnorm(0,0,1, lower.tail=F)

Histogramme pour 100 valeurs échantillonnées aléatoirement suivant une N(0,1) :

hist(rnorm(100,0,1))


Application à la forêt de séquoias

Utilisez la moyenne et l'écart-type des tailles des séquoias chargées précédemment pour calculer la probabilité qu'un arbre mesure 80m selon la loi normale paramétrée ainsi.

Quelle est la probabilité qu'un séquoia mesure plus de 110m ?


Calculs basés sur une loi Binomiale B(n,p)

Exemple: pour un lot de n=100 graines d'Arabidopsis ayant chacune une probabilité de germination p=0.8 donc selon B(100,0.8)

La probabilité que k=80 graines germent (k=nombre de succès)

dbinom(80,100,prob=0.8)

La probabilité que au maximum 80 graines germent

pbinom(80,100,prob=0.8)

La probabilité que plus de 80 graines germent

pbinom(80,100,prob=0.8, lower.tail=F)

En moyenne, combien de graines devraient germer ?

La distribution Binomiale B(100,0.8):

k=seq(0,100,1)    # k=nombre de graines qui germent
plot(dbinom(k,100,prob=0.8),xlab="k",ylab="Prob(X=k)")


Ajoutez ces parties à votre compte rendu.

Estimer un intervalle de confiance (IC) d'un paramètre populationnel (µ, p,...) à partir d'un échantillon

Intervalle de confiance de la moyenne théorique µ de la population, à partir de la moyenne m calculée sur un échantillon

Calculer l'intervalle de confiance (IC) du poids moyen théorique d'une récolte de tomates cerises, à partir d'un échantillon. Téléchargez et sauvegardez le fichier le fichier tomates_cerises.txt (click droit de la souris -- enregistrer la cible sous...).

Chargez le fichier avec la commande read.table dans une variable tomates, puis utilisez la commande attach pour accéder directement aux noms des colonnes. Utilisez également la commande names pour connaître le nom des colonnes dans le fichier. Affichez le contenu de la colonne poids_tomate.

Script pour calculer l'IC:

m = mean(poids_tomate)  # moyenne calculée sur l'échantillon
s = sd(poids_tomate)    # écart-type calculé sur l'échantillon
n = length(poids_tomate)# taille de l'échantillon

Si n≤100, on utilise la distribution de Student pour calculer la marge d'erreur. Ici, nous allons calculer l'intervalle ayant 95% de chance de contenir la vraie moyenne de la population (et pas celle de l'échantillon).

error <- qt(0.975,df=n-1)*s/sqrt(n) # qt(0.975) = quantile 97,5% de la distribution de Student ; s/sqrt(n) = erreur standard de la moyenne
left  <- m-error
right <- m+error

L'intervalle ayant 95% de chance de contenir la vraie moyenne µ de la population est:

left;right

Si n>100 on utilise la distribution normale pour calculer l'erreur.

error <- qnorm(0.975)*s/sqrt(n) # qnorm(0.975) = quantile 97,5% de la distribution normale
left  <- m-error
right <- m+error

L'intervalle ayant 95% de chance de contenir la vraie moyenne µ de la population est alors:

left;right

REMARQUE IMPORTANTE : si l'on veut un IC à 99% ==> remplacer 0.975 par 0.995.

IC de la proportion théorique p d'un caractère dans la population, à partir de la fréquence f calculée sur un échantillon

Calcul de l'IC de la proportion théorique de pommes rouges dans une récolte comportant des pommes rouges et vertes, à partir d'un échantillon de n=125 pommes contenant f=0.4 de pommes rouges (40%). On utilise une approximation normale de la distribution d'une proportion :

f <- 0.4	# fréquence de pommes rouges dans l'échantillon
n <- 125	# taille de l'échantillon
error <- qnorm(0.975)*sqrt(f*(1-f)/n)
left  <- f-error
right <- f+error

L'intervalle ayant 95% de chance de contenir la vraie proportion p de la population est:

left;right

Faites varier la taille de l'échantillon ainsi que la confiance.


Ajoutez tout cela au compte rendu de TP avant de l'envoyer à votre enseignant par mail (bonhomme@lrsv.ups-tlse.fr ou roland.barriot@univ-tlse3.fr). Le compte rendu est à envoyer avant de commencer le TP3. Envoyez les 2 fichiers (.Rmd et .html). Envoyez-vous aussi le mail en copie pour pouvoir vérifier que tout est bien passé. Mettez un titre tel que "Compte rendu TP2 TDB de -et votre Nom et Prénom-".

Liens

Chargement des données avec l'adresse des fichiers

genomes=read.table("http://silico.biotoul.fr/site/images/d/de/Bacterial_genomes.txt", sep="\t", header=TRUE)
sequoia=read.table("http://silico.biotoul.fr/site/images/e/e1/Sequoia.txt", sep="\t", header=TRUE)
tomates=read.table("http://silico.biotoul.fr/site/images/3/3d/Tomates_cerises.txt", sep="\t", header=TRUE)