PROGRAMMATION R AVANCEE : 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.   Table des matières cliquable   1. Comptages de sites en génomique   2. Production automatique de graphiques en protéinique   Il est possible d'afficher toutes les solutions via ?solutions=1.

 

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 :     masquer 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 :     masquer 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.

 

 

   Retour à la page principale de   (gH)