Exercices longs

EXERCICES LONGS

Les exercices proposés ici demandent nettement plus de réflexion et d'investissement que les petits exercices présentés précédemment. L'exécution des scripts R correspondant dure parfois plusieurs minutes. Il sera donc indispensable de prévoir des paramétrages des fonctions pour tester des petits jeux d'essais, voire d'effectuer un "profilage" des temps passés dans chaque fonction afin d'obtenir un script réutilisable et "assez rapide".

Il s'agit bien sûr de données réelles traitées avec des collègues ou des doctorant(e)s.

1. Comptages de sites en génomique

Une bioinformaticienne vient d'utiliser le site PLACE. Elle se retrouve avec une quatre-vingtaine de fichiers dans le dossier de longévité et un peu plus de mille six cent fichiers dans le répertoire de réseau. On trouvera ces fichiers dans l'archive exercicePLACE.zip (8 Mo en compressés, 60 Mo décompressés). Il faudrait produire un fichier de comptage de chaque site précédé de la lettre L pour longévité et de la lettre N pour réseau (network) sachant que tous les sites indiqués dans les fichiers ne sont pas tous présents de la même façon.

On trouvera dans placexmp_long.txt un exemple de fichier à traiter (le fichier 20.txt), dont voici une présentation simplifiée :


     Web Signal Scan Program
     
     Database Searched: PLACE
     
     This is the sequence you submitted
     >nc|dna|Medtr4g055610_upstream parent=4 chrom=4 range=17080451-17081950 length=1500 strand=1  sp=Medicago_truncatula, 1500 bases, A9ED66C3 checksum.
     AAACAAAATAAATTTACAGTCAGCGTGCCACATTAGCGAAAATGCGCATT
     CAGATGACCTATGGGGTATTTTGAAACAATTTTCTTTTACAAGAACCGAA
     TTGAATTTTTTTTATAGGAGGCAAACCAGGAAAAAAACTATATATTATAG
     GAGGG...
     
     RESULTS OF YOUR SIGNAL SCAN SEARCH REQUEST
     
     ../../tmp/sigscan//signaldone.10566: 1500 base pairs
     Signal Database File: user.dat
     
          Factor or Site Name          Loc.(Str.)       Signal Sequence             SITE #
     _____________________________________________________________________________________
     -300ELEMENT               site     84 (-) TGHAAARK                             S000122
     2SSEEDPROTBANAPA          site   1149 (+) CAAACAC                              S000143
     2SSEEDPROTBANAPA          site    589 (-) CAAACAC                              S000143
     ...
     WRKY71OS                  site     55 (+) TGAC                                 S000447
     WRKY71OS                  site    201 (+) TGAC                                 S000447
     WRKY71OS                  site    835 (+) TGAC                                 S000447
     WRKY71OS                  site   1075 (+) TGAC                                 S000447
     WRKY71OS                  site   1206 (+) TGAC                                 S000447
     WRKY71OS                  site   1437 (+) TGAC                                 S000447
     WRKY71OS                  site     19 (-) TGAC                                 S000447
     WRKY71OS                  site    287 (-) TGAC                                 S000447
     WRKY71OS                  site    322 (-) TGAC                                 S000447
     WRKY71OS                  site    814 (-) TGAC                                 S000447
     WRKY71OS                  site    987 (-) TGAC                                 S000447
     WRKY71OS                  site   1003 (-) TGAC                                 S000447
     WRKY71OS                  site   1040 (-) TGAC                                 S000447
     WRKY71OS                  site   1058 (-) TGAC                                 S000447
     WRKY71OS                  site   1319 (-) TGAC                                 S000447
     
     
     
     -------------------------------------------
     
     o If you use this program in published research, please cite:
       - Higo, K., Y. Ugawa, M. Iwamoto and T. Korenaga (1999) Plant cis-acting
         regulatory DNA elements (PLACE) database:1999. Nucleic Acids Research
         Vol.27 No.1 pp. 297-300.
       - Prestridge, D.S. (1991) SIGNAL SCAN: A computer program that scans DNA
         sequences for eukaryotic transcriptional elements. CABIOS 7, 203-206.

On doit ici comptabiliser 1 occurence de -300ELEMENT, 2 occurences de 2SSEEDPROTBANAPA et 15 occurences de WRKY71OS.

On trouvera dans exoLong01.csv (1,2 Mo) le résultat à obtenir... soit 1708 lignes et 344 colonnes.

Avec un peu de patience, votre navigateur en affichera la version texte. En voici un extrait :


     Fichier     -10PEHVPSBD   2SSEEDPROTBANAPA   -300CORE   -300ELEMENT   -300MOTIFZMZEIN   5256BOXLELAT5256   5659BOXLELAT5659   AACACOREOSGLUB1
     L_2.txt      1            1                   0         1   0   0   0   3
     L_3.txt      4            0                   0         0   0   0   0   2
     L_4.txt      1            2                   0         2   0   0   0   0
     L_5.txt      2            0                   0         4   0   0   0   2
     L_6.txt      2            0                   1         3   0   0   0   0
     L_7.txt      2            0                   0         2   0   0   0   1
     L_8.txt      5            0                   0         3   0   0   0   0
     ...
     L_82.txt     1            0                   0         0   0   0   0   0
     L_83.txt     2            0                   0         0   0   0   0   0
     L_84.txt     0            0                   0         1   0   0   0   0
     L_85.txt     1            0                   1         2   0   0   0   0
     N_1.txt      3            0                   1         4   0   0   0   3
     N_2.txt      2            1                   1         0   0   0   0   0
     N_3.txt      2            1                   0         0   0   0   0   2
     N_4.txt      0            3                   0         4   0   0   0   0
     N_5.txt      3            1                   0         1   0   0   0   2
     N_6.txt      4            0                   1         3   0   0   0   0
     ...
     N_1621.txt   3            0                   0         5   0   0   0   1
     N_1622.txt   4            2                   1         3   0   0   0   2
     N_1623.txt   0            0                   0         4   0   0   0   1
     N_1624.txt   2            1                   1         2   0   0   0   2

Données gentiment fournies par Karima RIGHETTI dans le cadre d'un travail de post doctorat encadré par Julia BUITINK sur Medicago truncatula.

Solution : afficher la solution

Analyse et explications

Etant donné le nombre de fichiers à traiter, il est certainement prudent de construire une maquette de programme avec une dizaine de fichiers seulement. C'est ce que nous avons fait, avec un répertoire miniplace qui contient 4 fichiers dans le répertoire miniplace/lon et 8 fichiers dans le répertoire miniplace/net.

Déterminer le nomber de fichiers à passer en revue n'est pas un problème, puisque R dispose de la fonction list.files().

Si on arrive à lire chaque fichier et à renvoyer chaque site vu avec son nombre d'occurences, il suffira d'utiliser la fonction merge() pour tout comptabiliser.

Vu le nombre de fichiers, il faudra sans doute transposer les résultats avant de les exporter sous Excel.

Pour trouver les sites, il suffit de répérer la ligne de tiret et la ligne de pointillés. Entre les deux, on a des sites comme mot numéro 1.

Comme premier jet de programme, on peut donc écrire quelque chose comme ce qui suit :


     0001    # version maquette :
     0002
     0003    #  4 fichiers dans miniplace/lon
     0004    #    "2.txt" "3.txt" "4.txt" "5.txt"
     0005
     0006    #  8 fichiers dans miniplace/net
     0007    #    "1.txt" "2.txt" "3.txt" "4.txt" "5.txt" "7.txt" "8.txt" "9.txt"
     0008
     
     0009    #  étape 1 : liste des fichiers longévité
     0010
     0011    path1 <- "miniplace/lon"
     0012    ldf1  <- list.files(pattern="^[0-9]*\\.txt",path=path1)
     0013
     
     0014    # étape 2 : liste des fichiers networks
     0015
     0016    path2 <- "miniplace/net"
     0017    ldf2 <- list.files(pattern="^[0-9]*\\.txt",path=path2)
     0018
     0019    nbf1 <- length(ldf1)
     0020    nbf2 <- length(ldf2)
     0021    cat(" on doit traiter ",nbf1," fichiers longévité\n")
     0022    cat(" on doit traiter ",nbf2," fichiers network\n")
     0023
     
     0024    # étape 3 : on passe en revue les fichiers longévité
     0025
     0026    cats("longévité")
     0027    numf <- 0
     0028
     0029    for (idf in ldf1) {
     0030
     0031      numf <- numf + 1
     0032      nomf <- paste(path1,"/",idf,sep="")
     0033      cat("fichier ",sprintf("%4d",numf)," / ",nbf1," : ",nomf,"\n")
     0034
     0035      texte  <- ""
     0036      lignes <- readLines(nomf)
     0037      nbligs <- length(lignes)
     0038      sites  <- 0 # indique si on est en mode lecture de site ou non
     0039
     0040      for (idl in (1:nbligs)) {
     0041        laLigne <- lignes[idl]
     0042        if (nchar(laLigne)>0) {
     0043           if (substr(laLigne,1,5)=="-----") { sites <- 0 }
     0044           if (sites==1) {
     0045             mots <- strsplit(laLigne," ")
     0046             mot1 <- mots[[1]][1]
     0047             # pour debug :
     0048             ## cat(idl," : on retient juste ",mot1," dans ",laLigne,"\n")
     0049             texte <- paste(texte,mot1)
     0050           } # fin si
     0051           if (substr(laLigne,1,5)=="_____") { sites <- 1 }
     0052        } # fin si
     0053      } # fin pour idl
     0054
     0055      # pour debug :
     0056      ## print(texte) # tous les mots rencontrés
     0057
     0058      # comptages avec notre fonction analexies
     0059
     0060      vus <- analexies(texte,FALSE,FALSE,FALSE)$tmots
     0061
     0062      # mise en forme en long avec nom de fichier
     0063
     0064      dfm <- as.data.frame(matrix(nrow=nrow(vus),ncol=2))
     0065      colnames(dfm) <- c("site",nomf)
     0066      dfm[,1] <- row.names(vus)
     0067      dfm[,2] <- vus[1] # attention !
     0068
     0069      # pour debug
     0070      ## print(dfm)
     0071
     0072      # si c'est le premier fichier, on initialise, sinon on fusionne
     0073
     0074      if (numf==1) {
     0075        cntSites1 <- dfm
     0076      } else {
     0077        cntSites1 <- merge(cntSites,dfm,by="site",all.x=TRUE,all.y=TRUE)
     0078      } # fin si
     0079
     0080    } # # fin pour sur idf
     0081
     0082    # pour voir les résultats, mieux vaut uiliser un navigateur
     0083
     0084    ## hprint(cntSites1)
     0085
     
     0086    # étape 4 : on passe en revue les fichiers network
     0087
     0088    cats("network")
     0089    numf <- 0
     0090
     0091    for (idf in ldf2) {
     0092
     0093      numf <- numf + 1
     0094      nomf <- paste(path2,"/",idf,sep="")
     0095      cat("fichier ",sprintf("%4d",numf)," / ",nbf2," : ",nomf,"\n")
     0096
     0097      texte  <- ""
     0098      lignes <- readLines(nomf)
     0099      nbligs <- length(lignes)
     0100      sites  <- 0 # indique si on est en mode lecture de site ou non
     0101
     0102      for (idl in (1:nbligs)) {
     0103        laLigne <- lignes[idl]
     0104        if (nchar(laLigne)>0) {
     0105           if (substr(laLigne,1,5)=="-----") { sites <- 0 }
     0106           if (sites==1) {
     0107             mots <- strsplit(laLigne," ")
     0108             mot1 <- mots[[1]][1]
     0109             # pour debug :
     0110             ## cat(idl," : on retient juste ",mot1," dans ",laLigne,"\n")
     0111             texte <- paste(texte,mot1)
     0112           } # fin si
     0113           if (substr(laLigne,1,5)=="_____") { sites <- 1 }
     0114        } # fin si
     0115      } # fin pour idl
     0116
     0117      # pour debug :
     0118      ## print(texte) # tous les mots rencontrés
     0119
     0120      # comptages avec notre fonction analexies
     0121
     0122      vus <- analexies(texte,FALSE,FALSE,FALSE)$tmots
     0123
     0124      # mise en forme en long avec nom de fichier
     0125
     0126      dfm <- as.data.frame(matrix(nrow=nrow(vus),ncol=2))
     0127      colnames(dfm) <- c("site",nomf)
     0128      dfm[,1] <- row.names(vus)
     0129      dfm[,2] <- vus[1] # attention !
     0130
     0131      # pour debug
     0132      ## print(dfm)
     0133
     0134      # si c'est le premier fichier, on initialise, sinon on fusionne
     0135
     0136      if (numf==1) {
     0137        cntSites2 <- dfm
     0138      } else {
     0139        cntSites2 <- merge(cntSites,dfm,by="site",all.x=TRUE,all.y=TRUE)
     0140      } # fin si
     0141
     0142    } # # fin pour sur idf
     0143
     0144    # pour voir les résultats, mieux vaut uiliser un navigateur
     0145
     0146    hprint(cntSites2)
     0147
     
     0148    # étape 5 : il faut fusionner tous les résultats
     0149
     0150    cnt3  <- merge(cntSites1,cntSites2,by="site",all.x=TRUE,all.y=TRUE)
     0151    cnt3[ is.na(cnt3) ] <- 0
     0152    cnt4  <- t(cnt3[,-1])
     0153    colnames(cnt4)  <- cnt3[,1]
     0154
     0155    # on vérifie :
     0156
     0157    bloc(cnt4)
     0158
     0159    # il reste à exporter en .csv pour excel...
     0160

Le résutat d'exécution simplifié est alors


      on doit traiter  4  fichiers longévité
      on doit traiter  8  fichiers network
     
     longévité
     =========
     
     fichier     1  /  4  :  miniplace/lon/2.txt
     fichier     2  /  4  :  miniplace/lon/3.txt
     fichier     3  /  4  :  miniplace/lon/4.txt
     fichier     4  /  4  :  miniplace/lon/5.txt
     
     network
     =======
     
     fichier     1  /  8  :  miniplace/net/1.txt
     fichier     2  /  8  :  miniplace/net/2.txt
     fichier     3  /  8  :  miniplace/net/3.txt
     fichier     4  /  8  :  miniplace/net/4.txt
     fichier     5  /  8  :  miniplace/net/5.txt
     fichier     6  /  8  :  miniplace/net/7.txt
     fichier     7  /  8  :  miniplace/net/8.txt
     fichier     8  /  8  :  miniplace/net/9.txt
     
                           -10PEHVPSBD 2SSEEDPROTBANAPA -300ELEMENT AACACOREOSGLUB1
     miniplace/lon/2.txt.x           1                1           1               3
     miniplace/lon/3.txt.x           4                0           0               2
     miniplace/lon/4.txt.x           1                2           2               0
     miniplace/lon/5.txt.x           2                0           4               2
     miniplace/lon/5.txt.y           2                0           2               2
     miniplace/lon/2.txt.y           1                1           1               3
     miniplace/lon/3.txt.y           4                0           0               2
     miniplace/lon/4.txt.y           1                2           2               0
     miniplace/lon/5.txt             2                0           4               2
     miniplace/net/9.txt             2                2           2               0

Même si cela semble fonctionner, c'est-à-dire produire un résultat, celui-ci est faux (nous y reviendrons). Il est clair aussi que les étapes 3 et 4 sont identiques, au nom prés de la liste des fichiers. Il est donc plus court d'écrire un sous-programme pour automatiser le même traitement. La lecture du fichier place peut aussi se mettre sous forme de sous-programme. Enfin, au lieu d'utiliser les noms de répertoire, une initiale de dossier serait plus courte à lire dans les résultats. Si on tient compte des ces remarques, on peut alors une version plus propre comme dans l'encadré ci-dessous


     #  (gH)   -_-  countsites.r  ;  TimeStamp (unix) : 26 Juin 2013 vers 16:57
     
     # source("statgh.r")
     
     library(hash)
     library(lubridate)
     
     ##############################################################################
     
     countSites <- function(rep1="",rep2="") {
     
     ##############################################################################
     
        if (missing(rep1) | (rep1=="")) {
          cat(" counSites() : fusion de fichiers txt pour compter les sites \n\n")
          cat("syntaxe : countSites(rep1_longevity,rep2_network) \n")
          cat("exemple : countSites(\"results_longevity_module\",\"results_network\") \n")
          cat("\n")
          cat("le mieux est d'exécuter cette fonction dans le répertoire parent de ces deux répertoires\n\n")
          return("")
        } # fin d'aide
     
        tdeb <- Sys.time()
     
        # ---------------------------------------------------------------------------
     
        cat("\n\n")
        cat("COMPTAGE DES SITES dans les répertoires ",rep1," et ",rep2)
        cat("\n\n")
     
        # ---------------------------------------------------------------------------
     
        nomf1 <- "longevity.csv"
        nomf2 <- "network.csv"
        nomf3 <- "global.csv"
        nomf4 <- "tglobal.csv"
     
        tdf1  <- compteFichiers(rep1)
        tdf2  <- compteFichiers(rep2)
     
        # ---------------------------------------------------------------------------
     
        cnt1  <- compteSite(tdf1,rep1,"L")
        cnt2  <- compteSite(tdf2,rep2,"N")
     
        # ---------------------------------------------------------------------------
     
        cat("\n",format(Sys.time(),"%X")," Fusion...\n")
        cnt3  <- merge(cnt1,cnt2,by="site",all.x=TRUE,all.y=TRUE)
        cnt3[ is.na(cnt3) ] <- 0
        cnt4  <- t(cnt3[,-1])
        colnames(cnt4)  <- cnt3[,1]
     
        # ---------------------------------------------------------------------------
     
        write.csv(x=cnt1,file=nomf1,row.names=FALSE)
        write.csv(x=cnt2,file=nomf2,row.names=FALSE)
        write.csv(x=cnt3,file=nomf3,row.names=FALSE)
        write.csv(x=cnt4,file=nomf4,row.names=TRUE)
     
        cat("\n")
        cat("FIN DE COMPTAGE DES SITES dans les répertoires ",rep1," et ",rep2,"\n\n")
        tfin <- Sys.time()
        cat(" durée : ",as.duration(tfin-tdeb)," secondes entre ",format(tdeb,"%X")," et ",format(tfin, "%X"))
        cat("\n")
        cat("\n")
        cat(" Fichiers résultats : ",nomf1,dims(cnt1),", ",nomf2,dims(cnt2),",\n ",nomf3,dims(cnt3)," et ",nomf4,dims(cnt4),"\n")
     
     } # fin de fonction countSites
     
     ##############################################################################
     ##############################################################################
     
     fichiersTxt <- function(rep=".") {
     
     ##############################################################################
     
        fs  <- paste("/tmp/liste.",runif(1,1,10**8),sep="")
        cmd <- paste("(cd ",rep," ; ls *.txt | sort -n > ",fs," ) ",sep="")
        # cat(cmd,"\n")
        system(cmd)
     
        return(as.data.frame(read.table(fs,as.is=c(1))))
     
     } # fin de fonction fichiersTxt
     
     ##############################################################################
     
     compteFichiers <- function(rep) {
     
     ##############################################################################
     
        tdf <- fichiersTxt(rep)
        nbf <- nrow(tdf)
        cat(" il y a ",sprintf("%4d",nbf)," fichiers .txt dans ",rep,"\n")
     
        return(tdf)
     
     } # fin de fonction compteSite
     
     ##############################################################################
     
     tdm2df <- function(tdm,nomcol) {
     
     ##############################################################################
     
       leDf           <- as.data.frame(matrix(nrow=nrow(tdm),ncol=2))
       colnames(leDf) <- c("site",nomcol)
       leDf[,1]       <- row.names(tdm)
       leDf[,2]       <- tdm[1]
     
       return(leDf)
     
     } # fin de fonction tdm2df
     
     ##############################################################################
     
     compteSite <- function(tdf,rep,inirep) {
     
     ##############################################################################
     
        cat("\n")
        nbf <- nrow(tdf)
        cat(" Traitement des ",sprintf("%4d",nbf)," fichiers .txt dans ",rep,"\n")
        for (idf in (1:nbf)) {
            leFichier  <- tdf[idf,1]
            nomFichier <- paste(rep,"/",leFichier,sep="")
            cat("   fichier ",sprintf("%5d",idf)," / ",nbf," : ", leFichier,"\n")
            lesmots <- sitesVus(nomFichier)
            if (lesmots=="") {
               cat("           pas de site vu\n")
            } else {
               # cat(" on comptabilise les mots de ",lesmots,"\n")
               tdmc      <- analexies(lesmots,FALSE,FALSE,FALSE)$tmots
               idFichier <- paste(inirep,"_",leFichier,sep="")
               tdmcDf    <- tdm2df(tdmc,idFichier)
               if (idf==1) {
                 cntSites <- tdmcDf
               } else {
                 cntSites <- merge(cntSites,tdmcDf,by="site",all.x=TRUE,all.y=TRUE)
               } # fin si
            } # fin si
        } # fin pour idf
        cat("\n")
     
        # merge produit des NA au lieu de 0, donc :
     
        cntSites[ is.na(cntSites) ] <- 0
     
        return(cntSites)
     
     } # fin de fonction compteSite
     
     ##############################################################################
     
     sitesVus <- function(unFichier) {
     
     ##############################################################################
     
     # on analyse unFichier, produit par PLACE (http://www.dna.affrc.go.jp/PLACE/index.html)
     # le premier site est indiqué après une ligne de "_"
     # lorsqu'on voit plusieurs "--" qui se suivent, il n'y a plus de site à lire
     # chaque site correspond au premier mot de la ligne
     
     texte  <- ""
     lignes <- readLines(unFichier)
     nbligs <- length(lignes)
     sites  <- 0
     
     for (idl in (1:nbligs)) {
       laLigne <- lignes[idl]
       if (nchar(laLigne)>0) {
          if (substr(laLigne,1,5)=="-----") { sites <- 0 }
          if (sites==1) {
            mots <- strsplit(laLigne," ")
            mot1 <- mots[[1]][1]
            # cat(idl," : on retient juste ",mot1," dans ",laLigne,"\n")
            texte <- paste(texte,mot1)
          } # fin si
          if (substr(laLigne,1,5)=="_____") { sites <- 1 }
      } # fin si
     } # fin pour idl
     
     # print(texte)
     
     return(texte)
     
     } # fin de fonction sitesVus
     
     ##############################################################################
     
     ##############################################################################
     ##############################################################################
     
     # countSites("testmo","testne") ## (maquette)
     # cat("\n")
     countSites("results_longevity_module","results_network")

Le résutat d'exécution simplifié est lisible à l'adresse countsites_v2sor.txt et le listing numéroté associé est countsite_V1.

La programmation est un art. Si le programme précédent semble bien remplir le tableau demandé, une vérification "au hasard" sur une dizaine de fichiers montre que le nombre de mots renvoyés est faux. Après une petite séance de "debug", le coupable est trouvé : il s'agit de la fonction tdm2df (tableau des mots vers data frame) qui renvoie tdm[1] au lieu de tdm[,1]. Moral de l'histoire : il faut toujours vérifier ses résultats, même si sur la maquette on trouve les bons résultats.

On pourrait se demander si construire deux tableaux des sites et les fusionner (tableaux cnt1, cnt2, cnt3) est plus long que de garder un seul tableau cnt dès le départ. Si on teste ces différentes versions, le profilage montre qu'on n'y gagne rien en performance, la partie longue du traitement étant dues à la fonction analexies().

Version correcte et optimisée : countsite_V2.

2. Production automatique de graphiques en protéinique

Un bioinformaticien dispose de 7 groupes de séquences Fasta de protéines. Grâce au site protpc il arrive à produire 7 fichiers .CSV de propriétés physico-chimiques relatives aux acides aminés. Ces fichiers sont disponibles dans l'archive exercicePROTPC.zip. Il faudrait construire un fichier global protpc7grps.Rdata avec indication du numéro de groupe sachant que les groupes doivent être dans l'ordre Pool1, Pool2, Pool3, WHy, Class8, IDP, FS avant de produire des graphiques systématiquement en boxplot et, quand c'est possible, en beanplot, de façon à réaliser un manuel des graphiques comme protpc7grps.pdf.

On voudrait également un tableau résumé des analyses par groupe paramétrique (ANOVA) et non paramétrique (KRUSKAL-WALLIS) sous la forme d'un tableau de p-values et un tableau d'interprétation de médianes (+1 si la médiane du groupe est supérieure à la médiane globale, -1 si inférieure, avec une tolérance de 5 %).

Voici un extrait des données issues de PROTPC :


     PROTID               ; LENGTH   ;  PI        ;  MW        ;  FI        ;  GRAVY     ;  CHARGE    ;  HYDROPHI
     AAA34137             ;      115 ;     6.9400 ; 13131.4200 ;    -0.1392 ;    -1.1739 ...
     AAA86052             ;      263 ;     4.7100 ; 29054.3800 ;    -0.3565 ;    -1.6008 ...
     AAB72175             ;      114 ;     9.7800 ; 12395.7300 ;    -0.0245 ;    -0.6895 ...
     AAC41651             ;      282 ;    11.1400 ; 29984.7300 ;    -0.1937 ;    -0.7762 ...
     AAP37981             ;      112 ;     6.7000 ; 12607.6100 ;    -0.1874 ;    -1.2705 ...
     AAT38812             ;      800 ;     6.0900 ; 85093.8400 ;    -0.2181 ;    -1.4731 ...
     AAT42177             ;      224 ;     8.3900 ; 25340.7000 ;    -0.0740 ;    -0.9473 ...
     AAV76001             ;      234 ;     8.1700 ; 26084.3900 ;    -0.0820 ;    -1.0316 ...
     AAX96480             ;      138 ;     6.2600 ; 15466.9100 ;    -0.1940 ;    -1.2435 ...
     AAY97997             ;      111 ;     6.5800 ; 12693.7400 ;    -0.2028 ;    -1.2901 ...
     ABD28662             ;      320 ;     7.0300 ; 35549.3400 ;    -0.0260 ;    -0.8544 ...
     ABN09808             ;      124 ;     6.0800 ; 13977.0900 ;    -0.1533 ;    -1.1976 ...
     ACG35620             ;      103 ;     6.9400 ; 11689.4800 ;    -0.1910 ;    -1.3350 ...
     ACG48760             ;      192 ;    10.5700 ; 20260.9600 ;    -0.1944 ;    -1.2234 ...
     ACJ37937             ;       95 ;     6.1900 ; 10630.0800 ;    -0.2549 ;    -1.4000 ...
     ACP28171             ;      322 ;     7.0100 ; 35062.2100 ;    -0.0833 ;    -1.0497 ...
     ACZ60123             ;      143 ;     6.1400 ; 16202.7800 ;    -0.2028 ;    -1.2098 ...
     ACZ60133             ;      143 ;     6.1900 ; 16179.7800 ;    -0.1975 ;    -1.2154 ...
     ADG57880             ;      185 ;     9.3700 ; 20486.4900 ;    -0.1327 ;    -1.0870 ...
     ADJ67682             ;      732 ;     6.0800 ; 79198.5900 ;    -0.2202 ;    -1.4656 ...
     AED91653             ;      168 ;    10.0100 ; 18465.6200 ;    -0.1568 ;    -1.2101 ...
     AED94821             ;      338 ;     9.5700 ; 37685.5000 ;    -0.0562 ;    -0.8568 ...
     AED96183             ;      211 ;     4.2900 ; 23351.2000 ;    -0.1877 ;    -1.0654 ...
     ...

le début du tableau résumé des analyses que l'on veut obtenir


     Variable        p ANOVA    p KRUSKAL
     LENGTH     5.391201e-24   1.779741e-89
     PI         1.320983e-16   1.779741e-89
     MW         3.782924e-24   1.779741e-89
     FI         1.273242e-198  1.779741e-89
     GRAVY      2.158512e-231  1.779741e-89
     CHARGE     2.651195e-11   1.779741e-89
     HYDROPHI   2.368096e-157  1.779741e-89
     HYDROPHO   1.145593e-210  1.779741e-89
     FLEXI      2.719533e-54   1.779741e-89
     BULKI      7.618586e-72   1.779741e-89
     BURIED     5.963329e-120  1.779741e-89
     ACCESS     9.461772e-93   1.779741e-89
     ...

et le tableau désiré des médianes avec leur interprétations :


            LENGTH    PI       MW       FI    GRAVY   CHARGE HYDROPHI HYDROPHO   FLEXI    BULKI  BURIED  ACCESS   TRANSM MW.LENGTH
     Pool1   209.5 5.185 23737.33 -0.32370 -1.59640 -0.03640  1.12865 -0.35960 0.46365 14.25390 4.58020 6.64730 -1.50985  113.2551
     Pool2   151.0 7.250 16497.14  0.05990 -0.46520  0.00200  0.34690 -0.10550 0.43980 14.39760 6.23110 6.26940 -0.74830  107.0611
     Pool3   168.0 6.950 19028.07 -0.23160 -1.25980  0.00030  0.52810 -0.28750 0.46110 13.80100 4.99080 6.30480 -1.09600  112.5941
     WHy     101.0 5.040 10934.49  0.24415  0.07425 -0.01880 -0.06240  0.12275 0.44665 15.13515 6.81935 5.70590 -0.31495  108.6133
     Class8  183.0 8.190 20858.55  0.24800  0.07255  0.01080 -0.05565  0.08750 0.44245 15.17150 6.54400 5.71060 -0.33710  109.9836
     IDP     192.5 6.915 21291.82 -0.14625 -1.10885 -0.00335  0.53445 -0.17925 0.44725 14.12460 5.32515 6.02880 -0.99945  111.2180
     FS      161.0 6.110 17772.35  0.12730 -0.28835 -0.01195 -0.00600  0.03350 0.44115 14.88665 6.12510 5.65365 -0.47205  110.5061
     global  157.5 6.380 17013.33  0.15535 -0.16695 -0.00710  0.00465  0.05410 0.44380 14.86335 6.23540 5.74650 -0.44615  110.0492
     
            LENGTH PI MW FI GRAVY CHARGE HYDROPHI HYDROPHO FLEXI BULKI BURIED ACCESS TRANSM MW.LENGTH
     Pool1  +1     -1 +1 -1 -1    -1     +1       -1        0     0    -1     +1     -1      0
     Pool2   0     +1  0 -1 -1    +1     +1       -1        0     0     0     +1     -1      0
     Pool3  +1     +1 +1 -1 -1    +1     +1       -1        0    -1    -1     +1     -1      0
     WHy    -1     -1 -1 +1 +1    -1     -1       +1        0     0    +1      0     +1      0
     Class8 +1     +1 +1 +1 +1    +1     -1       +1        0     0     0      0     +1      0
     IDP    +1     +1 +1 -1 -1    +1     +1       -1        0     0    -1      0     -1      0
     FS      0      0  0 -1 -1    -1     -1       -1        0     0     0      0     -1      0
     relative tolerance:  5  %

Données fournies par Emmanuel JASPARD pour un article sur des protéines LEA.

Solution : afficher la solution

Analyse et explications

Il y a clairement deux parties distinctes la production de graphiques et la production de tableaux statistiques résumés. Vouloir les mélanger dans un seul programme ou fonction n'a de sens que si cela ne dure pas longtemps (ce qui est le cas ici). La production de graphiques PNG et PS pour intégration dans LaTeX ne pose pas de problème particulier, si ce n'est que la faible dispersion de certaines variables empêche de tracer certains beanplot. Il faut donc construire une liste de ces variables pour éviter que le programme ne s'arrête en plein milieu d'exécution.

Une seule boucle de traitement est donc certainement suffisante pour tout gérer, à condition de préparer des matrices de résultats qu'on viendra remplir au fur et à mesure.


     #  (gH)   -_-  art07.r  ;  TimeStamp (unix) : 01 Juillet 2013 vers 14:28
     
     library(beanplot)      # pour beanplot, forcément !
     library(stringr)       # pour str_locate
     load("article3.Rdata") # le df est art3
     
     source("art07_inc.r",encoding="latin1") # qui contient la fonction decritQTparFacteurArt3
     
     lesGroupes <- c("Pool1","Pool2","Pool3","WHy","Class8","IDP","FS")
     nbg        <- length(lesGroupes)
     
     pdv  <- 2:86  # col 1 : identifiant de la protéine ; col 87 : numéro de groupe
     nomc <- colnames(art3)
     
     vars              <- 2:15   # tableau des médianes,
     nblm              <- 1+nbg  # seulement 14 variables
     
     #################################################################
     
     # Propriétés des séquences d'acides aminés issues du site PROTPC
     #
     #    http://forge.info.univ-angers.fr/~gh/wstat/Protpc/
     #
     
     # on analyse les colonnes 02 à 15, de LENGTH à MW.LENGTH
     # on analyse les colonnes 16 à 35, de CNTA à CNTY
     # on analyse les colonnes 36 à 55, de pctA à pctY
     # on analyse les colonnes 56 à 75, de pctA.Unip à pctY.Unip
     # on analyse les colonnes 76 à 86, de D.E à C.F.Y.W
     
     #################################################################
     
     # liste des variables pour lesquelles beanplot n'est pas possible
     
     nobp <- c("LENGTH","cntC","cntW","pctC","pctW","pctC.Unip","pctW.Unip")
     
     # table de correspondance entre nom de colonne R et vrai nom de colonne
     
     vraiNom <-c(
        "MW.LENGTH",    "MW/LENGTH",
        "C.W",          "C+W",
        "D.E",          "D+E",
        "K.R",          "K+R",
        "N.Q",          "N+Q",
        "S.T",          "S+T",
        "F.W.Y",        "F+W+Y",
        "D.E.K.R",      "D+E+K+R",
        "D.E.K.R.1",    "D+E-K-R",
        "A.I.L.V",      "A+I+L+V",
        "R.E.S.P",      "R+E+S+P",
        "C.F.Y.W",      "C+F+Y+W"
     ) # fin de vraiNom
     
     # tableau des p-values par variable
     
     matRes            <- matrix(nrow=length(pdv),ncol=2) # matrice des résultats
     row.names(matRes) <- nomc[pdv]
     colnames(matRes)  <- c("p ANOVA","p KRUSKAL")
     
     # tableau des valeurs médianes globale et par groupe
     
     matMed            <- matrix(nrow=nblm,ncol=length(vars))
     row.names(matMed) <- c(lesGroupes,"global")
     colnames(matMed)  <- nomv[vars]
     
     # tableau des interprétations des médianes en +1/ 0/-1
     # on tolère une différence de eps=5  %
     
     matMedInt <- matMed[(1:(nblm-1)),]
     
     # boucle de traitement
     
     for (idv in pdv) {
     
       jdv   <- idv -1 # car on commence à 2
       laCol <- nomc[idv]
       laVar <- art3[,idv]
       bp    <- TRUE
       nbp   <- " "
       if (laCol %in% nobp) {
         bp <- FALSE
         nbp <- " (sans beanplot) "
       } # fin si
     
       nomOk <- laCol
       if (!is.na(str_locate(laCol,"Unip")[1])) { nomOk <- sub("\\.","/",laCol) }
       pdc <- which(laCol==vraiNom)
       if (length(pdc)>0) { nomOk <- vraiNom[ pdc + 1 ] }
       cats(paste("ANALYSE DE LA VARIABLE ",nomOk,nbp,sep=""))
     
       ba <- (-1)
       if (laCol=="FI") { ba <- 0 }
     
       nomgr <- paste("art3_",laCol,".png",sep="")
       gr(nomgr)
       res <- decritQTparFacteurArt3(laCol,laVar," ","group",art3$group,lesGroupes,TRUE,beanp=bp,barre=ba,vraiNom=nomOk)
       dev.off()
     
       nomps <- paste("art3_",laCol,".eps",sep="")
       postscript(nomps)
       res <- decritQTparFacteurArt3(laCol,laVar," ","group",art3$group,lesGroupes,TRUE,beanp=bp,barre=ba,vraiNom=nomOk,calc=FALSE)
       dev.off()
     
       matRes[jdv,] <- res # stockage des p-values
       cat(" ok pour ",nomgr," et ",nomps," variable ",nomOk,"\n")
     
       # gestion des médianes pour les premières variables seulement
     
       if (idv<=15) {
     
         globMed              <- median(laVar)
         matMed[nblm,jdv]     <- globMed
         meds                 <- tapply(X=laVar,INDEX=art3$group,FUN=median)
         matMed[(1:nbg),jdv]  <- meds
     
         for (idg in (1:nbg)) {
            dif <- (matMed[idg,jdv]-globMed)/abs(globMed)
            if (abs(dif)<eps) {
               matMedInt[idg,jdv] <- " 0" # noter l'espace devant zéro
            } else {
              if ((dif>= eps) & (dif>0)) { matMedInt[idg,jdv] <- "+1" }
              if ((dif<=-eps) & (dif<0)) { matMedInt[idg,jdv] <- "-1" }
            } # fin si
            # cat("groupe ",idg,"global med ",matMed[nblm,idv]," groupe med ",matMed[idg,idv]," rel dif ",dif,"\n")
         } # fin pour idg
     
       } # fin si
     
     } # fin pour idv
     
     # affichage de la matrice des résultats pour les p-values
     
     print(matRes)
     
     # affichage de la matrice des médianes
     
     print(matMed)
     
     # affichage de la matrice des interprétations des médianes
     
     print(matMedInt,quote=FALSE)
     cat("relative tolerance: ",eps*100," % \n")

Listing numéroté : protppc_7grps dont le résutat d'exécution simplifié est protpc_v2sor.txt.

Listing numéroté du fichier inclus : protpcv2_inc.

PROGRAMMATION R AVANCEE :

EXERCICES LONGS

Table des matières cliquable

1. Comptages de sites en génomique

Solution : afficher la solution

Analyse et explications

2. Production automatique de graphiques en protéinique

Solution : afficher la solution

Analyse et explications