silico.biotoul.fr
 

M1 Traitement de Donnees Biologiques - Rice Expression Atlas Guided Tour

From silico.biotoul.fr

(Difference between revisions)
Jump to: navigation, search
m (Analyse en composantes principales)
m (Analyse différentielle)
(11 intermediate revisions not shown)
Line 136: Line 136:
-
Souvent lors d'une ACP, on affiche aussi les variables sur un ''biplot''. Vous pouvez essayer mais on peut s'attendre à ce que cela ne donne rien d'exploitable avec plus de 57k variables. Une solution alternative est de visualiser les variables qui contribuent le plus à une composante. Pour cela, nous allons triés les variables de manière décroissante (en valeur absolue) et faire un diagramme en barre :
+
Souvent lors d'une ACP, on affiche aussi les variables sur un ''biplot''. Vous pouvez essayer mais on peut s'attendre à ce que cela ne donne rien d'exploitable avec plus de 57k variables. Une solution alternative est de visualiser les variables qui contribuent le plus à une composante. Pour cela, nous allons trier les variables de manière décroissante (en valeur absolue) et faire un diagramme en barre :
<source lang='rsplus'>
<source lang='rsplus'>
probes.top = tibble(probeset=rownames(acp$rotation), PC1=acp$rotation[,1]) %>%    # création d'un tibble avec 2 colonnes : les identifiants de probeset et le premier vecteur propre
probes.top = tibble(probeset=rownames(acp$rotation), PC1=acp$rotation[,1]) %>%    # création d'un tibble avec 2 colonnes : les identifiants de probeset et le premier vecteur propre
-
   arrange(desc(abs(PC1))) %>%                                                    # trie décroissant par écart-type/variance
+
   arrange(dplyr::desc(abs(PC1))) %>%                                                    # trie décroissant par écart-type/variance
-
   mutate(probeset=fct_reorder(probeset, desc(abs(PC1)))) %>%                      # trie des identifiants de proebset avec le même ordre
+
   mutate(probeset=fct_reorder(probeset, dplyr::desc(abs(PC1)))) %>%                      # trie des identifiants de probeset avec le même ordre
   head(100)                                                                      # les 100 meilleures
   head(100)                                                                      # les 100 meilleures
probes.top
probes.top
Line 159: Line 159:
Dans cette partie, nous allons effectuer une analyse différentielle afin d'identifier les gènes dont l'expression est particulière dans un organe. Pour cela, nous allons utiliser la librairie <tt>limma</tt> qui a été développée tout d'abord pour l'analyse de données de microarray mais qui a été depuis adaptée pour des données RNA-Seq. Elle est très largement utilisée dans les pipelines d'analyse de données de transcriptome. Nous n'en verrons qu'une petite partie ensemble.
Dans cette partie, nous allons effectuer une analyse différentielle afin d'identifier les gènes dont l'expression est particulière dans un organe. Pour cela, nous allons utiliser la librairie <tt>limma</tt> qui a été développée tout d'abord pour l'analyse de données de microarray mais qui a été depuis adaptée pour des données RNA-Seq. Elle est très largement utilisée dans les pipelines d'analyse de données de transcriptome. Nous n'en verrons qu'une petite partie ensemble.
-
Il s'agit, pour cette analyse, de spécifier quelle colonnes de la matrice correspondent à quelles conditions expérimentales (design), puis de spécifier les conditions que l'on veut comparer (contasts). Ensuite, le traitement consiste à utiliser des modèles linéaires suivis de modèles bayésiens pour identifier ce qui diffère de manière statistiquement significative.
+
Il s'agit, pour cette analyse, de spécifier quelles colonnes de la matrice correspondent à quelles conditions expérimentales (''design''), puis de spécifier les conditions que l'on veut comparer (''contrasts''). Ensuite, le traitement consiste à utiliser des modèles linéaires suivis de modèles bayésiens pour identifier ce qui diffère de manière statistiquement significative.
Nous avons choisi pour ce TP d'étudier en particulier les conditions YL (Young Leaves). Chargement de la librairie et création de la variable qui contient le <tt>design</tt> :
Nous avons choisi pour ce TP d'étudier en particulier les conditions YL (Young Leaves). Chargement de la librairie et création de la variable qui contient le <tt>design</tt> :
Line 180: Line 180:
Analyse différentielle :
Analyse différentielle :
<source lang='rsplus'>
<source lang='rsplus'>
-
fit  = lmFit(t(dm), design)          # modèle linéaire: calcule la moyenne de P1 et P3 pour chaque probeset
+
fit  = lmFit(t(dm), design)          # modèle linéaire: calcule la moyenne de YL et other pour chaque probeset
-
fit2 = contrasts.fit(fit, contrasts) # calcul des contrastes (différence des moyennes P1 et P3)
+
fit2 = contrasts.fit(fit, contrasts) # calcul des contrastes (différence des moyennes YL et other)
efit = eBayes(fit2)                # modèle bayésien
efit = eBayes(fit2)                # modèle bayésien
all  = topTable(efit, adjust='BH', p.value = 1, number=60000)  # tous les résultats (60000 les plus significatifs) avec FDR (adjust='BH') sans seuil alpha (p.value=1)
all  = topTable(efit, adjust='BH', p.value = 1, number=60000)  # tous les résultats (60000 les plus significatifs) avec FDR (adjust='BH') sans seuil alpha (p.value=1)
Line 194: Line 194:
   ggplot(aes(x=logFC, y=-log10(adj.P.Val))) +  
   ggplot(aes(x=logFC, y=-log10(adj.P.Val))) +  
   geom_point(alpha=.2, size=.5) +
   geom_point(alpha=.2, size=.5) +
-
   annotate(geom='rect', xmin=min(all$logFC, na.rm=T), xmax=max(all$logFC, na.rm=T), ymin=0, ymax=-log10(0.05), fill='red', alpha=0.25) + # bellow p-value threshold
+
   annotate(geom='rect', xmin=min(all$logFC, na.rm=T), xmax=max(all$logFC, na.rm=T), ymin=0, ymax=-log10(0.05), fill='red', alpha=0.25) + # below p-value threshold
-
   annotate(geom='rect', xmin=-log2(10), xmax=log2(10), ymin=-log10(0.05), ymax=max(-log10(all$adj.P.Val)), fill='red', alpha=0.25) +    # bellow logFC threshold
+
   annotate(geom='rect', xmin=-log2(10), xmax=log2(10), ymin=-log10(0.05), ymax=max(-log10(all$adj.P.Val)), fill='red', alpha=0.25) +    # below logFC threshold
   ggtitle("Volcano plot YL vs. other")
   ggtitle("Volcano plot YL vs. other")
</source>
</source>
Line 206: Line 206:
             adjust='BH',  
             adjust='BH',  
             p.value=0.05,  
             p.value=0.05,  
-
             number=10000,  
+
             number=60000,  
             resort.by="logFC",  
             resort.by="logFC",  
             lfc=log2(10))
             lfc=log2(10))
Line 226: Line 226:
Création d'une matrice plus petite ne contenant que les valeurs d'expression des gènes différentiellement exprimés identifiés précédemment :
Création d'une matrice plus petite ne contenant que les valeurs d'expression des gènes différentiellement exprimés identifiés précédemment :
<source lang='rsplus'>
<source lang='rsplus'>
-
YL.probesets = rownames(deg.YL)
+
YL.dm = dm[ , deg.YL$probeset ]
-
YL.dm = dm[ , YL.probesets ]
+
rownames(YL.dm) = info$id
-
rownames(YL.dm) = YL.probesets
+
colnames(YL.dm) = deg.YL$probeset
 +
dim(YL.dm)
</source>
</source>
Line 237: Line 238:
<span style='color: #990000'>Commenter la heatmap obtenue.</span>
<span style='color: #990000'>Commenter la heatmap obtenue.</span>
 +
 +
= Caractérisation d'une liste de gènes =
 +
 +
Pour cette partie, nous allons nous concentrer sur les gènes sur-exprimés dans la condition YL afin d'étudier s'ils sont annotés comme impliqués dans un ou des processus biologiques particuliers.
 +
 +
Première étape : récupération des identifiants des probesets sur-exprimés :
 +
<source lang='rsplus'>
 +
deg.YL %>%
 +
  filter(logFC>0) %>%
 +
  write_tsv("deg.YL.up.tsv")
 +
</source>
 +
 +
La fonction <tt>filter</tt> permet de ne garder que les lignes qui satisfont les critères de sélection (ici un log fold change positif). Le résultat est enregistré dans un fichier TSV.
 +
 +
<span style='color: #990000'>Ouvrez le fichier avec LibreOffice afin de copier les identifiants des probesets pour les coller dans l'interface d'analyse d'enrichissement du site Web http://silico.biotoul.fr/enrichment/ pour réaliser ce type d'analyse.</span>
 +
 +
<span style='color: #990000'>Commenter les résultats obtenus.</span>
 +
 +
<span style='color: #990000'>Utiliser le bouton "Start Revigo" pour vous y aider.</span>
 +
 +
= Compte rendu de TP =
 +
 +
 +
<span style='color: #990000;'>Envoyer le compte rendu à vos enseignants par mail ([mailto:maxime.bonhomme@univ-tlse3.fr maxime.bonhomme@univ-tlse3.fr] et [mailto:roland.barriot@univ-tlse3.fr roland.barriot@univ-tlse3.fr]).
 +
 +
<span style='color: #990000;'>Le compte rendu est à envoyer '''avant une semaine'''.
 +
 +
<span style='color: #990000;'>Envoyez les 2 fichiers (.Rmd et .html). Envoyez-vous aussi le mail en copie pour pouvoir vérifier que tout est bien passé.
 +
 +
<span style='color: #990000;'>Mettez un titre tel que "Compte rendu TP5 TDB de -et vos Nom et Prénom-".
 +
 +
<span style='color: #990000;'>'''Votre compte rendu ne doit pas être un simple "copié-collé" des commandes. Nous attendons une mise en forme correcte (titres, sous-titres, paragraphes, ...), des commentaires (notamment sur les résultats des commandes), des essais de votre part pour améliorer par exemple la visualisation d'un graphique, ou pour rechercher des fonctions R alternatives ou complémentaires à ce qui est proposé ; bref une personnalisation de votre compte-rendu.'''</span>
 +
 +
= Annexes =
 +
 +
'''Liens :'''
 +
* Un livre entier (libre et gratuit) consacré à [https://r4ds.had.co.nz/ R for Data Science] avec des explications complètes depuis l'utilisation de RStudio jusqu'à des traitement plus élaborés, et entre autres les tibbles, l'import de données, les types de données et leur manipulation, les chaînes de caractères, les facteurs, ...
 +
 +
 +
Les <tt>cheat sheets</tt> faisant '''aides-mémoires''' :
 +
* [https://github.com/rstudio/cheatsheets/blob/master/data-visualization-2.1.pdf ggplot2]
 +
* [https://raw.githubusercontent.com/rstudio/cheatsheets/main/data-import.pdf readr] pour lire/écrire des fichiers
 +
* [https://github.com/rstudio/cheatsheets/blob/main/data-transformation.pdf dplyr]
 +
* [https://github.com/rstudio/cheatsheets/blob/master/tidyr.pdf tidyr]
 +
* [https://raw.githubusercontent.com/rstudio/cheatsheets/main/factors.pdf forcats] pour les facteurs
 +
* [https://github.com/rstudio/cheatsheets/blob/main/strings.pdf stringr] manipulation de texte
 +
 +
== Installation des librairies seulement si nécessaire ==
 +
<source lang='rsplus'>
 +
installed = rownames(installed.packages())
 +
cran.libs = c('tidyverse', 'factoextra','BiocManager')
 +
for (lib.name in cran.libs)
 +
  if (!lib.name %in% installed)
 +
    install.packages(lib.name)
 +
 +
bioconductor.libs = c('limma','made4')
 +
for (lib.name in bioconductor.libs)
 +
  if (!lib.name %in% installed)
 +
    BiocManager::install(lib.name)
 +
</source>
 +
 +
 +
== Variables contribuant le plus à un axe ==
 +
 +
Création du tibble et ajout des valeurs absolues
 +
<source lang='rsplus'>
 +
tmp1 = tibble(probeset=rownames(acp$rotation), PC1=acp$rotation[,1]) %>%
 +
  mutate(PC1.abs = abs(PC1))
 +
tmp1
 +
</source>
 +
 +
Tri par valeur décroissante en valeur absolue et ajout de l'indice
 +
<source lang='rsplus'>
 +
tmp2 = tmp1 %>%
 +
  arrange(desc(PC1.abs)) %>%
 +
  rowid_to_column
 +
tmp2
 +
</source>
 +
 +
Les 100 meilleures/premiers
 +
<source lang='rsplus'>
 +
tmp3 = tmp2 %>%
 +
  head(100)
 +
tmp3
 +
</source>
 +
 +
plot
 +
<source lang='rsplus'>
 +
tmp3 %>%
 +
  ggplot(aes(x=rowid, y=PC1)) +
 +
  geom_bar(stat = 'identity')
 +
</source>

Revision as of 13:57, 5 October 2022

Contents

Préparation de l'environnement

Au cours de ce TP, vous allez travailler sur un jeu de données d'expression du génome du riz dans différents tissus et organes au cours du développement. Il s'agit de données obtenues par hybridation sur des microarrays pour chacun des 57k probesets présents sur le microarray.

Il s'agit pour commencer de créer l'environnement de travail.

Créez tout d'abord un répertoire de travail sur le bureau (par exemple TDB-TP5) et commencez par télécharger le fichier source que vous allez utiliser et compléter pour générer le compte rendu de TP : M1.TDB.TP_Rice_Expression_Atlas.Rmd (click droit de la souris -- enregistrer la cible sous...). Ouvrez le logiciel RStudio et chargez ce fichier puis lancez sa compilation pour voir le compte rendu. Pour cela cliquez sur le bouton Knit HTML ou bien utilisez la combinaison de touches Ctrl + shift + K.


Récupération des données. Il s'agit ensuite de récupérer les fichiers de données. Téléchargez les fichiers suivant dans votre répertoire de travail :

Contexte

Dans ce TP, nous allons analyser les données de transciptome obtenues sur Oryza sativa par microarray. Le déroulement des analyses suit les travaux de Fujita et al. publiés en 2010 (pmid:21062870). Ces données sont publiques et ont été mises à disposition sur la banque GEO et sont accessibles avec l'identifiant GSE14304. Il s'agit de 98 hybridations réalisées avec un microarray de la société Affymetrix ayant l'identifiant GPL2025 comportant 57 381 spots (ou probesets) correspondant à 51 279 transcrits d'Oryza sativa japonica et indica.

Pour réaliser ce TP, les données brutes ont été prétraitées afin d'obtenir pour chaque probeset une valeur d'expression normalisée (= comparable entre une hybridation et une autre). Le traitement effectué suit les étapes décrites dans le paragraphe Microarray data extraction, processing and cluster analysis de la partie Materials and Methods :

"For Affymetrix array data, CEL files produced by GCOS 1.3 (Affymetrix, Inc.) were analyzed using the statistical software R with bioconductor package ‘affy’. Signal intensities were extracted by expresso algorithm with parameters: bgcorrect.method = ‘mas’, normalize = ‘F’, pmcorrect.method = ‘pmonly’, summary.method = ‘mas’. Extracted signal intensities were introduced into GeneSpring 7.3.1 (Agilent Technologies, Inc.) and scaled to the 75th percentile per chip."

Les objectifs de ce TP sont multiples et vont être :

  • Prise en main des librairies et fonctions nécessaires,
  • Déterminer si les résultats sont reproductibles : ACP sur les 98 conditions expérimentales → est-ce que les réplicats se regroupent ?
  • Réaliser une analyse d'expression différentielle afin de déterminer les gènes différentiellement exprimés dans une des conditions expérimentales
  • Faire une analyse de clustering sur les profils d'expression de ces gènes et visualiser le résultat
  • Caractériser la liste des gènes sur-exprimés dans la condition expérimentale étudiée

Librairies tidyverse, chargement et préparation des données

Pour ce TP, vous allez commencer à utiliser un ensemble de librairies de plus en plus utilisées en traitement de données (et data science). Les 2 principaux avantages sont que :

  • il y a énormément de documentation en ligne, dans différentes langues, sous forme de tutoriels, de forums, ..., notamment stackoverflow, et que même si la prise en main peut-être un peu ardue au tout début, il y a, le plus généralement, quand on cherche à faire un traitement, quelqu'un qui a eu les mêmes problèmes et donc, la solution disponible mais à adapter à votre jeu de données,
  • le code gagne en lisibilité, et donc plus facile à reprendre, améliorer et maintenir.

Pour ce TP, nous allons utiliser

  • la lecture de fichier,
  • les tibbles : un data.frame amélioré,
  • un opérateur particulier : le pipe %>% qui permet plus de lisibilité dans l’enchaînement des traitements effectué, le résultat d'une fonction étant passé à la fonction suivante,
  • les fonctions graphiques améliorés de la librairie ggplot2.

Les fichiers fournis sont au format TSV, c'est-à-dire que la 1ère ligne contient le nom des colonnes, et les colonnes sont séparées par des tabulations. Chargez le fichier contenant les informations relatives aux 98 hybridations, et affichez le contenu comme suit :

library(tidyverse)
info = read_tsv("rice.atlas.info.tsv")
info

La fonction read_tsv renvoie un tibble qui est stocké dans la variable info. L'affichage du contenu est plus lisible que pour un data.frame : on a la taille, les types des colonnes, et s'il y a beaucoup de lignes et de colonnes, la console n'est pas submergée par des centaines de lignes de données.

Adaptez la commande pour charger les données d'expression du fichier rice.atlas.expr.normalized.tsv dans une variable que vous pourrez appeler tb, et affichez son contenu. Combien de lignes et de colonnes contient-elle ?

La première colonne des 2 tibbles correspond à l'identifiant des hybridations et permet de faire le lien entre les 2 tables. Pour réaliser l'ACP et les autres traitements du TP, nous allons isoler les valeurs numériques d'expression dans une matrice (numérique) comme suit :

dm = tb %>% 
  select(-id) %>% 
  as.matrix

Cela nous permet d'utiliser le pipe %>% pour la première fois, explication :

  • tb %>% : a pour effet que le résultat de tb (normalement l'affichage) est transmis/redirigé à la commande suivante (de la ligne suivante)
  • select(-id) %>% : sur le tibble renvoyé par la ligne précédente (la totalité des données), on utilise la fonction select qui permet de sélectionner certaines colonnes dans un tibble (ou data.frame) ; ici, on sélectionne toutes les colonnes sauf id, et le résultat est passé à la ligne suivante avec le %>%
  • as.matrix : sur le tibble renvoyé par la ligne précedent, on appelle la fonction as.matrix qui a pour effet de transformer le tibble en matrix (tableau dont toutes les cases ont le même type).
  • Le résultat est donc une variable de type matrix et est stocké dans la variable dm (pour data matrix).


Premières explorations. Avec de telles dimensions, notamment le nombre de colonnes, il ne va pas être possible d'utiliser la fonction summary sur les données d'expression telles quelles. Il est néanmoins possible de transposer la matrice (consiste à inverser les lignes et les colonnes), afin d'avoir un summary par hybridation (et non par probeset) :

dm %>%
 t %>% # fonction transpose
 summary

Cela laisse néanmoins beaucoup d'informations à lire...

Faire un histogramme afin d'apprécier la distribution des valeurs d'expression.

Quel problème rencontrez-vous ?

Faire l'histogramme sur les données transformées en log2 (fonction log2)..


Distributions des valeurs d'expression par hybridation. Le résultat de l'effet de la normalisation effectuée pourra se visualiser à l'aide de boîtes à moustaches :

dm %>% 
  t %>% 
  log2 %>% 
  boxplot

Cela semble-t-il correspondre à ce qui est décrit dans le Materials and methods ?


Un premier plot avec ggplot. Pour finir cette partie, nous allons visualiser le niveau d'expression d'un probeset en particulier, rRNA 18S, à travers les différentes conditions et réplicats. ggplot prend en paramètre un data.frame ou tibble et permet de construire petit à petit la visualisation : il s'agit de préciser ce qui est affiché en x, y, ..., quelles transformations sont à effectuer (pour les intervalles d'un histogramme par exemple), sur quels axes avec éventuellement leur échelle (logarithmique, exponentielle), les couleurs, les formes et autres aspects esthétiques :

tb %>%                                                                 # le tibble
  ggplot(aes(y=AFFX.OS.18SrRNA_s_at, x=id)) +                          # l'appel à ggplot en précisant ce que l'on va "plotter" via aes (pour aesthetics) en x et en y
  geom_point(aes(color=info$stage)) +                                  # ajout de la couche avec les "points" en précisant la couleur
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))  # modification des légendes de l'axe des x pour les afficher verticalement

Analyse en composantes principales

La fonction prcomp fait partie de la librairie stats qui est normalement chargée au démarrage de R.

Pour réaliser une ACP, il suffit d'appeler cette fonction sur la matrice de données. Il faudra préciser s'il s'agit d'une ACP normée ou non (paramètre scale=TRUE ou scale=FALSE, cf. l'aide ?prcomp). Ici, les données ont été préalablement normalisées, donc nous allons nous orienter vers une ACP non normée :

acp = dm %>% 
  log2 %>% 
  prcomp

Commenter ce qu'affiche la fonction summary qui nous permet d'avoir la variance portée par les axes :

acp %>% summary

Pour visualiser cette information, il est de coutume de faire un screeplot. Pour cela, la librairie FactoExtra est toute adaptée :

library(factoextra)
fviz_screeplot(acp, addlabels=T)

Affichage (avec ggplot) des individus projetés (acp$x) sur les 2 premières composantes :

acp$x %>%
  as_tibble %>%
  ggplot(aes(x=PC1, y=PC2, label=info$stage)) + 
  ggtitle("ACP sur les différents stades de développement avec répétitions\nà partir des niveaux d'expression des 57k probesets") +
  geom_text(color=info$color) + 
  theme_light()

Commentez le résultat obtenu.

Effectuez la même chose avec les composantes PC1 et PC3.


Souvent lors d'une ACP, on affiche aussi les variables sur un biplot. Vous pouvez essayer mais on peut s'attendre à ce que cela ne donne rien d'exploitable avec plus de 57k variables. Une solution alternative est de visualiser les variables qui contribuent le plus à une composante. Pour cela, nous allons trier les variables de manière décroissante (en valeur absolue) et faire un diagramme en barre :

probes.top = tibble(probeset=rownames(acp$rotation), PC1=acp$rotation[,1]) %>%    # création d'un tibble avec 2 colonnes : les identifiants de probeset et le premier vecteur propre
  arrange(dplyr::desc(abs(PC1))) %>%                                                     # trie décroissant par écart-type/variance
  mutate(probeset=fct_reorder(probeset, dplyr::desc(abs(PC1)))) %>%                      # trie des identifiants de probeset avec le même ordre
  head(100)                                                                       # les 100 meilleures
probes.top

Diagramme en barres :

probes.top %>% 
  ggplot(aes(x=probeset, y=PC1)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=1))

Dans la partie suivante, nous allons sélectionner les gènes différentiellement exprimés dans un tissu particulier afin d'avoir un peu moins de données à manipuler.

Analyse différentielle

Dans cette partie, nous allons effectuer une analyse différentielle afin d'identifier les gènes dont l'expression est particulière dans un organe. Pour cela, nous allons utiliser la librairie limma qui a été développée tout d'abord pour l'analyse de données de microarray mais qui a été depuis adaptée pour des données RNA-Seq. Elle est très largement utilisée dans les pipelines d'analyse de données de transcriptome. Nous n'en verrons qu'une petite partie ensemble.

Il s'agit, pour cette analyse, de spécifier quelles colonnes de la matrice correspondent à quelles conditions expérimentales (design), puis de spécifier les conditions que l'on veut comparer (contrasts). Ensuite, le traitement consiste à utiliser des modèles linéaires suivis de modèles bayésiens pour identifier ce qui diffère de manière statistiquement significative.

Nous avons choisi pour ce TP d'étudier en particulier les conditions YL (Young Leaves). Chargement de la librairie et création de la variable qui contient le design :

library(limma)
design = tibble(YL=as.numeric(info$stage=='YL'), other=1-YL) %>% as.matrix
design %>% tail

On a des 1 pour la colonne YL seulement pour les 4 dernière lignes.

Remarque : on peut utiliser n'importe quelle méthode pour créer cette matrice, du moment qu'on aboutit à des 1 et des 0 qui indiquent quelles colonnes de la matrice qui va être analysée correspondent à quelle conditions expérimentales (ici, soit YL, soit other).

Définition des contrasts → qu'est-ce que l'on compare (ici YL - other) :

contrasts=makeContrasts(controlVSah=YL-other, levels=design)
contrasts

Analyse différentielle :

fit  = lmFit(t(dm), design)          # modèle linéaire: calcule la moyenne de YL et other pour chaque probeset
fit2 = contrasts.fit(fit, contrasts) # calcul des contrastes (différence des moyennes YL et other)
efit = eBayes(fit2)                 # modèle bayésien
all  = topTable(efit, adjust='BH', p.value = 1, number=60000)  # tous les résultats (60000 les plus significatifs) avec FDR (adjust='BH') sans seuil alpha (p.value=1)
all %>% head

Commenter la sortie de la dernière ligne.

L'ensemble des résultats obtenus est généralement illustré par un Volcano plot :

all %>% 
  ggplot(aes(x=logFC, y=-log10(adj.P.Val))) + 
  geom_point(alpha=.2, size=.5) +
  annotate(geom='rect', xmin=min(all$logFC, na.rm=T), xmax=max(all$logFC, na.rm=T), ymin=0, ymax=-log10(0.05), fill='red', alpha=0.25) + # below p-value threshold
  annotate(geom='rect', xmin=-log2(10), xmax=log2(10), ymin=-log10(0.05), ymax=max(-log10(all$adj.P.Val)), fill='red', alpha=0.25) +     # below logFC threshold
  ggtitle("Volcano plot YL vs. other")

Peut-on s'attendre à ce qu'il y ait beaucoup de gènes différentiellement exprimés ?

Pour cela, nous allons extraire avec la fonction topTable précédente les gènes qui ont une p-valeur ajustée (FDR) inférieure à 0.05 et dont l'expression varie beaucoup (au minimum +10x ou -10x) :

deg.YL=topTable(efit, 
             adjust='BH', 
             p.value=0.05, 
             number=60000, 
             resort.by="logFC", 
             lfc=log2(10))
deg.YL$probeset = rownames(deg.YL)
dim(deg.YL)

Combien de gènes différentiellement exprimés obtenez-vous ?

Dans la partie suivante, nous allons visualiser les profils de ces gènes à l'aide d'une heatmap.

Clustering

Il existe différentes librairies permettant de tracer des heatmap/heatplot, nous allons utiliser made4 :

library(made4)

Création d'une matrice plus petite ne contenant que les valeurs d'expression des gènes différentiellement exprimés identifiés précédemment :

YL.dm = dm[ , deg.YL$probeset ]
rownames(YL.dm) = info$id
colnames(YL.dm) = deg.YL$probeset
dim(YL.dm)

A présent, nous allons utiliser la fonction heatplot pour générer la heatmap. La fonction va réaliser un clustering hiérarchique des lignes et/ou des colonnes avec éventuellement une normalisation des lignes et/ou des colonnes. Nous allons utiliser la distance Ward entre clusters (cf. cours) et ajouter de la couleur pour les différentes conditions expérimentales (YL, ...) :

heatplot(log2(YL.dm), cols.default = F, lowcol = 'blue', highcol='yellow', dualScale = F, scale='column', method='ward.D2', RowSideColors=info$color)

Commenter la heatmap obtenue.

Caractérisation d'une liste de gènes

Pour cette partie, nous allons nous concentrer sur les gènes sur-exprimés dans la condition YL afin d'étudier s'ils sont annotés comme impliqués dans un ou des processus biologiques particuliers.

Première étape : récupération des identifiants des probesets sur-exprimés :

deg.YL %>% 
  filter(logFC>0) %>%
  write_tsv("deg.YL.up.tsv")

La fonction filter permet de ne garder que les lignes qui satisfont les critères de sélection (ici un log fold change positif). Le résultat est enregistré dans un fichier TSV.

Ouvrez le fichier avec LibreOffice afin de copier les identifiants des probesets pour les coller dans l'interface d'analyse d'enrichissement du site Web http://silico.biotoul.fr/enrichment/ pour réaliser ce type d'analyse.

Commenter les résultats obtenus.

Utiliser le bouton "Start Revigo" pour vous y aider.

Compte rendu de TP

Envoyer le compte rendu à vos enseignants par mail (maxime.bonhomme@univ-tlse3.fr et roland.barriot@univ-tlse3.fr).

Le compte rendu est à envoyer avant une semaine.

Envoyez les 2 fichiers (.Rmd et .html). Envoyez-vous aussi le mail en copie pour pouvoir vérifier que tout est bien passé.

Mettez un titre tel que "Compte rendu TP5 TDB de -et vos Nom et Prénom-".

Votre compte rendu ne doit pas être un simple "copié-collé" des commandes. Nous attendons une mise en forme correcte (titres, sous-titres, paragraphes, ...), des commentaires (notamment sur les résultats des commandes), des essais de votre part pour améliorer par exemple la visualisation d'un graphique, ou pour rechercher des fonctions R alternatives ou complémentaires à ce qui est proposé ; bref une personnalisation de votre compte-rendu.

Annexes

Liens :

  • Un livre entier (libre et gratuit) consacré à R for Data Science avec des explications complètes depuis l'utilisation de RStudio jusqu'à des traitement plus élaborés, et entre autres les tibbles, l'import de données, les types de données et leur manipulation, les chaînes de caractères, les facteurs, ...


Les cheat sheets faisant aides-mémoires :

Installation des librairies seulement si nécessaire

installed = rownames(installed.packages())
cran.libs = c('tidyverse', 'factoextra','BiocManager')
for (lib.name in cran.libs)
  if (!lib.name %in% installed)
    install.packages(lib.name)
 
bioconductor.libs = c('limma','made4')
for (lib.name in bioconductor.libs)
  if (!lib.name %in% installed)
    BiocManager::install(lib.name)


Variables contribuant le plus à un axe

Création du tibble et ajout des valeurs absolues

tmp1 = tibble(probeset=rownames(acp$rotation), PC1=acp$rotation[,1]) %>% 
  mutate(PC1.abs = abs(PC1))
tmp1

Tri par valeur décroissante en valeur absolue et ajout de l'indice

tmp2 = tmp1 %>% 
  arrange(desc(PC1.abs)) %>% 
  rowid_to_column
tmp2

Les 100 meilleures/premiers

tmp3 = tmp2 %>% 
  head(100)
tmp3

plot

tmp3 %>% 
  ggplot(aes(x=rowid, y=PC1)) +
  geom_bar(stat = 'identity')