Introduction au logiciel d’analyse statistique R

Denis Puthier
'2015-10-07'

Installer 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.

Création d'un objet.

  • 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

Opérateur d'assignation

  • 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

Plusieurs instructions sur une ligne.

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 sont stockés en mémoire

  • 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())

Fonctions. Informations de base sur les objets

  • 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,...)

Exemple d'utilisation de fonctions

  • 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" 

Exemple d'utilisation de fonctions

  • 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

Les modes

  • 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"

L’aide dans R

  • 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 objets de type vecteur : création

  • 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. 

Création de vecteurs avec des valeurs pseudo-aléatoires.

  • On peut utiliser par exemple:

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

Exemple: générer des distributions

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))

Exemple: générer des distributions

plot of chunk unnamed-chunk-11

Opérations de tri et de randomisation sur les vecteurs

  • 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

Indexation des vecteurs

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.

Indexation des vecteurs (position)

  • 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

Indexation des vecteurs (logique)

  • 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

Indexation des vecteurs (noms)

  • 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 

Lire un fichier

  • 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.

Lire un fichier

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.

Lire un fichier

Lire un fichier

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"

L'objet 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 ($).

Indexation de L'objet data.frame

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

Notation pour la suite.

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

Calculer le nombre d'occurrences des k-mers

  • 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.

Calculer le nombre d'occurrences des k-mers

sum_met <- sum(yeast.df$MET_occ)
print(sum_met)
[1] 9350
sum_all <-  sum(yeast.df$all_occ)
print(sum_all)
[1] 2421832

Boucles for et vectorisation

  • 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]
  } 

La vectorisation

  • 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

Calcul de la fréquence de chaque mots dans la classe MET et la class ALL

  • 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”

Calcul de la fréquence de chaque mots dans la classe MET et la class ALL

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

Calcul du nombre attendu d'occurrences dans les promoteurs MET

  • 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.

Calcul du nombre attendu d'occurrences dans les promoteurs MET

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

Diagramme avec R

  • 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

Exemple

  • Representons le nombre d'occurence dans all et le nombre d'occurrence dans MET.
  plot(yeast.df$all_occ, yeast.df$MET_occ)

plot of chunk unnamed-chunk-23

Exemple (amélioré)

  • Une version plus élégante
  plot(yeast.df$all_occ, yeast.df$MET_occ,
       pch=16, col="blue", panel.first=grid())

Exemple (amélioré)

  • Une version plus élégante

plot of chunk unnamed-chunk-25

Exemple (amélioré encore)

  • Une version un peu plus élégante
  col.palette <- densCols(yeast.df$all_occ, yeast.df$MET_occ)
  plot(x=yeast.df$all_occ, 
       y=yeast.df$MET_occ,
       pch=16,
       col=col.palette,
       panel.first=grid(),
       xlab="Nombre d'occurence dans l'ensemble des promoters",
       ylab="Nombre d'occurence dans les promoters MET",
       main="Nombre d'occurrences des k-mers"   
  )
  lm <- lm(yeast.df$MET_occ ~ yeast.df$all_occ) 
  abline(lm, col="red", lwd=2)