|
Partie Statistiques du cours de BioInformatiqueMaster BTV, UFR Sciences - Université d'AngersSolutions du TD numéro 4 (énoncés)
|
Retour à la page principale de
(gH)
Partie Statistiques du cours de BioInformatique
Master BTV, UFR Sciences - Université d'Angers
Solutions du TD numéro 4 (énoncés)
Au lieu de calculer "bêtement" mean(V1), sd(V1), mean(V2)... nous profitons des fonctions de R pour écrire des instructions plus concises et matricielles :
md <- cbind(V1,V2,V3,V4,V5) fr <- function(x) { sprintf("%8.3f",x) } mr <- rbind( fr(colMeans(md)), fr(sd(md))) # #### ou mr <- apply( rbind( colMeans(md), sd(md) ) , 2 , fr ) row.names(mr) <- c("Moyenne","Ecart-type") colnames(mr) <- colnames(md) print(mr,quote=FALSE)L'exécution de ces instructions montre que V1 a une moyenne et un écart-type quelconque, que V2 a une moyenne nulle, que V3 a un écart-type de 1. V4 et V5 combinent le centrage (soustraction de la moyenne) et la réduction (division par l'écart-type) mais seul V5 aboutit à une variable centrée réduite (moyenne=0, écart-type=1).
V1 V2 V3 V4 V5 Moyenne 204.797 -0.000 1.389 -203.408 -0.000 Ecart-type 147.399 147.399 1.000 1.000 1.000Voici des fonctions qui implémentent ces opérations :
centrer <- function(v) { return( v-mean(v) ) } reduire <- function(v) { return( v/sd(v) ) } bien1 <- function(v) { return( (v -mean(v))/sd(v) ) } bien2 <- function(v) { return( (1/sd(v)*(v - mean(v)) ) ) } mal <- function(v) { return( v/sd(v) - mean(v) ) } msf <- function(tit,v) { cat(sprintf("Variable %-15s",tit) ) cat(sprintf(" moyenne %8.3f",mean(v) ) ) cat(sprintf(" écart-type %8.3f",sd(v) ) ) cat("\n") } # fin de fonction msf V <- V1 msf("V",V) msf("c(V)" ,centrer( V ) ) msf("r(V)" ,reduire( V ) ) msf("cr(V)" ,centrer( reduire(V) ) ) msf("rc(V)" ,reduire( centrer(V) ) ) msf("bien1" ,bien1( V ) ) msf("bien2" ,bien2( V ) ) msf("mal" ,mal( V ) )Et leurs résultats :
Variable V moyenne 204.797 écart-type 147.399 Variable c(V) moyenne -0.000 écart-type 147.399 Variable r(V) moyenne 1.389 écart-type 1.000 Variable cr(V) moyenne 0.000 écart-type 1.000 Variable rc(V) moyenne -0.000 écart-type 1.000 Variable bien1 moyenne -0.000 écart-type 1.000 Variable bien2 moyenne -0.000 écart-type 1.000 Variable mal moyenne -203.408 écart-type 1.000L'article traite de données issues de puces à ADN pour des cellules de levure (yeast en anglais), soumises à un stress chimique dans le but d'identifier les gènes impliqués dans la réponse à la drogue utilisée (le Bénomyl, en anglais Benomyl). Voir l'introduction à cette étude par Gaëlle Lelandais que remercions ici pour son autorisation à utiliser ces textes et ces données.
Le suffixe gpr correspond à un des formats de sortie de Genepix, qui est un logiciel de la société Molecular Devices et qui permet d'exploiter les résultats d'une expérimentation avec une puce à ADN. En R, on peut lire un fichier gpr avec read.GenePix ; la fonction read.GenePix fait partie du package marray dont qui sert à l'Exploratory analysis for two-color spotted microarray data. En R, il y a aussi une fonction read.genepix (sans majuscules) qui fait partie du package sma dont le nom signifie Statistical Microarray Analysis mais elle ne s'applique pas directement ici (de plus, le package semble abandonné) :
# Fichier gprnon.r library(sma) resultatPuce <- read.genepix("20min_Beno_Cy5_DMSO_Cy3.gpr",skip=36)> source("gprnon.r") Erreur dans scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings, : la ligne 5780 n'avait pas 56 élémentsVoici la lecture du fichier et son résultat avec read.GenePix du package marray :
options(width=450) library(marray) puce <- read.GenePix("20min_Beno_Cy5_DMSO_Cy3.gpr") print( dim(puce) ) print( class(puce) ) print( names(puce) ) print( head(puce) )############################################### ## ## listing nommé demogpr.sortie obtenu via ## > sink("demogpr.sortie") ## > source("demogpr.r",echo=TRUE) ## > sink() ## ############################################### # # les sorties sont "aménagées" pour prendre # moins de place à l'affichage # ############################################### > options(width=450) > library(marray) > puce <- read.GenePix("20min_Beno_Cy5_DMSO_Cy3.gpr") Reading ... 20min_Beno_Cy5_DMSO_Cy3.gpr > print( dim(puce) ) ==> [1] 15552 1 > print( class(puce) ) ==> [1] "marrayRaw" ==> attr(,"package") [1] "marray" > print( names(puce) ) ==> NULL > print( head(puce) ) ==> An object of class "marrayRaw" @maRf 20min_Beno_Cy5_DMSO_Cy3.gpr [1,] 992 @maGf 20min_Beno_Cy5_DMSO_Cy3.gpr [1,] 2561 @maRb 20min_Beno_Cy5_DMSO_Cy3.gpr [1,] 495 @maGb 20min_Beno_Cy5_DMSO_Cy3.gpr [1,] 2358 @maW 20min_Beno_Cy5_DMSO_Cy3.gpr [1,] 0 @maLayout An object of class "marrayLayout" @maNgr [1] 12 @maNgc [1] 4 @maNsr [1] 18 @maNsc [1] 18 @maNspots [1] 15552 @maSub [1] TRUE FALSE FALSE FALSE FALSE 15547 more elements ... @maPlate [1] 1 Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20... @maControls [1] probes Levels: probes @maNotes character(0) @maGnames An object of class "marrayInfo" @maLabels [1] "295078" @maInfo ID X295078 Name 295078 YDR060W_01 |YDR060W|MAK21|essential for 60s ribosome biogenesis; involved in nuclear export of pre-ribosomes @maNotes [1] "" @maTargets An object of class "marrayInfo" @maLabels character(0) @maInfo tableau de données (data frame) avec 0 colonnes et 0 lignes @maNotes character(0) @maNotes [1] "GenePix Data"L'exploitation de ce fichier sous forme de TP de BioInformatique peut se trouver ici. A nouveau, merci, Gaëlle !
R signifie "rouge" (Red), G signifie "vert" (Green), b signifie "fond (ou arrière-plan)" (Background) et f signifie "avant-plan" (Foreground). Comme indiqué dans la présentation des données, Cy5 (rouge) correspond à la culture avec drogue et Cy3 (vert) correspond à la culture sans drogue.
La page d'aide fait référence à une expérience avec une puce à ADN sur le poisson zèbre ; plus de détails ici (section 3, page 2) avec au passage quelques compléments sur Bioconductor et le package marray.
Pour le détail sur préparation statistique des données, nous renvoyons à nouveau aux transparents de Gaëlle Lelandais ; en résumé, il faut identifier et éliminer les spots saturants ("non conformes") puis corriger les biais expérimentaux, effectuer une normalisation inter-canaux des deux intensités de fluorescence.
On peut aussi consulter ce texte de S. Le Crom qui explique qu'à partir des variables R et G (pour rouge et vert) on calcule LG et LR qui en sont les logarithmes en base deux pour tracer ensuite M (différence des intensités) en fonction de A (moyenne des intensités). Au passage, savez-vous pourquoi on utilise les lettres A et M ?
Au vu des calculs effectués et des graphiques réalisés dans les séances, Excel est très incomplet pour les statistiques de base : pas de fonction t.test() ni boxplot() par exemple... De plus il ne dispose d'aucune fonction pour réaliser des ACP ou traiter des données de puces à ADN. Consultez aussi la section Pourquoi ne faut-il pas utiliser Excel de ma page statgen.
Les "grands" logiciels statistiques classiques sont SAS, SPSS, Statistica, Stata. Le plus grand avantage de R est d'être gratuit. On peut le télécharger pour Windows ici. Les autres grands points forts de R sont la diversité des structures de données (vecteurs, listes, "dataframes", objets...) alliée à la facilité de programmation et à la disponibilité d'un grand nombre de packages (1726 en mars 2009). On pourra consulter ici la liste des packages.
On pourra consulter les questions et réponses de la séance 5 de mes cours de biostatistiques (module 1) pour l'Ecole Doctorale afin d'obtenir des compléments sur la comparaison entre "grands" logiciels statistiques.
Pour la seule partie analyse statistique des données de puces à ADN, il existe de nombreux logiciels, dont J-express Pro (brochure ici), Mev (MultiExperiment Viewer) et Genesis Server (brochure ici) et bien sûr GenePix Pro.
Pour les analyses statistiques en général, il y a de nombreux packages qui proposent des menus et des interfaces et il est même possible d'interfacer R avec Excel. Voir http://www.sciviews.org/_rgui/. Les images ci-dessous sont cliquables...
Les taskviews permettent de visualiser les 3700 (avril 2012) packages regroupés par tâche, une tâche étant un thème particulier, comme la statistique génétique ou la phylogénétique génétique.
La liste des tâches est disponible à l'adresse http://cran.r-project.org/web/views/.
BioConductoR (wiki_en) est un ensemble d'outils pour traiter les données génomiques qui utilise le logiciel R. Pour plus d'infos, consulter le site officiel à l'adresse http://www.bioconductor.org/.
En particulier, la liste des packages R liés à BioConductoR est ici : http://www.bioconductor.org/help/bioc-views/release/bioc/. Un manuel de R et de bioconductor est disponible à l'adresse http://manuals.bioinformatics.ucr.edu/home/R_BioCondManual.
Un exemple en français d'analyse de transcriptome avec GNU R et bioconductor est à l'URL TD_transcriptome.html et une liste de documents associés se trouve chez biotoul mais vous pouvez aussi consulter le TP nommé Tp_Affy.
Comme son nom l'indique (seqinr = sequences in R), seqinr permet de récupérer et d'analyser des séquences biologiques. Ce package est écrit et géré par des français, il a sa page web en français : http://pbil.univ-lyon1.fr/software/seqinr/. Vous pouvez consulter ici l'index des fonctions de ce package.
Celui d'E. Paradis chez Springer :
Au vu des exemples proposés mais aussi de l'index des fonctions du package ape, il est clair qu'il est assez facile de générer des arbres phylogénétiques en R.
Surtout pas en croyant que vous savez tout faire ! Le but de ces séances était de montrer ce qu'on fait classiquement en statistiques, en biostatistiques et en statistiques pour la bioinformatique. Il reste à progresser en biologie, en biochimie, en informatique, en statistiques etc. avant de commencer à être autonome.
Par contre tout ce qui a été présenté permet de se rendre compte que tout est à votre portée et qu'il suffit d'approfondir. En ce qui concerne R, il faut certainement
Aller plus loin avec les manipulations de base, comme par exemple dans notre cours nommé RMFE ;
Progresser au niveau programmation en général, connaissance des bases de données et programmation en R, tant au niveau des structures de données, des fonctions que des calculs statistiques, par exemple en lisant R-intro ;
Progresser au niveau bioinformatique à l'aide du support de cours avec exercices corrigés de Wim P. Krijnen nommé Applied Statistics for Bioinformatics Using R disponible dans les manuels du CRAN -- contributed documentation (copie locale ici) .
Envisager de passer à une utilisation soutenue de la programmation R, avec l'écriture de fonctions efficaces, voire même de classes d'objets et de packages. Consultez par exemple nos séances R approfondi.