Valid XHTML     Valid CSS2    

 

R et la bioinformatique en 15 minutes

                gilles.hunault "at" univ-angers.fr

 

non su          non su          non su

 

Il y a une version anglaise disponible ici.

 

Décembre 2015, ultra-courte présentation pour le GBSA

(Groupe Bioinformatique Santé Angevin).

 

Table des matières cliquable

  1. Mini-présentation de R

  2. R et la bioinformatique dont la métagénomique

  3. Bibliographie R volontairement restreinte et orientée bioinformatique

  4. Démonstrations par l'exemple via R dans RStudio

  5. Bilan (rapide) : avantages et inconvénients de R

 

       En guise de conclusion : quelques bonnes nouvelles (pour survivre à R !)

 

Introduction

Après «ma thèse en 3 secondes», «toute l'Encyclopedia Universalis en 20 lignes» et «les fondements de la rationalité de l'être chez Schopenhauer et Nietzsche en moins de 120 signes», voici une présentation ultra-courte de 10 minutes plus 5 minutes de démonstration interactive sur R avec une forte orientation bioinformatique.

Ce n'est, bien sûr, ni pédagogique ni faussement inutile, même si cela ressemble un peu à de la mise en scène.

rideauxR15.jpg

1. Mini-présentation de R

R est un logiciel scriptable et programmable gratuit, disponible pour Windows, MacOs et Linux. Il est exhaustivement complet pour tout ce qui touche aux statistiques, ce qui inclut les biostatistiques et la bioinformatique.

R comprend :

  • des centaines de commandes de base ;

  • la possibilité d'enchainer les commandes et d'automatiser les traitements dans des scripts ;

  • des milliers de données prêtes à l'emploi, y compris en génomique, métagénomique ;

  • des structures de données adaptées aux statistiques, aux biostatistiques et à la bioinformatique : vecteurs, matrices, listes, "data frames", séries chronologiques, modèles statistiques... ;

  • des environnements d'exécution avec ou sans menu, avec production automatisée ou non de documents HTML, Word, LaTeX ou PDF ;

  • des milliers de fonctions déjà disponibles et la possibilité d'en programmer de nouvelles ;

  • des graphiques de qualité prévus pour les articles et les posters ;

  • une organisation en packages et en tâches pour [essayer de] s'y retrouver ;

  • des manuels de références, des documentations PDF nommées vignettes et des "vrais" livres ;

  • des millions d'utilisateurs et de professionnels, des blogs, des sites de référence, des conférences...

2. R et la bioinformatique dont la métagénomique

Les analyses bioinformatiques couvrent un large spectre de calculs et de techniques, que ce soit les constructions d'alignements de séquences, la réalisation d'arbres phylogénétiques, la visualisation sous forme de heatmaps, l'expression différentielle de gènes, l'analyse de réseaux d'interactions ou de réseaux de gènes, la prédiction de structures, la génétique des populations, la représentation en bigraphes de gènes et d'articles PUBMED, le calcul de distances génomiques euclidiennes ou non euclidiennes.. .

           ../Introduction_R/heatmap01.png ../Introduction_R/bigraph.png http://mixomics.org/wp-content/uploads/2015/03/network_spls.png
           Heatmap (code) Bigraph (code) Network (code)

R est particulièrement bien adapté à ce genre de calculs et d'analyses parce que R est un langage assez ouvert et qui s'interface bien avec des exécutables utilisables en ligne de commandes. De plus les possibilités graphiques de R et sa lecture de fichiers sur Internet lui permettent d'accéder dynamiquement aux grandes bases de données biologiques du NCBI, de l'EBI, d'UNIPROT, sans compter son site dédié nommé Bioconductor...

            ncbi  pubmed workflows clinical trials

De plus R n'est pas en reste pour tout ce qui concerne les données massives, que ce soit le faux Big Data, nommé aussi Big Data Statistique (disons moins de 100 Go de données, mais des tableaux de dimension (n,p) avec p>>n) ou le vrai Big Data, pour lequel les données ne tiennent pas en mémoire centrale, et qui nécessite plusieurs ordinateurs, plusieurs serveurs via un cluster ou un grid.

                     Disease Model           R Hadoop
          SRAdb metadata 15 Go            non su data 90 To     MGRASTer

Enfin, la centralisation et la mise à disposition des calculs, données, exemples et fichiers d'aides font de R un laboratoire interactif d'outils bioinformatiques inégalable, même si le langage Python, plus généraliste, se révèle un concurrent sérieux et efficace, notamment quand on sait déjà programmer...

  cRan 7 400 packages 150 000 fonctions 52 000 graphiques

3. Bibliographie R volontairement restreinte et orientée bioinformatique

Des recherches sur le web avec les expressions bioinformatics, "with R", "using R" et books ou filetype:pdf aboutissent à de nombreux ouvrages, parfois publiés, parfois simplement distribués sur le Web, comme on peut le vérifier ci-dessous.

   

   

Voici mes favoris :

non su non su non su
Bioconductor Case Studies Bioinformatics[...] using R[...] Statistical Bioinformatics with R
(2008, 284 p.) (2005, 474 p.) (2010, 336 p.)
non su non su non su
Bioinformatics with R Cookbook Genomic Data with R Microarays with R
(2014, 300 p.) (2015, 280 p.) (2011, 1036 p.)

4. Démonstrations par l'exemple via R dans RStudio

Voici un texte R simple à lire pour compter puis analyser le nombre de nucléotides du gène X94991. C'est bien sûr R qui fait la lecture au NCBI de la séquence du gène via Internet (tout ce qui commence par # est un commentaire, ignoré par R) :


     #############################################
     #                                           #
     # comptage des nucléotides du gène X94991   #
     #                                           #
     #############################################
     #
     # pour le GBSA          (gH). novembre 2015 #
     
     library(ape)                                        # chargement de la librairie
     geneX <- read.GenBank("X94991.1",as.character=TRUE) # lecture du gène demandé
     cnt <- table(geneX)                                 # comptages par type de nucléotide
     print( cnt )                                        # affichage
     print(chisq.test(cnt))                              # calcul du chi-deux d'adéquation
     
     

Voici le résultat de l'exécution 


     geneX
       a   c   g   t
     410 789 573 394
     
     Chi-squared test for given probabilities
     data:  cnt X-squared = 187.07, df = 3, p-value < 2.2e-16
     
     

Moyennant quelques autres commandes on peut obtenir un affichage plus propre et graphique des comptages 


     #############################################
     #                                           #
     # comptage des nucléotides du gène X94991   #
     # puis pourcentages et graphiques           #
     #                                           #
     #############################################
     #
     # pour le GBSA          (gH). novembre 2015 #
     
     # tri à plat des nucléotides du gène
     
     tap <- table(geneX)
     
     # répartition détaillée des nucléotides
     
     props        <- as.numeric(prop.table(tap))
     pcts         <- paste(round(100*props,1),"%")
     dfRes        <- data.frame(tap,props,pcts)
     names(dfRes) <- c("Nucléotides","Comptages","Proportions","Pourcentages")
     
     print(dfRes)
     
     # visualisation graphique
     
     lng   <- sum(tap)
     coul  <- rainbow(dim(tap)) # mieux que c("red","green","blue","yellow")
     titre <- paste("Nucléotides du gène X94991\n (longueur=",lng,")",sep="")
     
     barplot(tap,main=titre,col=coul,legend.text=pcts,xlim=c(0,6))
     

On obtient alors 


       Nucléotides Comptages Proportions Pourcentages
     1           a       410   0.1892890       18.9 %
     2           c       789   0.3642659       36.4 %
     3           g       573   0.2645429       26.5 %
     4           t       394   0.1819021       18.2 %
     

          genex02.png

Ensuite, si on programme «un peu», on peut fournir des sorties en français et automatiser des traitements pour les utilisateurs de base :


     #############################################
     #                                           #
     # comptage des nucléotides du gène X94991   #
     # calcu du chi2 d'équirépartition           #
     #                                           #
     #############################################
     #
     # pour le GBSA          (gH). novembre 2015 #
     
     # un chi-deux détaillé en français avec conclusion automatique
     # issu des fonctions du "p'titgros"
     # disponibles à l'adresse
     #   http://www.info.univ-angers.fr/~gh/wstat/statghfns.php
     
     source("http://forge.info.univ-angers.fr/~gh/wstat/statgh.r",encoding="latin1")
     
     # réalisation du chi-deux et explications en français
     
     # il s'agit ici du chi-deux d'adéquation à une distribution théorique
     # (de même effectif global) soit ici la loi uniforme discrète,
     # autrement dit, la loi sous hypothèse d'équirépartition
     
     chi2Adeq(vth=rep(length(geneX[[1]])/4,4),vobs=tap,graph=TRUE)
     
     

      CALCUL DU CHI-DEUX D'ADEQUATION
     
      Valeurs théoriques   541.5 541.5 541.5 541.5
      Valeurs observées    410   789   573   394
      Valeur du chi-deux   187.0674
     
      Chi-deux max (table) à 5 %  7.814728 pour   3  degrés de liberté ; p-value  2.62465e-40
     
      décision : au seuil de  5 % on peut rejeter l'hypothèse
      que les données observées correspondent aux valeurs théoriques.
     
     

Enfin, avec quelques commandes Markdown bien choisies, par exemple si on ajoute en début de fichier  


     ---
     title: bRavo !
     output: all
     ---
     # un exemple vite fait
     ```{r}
     
     
     [...] ici le code R précédent
     

alors on produit en un seul clic les fichiers suivants :

Word genex04.docx
HTML genex04.html
PDF genex04.pdf

Vous pouvez vérifier tout ceci en exécutant le fichier genex04.rmd suivant 


     ---
     title: "demoR15"
     output: all
     ---
     
     # Gène X94991
     
     ```{r}
     library(ape)                                        # chargement de la librairie
     geneX <- read.GenBank("X94991.1",as.character=TRUE) # lecture du gène demandé
     cnt <- table(geneX)                                 # comptages par type de nucléotide
     print( cnt )                                        # affichage
     print( cnt )                                        # affichage
     ```
     
     # Histogramme des fréquences
     
     ```{r, echo=FALSE}
     # répartition détaillée des nucléotides
     
     tap <- table(geneX)     
     props        <- as.numeric(prop.table(tap))
     pcts         <- paste(round(100*props,1),"%")
     dfRes        <- data.frame(tap,props,pcts)
     names(dfRes) <- c("Nucléotides","Comptages","Proportions","Pourcentages")
     print(dfRes)
          
     # visualisation graphique
          
     lng   <- sum(tap)
     coul  <- rainbow(dim(tap)) # mieux que c("red","green","blue","yellow")
     titre <- paste("Nucléotides du gène X94991\n (longueur=",lng,")",sep="")
          
     barplot(tap,main=titre,col=coul,legend.text=pcts,xlim=c(0,6))
     ```
     # Avec les fonctions (gH)
     ```{r,echo=FALSE,print=FALSE}
     source("http://forge.info.univ-angers.fr/~gh/wstat/statgh.r",encoding="latin1")
     ```
     ```{r}
     # réalisation du chi-deux et explications en français
          
     # il s'agit ici du chi-deux d'adéquation à une distribution théorique
     # (de même effectif global) soit ici la loi uniforme discrète,
     # autrement dit, la loi sous hypothèse d'équirépartition
          
     chi2Adeq(vth=rep(length(geneX[[1]])/4,4),vobs=tap,graph=TRUE)
     ```
     

5. Bilan (rapide) : avantages et inconvénients de R

Avantages de R

Très bel outil statistique, gratuit, complet, bien documenté et bien adapté à la bioinformatique, aux publications de recherche et aux études reproductibles car il permet de s'affranchir des copier/coller pour produire les rapports et articles.

Inconvénients de R

Difficile à maitriser rapidement sans y passer du temps. Oblige (mais est-ce vraiment un défaut ?) à programmer pour être très efficace.

Quelques bonnes nouvelles
(pour survivre à R !)

Pour finir cet exposé, quatre bonnes nouvelles pour éviter qu'utiliser R ne soit une galère :

                         ../../Phenotic/galere.jpg

 

 

retour gH    Retour à la page principale de   (gH)