InfoBio TD transcriptome
From silico.biotoul.fr
Contents |
Aperçu de la séance
- Récupération d'un jeu de données brutes
- Premier pas avec le logiciel R
- Chargement des données d'expression
- Filtrage et normalisation
- Identification des gènes différentiellement exprimés (au moins 5x)
- Analyse de la liste de gènes obtenue
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).
Combien d'hybridations sont disponibles sur la banque GEO ? Pour combien de puces différentes ?
Au cours de cette séance, nous allons analyser les données telles que fournies par le logiciel d'analyse d'image. Il s'agit de séries d'hybridations obtenues dans le cadre d'une étude d'hépatites sévères liées à la consommation d'alcool. Commencez par lire le résumé de la publication associée ; son identifiant PubMed est le 22637703.
Combien d'hybridations ont été utilisées dans cette étude ? Combien de gènes différentiellement exprimés (>x5 fois) ont été identifiés ?
A partir de GEO du NCBI, retrouvez et téléchargez la totalité des données associées à la publication (il s'agit de la série GSE28619). Ensuite, décompressez le fichier (avec les commandes tar/gunzip sous linux par exemple).
Avec quelle version de microarray ont été obtenues les données ? Combien y a-t-il de spots sur la puce ? Combien d'hybridations sont disponibles pour cette version de microarray ? Regroupées en combien de séries d'hybridations ?
Premiers pas avec R
R est un environnement de statistiques utilisé pour des recherches et applications variées (physique, astrophysique, écologie, biologie moléculaire, génomique, ...). Pour plus d'info : http://www.r-project.org et Aide mémoire des commandes R
Type de données :
- scalar : une valeur de type numérique, chaîne de caractères, booléenne, nominale.
- vector : tableau à 1 dimension de scalaires de même type.
- matrix : tableau à 2 dimensions de scalaires de même type.
- array : tableau à n dimensions de scalaires de même type.
- data frame : tableau à 2 dimensions, les colonnes peuvent avoir des types différents.
- list : tableau à une dimension pouvant contenir des valeurs de types différents.
Pour se familiariser avec la syntaxe, nous allons :
- générer aléatoirement 2 séries de valeurs suivant une loi normale avec des paramètres différents
- tracer des histogrammes et courbes de densité
- faire un test de comparaison de moyennes
Pour lancer le logiciel, lancer la commande R dans un terminal.
Génération des 2 séries de nombres aléatoires :
# Première série serieA=rnorm(20, 50, 10) # échantillon de 20 valeurs d'une population ayant pour moyenne 50 et écart-type 10 serieA # affichage du contenu mean(serieA) ; sd(serieA) # moyenne, écart-type summary(serieA) boxplot(serieA) hist(serieA, breaks=10) # histogramme de la première série plot(density(serieA), col='orange', xlim=(c(0,100))) abline(v=mean(serieA), col='orange') # moyenne observée abline(v=50, col='red') # moyenne théorique xi=1:100 ; lines(xi, dnorm(xi, 50, 10), col='red') # distribution théorique # Deuxième série serieB=rnorm(25, 60, 20) lines(density(serieB), col='blue') abline(v=mean(serieB), col='blue') lines(xi, dnorm(xi, 60, 20), col='pink') abline(v=60, col='pink') # comparaison des moyennes : # test de normalité des échantillons shapiro.test(serieA) shapiro.test(serieB) # test d'homoscédasticité (comparaison des variances) var.test(serieA, serieB) # test de student t.test(serieA, serieB) # faire t.test(serieA, serieB, var.equal=TRUE) si les échantillons alétaoires ont même variance
Chargement et visualisation des données
Vous avez normalement téléchargé un fichier tar contenant l'ensemble des hybridations au format CEL (format développé par la société Affymetrix qui commercialise des microarray). Après décompression, vous devriez obtenir 22 fichier avec l'extension .CEL.gz
Pour charger et analyser ces données dans R, nous allons utiliser la librairie affy de la suite R/Bioconductor.
library(affy) # chargement de la librairie # si la librairie n'est pas installée, faire : source("http://bioconductor.org/biocLite.R") ; biocLite("affy") ; library(affy) cel=ReadAffy() # lecture des fichiers .CEL.gz contenus dans le répertoire de travail
Données concernant chaque hybridation :
summary(cel) cel boxplot(cel) # distributions sous forme de boite à moustache hist(cel) # sous forme de courbes de densité
Sous l'hypothèse que la grande majorité des gènes ont le même niveau d'expression dans les 2 conditions (foie sain vs. foie malade), les distributions devraient ce superposer. Or, ce n'est pas le cas car les données ont besoin d'être normalisées.
Pour plus d'information sur Biocondutor : http://bioconductor.org et la librairie affy : http://bioconductor.org/packages/release/bioc/html/affy.html notamment les PDF "Primer"
Les 7 premières hybridations correspondent aux foies sains (réplicats biologiques). On peut tracer des MA-plot de la manière suivantes :
MAplot(cel[,1:3], pairs=TRUE, plot.method="smoothScatter") # on ne trace ici que pour les 3 premières hybridations
Contrôle qualité : Avant de réellement utiliser les données, il faut s'assurer que les hybridations se sont bien déroulées (pas d'erreur systèmatiques de type microarray défectueux par exemple, ou de corrélation spatiale des intensités). Affichage des intensités sur la puce :
image(cel[,1]) # Première hybridation, à vous de tester sur les autres
Normalisation des données
Pour la normalisation, nous avons plusieurs méthodes à notre disposition. Il existe aussi différentes manières de prendre en compte le bruit de fond.
Etudier l'aide de la fonction et des méthodes :
# pour les méthodes disponibles des différentes étapes : bgcorrect.methods() normalize.methods(cel) # puis essayer, par exemple : expr=bg.correct.mas(cel) hist(expr) expr=normalize(expr, method='constant') hist(expr)
Identification des gènes différentiellement exprimés
Cette étape sera réalisée avec une autre librairie de R/Bioconductor : limma (pour Linear Model for Microarray Analysis).
Chargement de la librairie :
library(limma) # si elle n'est pas installée : source("http://bioconductor.org/biocLite.R") ; biocLite("limma") ; library(limma)
# conditions à comparer design=matrix( c(rep(0:1,7), rep(1:0,15)), ncol=2, byrow=TRUE) colnames(design)=c('control','AH') design cont.matrix=makeContrasts(controlVSah=control-AH, levels=design) # extraction et transformation (log2) des valeurs d'expression eset <- computeExprSet(expr, pmcorrect.method='pmonly', summary.method='avgdiff') # 1000 ! eset=log2(exprs(eset)) # 37 deg # identification des gènes différentiellement exprimés (>5 fold, p-valeurs<=0.05) fit = lmFit(eset, design) fit2=contrasts.fit(fit, cont.matrix) efit = eBayes(fit2) deg=topTable(efit, adjust='BH', p.value=0.05, number=1000, sort.by="logFC", lfc=log2(5))
