From silico.biotoul.fr

Revision as of 20:19, 23 October 2011

Analyse de transcriptome de souches de Pseudomonas aeruginosa ayant infecté des patients atteints de mucoviscidose

Etude basée sur la publicaton de Lei Yang et al. Bacterial adaptation during chronic infection revealed by independent component analysis of transcriptomic data.

Récupération des données

Récupérer les données relatives à cette publication sur GEO. Pour cela, il vous faudra retrouver dans la publication les références du jeu de données qu'ils ont utilisé.

Normalisation

Le format est celui d'Affymetrix (fichiers .CEL). Nous utiliserons donc la librarie affy de R/Bioconductor.

Après avoir désarchivé les données dans un répertoire, changer le répertoire courant dans R pour se placer dans le répertoire contenant les fichiers .CEL. L'ensemble des données contenues dans les fichiers GSMxxx.CEL est chargé sous R avec les commandes suivantes :

# chargement de la librairie
library(affy)
raw=ReadAffy()

En effet, par défaut, la fonction ReadAffy() charge tous les fichiers .CEL du répertoire courant.

On peut s'assurer du besoin de normaliser avec une boite à moustache :

boxplot(raw)

Pour la normalisation, nous utiliserons la fonction rma() (Robust Multiarray Averaging) de la librairie affy :

raw.rma=rma(raw)

L'objet raw.rma est de type ExpressionSet et contient l'ensemble des données normalisées. Nous nous intéresserons uniquement aux intensités que nous récupérons avec la commande suivante :

rma.expr=exprs(raw.rma)

Pour plus d'information sur le type ExpressionSet vous pouvez consulter l'aide :

?ExpressionSet

Le résultat de la normalisation peut aussi se visualiser avec une boite à moustache :

par(mfrow=c(2,1))
boxplot(raw)
boxplot(rma.expr)

Données annexes sur les hybridations

Afin de plus facilement analyser les données, les informations de la publication concernant le stade ou phénotype (early, late, mucoid) des bactéries a été compilé dans le fichier suivant :

Charger ces données dans un objet que l'on appellera info.

Clustering des profils

Pour effectuer des analyses de clustering, nous utiliserons la fonction heatplot de la librairie made4. Cette fonction reprend les fonctionnalités de la fonction heatmap vue dans un TD précédent. Elle permet notamment d'utiliser une méthode de distance entre profils d'expression basé sur le coefficient de corrélation de Pearson.

Effectuer un premier clustering avec les paramètres par défaut :

# chargement de la librairie
library(made4)
# on revient à un seul graphique
par(mfrow=c(1,1))
# affichage du clustering
heatplot(rma.expr)

Si l'affichage est lent, vous pouvez n'afficher que le dendrogramme sur les colonnes (hybridations) :

heatplot(rma.expr, dend='column')

Nous allons ensuite changer les couleurs pour un meilleur contraste :

heatplot(rma.expr, dend='column', lowcol="blue",highcol="yellow")

Puis, ajouter les informations concernant les puces que vous avez normalement chargé dans l'objet info :

heatplot(rma.expr, dend='column', lowcol="blue",highcol="yellow", ColSideColors=as.character(info$color),cols.default=FALSE)

• Est-ce que les réplicats sont regroupés dans les mêmes clusters ?

• Est-ce que les différents stades/phénotypes forment de clusters ?

Par défaut, la méthode pour mesurer la distance entre 2 profils est basée sur le coefficient de corrélation de Pearson (1 - r). Essayer avec une autre méthode de distance : la distance euclidienne

heatplot(rma.expr, dend="column", lowcol="blue",highcol="yellow", ColSideColors=as.character(info$color),cols.default=FALSE, distfun="euclidean")

• Que pensez vous du résultat : vous paraît-il meilleur ?

Par défaut, la méthode pour identifier les clusters à fusionner est average linkage. Essayer avec une autre méthode de distance entre clusters : single linkage

heatplot(rma.expr, dend="column", lowcol="blue",highcol="yellow", ColSideColors=as.character(info$color),cols.default=FALSE, method="single")

• Que pensez vous du résultat : vous paraît-il meilleur ?

Essayer également avec la méthode complete linkage :

heatplot(rma.expr, dend="column", lowcol="blue",highcol="yellow", ColSideColors=as.character(info$color),cols.default=FALSE, method="complete")
heatplot(rma.expr, dend="column", lowcol="blue",highcol="yellow", ColSideColors=as.character(info$color),cols.default=FALSE, method="ward")

Consultez l'aide de la fonction heatplot. Y a-t-il d'autres méthodes disponibles ? Essayez-les et notez le résultat qui vous semble le meilleur.

Préparation des données

Chaque expérience a été menée en triplicat. Afin de réaliser le clustering, nous allons faire la moyenne des intensités obtenues pour chaque réplicat :

rma.mean=by(t(rma.expr), INDICES=info$Experience, FUN=mean)
mat=matrix(unlist(rma.mean),nrow=nrow(rma.expr),byrow=F,dimnames=list(rownames(rma.expr),names(rma.mean)))

La fonction by() de R permet d'appliquer une fonction (ici mean) par modalité d'un facteur (ici correspondant à info$Experience dans lequel on a attribué le même nom aux réplicats). Le résultat est un objet un peu commode : nous avons pour chaque expérience, un vecteur contenant le niveau d'expression pour chaque spot. Nous re-transformons cette objet en matrice spots x conditions avec la fonction matrix.

Afin de pouvoir utiliser les couleurs affectées à chaque stade/phénotype, il nous faut sélectionner une ligne sur trois de l'objet info :

# récupération d'une ligne sur trois
mat.info=info[ seq(1,78,by=3), 2:3]
# affectation du nom de puce aux ligne
rownames(mat.info)=mat.info$Experience
# affectation de la couleur en fonction du nom de la puce 
mat.colors=mat.info[ colnames(mat), "color" ]

Vous pouvez maintenant réeffectuer des clustering en faisant varier le paramètre method et le paramètre distfun. Lesquels vous donnent le meilleur résultat ?

Avant de passer à un autre logiciel, nous avons une dernière manipulation à effectuer : notre matrice d'expression contient les noms des probesets et non ceux des gènes. Pour récupérer les noms des gènes nous aurons besoins du fichier de GEO décrivant le microarray (ici GPL84). Récupérer ce fichier au format soft sur GEO. Puis charger le de la manière suivante :

# on ne tient pas compte des 27 premières constituant une en-tête (paramètre skip)
# on ne garde que les colonnes 1 (petite descr), 2 (nom le + commun), et 11 (identifiant systématique, celui qui nous intéresse)
GPL=read.table("GPL84.annot", header=TRUE, sep="\t", skip=27, fill=TRUE, row.names=1, quote="")[,c(1,2,11)]
# on se débarrase de la dernière ligne qui contient !platform_table_end                                    
GPL=GPL[-5901,]

La colonne 11 (platform_ORF) que nous avons gardé contient soit l'identifiant du gène soit rien. Nous allons donc l'utiliser pour ne garder que les gènes :

rownames(mat) = GPL$Platform_ORF
mat.gene = mat[ rownames(mat)!='', ]
# on sauvegarde pour l'ouvrir dans MeV
write.table(mat.gene,"GSE31227_mean.txt", col.names=TRUE, row.names=TRUE, quote=FALSE, sep="\t")

Clustering des profils avec MeV

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.

Charger ensuite le jeu de données que vous avez créé à l'étape précédente (GSE31227_mean.txt).

Essayer de transformer les données (menu Adjust data et d'ajuster les couleurs (menu Display) pour un meilleur rendu.

Effectuer des clustering (k-means et hiérarchiques avec différents paramètres).

Sélectionner des clusters, notamment ceux pour lesquels l'expression des gènes diffèrent entre les stades/phénotypes early, late et mucoid.

Analyser ensuite ces clusters à l'aide du programme sook que vous trouverez sur l'environnement mobyle hébergé au LMGM.

Analyse d'un compendium de données de transcriptome de levure Saccharomyces cerevisiae