Introduction au logiciel d’analyse statistique R

• Pour installer R aller sur le CRAN.
• Pour lancer R sous Windows, cliquez, sur le raccourci présent sur le bureau.
• Sous Linux/Mac OSX, tapez “R” dans un terminal.
• On conseille l'utilisation de RStudio pour un environnement plus convivial .

Premiers pas

• R est un langage basé sur des objets
  • vector, matrix, data.frame, list, factor…
• Nous reviendrons sur cette notion importante plus tard.

• Pour commencer, créons un objet simple de type “vecteur” qui contiendra une valeur numérique.

x <- 15
print(x) 

[1] 15

show(x)

[1] 15

x

[1] 15

• L’opérateur d’assignation = peut se substituer à l’operateur <-
• Cependant l’opérateur <- est très recommandé.
• Le code pourra être commenté a l’aide du caractère #.
  • Toute commande à la suite de ce caractère ne sera pas interprétée.

x = 23
# x <- 32
x

[1] 23

Les instructions peuvent être séparées par un retour à la ligne ou par le caractère ;.

x <- 12; y <- 13
cat(x) # print(x)

12

cat(y) # print(y)

13

• Les objets créés sont stockés dans la mémoire vive de l’ordinateur
  • Ils ne sont pas, à ce stade, stockés sur le disque).
• On peut lister les objets disponibles en mémoire avec la fonction ls.
• Si vous quittez R ces objets seront détruits.
• Une autre méthode pour les détruire consiste à utiliser la fonction rm (remove).

ls()

[1] "x" "y"

rm(list=ls())

• Nous venons de construire des vecteurs (“vector” pour R) contenant des données numériques.
• Nous avons, par ailleurs utilisé les fonctions ls et rm.
• Dans R, on peut appeler des fonctions pour interroger des objets et réaliser des actions à partir de ceux-ci.
  • Les fonctions se présentent sous la forme suivante:

# Arg1 et arg2 (...) correspondent au noms des arguments
# a et b correspondent aux objets que l’on souhaite passer à la fonction

nomdelafonction(arg1= a, arg2= b,...)

• Ex: utilisation de la fonction c() (combiner).
  • Combiner des valeurs (e.g numériques) dans un vecteur (vector)

y <- c(-50, 1:8, 70) ; print(y)

 [1] -50   1   2   3   4   5   6   7   8  70

length(y)

[1] 10

is(y)

[1] "numeric" "vector" 

• Ex: utilisation de la fonction mean
  • Calculer d'une moyenne tronquée
  • Notez que les noms des l’arguments peuvent être abbrégés (si abréviation non ambigue) voire omis (si dans l'ordre attendus).

mean(y)

[1] 5.6

mean(y, trim = 0.2) # 10% deleted at each ends

[1] 4.5

mean(y, tr = 0.2)

[1] 4.5

mean(y, 0.2)

[1] 4.5

• Quatre modes: “numérique”“, "caractère”“, "logique”“,"complexe”.
• On peut stocker ces variables dans des vecteurs en utilisant la fonction c (combine) .
  • Attention, un vecteur n’accepte qu’un seul type de mode.

noms <-  c("celine", "alain", "robert") 
class(noms) # is(noms) 

[1] "character"

logic <- c(TRUE, FALSE, TRUE) # logic <- c(T, F, T)
class(logic) # is(logic)

[1] "logical"

• La fonction help().
  • Elle peut aussi être appelée à l’aide du caractère “?”.
  • Elle permet d’appeler l’aide sur une fonction donnée.

help(mean) # ?mean

• La fonction help.search().
  • Elle permet de chercher un terme donné dans les fichiers d’aide des fonctions.
• La fonction apropos().
  • Renvoie l’ensemble des fonctions dont le nom contient la chaîne de caractères recherchée

• Les fonctions c(), rep() et seq().

x <- 1:10 # une suite d'entiers
y <- c(2, 5, 7) # combiner 2 5 et 7 dans un vecteur
z <- rep(x, times = 5) # répéter x cinq fois. 

• On peut utiliser par exemple:

  • rnorm(): random normal.
  • runif(): random uniform.
  • rpois: random poisson.
  • rhyper: random hypergeometric.
  • …

par(mfrow=c(2, 2))

x <- rnorm(n=1000, mean=10, sd=2)
hist(x, main="Distribution normale (m=10, sd=2)")

x <- runif(n=1000, min = 10, max=20)
hist(x, main="Distribution uniforme (min=10, max=20)")

x <- rpois(1000, lambda=10)
hist(x, main="distribution de Poisson (lambda=10)")

x <- rhyper(nn=1000, m=100, n=100, k = 20)
hist(x, main="Distribution hypergéométrique (m=100, n=100, k = 20)")

par(mfrow=c(1, 1))

• On pourra faire des opérations simples de tri ou de randomisation sur les vecteurs.
• Notez que, dans l’exemple suivant, quand l’argument replace est positionné sur TRUE, le tirage aléatoire est effectué avec remise.

x <- 1:10 
sort(x, decreasing = TRUE)

 [1] 10  9  8  7  6  5  4  3  2  1

sample(x)

 [1] 10  6  8  3  2  4  5  7  1  9

sample(x, replace = TRUE)

 [1]  6  8  4  3 10  2  8  2  5  3

Pour interroger les éléments d’un vecteur on utilisera:

1. Un vecteurs d’indices (ie; de positions)
  * les positions sont notées de 1 à n (taille du vecteur).
2. Un vecteur logique 
3. les noms des éléments du vecteur.

• Exemple:

x <- c(-2, 4, -5, 1, 7)
x[3]  # position 3

[1] -5

x[3:5] # position 3 à 5

[1] -5  1  7

x[c(2,5)] # position 2 et 5

[1] 4 7

x[-3] # x sans la position 4 (-5)

[1] -2  4  1  7

• Exemple:

x <- c(-2, 4, -5, 1, 7)
x < 4

[1]  TRUE FALSE  TRUE  TRUE FALSE

x[x < 4] # renvoie les positions pour lesquelles le test est vrai

[1] -2 -5  1

x[x < 4 & x > 0] # & correspond au et logique

[1] 1

x[x < 4 | x ==  7] # | correspond au ou logique, == à l'opérateur d'égalité

[1] -2 -5  1  7

• Exemple:

x <- c(-2, 4, -5, 1, 7)
names(x) <- letters[1:length(x)]
x

 a  b  c  d  e 
-2  4 -5  1  7 

x["a"]

 a 
-2 

x[c("a", "d")]

 a  d 
-2  1 

• Fonction read.table

  • Cette commande est incontournable
  • Permettra la lecture de tableaux de données.
    • locaux ou distants (htttp, ftp)
  • Cette fonction renvoie un objet de type data.frame.
  • Contient de nombreux arguments qui permettent de paramètrer la lecture des données.
  • Tapez ?read.table pour avoir accès à l’aide sur cette fonction.

Fonction read.table, arguments importants

• file: le nom du fichier
• header: la première ligne correspond aux noms des colonnes.
• skip: Passer les n premières lignes avant la lecture.
• sep: le type de séparateurs de colonnes (e.g “\t”)
• row.names: la colonne contenant les noms des lignes (e.g, 1)
• quote: le délimiter de champs (à positionner plutôt sur “”)
• comment.char: par défaut “#”. Le texte précédé de ce caractère n'est pas lu. A utiliser pour des lignes de commentaires.

• Exercice:

  • Lisez le fichier à l'adresse suivante à l'aide de la fonction read.table. http://dputhier.github.io/STAT1_R/data/yeast_up500_6nt-1str-ovlp_MET_vs_all_sub.tsv

base.dir <- "http://dputhier.github.io/STAT1_R/data"
file.name <- "yeast_up500_6nt-1str-ovlp_MET_vs_all_sub.tsv"
yeast.df <- read.table(file.path(base.dir,file.name), 
                sep="\t", 
                header=T, 
                row=1)
class(yeast.df)

[1] "data.frame"

• C'est une matrice dont chacune des colonnes accepte un mode différent (character, numeric, boolean,…)
• On peut indexer un data.frame en interrogeant les lignes et colonnes

mon.data.frame[ligne, colonne]

• On pourra utiliser l'indexation par nom, par position ou l'indexation logique.
• Pour les noms de colonne, on pourra utiliser l'indexation à l'aide de l'opérateur dollar ($).

yeast.df[1,] # première ligne
yeast.df[1:10, 1:2] # 10 premières lignes des colonnes 1 et 2
yeast.df[c(1,3), 1] # lignes 1 et 3 de la colonne 1
yeast.df[ ,"MET_occ"] # la colonne MET_occ
yeast.df$MET_occ # la colonne MET_occ
yeast.df$MET_occ[1:3] # les trois premiers éléments de la colonne MET_occ

Exemple: le k-mère W = TTCACG

• \( N_{W,all}=424 \) Nombre d'occurrences d'un k-mère dans l'ensemble des séquences promotrices (colonne “all_occ”).

• \( N_{W,MET}=7 \) Nombre d'occurrences dans les gènes de la méthionine (colonne “MET_occ”“)

• \( N_{tous,all}=2421832 \) Nombre total de mots comptés dans tous les promoteurs

• \( N_{tous,MET}=9350 \) Nombre total de mots comptés dans les promoteurs de méthionine

• Demander le l'aide sur la fonction sum().
• Créez une variable sum_all et sum_met contenant le nombre total d'occurrences des k-mers dans les promoteurs de classe MET ou ALL.

sum_met <- sum(yeast.df$MET_occ)
print(sum_met)

[1] 9350

sum_all <-  sum(yeast.df$all_occ)
print(sum_all)

[1] 2421832

• Les boucles for
• La syntaxe des boucles for est la suivante:

for( i in vec){
  faireQuelqueChose(i)
  ...
} 

• Considérons l’addition de deux vecteurs

  x <- 1:10 
  y <- -x 
  z <- vector() 

  for (i in 1:length(x)){
    z[i] <- x[i] + y[i]
  } 

• Dans R, les boucles peuvent de manière générale être évitées
• Caractéristique du langage R
  • la vectorisation. La structure vectorielle rend les boucles implicites dans les expressions.

z <- x + y
z

 [1] 0 0 0 0 0 0 0 0 0 0

# Autres exemples
x/2

 [1] 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0

x^2

 [1]   1   4   9  16  25  36  49  64  81 100

• Définissons par \( P(W|all) \), ou \( p_W \), la probabilité d'observer un mot (\( W \)) donné à une position donnée, dans l'ensemble des promoteurs (\( all \)). On appelle cette probabilité probabilité a priori (prior).

• Pour le mot “TTCACG”, on obtient

  \( P(W|all) = p_W = \frac{N_{W,all}}{N_{tous,all}} = \frac{424}{2421832} = 0.000175 \)

• En vous basant sur ce que nous venons de voir, calculez la fréquence de chaque mots dans la classe MET et la class ALL (variables MET_freq et all_freq).

• Ajouter ces variables comme nouvelles colonnes au data.frame

• Testez que le résultat obtenu est conforme pour le k-mer “TTCACG”

MET_freq <- yeast.df$MET_occ / sum_met
all_freq <- yeast.df$all_occ / sum_all
yeast.df$MET_freq <- MET_freq
yeast.df$all_freq <- all_freq
head(yeast.df, n=2)

       MET_occ MET_mseq all_occ all_mseq     MET_freq     all_freq
CACGTG      16       14     317      285 0.0017112299 1.308926e-04
CGTGAC       6        6     197      192 0.0006417112 8.134338e-05

yeast.df["TTCACG","all_freq"]

[1] 0.0001750741

round(yeast.df["TTCACG","all_freq"], 6) == 0.000175

[1] TRUE

• On peut calculer le nombre d'occurrences attendues de la manière suivante:

  • $x = P(W|all) * $N_{tous,MET}$
  • I.e probabilité d'occurrence du mot dans all * Nombre de mots dans les promoteurs MET.
  • Créez une nouvelle colonne (MET_exp_occ) dans le data.frame contenant le nombre attendu d'occurrences.

yeast.df$MET_exp_occ <- yeast.df$all_freq * sum_met

• R dispose de librairies graphiques puissantes.
• On distingue les fonctions graphiques
  • primaires
    • permettant de créer un diagramme
    • e.g plot, hist, barplot, boxplot, …
  • secondaires
    • e.g points, text, identify, legend, …
    • permettant d'ajouter des éléments au diagramme

• Une version plus élégante

  plot(yeast.df$all_occ, yeast.df$MET_occ,
       pch=16, col="blue", panel.first=grid())

• On peut déjà identifier quelques k-mers remarquables.
• Afin d'ajouter leurs noms au graphique on peut utiliser text() ou identify()

# C'est une fonction graphique secondaire
# le graphique doit être généré au préalable

identify(x=yeast.df$all_occ, 
        y=yeast.df$MET_occ,
        labels=rownames(yeast.df),
        cex=0.7,
        col="darkgray")

• Faites un graphique dans lequel vous représenterez les occurrences observées (MET_occ) versus les occurrences attendues (MET_exp_occ) dans les promoteurs MET.

• Calculez le ratio des fréquences observées sur attendues.
• Ajoutez les ratio au data.frame dans une colonne occ_ratio.
• Transformez ces ratio en logarithme base 2.
• Ajoutez une colonne log2_occ_ratio dans le data.frame.
• représentez sur le graphique les points ayant une valeur absolue de log ratio supérieur à 1.

# Attention il y a des valeurs à Inf...
yeast.df$occ_ratio <- yeast.df$MET_occ / yeast.df$MET_exp_occ
yeast.df$log_occ_ratio <- log2(yeast.df$occ_ratio)

# Les points à mettre en exergue
to_point <- abs(yeast.df$log_occ_ratio) > 1

# Eventuellement
# to_point <- abs(yeast.df$log_occ_ratio) > 1 & is.finite(yeast.df$log_occ_ratio)
points(x=yeast.df$MET_exp_occ[to_point], 
       y=yeast.df$MET_occ[to_point], pch=1, lwd=2, col="green"
)
abline(a=0, b=0.5, col="red")
abline(a=0, b=2, col="red" )

• La loi de Poisson est une une loi de probabilité discrète qui ne dépend que d'un seul paramètre, \( \lambda \).

• Dans notre cas, \( \lambda \) correspond au d’occurrences attendues au hasard dans les promoteurs de la méthionine.

• On va chercher ici à apprécier le caractère remarquable du nombre d'occurrences observées.

• On dispose des informations suivantes sur ce k-mer:

# On index le data.frame par le nom du k-mer 
# On transpose avec t() pour faciliter la lecture
t(yeast.df["TTCACG",])

                    TTCACG
MET_occ       7.000000e+00
MET_mseq      6.000000e+00
all_occ       4.240000e+02
all_mseq      4.070000e+02
MET_freq      7.486631e-04
all_freq      1.750741e-04
MET_exp_occ   1.636943e+00
occ_ratio     4.276265e+00
log_occ_ratio 2.096351e+00

• Calculons la probabilité d'observer 7 occurrences:

# Surtout pas, mais c'est l'idée
sum(dpois(7:10000000, lambda=lambda))

[1] 0.001518595

# une autre solution, plus satisfaisante
1 - sum(dpois(0:(7-1), lambda=lambda))

[1] 0.001518595

# Une dernière solution (la plus satisfaisante)
ppois(q=7-1, lamb = lambda, lower = FALSE)

[1] 0.001518595

* Calculez la probabilité d'observer 4 occurrences
* Représentez le résultat graphiquement

ppois(q=4-1, lamb = lambda, lower = FALSE)
plot(0:15, dpois(0:15, lambda = lambda), type="h", col="blue", lwd=4)
points(4:15, dpois(4:15, lambda=lambda), type="h", lwd=4, col="red")

• Calculez la p-valeur de la Poisson pour tous les mots en utilisant une boucle for.
  • Notez qu'on pourrait aussi utiliser la fonction apply() en créant une petite fonction.

pois.p.val <- vector()
for(i in 1:nrow(yeast.df)){
  lambda <- yeast.df[i, "MET_exp_occ"]
  obs <- yeast.df[i, "MET_occ"]
  pois.p.val[i] <- ppois(q=obs-1, lamb = lambda, lower = FALSE)
}

yeast.df$pois.p.val <- pois.p.val
head(yeast.df[order(yeast.df$pois.p.val), c("log_occ_ratio","pois.p.val")], n=3)

       log_occ_ratio   pois.p.val
CACGTG      3.708578 3.835354e-13
TATATA      1.348107 1.959099e-08
ATATAT      1.306190 1.633240e-07

• Representez le log-ratio en fonction de la p-valeur calculée.

• On peut modéliser le problème à l'aide d'un test hypergéométrique.
• Pour chaque k-mer on considère:
  • \( N \) le nombre de mots (la taille de l'urne, le nombre de promoteurs).
  • Un tirage de taille \( k=19 \) (gènes régulés par la méthionine) éléments dans une urne contenant:
    • m boules blanches (marquées).
    • i.e : promoteurs contenant l'hexamère.
    • n boules noires (non marquées).
    • i.e : promoteurs ne contenant le l'hexamère.
• On cherche à connaître la probabilité d'obtenir \( X=x \) boules marquées.

# Pour le hexamer TTCACG
N <- 5880
m <- 407 # yeast$all_mseq 
n <- N -m
k <- 19
x <- 6 # yeast$MET_mseq

# La p-valeur est donnée par:
phyper(q=x-1, m = m, n=n, k=k, lower.tail = FALSE)

[1] 0.001324131

• Calculez la p-valeur du test hypergéométrique pour tous les gènes.

N <- 5880
k <- 19

hyper.p.val <- vector()
for(i in 1:nrow(yeast.df)){
  m <- yeast.df[i, "all_mseq"]
  n <- N - m
  x <-  yeast.df[i, "MET_mseq"]
  hyper.p.val[i] <- phyper(q=x-1, m = m, n=n, k=k, lower.tail = FALSE)
}

yeast.df$hyper.p.val <- hyper.p.val