silico.biotoul.fr
 

M1 Traitement de Donnees Biologiques - Transcriptome - Oryza sativa

From silico.biotoul.fr

(Difference between revisions)
Jump to: navigation, search
m (Clustering hiérarchique)
Line 234: Line 234:
</source>
</source>
-
Afin de visualiser la similarité/dissemblance entre profils, il est possible de les projeter sur 2 dimensions (à partir des 98) avec une analyse en décomposition des valeurs singulière de la matrice de distance entre chaque paire de profils (ici <tt>d</tt>)
+
Afin de visualiser la similarité/dissemblance entre profils, il est possible de les projeter sur 2 dimensions (à partir des 98) avec une analyse en décomposition des valeurs singulières de la matrice de distance entre chaque paire de profils (ici <tt>d</tt>)
<source lang='rsplus'>
<source lang='rsplus'>
fit=cmdscale(d, eig=TRUE,k=2)
fit=cmdscale(d, eig=TRUE,k=2)

Revision as of 15:19, 17 October 2016

Contents

Préparation de l'environnement

Dans ce TP, nous allons analyser les données de transciptome obtenues sur Oryza sativa par microarray. Le déroulement des analyses suit les travaux de Jujita et al. publiés en 2010 (pmid:21062870). Il vous est conseillé de prendre connaissance de l'article, tout du moins son résumé avant de commencer le TP.

Créez tout d'abord un répertoire de travail sur le bureau (par exemple TDB-TP5) et commencez par télécharger le fichier source que vous allez utiliser et compléter pour générer le compte rendu de TP : M1.TDB.TP_Transcriptome.Rmd (click droit de la souris -- enregistrer la cible sous...). Ouvrez le logiciel RStudio et chargez ce fichier puis lancez sa compilation pour voir le compte rendu. Pour cela cliquez sur le bouton Knit HTML ou bien utilisez la combinaison de touches Ctrl + shift + K.


Récupération des données

Les données relatives à l'analyse ont été mises à disposition sur la banque GEO et sont accessibles avec l'identifiant GSE14304.

Il s'agit de 98 hybridations réalisées avec un microarray de la société Affymetrix ayant l'identifiant GPL2025 comportant 57 381 spots correspondant à 51 279 transcrits d'Oryza sativa japonica et indica.

La quantité de données étant un peu volumineuse pour ce TP, nous allons travailler sur certaines conditions en particulier ; celles spécifiées dans la publication : les gènes différentiellement exprimés entre les stades P1 (Anther, uni-nuclear pollen) et P3 (Anther, tri-cellular pollen).

Oryza.sative.P1-P3.CEL.zip

Chargement et exploration des données

Librairies R/BioConductor nécessaires au TP :

library(affy)
library(limma)

Chargement des données de microarray après analyse d'image :

P1P3.files = c('GSM351444.CEL.gz', 'GSM351445.CEL.gz', 'GSM351446.CEL.gz', 'GSM351450.CEL.gz', 'GSM351451.CEL.gz', 'GSM351452.CEL.gz')
P1P3.raw=ReadAffy(filenames=P1P3.files)

Visualisation et exploration du contenu

P1P3.raw

Distribution :

boxplot(P1P3.raw)

Affichage en boîtes à moustache (y aurait-il besoin d'une normalisation ?)

La commande hist, comme son nom ne l'indique pas, permet de tracer des courbes de densités des distributions des intensités sur la puce.

hist(P1P3.raw)


Une pensée ?

Prétraitement des données

Nous allons utiliser les mêmes paramètres que dans la publication des auteurs afin d'essayer d'obtenir les même résultats : reportez-vous à la section Materials & Methods dans la partie Microarray data extraction, processing and cluster analysis.


Prétraitement : correction du bruit de fond sans normalisation avec les paramètres décrits dans la publi :

P1P3.bgcorr = expresso(P1P3.raw, bgcorrect.method='mas', normalize=F, pmcorrect.method='pmonly', summary.method='mas')

Le paramètre pmcorrect est spécifique à Affymetrix : il ne tient compte que des oligos (probesets) qui ont un perfect match.

Mise à l'échelle au 75ème percentile

P1P3.scaled = exprs(P1P3.bgcorr)
P1P3.scaled[,1] = P1P3.scaled[,1]/ quantile(P1P3.scaled[,1], probs=seq(0,1,0.01))[76]
P1P3.scaled[,2] = P1P3.scaled[,2]/ quantile(P1P3.scaled[,2], probs=seq(0,1,0.01))[76]
P1P3.scaled[,3] = P1P3.scaled[,3]/ quantile(P1P3.scaled[,3], probs=seq(0,1,0.01))[76]
P1P3.scaled[,4] = P1P3.scaled[,4]/ quantile(P1P3.scaled[,4], probs=seq(0,1,0.01))[76]
P1P3.scaled[,5] = P1P3.scaled[,5]/ quantile(P1P3.scaled[,5], probs=seq(0,1,0.01))[76]
P1P3.scaled[,6] = P1P3.scaled[,6]/ quantile(P1P3.scaled[,6], probs=seq(0,1,0.01))[76]

Visualisation du résultat du traitement (la fonction exprs extrait les valeurs d'expression de chaque probeset) :

boxplot(log2(P1P3.scaled), las=3)

ou bien, comme d'habitude, sous forme d'histogramme (des logarithmes des valeurs ; sinon on n'y voit rien) :

hist(log2(P1P3.scaled))

Et enfin, sous forme de courbes :

plot(density(log2(P1P3.scaled)[,5]), col=5) # on commence par celle qui a la valeur maximale (ici, la 5)
lines(density(log2(P1P3.scaled)[,1]),col=1)
lines(density(log2(P1P3.scaled)[,2]),col=2)
lines(density(log2(P1P3.scaled)[,3]),col=3)
lines(density(log2(P1P3.scaled)[,4]),col=4)
lines(density(log2(P1P3.scaled)[,6]),col=6)

Un commentaire ?

  • On observe que la 5 et la 6 se distinguent des autres,
  • ainsi que des valeurs négatives. Comme on a divisé par le 75ème percentile, 75% des valeurs d’expression sont <1, et donc on obtient un log négatif. Sur le boxplot, c’est très net : les boites sont alignées sur le 3ème quartile.

Analyse différentielle

Utilisation de la librairie limma et préparation à l'analyse différentielle : on indique quelle hybridation correspond à quelle condition :

design=matrix( c(rep(1:0,3), rep(0:1,3)), ncol=2, byrow=TRUE)
colnames(design)=c('P1','P3')
design

Analyse :

contrasts=makeContrasts(controlVSah=P1-P3, levels=design)
contrasts
fit = lmFit(P1P3.scaled, design)   # modèle linéaire: calcule la moyenne de P1 et P3 pour chaque probeset
fit2=contrasts.fit(fit, contrasts) # calcul des contrastes (différence des moyennes P1 et P3)
efit = eBayes(fit2)                # modèle bayésien (calcule les statistiques modifiées de t, F, puis p-value)

Puis extraction des probesets ayant une p-valeur significative (seuil de 0.05 sur p-valeur ajustée par la FDR (BH pour Benjamini et Hochberg)) et une variation d’expression supérieure à 10 :

deg=topTable(efit, adjust='BH', p.value=0.05, number=10000, resort.by="logFC", lfc=log2(10))
dim(deg)
head(deg)   # premières lignes du tableau (remarque: logFC correspond ici à la différence des moyennes P1 et P3)

Combien de gènes obtient-on ? pouvait-on s'y attendre ? dans quoi (quel processus biologique) sont-ils impliqués ? Que font-ils ?

Est-ce qu’on a le probeset de la figure 3-B et pas ceux de C et D (car non différentiellement exprimés entre P1-P3) ?

deg[deg$ID=='Os.5780.1.S1_at',]
deg[deg$ID=='Os.55400.2.S1_at',]
deg[deg$ID=='Os.52092.1.S1_at',]

Vérification de la figure 3

Comme le traitement des données prend un peu de temps et de mémoire, le résultat du traitement vous est fourni ici : P1-P3.induced.profiles.txt. Ce fichier contient les profils d'expression sur 33 organes et stades de développement (un total de 98 hybridations), de 150 gènes différentiellement exprimés (avec un ratio d'expression >10) entre P1 et P3.

induced.profiles=read.table('http://silico.biotoul.fr/site/images/c/c8/P1-P3.induced.profiles.txt')
dim(induced.profiles)

Il s'agit ensuite d'afficher les profils en spécifiant les 26 conditions étudiées dans la figure 3 de la publication:

plot(x=1:26, y=rep(0,26), ylim=c(0, max(induced.profiles[,1:26])), pch='', las=3, xaxt='n', xlab='', ylab='')
for (i in 1:150) lines(1:26, induced.profiles[i,1:26], col='green')
axis(1, at=1:26, labels=c('An1.1A','An1.1B','An1.1C','Mei1.2A','Mei1.2B','Mei1.2C','M1.3A','M1.3B','M1.3C','M1.3D','M2.4A','M2.4B','M2.4C','M3.5A','M3.5B','M3.5C','M3.5D','P1.6A','P1.6B','P1.6C','P2.7A','P2.7B','P2.7C','P3.8A','P3.8B','P3.8C'), las=3)

Analyse de clustering

Visualisation des profils d'expression des 150 gènes sur l'ensemble des hybridations (98):

plot(x=1:98, y=rep(0,98), ylim=c(0, max(induced.profiles)), pch='', las=3, xlab='', ylab='')
for (i in 1:150) lines(1:98, induced.profiles[i,], col='green')

Clustering hiérarchique

Chargement de la librairie pour le graphique

library(made4)

Clustering hiérarchique avec la fonction heatplot. Cette fonction effectue le clustering ascendant en fusionnant à chaque étape les 2 clusters les plus proches.

heatplot(induced.profiles, method='ave')

Utilisation du lien simple pour la fusion des clusters :

heatplot(induced.profiles, method='single')

Que remarquez-vous ?

Tentatives avec d'autres méthodes de fusion ascendante des clusters :

heatplot(induced.profiles, method='complete')
heatplot(induced.profiles, method='ward.D',scale='none')
heatplot(induced.profiles, method='ward.D2', scale='none')   # si "ward.D" ou "ward.D2" ne fonctionne pas: "ward"


Laquelle sélectionnerez-vous ?


Calcul de la matrice de distance entre chaque paire de gènes/probesets/profils

d=dist(induced.profiles)

Utilisation de hclust et cutree pour sélectionner les clusters

hc = hclust(d, method='ward.D')
plot(hc)
rect.hclust(hc, k=3) # k=3 : on veut 3 groupes

Récupération des groupes obtenus

membership=cutree(hc, k=3)
table(membership)


On peut, à présent, éventuellement visualiser les profils d'expression

par(mfrow=c(3,1))
c1 = which(membership==1)
plot(x=1:98, y=rep(0,98), ylim=c(0, max(induced.profiles)), pch='', las=3, xlab='', ylab='')
for (i in c1) lines(1:98, induced.profiles[i,], col=1)
 
c2 = which(membership==2)
plot(x=1:98, y=rep(0,98), ylim=c(0, max(induced.profiles)), pch='', las=3, xlab='', ylab='')
for (i in c2) lines(1:98, induced.profiles[i,], col=2)
 
c3 = which(membership==3)
plot(x=1:98, y=rep(0,98), ylim=c(0, max(induced.profiles)), pch='', las=3, xlab='', ylab='')
for (i in c3) lines(1:98, induced.profiles[i,], col=3)

Afin de visualiser la similarité/dissemblance entre profils, il est possible de les projeter sur 2 dimensions (à partir des 98) avec une analyse en décomposition des valeurs singulières de la matrice de distance entre chaque paire de profils (ici d)

fit=cmdscale(d, eig=TRUE,k=2)
plot(fit$points, pch=19)

Cette visualisation bénéficiera de la coloration par cluster/groupe attribué :

plot(fit$points, pch=19, col=membership)

k-means

Nous allons utiliser maintenant la méthode k-means pour détecter les clusters au sein des profils en spécifiant le nombre de groupes à détecter.

km = kmeans(induced.profiles, 3)

Pour visualiser le résultat, en utilisant le graphique précédent et en indiquant les couleurs à utiliser correspondant aux groupes détectés :

par(mfrow=c(1,2))
plot(fit$points, pch='+', col=membership, main='hierarchical')
plot(fit$points, col=km$cluster, pch='+', main='k-means')

Recherche d'enrichissement

Les "gènes" différentiellement exprimés détectés précédemment sont-ils impliqués dans un processus biologique particulier ?

Les probesets différentiellement exprimés entre les stades P1 et P3 sont disponibles dans le fichier précédent ou bien dans P1-P3.induced.probesets.txt


Allez sur le serveur silico.biotoul.fr afin d'effectuer cette recherche en enrichissement d'annotations de la Gene Ontology : http://silico.biotoul.fr/cgi-bin/mobyle/portal.py#forms::sook (s.v.p. mettez une adresse mail "valide", aucune utilisation n’en sera faite, c’est un serveur de l’équipe de R. Barriot. Cela évite de recevoir des mails indiquant que les correspondants ne sont pas joignables).


Réalisez toutes ces analyses et ajoutez-les au compte rendu de TP avant de l'envoyer à votre enseignant par mail (bonhomme@lrsv.ups-tlse.fr ou barriot@biotoul.fr). Le compte rendu est à envoyer avant une semaine. Envoyez les 2 fichiers (.Rmd et .html). Envoyez-vous aussi le mail en copie pour pouvoir vérifier que tout est bien passé. Mettez un titre tel que "Compte rendu TP5 TDB de -et votre Nom et Prénom-".