M1 Traitement de Donnees Biologiques - Transcriptome - Oryza sativa
From silico.biotoul.fr
m |
m (→Prétraitement des données) |
||
| Line 85: | Line 85: | ||
<span style='color: #990000'> | <span style='color: #990000'> | ||
Un commentaire ? | Un commentaire ? | ||
| + | </span> | ||
| + | |||
| + | * 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 <tt>limma</tt> et préparation à l'analyse différentielle : on indique quelle hybridation correspond à quelle condition : | ||
| + | <source lang='rsplus'> | ||
| + | design=matrix( c(rep(1:0,3), rep(0:1,3)), ncol=2, byrow=TRUE) | ||
| + | colnames(design)=c('P1','P3') | ||
| + | design | ||
| + | </source> | ||
| + | |||
| + | Analyse : | ||
| + | <source lang='rsplus'> | ||
| + | contrasts=makeContrasts(controlVSah=P1-P3, levels=design) | ||
| + | contrasts | ||
| + | fit = lmFit(P1P3.scaled, design) # modèle linéaire | ||
| + | fit2=contrasts.fit(fit, contrasts) # encore un | ||
| + | efit = eBayes(fit2) # modèle bayésien | ||
| + | </source> | ||
| + | |||
| + | 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 : | ||
| + | <source lang='rsplus'> | ||
| + | deg=topTable(efit, adjust='BH', p.value=0.05, number=10000, sort.by="logFC", lfc=log2(2)) | ||
| + | dim(deg) | ||
| + | </source> | ||
| + | |||
| + | <span style='color: #990000'> | ||
| + | Combien de gènes obtient-on ? pouvait-on s y attendre ? dans quoi (quel processus biologique) sont-ils impliqués ? Que font-ils ? | ||
</span> | </span> | ||
Revision as of 14:18, 16 October 2016
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.
Contents |
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 et P3.
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
summary(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 tien compte que des oligos (probesets) qui ont un perfect match.
Visualisation du résultat traitement (la fonction exprs extrait les valeurs d'expression de chaque probeset) :
boxplot(log2(exprs(P1P3.bgcorr)), las=3)
ou bien, comme d'habitude, sous forme d'histogramme (des logarithmes des valeurs ; sinon on n'y voit rien) :
hist(log2(exprs(P1P3.bgcorr)))
Et enfin, sous forme de courbes :
plot(density(log2(exprs(P1P3.bgcorr)[,5])), col=5) # on commence par celle qui a la valeur maximale (ici, la 5) lines(density(log2(exprs(P1P3.bgcorr)[,1])),col=1) lines(density(log2(exprs(P1P3.bgcorr)[,2])),col=2) lines(density(log2(exprs(P1P3.bgcorr)[,3])),col=3) lines(density(log2(exprs(P1P3.bgcorr)[,4])),col=4) lines(density(log2(exprs(P1P3.bgcorr)[,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 fit2=contrasts.fit(fit, contrasts) # encore un efit = eBayes(fit2) # modèle bayésien
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, sort.by="logFC", lfc=log2(2)) dim(deg)
Combien de gènes obtient-on ? pouvait-on s y attendre ? dans quoi (quel processus biologique) sont-ils impliqués ? Que font-ils ?
