silico.biotoul.fr
 

InfoBio TD transcriptome

From silico.biotoul.fr

Jump to: navigation, search

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 design de microarray différents ?

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 être identiques (se 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éplétitions 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()
 
# puis essayer, par exemple :
bg=bg.correct(cel, method='rma') 
# ou 
bg=bg.correct(cel, method='mas') 
#vérifier le résultat avec 
hist(bg)
 
# puis pour la normalisation
normalize.methods(bg)
# exemple 
no=normalize(bg, method='qspline')
hist(no)

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)
# limmaUserGuide() permet d'accéder au manuel utilisateur de la librairie
# conditions à comparer
design=matrix( c(rep(0:1,7), rep(1:0,15)), ncol=2, byrow=TRUE)
colnames(design)=c('control','AH')
# plan d'expérience
design
 
# conditions comparées
cont.matrix=makeContrasts(controlVSah=control-AH, levels=design)
 
# extraction et transformation (log2) des valeurs d'expression (autres méthodes avec pmcorrect.methods() et express.summary.stat.methods()
eset <- computeExprSet(no, pmcorrect.method='pmonly', summary.method='avgdiff')
# transformation en log2
expr=log2(exprs(eset))
 
# identification des gènes différentiellement exprimés (>5 fold, p-valeurs<=0.05)
fit = lmFit(expr, design)
fit2=contrasts.fit(fit, cont.matrix)
efit = eBayes(fit2)
# gènes différentiellement exprimés pour un seuil à 5%, au moins 5x de variation et ajustement FDR (BH) pour tests multiples
deg=topTable(efit, adjust='BH', p.value=0.05, number=1000, sort.by="logFC", lfc=log2(5))
pids=deg[,1] # probeset ids des spots

A ce stade, combien obtenez vous de gènes différentiellement exprimés ?

Nous allons maintenant convertir les identifiants des spots de la puce (probeset ids) en nom de gènes :

# récupération des informations sur les spots
cel # pour avoir le nom et design de la puce
cleancdfname("HG-U133_Plus_2") # pour connaître son identifiant bioconductor
source("http://bioconductor.org/biocLite.R") # script d'installation bioconductor
biocLite(c("hgu133plus2.db", "GO.db")) # installation des librairies 
library("hgu133plus2.db") # chargement de la librairie/données sur le microarray
 
cols(hgu133plus2.db) # liste des informations que l'on peut récupérer
 
keytypes(hgu133plus2.db) # champs qui peuvent servir de clé
 
# requête de sélection/conversion des identifiants de spot -> nom de gènes ou autres, cf. help(select)
unique(select(hgu133plus2.db, pids, "GENENAME", "PROBEID")[,2])
 
#conversion au format swissprot : cols(hgu133plus2.db)
uniprot=unique(select(hgu133plus2.db, pids, "UNIPROT", "PROBEID")[,2])
 
# écriture des identifiants dans un fichier
write.table(uniprot, 'deg.SP', quote=FALSE, row.name=FALSE)

Analyse et interprétation des relations entre les gènes et leurs produits sur http://string-db.org/ onglet multiple names pour l'organisme Homo sapiens.

Références et données

Installation de librairies

source("http://bioconductor.org/biocLite.R")

biocLite("affy")

biocLite('org.Hs.eg.db')

biocLite('AnnotationDbi')

biocLite('hgu133plus2.db')