# (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")
|
Retour à la page principale de
(gH)