R possède une aide pour chaque commande. Cette aide fournit des indications sur le calcul que fait R ainsi que sur les paramètres qui doivent figurer dans la commande pour qu'elle soit exécutée correctement. De plus, elle donne une description de la méthode et une liste de références. Il est donc très utile. Tapez R> help() pour obtenir l'aide de l'aide ou R> help(nom_de_la_commande), pour une commande spécifique.
Il est également possible de rechercher dans l'aide avec la fonction help.search('mots à rechercher').
La cryopréservation est une méthode largement utilisée pour le stockage à long terme de nombreuses cellules vivantes. Cette méthode implique des traitements de congélation et de décongélation qui causent des dommages aux cellules vivantes voire souvent à leur mort. Les puces à ADN permettent de mesurer les niveaux d'expression de milliers de gènes simultanément et de suivre les réponses biologiques à travers le niveau d'expression de presque tous les gènes de l'organisme. Nous allons analyser la réponse de cellules de levure après cryopréservation à partir de données de transcriptome publiées dans une étude de l'impact de la crypréservation sur la levure [Odani et al. 2003].
Il existe plusieurs entrepôts de données pour les données de transcriptome. Les principaux sont Gene Expression Omnibus (GEO) du NCBI, ArrayExpress de l'EBI, ainsi que le Stanford Microarray Database (SMD). A partir de GEO du NCBI, retrouvez les données associées à la publication [Odani et al. 2003].
Vous devriez trouver la série GSE9404. Lisez la page décrivant les données, puis téléchargez les données brutes (GSM239212..GSM239220) dans un répertoire que vous aurez créé pour nos analyses. Ensuite, décompressez ces fichiers (soit avec la commande gunzip sous linux, soit avec un logiciel de compression/décompression tel que 7-Zip.
Le package limma est spécialement conçu pour analyser les données de biopuces bi-couleurs. Il permet l'importation de fichiers de sortie des logiciels d'analyse d'image les plus courants dans le domaine des biopuces, la normalisation des données et l'analyse différentielle (pour une présentation détaillée, voir la page web http://bioinf.wehi.edu.au/limma/).
Le chargement du package est ensuite commandé par :
R> library(limma) # Chargement du package limmaRemarque : limma dispose d'un guide utilisateur accessible et pédagogique, incluant notamment de nombreux exemples d'utilisation. Pour y accéder :
R> limmaUsersGuide()Afin d'éviter d'avoir à indiquer de manière répétée le répertoire de travail au cours de la procédure d'importation des données, on peut le définir comme répertoire de travail pour l'ensemble de la session :
R> setwd("C:/My Documents/Transcriptome/") # sous Windows définit le répertoire C:\My Documents\Transcriptome comme répertoire de travail R> setwd("/home/myLogin/Transcriptome/") # sous linuxLe format des données de biopuces nécessaire à toute analyse statistique est un tableau dont les lignes sont associées à des puces et les colonnes à des gènes. Or, le plus souvent le point de départ de l'analyse repose sur autant de fichiers de sortie d'un logiciel d'analyse d'images qu'il y a de puces. Le package limma prévoit donc une procédure d'importation de ces fichiers conduisant à la création d'un tableau de données. La conversion de plusieurs fichiers en un seul tableau pucesxgènes nécessite quelques indications permettant de piloter l'ordonnancement des valeurs d'expressions géniques dans le jeu de données. Ces informations sont contenues dans un fichier supplémentaire nommé dans la suite targets.txt.
Le fichier targets.txt doit contenir les colonnes suivantes :
Les autres colonnes sont optionnelles. On donne le fichier targets_15min.txt suivant correspondant à la première série de réplicats effectués 15 minutes après décongélation.
| FileName | SlideNumber | Cy3 | Cy5 |
|---|---|---|---|
| GSM239212.gpr | 1 | ref | 15min |
| GSM239213.gpr | 2 | ref | 15min |
| GSM239214.gpr | 3 | ref | 15min |
La commande suivante permet de lire le fichier targets.txt et de créer un objet R, nommé ici targets, héritant de l'information contenue dans le fichier :
R> targets=readTargets("targets.txt",sep="\t") # Lecture de targets.txtL'argument sep définit le type de séparateur utilisé pour délimiter les colonnes dans le fichier targets.txt. Ici, on précise que le séparateur est la tabulation (\t).
Charger le fichier targets.
Remarque : pour obtenir de l'aide sur une commande ou fonction :
R> help(readTargets) # affiche l'aide de la fonction readTargetsLors de la lecture des fichiers de sortie de logiciels d'analyse d'image, certaines fonctionnalités de limma permettent de contrôler la qualité des données d'expressions en chaque spot. En effet, ces fichiers contiennent eux-mêmes différentes indications sur la qualité des spots. Par exemple, les fichiers Genepix sont structurés comme des tableaux de données dont les lignes correspondent aux spots et les colonnes à des caractéristiques des données en chacun de ces spots. On trouve parmi ces caractéristiques :
Il est donc possible de créer sa propre fonction de filtrage de spots : cette fonction, dont l'argument principal est le tableau spots x caractéristiques d'un fichier Genepix, prendra la valeur 1 si un spot est valide et 0 sinon. Par exemple, si on souhaite supprimer les spots dont les flags de Genepix sont inférieurs ou égaux à -49, il suffit de créer la fonction de filltrage suivante :
myFilter = function(X) { # X désigne le tableau Genepix okFLAG = X$Flags > -49; # okFLAG=TRUE si Flags>-49, FALSE sinon return (as.numeric(okFLAG)); # la fonction as.numeric convertit une valeur logique } # en valeur numérique selon la règle suivante : TRUE vaut 1, FALSE vaut 0La fonction read.maimages permet de lire les fichiers Genepix :
R> RG=read.maimages(files=targets$FileName,source="genepix",wt.fun=myFilter) # wt.fun [pour weight function] sert a spécifier un mode de filtrageChargez données des puces dans la variable RG.
Tel que défini ci-dessus, l'objet RG peut-être vu comme une armoire contenant toute l'information contenue dans les fichiers Genepix. Les noms des différents tiroirs de cette armoire sont les suivants :
R> names(RG) [1] "R" "G" "Rb" "Gb" "weights" "targets" "genes" "source" "printer"Plus précisément,
On peut utiliser RG$weights pour mesurer l'impact de la procédure de filtrage. Par exemple, la commande suivante calcule, pour chaque puce, le pourcentage de gènes non-filtrés :
R> unFiltered = apply(RG$weights,MARGIN=2,FUN=mean) # unFiltered reçoit les taux de spots non-filtrés pour chaque puce R> round(unFiltered,2) # Affichage des éléments de unFiltered arrondis à 2 décimalesLes noms des différents tiroirs de l'armoire nommée RG$genes sont accessibles par la commande suivante :
R> names(RG$genes) [1] "Block" "Row" "Column" "ID" "Name"Pour obtenir les identifiants des 10 premiers gènes, on peut utiliser le tiroir nommé RG$genes$ID :
R> RG$genes$ID[1:10]Dans les fichiers .gpr que vous avez téléchargés, les noms des gènes ne sont pas renseignés, ils sont en fait déjà indiqués au niveau des identifiants mais pas au niveau des noms de gène. Assurez-vous en et affectez les identifiants aux noms de gènes avec les commandes suivantes :
R> RG$genes$ID[1:10] # Affiche les identifiants des 10 premiers gènes (de 1 à 10) R> RG$genes$Name[1:10] # Affiche les noms des 10 premiers gènes R> RG$genes$Name=RG$genes$ID # Affecte/recopie les identifiants des gènes dans les noms des gènesDans l'optique d'une meilleure lisibilité des résultats de l'analyse, il peut être intéressant de tenir compte d'une typologie des spots à partir de leur statut. : gène, blanc, contrôle, ... On utilise pour cela un fichier texte (les séparateurs de colonnes sont des tabulations) appelé spotTypes.txt. Pour les données dont on dispose, ce fichier prend la forme du tableau suivant.
| SpotType | ID | Name | Color |
|---|---|---|---|
| gene | * | * | black |
| red_dye | Cy5 | * | red |
| green_dye | Cy3 | * | green |
| blank | - | * | yellow |
| ngt | NGT | * | blue |
La dernière colonne du fichier, nommée Color, permet d'associer une couleur à chaque type de spot, afin de mieux les distinguer dans les analyses à venir. Les commandes suivantes lisent le fichier spotTypes.txt et associent à chaque spot son statut dans RG$genes :
R> spottypes = readSpotTypes("spotTypes.txt") # spottypes reçoit le contenu du fichier spotTypes.txt R> RG$genes$Status = controlStatus(spottypes,RG) # RG$genes$Status est un nouveau tiroir de RG$genes contenant # les statuts de chaque spot Matching patterns for: ID Name Found 6272 gene Found 16 red_dye Found 16 green_dye Found 198 blank Found 151 ngt Setting attributes: values ColorSauvegarder le fichier spotTypes.txt dans votre répertoire de travail puis effectuer les commandes précédentes.
La normalisation des données est une étape importante dans la création du tableau de données sur lequel reposera l'analyse statistique. L'objectif de cette normalisation est d'identifier des sources parasites de variations des expressions. Parmi ces sources, on cible plus précisément :
De manière générale, il est important de s'assurer que les variations du bruit de fond des images analysées par Genepix ne sont pas structurées :
Que pensez-vous du bruit de fond des puces ?
La correction du signal par soustraction du bruit de fond peut être obtenue par la fonction backgroundCorrect.
L'hétérogénéité des signaux mesurés sur une puce est souvent attribuée à des variations uniquement liées aux fluorochromes. Cette différence entre les signaux est décelable sur le nuage dont les points, représentant des spots, ont pour abscisses les moyennes entre les logarithmes des signaux verts et rouges (dites valeurs A) et pour ordonnées les différences entre logarithmes des signaux verts et rouges (dites valeurs M). Ce graphique, dit graphe MA, est obtenu par la fonction plotMA :
R> plotMA(RG,status=RG[,1]$genes$Status) # L'argument status permet d'identifier les différents types de gènes. R> abline(0,0,col="blue") # ajoute la droite y = 0Afficher les graphes MA pour chacune des puces. Observez-vous des différences ?
Il est possible de visualiser les distributions des intensités des puces avec la commande suivante :
R> plotDensities(RG)
ou bien chaque puce individuellement :
R> plotDensities(RG[,1]) # pour la 1ère puce
Comment interprétez-vous ces graphiques ?
La fonction normalizeWithinArrays corrige cette hétérogénéité en ajustant un modèle des variations des valeurs M en fonction des valeurs A par une méthode de régression locale (dite méthode loess) :
R> MA = normalizeWithinArrays(RG,method="loess",bc.method="none") # MA contient désormais des données corrigées.Le résultat de la commande précédente, l'objet MA, ressemble en tous points à l'objet RG : seuls les tiroirs nommés R et G et contenant les signaux rouges et verts sont remplacés par d'autres, nommées M et A et contenant les valeurs du même nom. Le paramètre bc.method indique la méthode utilisée pour corriger le bruit de fond. La méthode par défaut est subtract qui soustrait les valeurs de bruit de fonds aux intensités. Il est parfois préférable, notamment pour l'analyse différentielle, d'utiliser la méthode normexp qui ajuste les intensités de manière adaptive au bruit de fond et mène à des valeurs d'intensité strictement positives. On peut en plus passer le paramètre offset qui sera ajouté à chaque intensité avant la transformation et qui permet de tirer vers 0 les log-ratios pour les faibles intensités. La fonction RG.MA permet d'obtenir les signaux verts et rouges normalisés. L'effet de la normalisation intra-puces est affiché par la commande suivante :
R> plotMA(MA,status=RG$genes$Status) R> abline(0,0,col="blue")Le graphique obtenu montre bien un autre effet de la normalisation décrite ci-dessus : le centrage par puces des valeurs M. La commande suivante permet de visualiser la dispersion des valeurs de M dans les différentes puces.
R> boxplot(data.frame(MA$M),main="log-ratios")La fonction plotDensities affiche les densités empiriques lissées pour les intensités rouges et vertes des puces. Sans normalisation, on observe une grande variabilité entre les canaux et entre les puces.
R> plotDensitities(RG) # affiche les fonctions de densité empiriques lissées des intensités rouges et vertes R> plotDensitities(MA) # après normalisationLa normalisation intra-puce vous semble-t-elle satisfaisante ? Essayer d'autres méthodes de normalisation (taper ?normalizeWithinArrays pour avoir l'aide et la liste des méthodes disponibles).
Après normalisation des valeurs M pour chaque puce, les distributions des intensités rouges et vertes devraient être essentiellement les mêmes pour chaque puce. Cependant, il peut être observé une grande variabilité entre les puces. La normalisaiton loess n'affecte pas les valeurs A. La fonction normalizeBetweenArrays permet de normaliser les distributions d'intensité entre les puces :
R> MA2=normalizeBetweenArrays(MA, method="Aquantile") R> plotDensities(MA2)Essayez différentes méthodes de correction du bruit et de fond, et de normalisation (intra-puce), jusqu'à obtenir des distributions d'intensité satisfaisantes.
limma permet d'identifier les gènes différentiellement exprimés au moyen de régression linéaire et de modèle Bayésiens (entre autres). Pour analyser vos puces avec limma, utilisez les commandes suivantes (correspondant au fichier targets_15min.txt du début du TD, dans lequel le contrôle (ref) correspond au cy3 cf. GSM239212 pour la 1ère puce) :
R> targets = readTargets("targets_15min.txt") R> design = modelMatrix(targets,ref="ref") R> RG = read.maimages(files=targets$FileName,source="genepix",wt.fun=myFilter) R> spottypes = readSpotTypes("spotTypes.txt") R> RG$genes$Name=RG$genes$ID R> RG$genes$Status = controlStatus(spottypes,RG) R> MA=normalizeWithinArrays(RG,method="robustspline",bc.method="none") R> MA2=normalizeBetweenArrays(MA,method="Aquantile") R> fit = lmFit(MA2, design) R> efit = eBayes(fit) R> results=topTable(efit,number=6000,adjust.method="BH",p.value=0.05,sort.by="logFC")Vous pouvez afficher les résultats de manière graphique avec un graphe volcano :
R> volcanoplot(efit)
Afficher les résultats contenus dans results. Vous remarquerez la présence d'identifiants tels que Cy5 ou NGT. Pour s'en débarrasser effectuer les commandes suivantes et enregistrer les résultats dans un fichier texte tabulé :
R> genes=results$Status=="gene" R> write.table(results[genes,],"deg_15min.txt",quote=FALSE,sep="\t",row.names=TRUE, col.names=TRUE)
Nous allons à présent analyser les gènes précédemment identifiés comme différentiellement exprimés. Pour cela, nous allons rechercher les annotations surreprésentées au sein de cet ensemble de gènes. Le principe est de tester si la fréquence d'une annotation est largement supérieure à la fréquence attendue (en générale la fréquence observée pour le génome) et d'effectuer ce test pour chaque annotation disponible.
Analysez les ensembles de gènes différentiellement exprimés à 15 minutes identifiés précédemment à l'aide du site YeastSpec/funSpec dédié aux annotations sur la levure.
Dans cette partie, vous allez utiliser des données de transcriptome de levure obtenues pour environ 300 mutants (approche compendium). L'idée est d'effectuer un clustering hiérarchique des profils d'expression des gènes dans toutes ces conditions expérimentales afin d'identifier des ensembles de gènes co-exprimés et donc potentiellement co-régulés et qui pourraient donc participer à un même processus biologique. Une fois les gènes regroupés en cluster, le but est de caractériser ces ensembles de gènes par des méthodes de surreprésentation statistique pour identifier dans quel(s) processus biologique(s) ils sont impliqués. L'analyse plus fine d'un cluster de gènes pourra permettre par la suite d'attribuer un processus biologique potentiel pour un gène de fonction inconnue dont le profil d'expression est fortement corrélé avec ceux d'un groupe de gènes bien connus.
Les données que vous allez utiliser ont été rassemblées et analysées et ont fait l'objet d'une publication dont l'identifiant PubMed es t 10929718. Retrouver et lire l'abstract sur PubMed. Avant de réellement analyser ces données, il est nécessaire d'appliquer un certain filtrage. Pour simplifier pour ce TD, on n'a retenu seulement les gènes au moins 2x plus ou 2x moins exprimés dans au moins 5 conditions ainsi que les conditions pour lesquelles on a au moins 5 gènes 2x plus ou 2x moins exprimés. Téléchargez les données sélectionnées et importez-les ensuite dans MeV (gratuit, normalement déjà installé sur votre machine) et effectuer un clustering hiérarchique en utilisant le coefficient de corrélation de Pearson. Remarque : vous pouvez ajuster la taille d'affichage des ratios dans le menu Display -> Set element size.
A présent, sélectionnez, enregistrez et analysez les clusters qui vous semblent pertinents (au moins 5). Pour l'analyse, utilisez l'outil YeastSpec/funSpec en utilisant les sources de données suivantes : MIPS Functional Classification, MIPS Phenotypes, MIPS Subcellular Localization, MIPS Protein Complexes, et GO Biological Process. Interpréter les résultats obtenus.