M1 MABS TDB TD Transcriptome - Analyse differentielle

From silico.biotoul.fr

Current revision as of 13:43, 14 October 2015

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.

Pour plus de facilités et éviter de laborieuses manipulations, vous trouverez l'ensemble des fichiers nécessaires au TP dans l'archive TP5_transcriptome-Analyse_differentielle.zip. Télécharger l'archive et extraire son contenu.

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

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

Le fichier targets.txt

Le fichier targets permet d'indiquer les fichiers de sortie de l'analyse d'image et les conditions expérimentales associées. En d'autres termes, cela permet de faire le lien entre les intensités de fluorescences (verte ou rouge) et les conditions (référence ou traitement).

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 :

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

Filtrage des spots et lecture des fichiers Genepix

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

Il est possible de créer sa propre fonction de filtrage de spots et d'appliquer les types de filtrage vu en cours. Aujourd'hui, nous utiliserons simplement la valeur de qualité (Flags) de chaque spot attribuée par genepix.

La fonction read.maimages permet de lire les fichiers Genepix :

RG = read.maimages(files = targets$FileName, source = "genepix", wt.fun=function(X) as.numeric(X$Flags> -49) )

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)

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

Intensités rouges pour les premiers spots des 3 hybridations :

head(RG$R)

Informations sur les premiers spots :

head(RG$genes)

Taux de filtrage

On peut utiliser RG$weights pour mesurer l'impact de la procédure de filtrage. RG$weights contient autant de colonnes que d'hybridations. Chaque ligne correspond à un spot et un 1 correspond à un spot conservé et un 0 à un spot filtré (c'est-à-dire non pris en compte pour la suite).

Par exemple, la commande suivante calcule, pour chaque puce, le pourcentage de gènes non-filtrés :

head(RG$weights)
# 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) 
# ou bien à l'aide de la fonction summary
summary(RG$weights)

Quels sont les pourcentages de spots éliminés par l'étape de filtrage ?

Différents types de spots

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. La colonne ID indique que l'attribution d'un type de spot se fera sur cette caractéristique (contenues dans l'objet RG : RG$genes$ID). 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 la colonne ID 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 = readSpotTypes("SpotTypes.txt")
RG$genes$Status = controlStatus(spottypes,RG)

Combien de spots correspondent à des gènes ?

Contrôle qualité visuel des hybridations

Les bruits de fond observés pour une même hybridation devraient être similaires. On ne devrait pas non plus observer de biais systèmatiques dans les données tels que rayures ou absence totale de bruit de fond/signal.

Pour visualiser ces différentes informations, on utilise la fonction imageplot qui permet de regarder la puce vue de dessus.

Pour le bruit de fond des 3 répétitions :

par(mfrow=c(2,3))
imageplot(RG$Rb[,1],layout=RG$printer, main='Red background rep1') 
imageplot(RG$Rb[,2],layout=RG$printer, main='Red background rep2') 
imageplot(RG$Rb[,3],layout=RG$printer, main='Red background rep3') 
imageplot(RG$Gb[,1],layout=RG$printer, main='Green background rep1') 
imageplot(RG$Gb[,2],layout=RG$printer, main='Green background rep2') 
imageplot(RG$Gb[,3],layout=RG$printer, main='Green background rep3')

Que pensez-vous du bruit de fond ?

Observe-t-on ici une corrélation spatiale des intensités du bruit de fond des spots ?

Pour les intensités :

par(mfrow=c(2,3))
imageplot(RG$R[,1],layout=RG$printer, main='Red signal rep1') 
imageplot(RG$R[,2],layout=RG$printer, main='Red signal rep2') 
imageplot(RG$R[,3],layout=RG$printer, main='Red signal rep3') 
imageplot(RG$G[,1],layout=RG$printer, main='Green signal rep1') 
imageplot(RG$G[,2],layout=RG$printer, main='Green signal rep2') 
imageplot(RG$G[,3],layout=RG$printer, main='Green signal rep3')

Observe-t-on ici une corrélation spatiale des intensités du signal des spots ?

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.

Selon les biopuces utilisées, la normalisation se fait en 1 ou 2 étapes :

• la normalisation intra-puce : pour les puces à plusieurs canaux (rouge, vert, …) a pour but de rendre comparables les intensités de fluorescence des différents canaux pour un même spot/gène
• la normalisation inter-puces : a pour but de rendre comparables les intensités issues d’hybridations différentes pour un même spot/gène

Pour s’en convaincre, on peut visualiser les distributions des bruits de fond ainsi que du signal des différentes hybridations :

par(mfrow=c(1,4))
boxplot(data.frame(log2(RG$Rb[,1:3])),main="Red background", las=3)
boxplot(data.frame(log2(RG$Gb[,1:3])),main="Green background", las=3)
boxplot(data.frame(log2(RG$R[,1:3])),main="Red", las=3)
boxplot(data.frame(log2(RG$G[,1:3])),main="Green", las=3)

Que peut-on observer ? par rapport aux bruis de fonds d'hybridations différentes ? au signal ?

La normalisation vous semble-t-elle nécessaire ?

Une autre manière de s’en convaincre est de visualiser ces distributions sous forme de courbes de densité.

Séparément :

par(mfrow=c(1,3))
plotDensities(RG[,1])
plotDensities(RG[,2])
plotDensities(RG[,3])

Ou sur le même plot :

plotDensities(RG)

Le détail de la fonction :

plotdensity = function(X) {
	plot(density( log2( X$R ) ), col='red' )
	lines(density( log2( X$G ) ), col='green' )
}

Normalisation intra-puces

La fonction normalizeWithinArrays corrige les valeurs afin de rendre comparables les différents canaux d'une même hybridation.

Différentes méthodes de normalisation sont disponibles dans limma que l'on spécifie avec le paramètre method qui peut prendre les valeurs none, median, loess, printtiploessn composite, control et robustspline.

Pour le détail des méthodes, il faut aller voir l'aide :

?normalizeWithinArrays

Lors de cette normalisation, il est possible de spécifier une méthode pour la prise en compte du bruit de fond avec le paramètre bc.method qui peut prendre les valeurs auto, none, subtract, ... Pour le détail des méthodes, il faut aller voir l'aide :

?backgroundCorrect

Il n'y a pas de méthode qui fonctionne mieux dans tous les cas, il faut donc les essayer et choisir celle qui est le mieux adaptée au jeu de données.

Ici on utilisera :

MA=normalizeWithinArrays(RG, method='robustspline', bc.method='subtract')

Pour visualiser le résultat, on visualise les distributions obtenues :

par(mfrow=c(1,3))
plotDensities(MA[,1])
plotDensities(MA[,2])
plotDensities(MA[,3])

Les distributions rouges et vertes se superposent à peu près sur les différentes hybridations.

Sur le même graphique :

plotDensities(MA)

Mais pas entre les hybridations : normlisation inter-puces nécéssaire

Normalisation inter-puces

Pour la 2ème étape, la fonction normalizeBetweenArrays permet de normaliser les distributions d'intensité entre les puces.

Là encore, plusieurs méthodes sont disponibles : none, scale, quantile, cyclicloess, Aquantile, ... Et il faut se référer à l'aide pour le détail des méthodes :

?normalizeBetweenArrays

Encore une fois, il n'y a pas une méthode qui donne toujours les meilleurs résultats et il faudra les essayer et choisir celle qui est le mieux adaptée au jeu de données traité.

Ici, on utilisera :

MA2=normalizeBetweenArrays(MA, method='Aquantile')

Visualisation des résultats

plotDensities(MA2)

A présent, les distributions se superposent : toutes les intensités se référant à un même spot sont comparables.

Gènes différentiellement exprimés

Pour chacune des répétitions, il est possible de tracer un MA-plot (ou RI-plot) :

par(mfrow=c(3,1))
plotMA(MA2[,1],status=RG$genes$Status, main="rep1") ; abline(h=0)
plotMA(MA2[,2],status=RG$genes$Status, main="rep2") ; abline(h=0)
plotMA(MA2[,3],status=RG$genes$Status, main="rep3") ; abline(h=0)

A quoi correspond la droite y=0 ?

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")
head(results)

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
deg = results[ results$Status=='gene', ]
head(deg)
# écriture dans un fichier des lignes correspondant à des gènes différentiellement exprimé
write.table(deg,"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
RG = read.maimages(files=targets$FileName,source="genepix",wt.fun=function(row) as.numeric(row$Flags>-49) ) 
# attribution des types de spots
spottypes = readSpotTypes("spotTypes.txt") 
RG$genes$Status = controlStatus(spottypes,RG) 
# normalisation intra-puce
MA=normalizeWithinArrays(RG,method="robustspline",bc.method="subtract") 
# 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.

Seconde analyse : plus de 2 conditions

Nous allons maintenant comparer le nombre de gènes différentiellement exprimés à chaque temps (15min, 30min et 60min après décongélation).

Le principe de l’analyse est le même :

• Création d’un fichier targets.txt pour l’ensemble des hybridations (de GSM239212 à GSM 239220)
• Filtrage et de normalisation (en vérifiant bien après chaque normalisation - MA et MA2 -, la qualité du résultat).
• Attribution des types de spots
• Identification des gène différentiellement exprimés

targets = readTargets("TargetsAll.txt")
design = modelMatrix(targets,ref="ref")
RG = read.maimages(files=targets$FileName,source="genepix",wt.fun=function(row) as.numeric(row$Flags>-49)) 
spottypes = readSpotTypes("SpotTypes.txt") 
RG$genes$Status = controlStatus(spottypes,RG) 
MA=normalizeWithinArrays(RG,method="robustspline",bc.method="subtract") 
MA2=normalizeBetweenArrays(MA,method="Aquantile") 
fit = lmFit(MA2, design) 
efit = eBayes(fit) 
results=topTable(efit,number=6000,adjust.method="BH",p.value=0.05)

La fonction decideTests permet d'accepter ou de rejeter l'hypothèse nulle (décider qu'un gène est différentiellement exprimé ou non) :

decisions = decideTests(efit, p.value=0.05, lfc=1)
# seuil de 0.05 sur la p-valeur
# gènes aux moins 2x plus ou 2x moins exprimés (lfc en log2)

Et l'on peut représenter graphiquement ces décisions à l'aide d'un diagramme de Venn :

vennDiagram(decisions)

Combien de gènes sont différentiellement exprimés à la fois 15 minutes et 30 minutes après décongélation ?

Combien de gènes sont différentiellement exprimés seulement 60 minutes après décongélations ?