Analyse de transcriptome

Partie 1. Premiers pas avec R

Introduction

R est un logiciel libre de calcul statistique basé sur un environnement orienté objet et sur le langage de programmation S. Il est disponible sous Microsoft Windows ou sur stations Unix ou Linux. R est utilisé en ligne de commandes. Cette approche permet de décomposer les différentes actions et ainsi de clarifier les étapes d'analyse du jeu de données d'intérêt. L'approche logiciel libre a permis à la communauté d'utilisateurs de proposer un large éventail de paquetages (packages), eux aussi suivant le principe du logiciel libre, faisant de R la plateforme de choix pour de nombreuses analyses telles que l'analyse de données de transcriptome (packages bioconductor et limma vus dans la deuxième partie), traitement de graphes, etc. Le site CRAN (Comprehensive R Archive Network) recense et héberge l'essentiel de ces packages.

Objectifs de cette partie

L'environnement R

Tout ce qu'on trouve sur R est un objet : les données importées, les graphiques, les fonctions, les résultats obtenus. Les objets ont pour avantage de pouvoir être manipulés de diverses manières : afficher (pour voir leur contenu), copier, coller, transformer. Il existe plusieurs types d'objets :

Les opérateurs

Tout d'abord, la fenêtre de commandes permet d'effectuer des calculs, avec les opérateurs mathématiques courants :

Notation :

Exemples :

R> 5*3 [1] 15 # le [1] correspond à l'élément 1 de la liste de résultats (un seul ici) R> 5*3+2 [1] 17

R vous donne immédiatement les résultats. Mais pour les conserver, vous devez leur donner un nom, sinon ils ne seront pas enregistrés. Pour cela, vous devez créer un objet.

Essayer les commandes précédentes, et en ajouter de nouvelles en utilisant les opérateurs vus précédemment.

Les variables

Sous R, il est possible de mémoriser les résultats d'un calcul ou d'un traitement (comme un graphique), une fonction, etc. Il s'agit en fait d'affecter un objet (résultat) à une variable. Différents opérateurs d'affectation sont disponible : =, <-, ou _. Exemple :

R> taille=1.78 R> poids=66 R> imc=poids/(taille*taille) R> imc [1] 20.83070

Remarque : l'indice de masse corporelle est une grandeur qui permet d'estimer la corpulence d'une personne. Selon la classification de l'OMS, une corpulence normale correspond à un IMC compris entre 18,5 et 25.

Pour pouvoir réutiliser le calcul à loisir, il peut être intéressant d'en faire une fonction :

R> imc = function(poids,taille) return (poids/(taille*taille)) # on définit la function imc # qui prend en paramètre poids et taille, calcule poids/(taille*taille) et retourne le résultat R> imc # on peut demander d'afficher le contenu d'une fonction function(poids,taille) return (poids/(taille*taille)) # on l'utilise ensuite comme suit : R> imc(66,1.78) [1] 20.83070

Définissez la fonction imc (copier/coller) et cherchez le poids que vous devriez faire pour votre taille pour atteindre un imc de 22.

Comme indiqué sur la page Wikipedia, il est important de garder à l'esprit que l'IMC n'est qu'un indicateur, non pas une donnée absolue. Et tout jugement doit également prendre en compte son indice de masse grasse (IMG) et la consultation d'un médecin nutritionniste ou d'un diététicien diplômé est recommandée.

Une des formules permettant de calculer l'IMG est : IMG = (1,2 x IMC) + (0,23 x âge) - (10,8 x S) - 5,4 avec (S = 0 pour la femme et S = 1 pour l'homme). On peut donc définir la fonction R :

img = function(poids, taille, sexe) { if (sexe=='M') s=1 else s=0 return ((1.2 * imc(poids, taille) ) + (0.23 * 33) - (10.8 * s) - 5.4) } R> img(66,1.78,'M') [1] 16.38684

Nous voyons au passage qu'il est possible d'effectuer certaines opérations (s=1) si une certaine condition est satisfaite (if), et sinon (else) une autre opération (s=0).

Nous remarquons également que nous avons réutilisé la fonction imc définie précédemment.

Définissez la fonction img et utilisez la pour calculer le votre ou celui de qui vous voulez.

Help

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

Avec les fonctions d'aide, trouver la description de la fonction qui permet d'effectuer un test de Student.

Remarque : Raccourcis pour la commande help('nom_de_commande') : R> ?nom_de_commande, et pour la commande help.search('mot') : R> ??mot.

Chargement d'un fichier

Nous allons nous appuyer sur un jeu de données associées à la publication [Ivshina et al.(2006) Cancer Research] consistant en des données cliniques sur des patientes atteintes d'un cancer du sein ainsi que de données de transcriptome.

Pour le choix du traitement d'un cancer du sein, les décisions sont guidées par la détermination du potentiel métastatique des tumeurs. Des mesures cliniques sont utilisées afin d'évaluer ce potentiel (état du ganglion lymphatique, taille de la tumeur) et prédire la réactivité endocrine (récepteurs à l'oestrogène et à la progestérone) résultant en une classification des tumeurs en sous-types associés à un pronostic. Dans cette étude, les auteurs soulignent les défauts (manque de précision, subjectivité des évaluations) de tels indicateurs et proposent d'identifier et d'utiliser une signature génétique, c'est-à-dire les profils d'expression de certains gènes marqueurs, pour améliorer la précision de la classification des sous-types de tumeurs et les pronostics associés.

Le système d'évaluation de Nottingham intègre des mesures de différentiation cellulaire et de potentiel réplicatif, et fournit un score qui quantifie l'agressivité d'une tumeur. Les tumeurs sont ainsi classées en type G1 (cellules bien différenciées à croissance lente), G2 (différenciation modérée) et G3 (cellules faiblement différenciées et très prolifératives). Des études ont montré que ce classement permet de bien prédire les risques de récidives et de décès indépendamment de l'état des ganglions lymphatiques et de la taille des tumeurs. En revanche pour la planification des traitement thérapeutique, cette classification a fait l'objet de controverse parmi les oncologues. De plus, les types G1 et G3 ont une pertinence clinique plus évidente que le type G2, plus hétérogène, et qui représente la moitié des cas de cancer du sein.

Nous allons nous intéresser aux données cliniques dont voici un extrait :

array_idgradesurvival_timealive_at_endpointestrogen_receptor_statuspopulationlymphnode_statusage_at_diagnosisp53_mutanttumor_size
NUH.BT114G2Singapore
NUH.BT119G2Singapore
X100B08G111.831Uppsala16819
X101B88G311.830Uppsala040012
X102B06G311.8301Uppsala151026
Table 1 : Extrait du fichier clinical_info.xls.

Description des colonnes :

Ouvrir le fichier clinical_info.xls, soit avec Ms Excel, soit avec OpenOffice, et enregistrer le fichier au format texte avec séparateur (vous pouvez utiliser des tabulations comme séparateur par exemple). Attention si vous utilisez un Excel français, il faudra au préalable remplacer les virgules (,) par des points (.) pour les colonnes qui contiennent des valeurs décimales.

Après conversion, vous pourrez charger le fichier dans R avec la commande suivante (le fichier converti s'appelle clinical_info.txt). Attention, avant de lancer la commande, il faut préciser à R où (dans quel dossier) trouver le fichier : pour cela, utiliser la fonction Changer de répertoire courant... depuis le menu Fichier de R, et sélectionner le répertoire dans lequel se trouve le fichier clinical_info.txt.

R> read.table("clinical_info.txt", sep="\t", header=TRUE) # la fonction prend en paramètres (entre autres) : # le nom du fichier # le type de séparateur, ici \t correspond à une tabulation # header=TRUE permet de préciser que la première ligne du fichier contient les noms des colonnes. # faire ?read.table pour avoir l'aide complète

Remarque : Il est nécessaire de stocker le résultat dans une variable, sinon le contenu du fichier est affiché puis "oublié". Charger le contenu du fichier dans la variable data comme suit :

R> data = read.table("clinical_info.txt", sep="\t", header=TRUE)

Faire afficher son contenu.

Statistiques descriptives

R étant un environnement statistique, il propose une large gamme de fonctions statistiques et graphiques. Nous allons utiliser quelques une d'entre elles pour appréhender notre jeu de données : summary (fonction générique produisant le résumé du contenu d'une variable), pie (camembert), hist (histogramme), boxplot (boite à moustaches), density (estimation d'une distribution), plot (fonction générique pour tracer un graphique), qplot (graphe des quantiles) et qqplot (graphe des quantiles-quantiles).

Pour la suite, nous allons avoir besoin d'accéder à certaines colonnes seulement, ou bien à certaines lignes de nos données. Ceci peut se faire de différentes manières, par exemple pour la colonne correspondant au type de tumeur :

Remarque : On peut obtenir le nom des colonnes avec colnames(variable), et les dimensions de la variables dim(variable).

Donner pour la colonne age_at_diagnosis : la valeur minimale (fonction min), maximale (max), la moyenne (mean), la médiane (median), les 4-quantiles ou quartiles (quantile), l'écart-type (fonction sd pour standard deviation en anglais). Cette colonne contient des valeurs manquantes (NA), il vous faudra ajouter le paramètre na.rm=TRUE lors de l'appel de ces fonctions. Essayer également la fonction summary sur cette colonne seulement, puis sur l'ensemble des données et commenter la sortie.

Camembert. Essayer maintenant la commande summary sur la colonne grade. Afficher le camembert en faisant : pie(summary(data$grade))

Distribution de l'âge des patientes lors du diagnostic : histogrammes, distributions et boites à moustaches. Utiliser la commande hist pour tracer un histogramme. Faire varier le nombre d'intervalles avec le paramètre breaks.

Utiliser maintenant la commande density sur age_at_diagnosis et afficher le résultat graphiquement. Il vous faudra appeler la fonction plot sur le résultat de density (avec le paramètre na.rm=TRUE comme précédemment).

Utiliser maintenant la fonction boxplot pour réaliser le tracer des boites à moustaches.

De la même manière, regarder les distributions pour les colonnes survival_time et tumor_size.

Nous allons maintenant sélectionner seulement ces colonnes (qui contiennent des valeurs numériques) et les sauvegarder dans une autre data frame afin de rechercher visuellement des corrélations :

R> vals = data.frame(data$survival_time, data$age_at_diagnosis, data$tumor_size) # data.frame() créé une frame avec les colonnes (vecteurs) passées en paramètres.

Interpréter le résultat de la commande pairs(vals,pch=19)

Il est possible de colorer les points en fonction d'un facteur. Nous choisirons ici le type de tumeur. Tout d'abord, il s'agit de replacer le type de tumeur par la couleur. Ceci s'efffectue de la manière suivante : c("green","yellow","red")[data$grade]. La fonction c() créé un vecteur. Ici, il contient 3 éléments. Ces éléments sont ensuite utilisés pour remplacer les valeurs de data$grade (s'il y a plus de 3 valeurs différentes, les couleurs seront réutilisées). On affiche donc le plot de la manière suivante : pairs(vals, col=c("green","yellow","red")[data$grade],pch=19). Effectuer cette commande. Observez-vous une corrélation entre certaines variables ?

Distributions et simulations

Nous allons maintenant nous intéresser à une distribution normale. Pour cela, R propose des fonctions pour la loi de probabilité (dnorm), la p-valeur (pnorm), ainsi qu'une fonction génératrice d'échantillons aléatoires suivant une loi normale (rnorm). Consulter l'aide de ces fonctions. Il est par exemple possible de tracer une loi normale de moyenne -0.5 et d'écart-type 1 de la manière suivante. Tout d'abord, nous allons définir les valeurs de x qui nous intéressent :

R> x=seq(-4,4,length.out=100)

Afficher le contenu de x.

Ensuite, nous calculons pour chaque valeur de x, sa probabilité :

R> t1=dnorm(x,-0.5)

Afficher le contenu de t1.

Tracer maintenant le graphique avec la commande suivante :

R> plot(x,t1,type="l",col=2)

Nous allons maintenant tirer un échantillon aléatoire issu de cette distribution avec la fonction rnorm. Effectuer la commande suivante :

R> d1=rnorm(100,-0.5)

Afficher le contenu de d1. Ensuite, ajouter la distribution estimée au graphique avec la commande suivante :

R> lines(density(d1),col=5)

Que constatez-vous ?

Nous allons maintenant générer un deuxième échantillon tiré d'une autre distribution (moyenne à 0.5 au lieu de -0.5) et ajouter cela au graphique précédent :

R> d2=rnorm(100,0.5) R> t2=dnorm(x,0.5) R> lines(density(d2),col=5) R> lines(x,t2,type="l",col=4)

Constater visuellement la différence entre les échantillons et les distributions théoriques.

A présent, nous allons utiliser le test de Student pour déterminer si les valeurs des deux échantillons proviennent de distributions ayant la même moyenne. Effectuer le test de Student (test t) de la manière suivante (le paramètre equal.var=TRUE permet d'indiquer que les deux échantillons ont même variance) :

R> t.test(d1,d2,var.equal=TRUE)

Que pouvez-vous conclure ?

Régression linéaire

Nous allons maintenant effectuer une régression linéaire à partir de 6 mesures de concentration prises par intervalle de 5 minutes.

Effectuer les commandes ci-dessous pour tracer les points.

R> x=1:6 R> y=c(1.2,1.8,2.9,4.1,5,5.1) R> plot(x,y)

La fonction lm permet d'inférer des modèles linéaires des données, et donc notamment, de faire une régression linéaire pour obtenir la droite qui va passer le plus près de chacun des points précédemment tracés. Effectuer les commandes suivantes pour faire la régression et ajouter la droite au graphique précédent.

R> m = lm(y~X) R> lines(x,m$fitted.values)

Pour une prise en main de R plus complète, vous trouverez un tutorial d'introduction (couvrant la plupart des points vus jusqu'ici) sur le site du CRAN ainsi que d'autres sources de documentation.

Partie 2. Analyse de transcriptome avec Bioconductor

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

Installation de bioconductor

Remarque: Bioconductor est normalement déjà installé sur l'ordinateur. Pour vérifier, dans R, taper la commande library(limma) si aucun message d'erreur n'apparaît, c'est que bioconductor a été chargé.

Pour la réalisation de ces analyses, nous utiliserons la suite bioconductor qui est un ensemble de packages pour le logiciel de traitements statistiques R. L'installation de bioconductor se fait très simplement aux moyens des deux commande s suivantes dans une session R :

> source("http://www.bioconductor.org/biocLite.R") # Chargement du fichier de commandes biocLite.R # disponible sur la page web http://www.bioconductor.org. > biocLite() # La fonction biocLite est définie dans biocLite.R. Elle commande le # téléchargement automatique des packages (version allégée de Bioconductor).

Récupération des données à partir d'un entrepôt de données de transcriptome

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.

Importation des données avec limma

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 limma

Remarque : 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 :

# sous Windows, remplacer "C:\My Documents\Transcriptome\" par l'endroit où vous avez placé vos fichiers R> setwd("C:\My Documents\Transcriptome\") # sous Windows définit le répertoire C:\My Documents\Transcriptome comme répertoire de travail # de même sous linux R> setwd("/home/myLogin/Transcriptome/") # sous linux

Remarque: Sous Windows, vous pouvez changer le répertoire de travail à partir du menu Fichier puis changer le répertoire courant...

Le 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

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.

FileNameSlideNumberCy3Cy5
GSM239212.gpr1ref15min
GSM239213.gpr2ref15min
GSM239214.gpr3ref15min
Table 1 : Contenu du fichier targets_15min.txt.

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_15min.txt",sep="\t") # Lecture de targets.txt

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

Chargez les trois fichiers targets dans trois variables (targets15, targets30 et targets60).

Remarque : pour obtenir de l'aide sur une commande ou fonction :

R> help(readTargets) # affiche l'aide de la fonction readTargets

Filtrage des spots

Lors 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 :

# X désigne une ligne du tableau du fichier Genepix (un spot) # okFLAG=TRUE si Flags>-49, FALSE sinon # la fonction as.numeric convertit une valeur logique # en valeur numérique selon la règle suivante : TRUE vaut 1, FALSE vaut 0 myFilter = function(X) { okFLAG = X$Flags > -49; as.numeric(okFLAG) }

Définissez cette fonction dans votre session R.

Lecture des fichiers Genepix

La 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 filtrage

Chargez l'ensemble des données correspondant aux trois réplicats à 15 minutes 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,

Taux de filtrage

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écimales

Nom et identifiant d'un gène

Les 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ènes

Statut d'un spot

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

SpotTypeIDNameColor
gene**black
red_dyeCy5*red
green_dyeCy3*green
blank-*yellow
ngtNGT*blue
Table 2 : Contenu du fichier spotTypes.txt.

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 Color

Sauvegarder le fichier spotTypes.txt dans votre répertoire de travail puis effectuer les commandes précédentes.

Normalisation des données

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 :

Hétérogénéité du bruit de fond

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 :

R> boxplot(data.frame(log2(RG$Rb[,1:3])),main="Red background") # Boîtes de dispersion des colonnes de log2(RG$Rb). # L'argument main ajoute un titre au graphique. R> imageplot(log2(RG$Rb[,1]),layout=RG$printer) # Image de log2(RG$Rb) # pour la 1ère puce. L'argument layout permet de tenir compte # des blocs d'aiguilles dans la puce. R> imageplot(log2(RG$Rb[,2]),layout=RG$printer) # Image de log2(RG$Rb) # pour la 2ème puce R> imageplot(log2(RG$Rb[,3]),layout=RG$printer) # Image de log2(RG$Rb) # pour la 3ème puce

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.

Hétérogénéité du signal

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[,1],status=RG$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 = 0

Afficher 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> plotDensities(RG) # affiche les fonctions de densité empiriques lissées des intensités rouges et vertes R> plotDensities(MA) # après normalisation

La 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> MA=normalizeBetweenArrays(MA, method="Aquantile") R> plotDensities(MA)

Essayez différentes méthodes de correction du bruit et de fond, et de normalisation (intra-puce), jusqu'à obtenir des distributions d'intensité satisfaisantes.

Gènes différentiellement exprimés

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

Afficher les résultats. 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=FALSE, col.names=TRUE)

Vous pouvez également afficher ces résultats de manière graphique avec un graphe volcano :

R> volcanoplot(efit)

Caractérisation d'un ensemble de gènes

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.

Partie 3 : Analyse d'un compendium de données de transcriptome

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. Télécharger les données sur les ratios d'expression pour les 300 conditions et les décompresser dans votre répertoire de travail. 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.

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écharger les données sélectionnées et refaites le clustering comme précédemment.

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.


Université Paul Sabatier © Roland Barriot 2009