From silico.biotoul.fr

Détermination des gènes différentiellement exprimés

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.

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

# Chargement du package limma
library(limma)

Remarque : limma dispose d'un guide utilisateur accessible et pédagogique, incluant notamment de nombreux exemples d'utilisation. Pour y accéder :

limmaUsersGuide()

Afin d'éviter d'avoir à indiquer de manière répétée le répertoire de travail contenant les fichiers de données 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 : changez 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 :

• FileName donnant le nom des fichiers d'images,
• Cy3 donnant la modalité du facteur marquée en vert,
• Cy5 donnant la modalité du facteur marquée en rouge.

Important : le marquage (Cy5 ou Cy3) pour chaque hybridation peut être différent. Il convient de vérifier dans la description des données quel canal correspond à quelle condition. Par exemple, pour l'hybridation GSM239312, dans la section Channel 1, on peut lire à la ligne Source Name : Recovery from freezing (15 min), Cy5 ; et dans la section Channel 2, Source name : control (before freezing treatment), Cy3 ; d'où le ref pour Cy3 et le 15min pour Cy5 pour cette puce.

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.

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 :

# Lecture de targets.txt
targets=readTargets("Targets15min.txt",sep="\t")

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

Créez ou récupérez le fichier targets puis chargez le dans un object R targets.

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

# affiche l'aide de la fonction readTargets
help(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 :

• Flags : type de Flag
• F635 Median : intensité médiane pour le signal rouge (F pour foreground, 635 pour la longueur d'onde 635nm)
• F532 Median : intensité médiane pour le signal vert
• F635 Mean : intensité moyenne pour le signal rouge
• F532 Mean : intensité moyenne pour le signal vert
• B635 Median : intensité médiane du bruit de fond pour le rouge (B pour background)
• ...

Il est 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 (correspond en fait au poids affecté à chacun des spots). 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 filtrage suivante :

# X désigne une ligne du tableau du fichier Genepix (un spot) 
# okFLAG vaut 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 :

RG = read.maimages(files = targets$FileName, source = "genepix", wt.fun=myFilter) # wt.fun pour ''weight function'' qui 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 :

names(RG) 
[1] "R" "G" "Rb" "Gb" "weights" "targets" "genes" "source" "printer"

Plus précisément,

• RG$R est un tableau spots x puces contenant les valeurs du signal rouge,
• RG$G est un tableau spots x puces contenant les valeurs du signal vert,
• RG$Rb est un tableau spots x puces contenant les valeurs du bruit de fond rouge,
• RG$Gb est un tableau spots x puces contenant les valeurs du bruit de fond vert,
• RG$weights est un tableau spots x puces contenant les valeurs 1 pour les spots conservés lors du filtrage et 0 sinon,
• RG$targets contient les information du fichier targets.txt,
• RG$genes est lui-même une armoire dont les différents tiroirs contiennent des informations sur chaque gène,
• RG$source rappelle le logiciel d'analyse d'images d'où viennent les données,
• RG$printer décrit la structure d'une puce (taille d'un bloc d'aiguilles, ...).

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 :

# unFiltered reçoit les taux de spots non-filtrés pour chaque puce
unFiltered = apply(RG$weights,MARGIN=2,FUN=mean) 
# Affichage des éléments de unFiltered arrondis à 2 décimales
round(unFiltered,2)

Nom et identifiant d'un spot/gène

Les noms des différents tiroirs de l'armoire nommée RG$genes sont accessibles par la commande suivante :

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 :

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 :

# Affiche les identifiants des 10 premiers spots (de 1 à 10) 
RG$genes$ID[1:10]
# Affiche les noms des 10 premiers spots
RG$genes$Name[1:10]
# Affecte/recopie les identifiants des spots dans les noms des spots 
RG$genes$Name=RG$genes$ID

Remarque : La manipulation précédente est spécifique à notre jeu de données et doit être adaptée à chaque analyse.

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) généralement appelé SpotTypes.txt. Pour les données dont on dispose, ce fichier prend la forme du tableau suivant.

La dernière colonne du fichier, nommée Color, permet d'associer une couleur à chaque type de spot (colonne SpotType), afin de mieux les distinguer dans les analyses à venir. Les colonnes ID et Name indiquent que l'attribution d'un type de spot se fera sur ces caractéristiques (contenues dans l'objet RG : RG$genes$ID et RG$genes$Name que nous avons manipulées dans la section précédente). Les * correspondent à un joker et indiquent que par défaut le type de spot sera gene. Ensuite, si la valeur de ID d'un spot correspond à Cy5, le type red_dye sera attribué au spot.

Important : chaque plateforme (version d'une puce pour un organisme) est spécifique et les valeurs dans les colonnes ID et Name du fichier Genepix sont propres à la plateforme. Il faut donc construire soit-même le fichier SpotTypes.txt afin de sélectionner les spots qui nous intéressent pour une analyse donnée.

Les commandes suivantes lisent le fichier SpotTypes.txt et associent à chaque spot son statut dans RG$genes :

# spottypes reçoit le contenu du fichier SpotTypes.txt
spottypes = readSpotTypes("SpotTypes.txt")
# RG$genes$Status est un nouveau tiroir de RG$genes contenant
# les statuts de chaque spot
RG$genes$Status = controlStatus(spottypes,RG) 
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

Créer ou 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 :

• l'hétérogénéité du bruit de fond.
  

Si le bruit de fond présente des variations très différentes d'une puce à une autre, ou très structurées spatialement sur une puce, alors on peut être amené à corriger le signal par soustraction du bruit de fond.

• l'hétérogénéité du signal.
  

De la même manière, le principe d'invariance d'une très grande majorité des expressions géniques d'une puce à une autre doit se traduire par une répartition comparable des valeurs des signaux entre les différentes puces. Si des différences marquées existent, il est judicieux de ramener les signaux moyens de chaque puce à la même valeur. D'autre part, il arrive que le fluorochrome soit lui-même source de variations du signal. Cet effet du fluorochrome dépend même parfois de la valeur d'expression et de la structure de la puce (partitionnement en blocs d'aiguilles).

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 :

• d'une part, par une comparaison des variations du bruit entre les différentes puces.

La commande suivante construit une boîte de dispersion du logarithme du bruit de fond pour les puces :

# Boîtes de dispersion des colonnes de log2(RG$Rb).
# L'argument main ajoute un titre au graphique.
boxplot(data.frame(log2(RG$Rb[,1:3])),main="Red background")

• d'autre part, par un examen de la répartition spatiale du bruit de fond de chaque puce. La commande suivante crée une image associant à chaque spot de la puce un pixel dont le niveau de couleur est proportionnel à la valeur du logarithme du bruit de fond :

# on va dessiner les 3 graphiques côte à côte
par(mfrow=c(1,3)) 
# pour la 1ère puce. L'argument layout permet de tenir compte des blocs d'aiguilles dans la puce.
imageplot(log2(RG$Rb[,1]),layout=RG$printer) # Image de log2(RG$Rb) 
# pour la 2ème puce
imageplot(log2(RG$Rb[,2]),layout=RG$printer) # Image de log2(RG$Rb)
# pour la 3ème puce  
imageplot(log2(RG$Rb[,3]),layout=RG$printer) # Image de log2(RG$Rb)
par(mfrow=c(1,1))

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 pour average dans MA-plot et I pour intensities dans RI-plot) et pour ordonnées les différences entre logarithmes des signaux verts et rouges (dites valeurs M pour minus et R pour ratio puisque log A - log B = Log A/B). Ce graphique, dit graphe MA, est obtenu par la fonction plotMA :

plotMA(RG[,1],status=RG$genes$Status) 
# remarque : le paramètre status permet d'identifier les différents types de spots
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 :

plotDensities(RG)

ou bien chaque puce individuellement :

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

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 :

par(mfrow=c(3,2)) 
# première puce
plotMA(RG[,1],status=RG$genes$Status, main="avant normalisation")
abline(0,0,col="blue")
plotMA(MA[,1],status=RG$genes$Status, main="apres normalisation intra-puces") 
abline(0,0,col="blue")
 
# deuxième puce
plotMA(RG[,2],status=RG$genes$Status, main="avant normalisation")
abline(0,0,col="blue")
plotMA(MA[,2],status=RG$genes$Status, main="apres normalisation intra-puces") 
abline(0,0,col="blue")
 
# troisième puce
plotMA(RG[,3],status=RG$genes$Status, main="avant normalisation")
abline(0,0,col="blue")
plotMA(MA[,3],status=RG$genes$Status, main="apres normalisation intra-puces") 
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.

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.

par(mfrow=c(1,2))
# affiche les fonctions de densité empiriques lissées des intensités rouges et vertes
plotDensities(RG)
# après normalisation
plotDensities(MA)

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 normalisation loess n'affecte pas les valeurs A. La fonction normalizeBetweenArrays permet de normaliser les distributions d'intensité entre les puces :

MA2=normalizeBetweenArrays(MA, method="Aquantile")
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.

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 Targets15min.txt du début du TD, dans lequel le contrôle (ref) correspond au cy3 cf. GSM239212 pour la 1ère puce) :

# indication de l'échantillon de référence
design = modelMatrix(targets,ref="ref")
# régression
fit = lmFit(MA2, design) 
# modèle bayésien
efit = eBayes(fit) 
# récupération des résultats
results=topTable(efit,number=6000,adjust.method="BH",p.value=0.05,sort.by="logFC")
# number=6000 permet de récupérer au maximum 6000 gènes différentiellement exprimés
# avec un seuil alpha de 0.05 (paramètre p.value=0.05)
# les résultats sont triés par log fold change (sort.by="logFC")

Selon la version de R/Bioconductor limma installée, il se peut que la fonction topTable ait été modifiée :

results=topTable(efit,number=6000,adjust.method="BH",p.value=0.05,resort.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écupération des lignes correspondant à des gènes
genes=results$Status=="gene"
# écriture dans un fichier des lignes correspondant à des gènes différentiellement exprimé
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 :

volcanoplot(efit)

Si l'on récapitule, l'ensemble de l'analyse tient en quelques lignes (après avoir déterminé les bons paramètres) :

# chargement du fichier décrivant les puces (images analysées par Genepix)
targets = readTargets("Targets15min.txt")
# indication de l'échantillon de référence
design = modelMatrix(targets,ref="ref")
# chargement des données
myFilter = function(X) { 
   okFLAG = X$Flags > -49; 
   as.numeric(okFLAG) 
}
RG = read.maimages(files=targets$FileName,source="genepix",wt.fun=myFilter) 
# attribution des types de spots
spottypes = readSpotTypes("spotTypes.txt") 
RG$genes$Name=RG$genes$ID 
RG$genes$Status = controlStatus(spottypes,RG) 
# normalisation intra-puce
MA=normalizeWithinArrays(RG,method="robustspline",bc.method="none") 
# normalisation inter-puces
MA2=normalizeBetweenArrays(MA,method="Aquantile") 
# identification des gènes différentiellement exprimés
fit = lmFit(MA2, design) 
efit = eBayes(fit) 
# résultat
results=topTable(efit,number=6000,adjust.method="BH",p.value=0.05,sort.by="logFC")

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 dans le génome ou le protéome) 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 FunSpec dédié aux annotations sur la levure. Vous utiliserez pour cela les sources de données MIPS.

Données

• GSE9404_RAW.tar

Détermination des gènes co-exprimés

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 est 10929718. Retrouver et lire l'abstract sur PubMed. Avant de réellement analyser ces données, il est préférable d'appliquer un certain filtrage. Pour simplifier pour ce TD, nous avons 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 les décompresser dans votre répertoire de travail. Nous allons maintenant utiliser un logiciel un peu plus convivial que la ligne de commande R. Trouver dans le menu démarrer le logiciel qui s'appelle MeV pour MultiExperiment Viewer. Importez-les ensuite dans MeV et effectuer un clustering hiérarchique. 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. 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.