Analyse de transcriptome

Partie 1. Premiers pas avec R et bioconductor, importation, filtrage, normalisation et analyse différentielle des données

Cette partie est largement inspirée d'une formation à l'analyse de données de microarray de David Causeur (Laboratoire de Mathématiques Appliquées, Rennes). Pour une prise en main de R, il existe un tutorial très bien fait.

La cryopréservation est une méthode largement utilisée pour le stockage à long terme de nombreuses cellules vivantes. Cette méthode implique des traitements de congélation et de décongélation qui causent des dommages aux cellules vivantes voire souvent à leur mort. Les puces à ADN permettent de mesurer les niveaux d'expression de milliers de gènes simultanément et de suivre les réponses biologiques à travers le niveau d'expression de presque tous les gènes de l'organisme. Nous allons analyser la réponse de cellules de levure après cryopréservation à partir de données de transcriptome publiées dans une étude de l'impact de la crypréservation sur la levure [Odani et al. 2003].

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). A partir de GEO du NCBI, retrouvez les données associées à la publication [Odani et al. 2003].

Vous devriez trouver la série GSE9404. Lisez la page décrivant les données, puis téléchargez les données brutes (GSM239212..GSM239220) dans un répertoire que vous aurez créé pour nos analyses. Ensuite, décompressez ces fichiers (soit avec la commande gunzip sous linux, soit avec un logiciel de compression/décompression tel que 7-Zip.

Installation de bioconductor

Pour la réalisation de ces analyses, nous utiliserons la suite bioconductor qui est un ensemble de packages pour le logiciel de traitements statistiques R. L'installation de bioconductor se fait très simplement aux moyens des deux commands suivantes dans une session R :

Importation des données avec limma

Le package limma est spécialement conçu pour analyser les données de biopuces bi-couleurs. Il permet l'importation de fichiers de sortie des logiciels d'analyse d'image les plus courants dans le domaine des biopuces, la normalisation des données et l'analyse différentielle (pour une présentation détaillée, voir la page web http://bioinf.wehi.edu.au/limma/).

Le chargement du package est ensuite commandé par :

Remarque : limma dispose d'un guide utilisateur accessible et pédagogique, incluant notamment de nombreux exemples d'utilisation. Pour y accéder :

Afin d'éviter d'avoir à indiquer de manière répétée le répertoire de travail au cours de la procédure d'importation des données, on peut le définir comme répertoire de travail pour l'ensemble de la session :

Le format des données de biopuces nécessaire à toute analyse statistique est un tableau dont les lignes sont associées à des puces et les colonnes à des gènes. Or, le plus souvent le point de départ de l'analyse repose sur autant de fichiers de sortie d'un logiciel d'analyse d'images qu'il y a de puces. Le package limma prévoit donc une procédure d'importation de ces fichiers conduisant à la création d'un tableau de données. La conversion de plusieurs fichiers en un seul tableau pucesxgènes nécessite quelques indications permettant de piloter l'ordonnancement des valeurs d'expressions géniques dans le jeu de données. Ces informations sont contenues dans un fichier supplémentaire nommé dans la suite targets.txt.

Le fichier targets.txt

Le fichier targets.txt doit contenir les colonnes suivantes :

• FileName donnant le nom des fichiers d'images,
• Cy3 donnant la modalité du facteur marquée en vert,
• Cy5 donnant la modalité du facteur marquée en rouge.

Les autres colonnes sont optionnelles. On donne le fichier targets_15min.txt suivant correspondant à la première série de réplicats effectués 15 minutes après décongélation.

Créez de la même manière les fichiers correspondant à 30 et 60 minutes après décongélation.

La commande suivante permet de lire le fichier targets.txt et de créer un objet R, nommé ici targets, héritant de l'information contenue dans le fichier :

L'argument sep définit le type de séparateur utilisé pour délimiter les colonnes dans le fichier targets.txt. Ici, on précise que le séparateur est la tabulation (\t).

Chargez les trois fichiers targets dans trois variables (targets15, targets30 et targets60).

Remarque : pour obtenir de l'aide sur une commande ou fonction :

Filtrage des spots

Lors de la lecture des fichiers de sortie de logiciels d'analyse d'image, certaines fonctionnalités de limma permettent de contrôler la qualité des données d'expressions en chaque spot. En effet, ces fichiers contiennent eux-mêmes différentes indications sur la qualité des spots. Par exemple, les fichiers Genepix sont structurés comme des tableaux de données dont les lignes correspondent aux spots et les colonnes à des caractéristiques des données en chacun de ces spots. On trouve parmi ces caractéristiques :

• Flags : type de Flag
• F635 Median : intensité médiane pour le signal rouge
• F532 Median : intensité médiane pour le signal vert
• F635 Mean : intensité moyenne pour le signal rouge
• F532 Mean : intensité moyenne pour le signal vert

Il est donc possible de créer sa propre fonction de filtrage de spots : cette fonction, dont l'argument principal est le tableau spots x caractéristiques d'un fichier Genepix, prendra la valeur 1 si un spot pest valide et 0 sinon. Par exemple, si on souhaite supprimer les spots dont les flags de Genepix sont inférieurs ou égaux à -49, il suffit de créer la fonction de filltrage suivante :

De la même manière, il est parfois souhaitable de filtrer les spots présentant une forte hétérogénéité du signal observé. Cette hétérogénéité se traduit par une grande différence entre le signal moyen et le signal médian. Un test permettant d'identifier simplement les cas d'hétérogénéité consiste donc à s'assurer que le rapport H suivant n'est pas trop grand :

H = ( |signal moyen - signal médian| ) / [ (signal moyen + signal médian)/2 ]

Si H_threshold désigne la valeur limite au-delà de laquelle le signal est considéré comme hétérogène, les commandes suivantes permettent de tester l'homogénéité des signaux rouge et vert :

Complétez la fonction myFilter afin d'intégrer le test d'hétérogénéité.

Lecture des fichiers Genepix

La fonction read.maimages permet de lire les fichiers Genepix :

Chargez les trois ensembles de puces dans trois variables (RG15, RG30 et RG60).

Tel que défini ci-dessus, l'objet RG peut-être vu comme une armoire contenant toute l'information contenue dans les fichiers Genepix. Les noms des différents tiroirs de cette armoire sont les suivants :

Plus précisément,

• RG$R est un tableau gènes x puces contenant les valeurs du signal rouge,
• RG$G est un tableau gènes x puces contenant les valeurs du signal vert,
• RG$Rb est un tableau gènes x puces contenant les valeurs du bruit de fond rouge,
• RG$Gb est un tableau gènes x puces contenant les valeurs du bruit de fond vert,
• RG$weights est un tableau gènes x puces contenant les valeurs 1 pour les spots conservés lors du filtrage et 0 sinon,
• RG$targets contient les information du fichier targets.txt,
• RG$genes est lui-même une armoire dont les différents tiroirs contiennent des informations sur chaque gène,
• RG$source rappelle le logiciel d'analyse d'images d'où viennent les données,
• RG$printer décrit la structure d'une puce (taille d'un bloc d'aiguilles, ...).

Taux de filtrage

On peut utiliser RG$weights pour mesurer l'impact de la procédure de filtrage. Par exemple, la commande suivante calcule, pour chaque puce, le pourcentage de gènes non-filtrés :

Nom et identifiant d'un gène

Les noms des différents tiroirs de l'armoire nommée RG$genes sont accessibles par la commande suivante :

Pour obtenir les identifiants des 10 premiers gènes, on peut utiliser le tiroir nommé RG$genes$ID :

Dans les fichiers .gpr que vous avez téléchargés, les noms des gènes ne sont pas renseignés, ils sont en fait déjà indiqués au niveau des identifiants. Affectez les identifiants aux noms de gènes.

Statut d'un spot

Dans l'optique d'une meilleure lisibilité des résultats de l'analyse, il peut être intéressant de tenir compte d'une typologie des spots à partir de leur statut. : gène, blanc, contrôle, ... On utilise pour cela un fichier texte (les séparateurs de colonnes sont des tabulations) appelé spotTypes.txt. Pour les données dont on dispose, ce fichier prend la forme du tableau suivant.

La dernière colonne du fichier, nommée Color, permet d'associer une couleur à chaque type de spot, afin de mieux les distinguer dans les analyses à venir. Les commandes suivantes lisent le fichier spotTypes.txt et associent à chaque spot son statut dans RG$genes :

Normalisation des données

La normalisation des données est une étape importante dans la création du tableau de données sur lequel reposera l'analyse statistique. L'objectif de cette normalisation est d'identifier des sources parasites de variations des expressions. Parmi ces sources, on cible plus précisément :

• l'hétérogénéité du bruit de fond.
  Si le bruit de fond présente des variations très différentes d'une puce à une autre, ou très structurées spatialement sur une puce, alors on peut être amené à corriger le signal par soustraction du bruit de fond.
• l'hétérogénéité du signal.
  De la même manière, le principe d'invariance d'une très grande majorité des expressions géniques d'une puce à une autre doit se traduire par une répartition comparable des valeurs des signaux entre les différentes puces. Si des différences marquées existent, il est judicieux de ramener les signaux moyens de chaque puce à la même valeur. D'autre part, il arrive que le fluorochrome soit lui-même source de variations du signal. Cet effet du fluorochrome dépend même parfois de la valeur d'expression et de la structure de la puce (partitionnement en blocs d'aiguilles).

Hétérogénéité du bruit de fond

De manière générale, il est important de s'assurer que les variations du bruit de fond des images analysées par Genepix ne sont pas structurées :

• d'une part, par une comparaison des variations du bruit entre les différentes puces.
  La commande suivante construit une boîte de dispersion du logarithme du bruit de fond pour les puces :

• d'autre part, par un examen de la répartition spatiale du bruit de fond de chaque puce. La commande suivante crée une image associant à chaque spot de la 3^ième puce un pixel dont le niveau de couleur est proportionnel à la valeur du logarithme du bruit de fond :

La correction du signal par soustraction du bruit de fond peut être obtenue par la fonction backgroundCorrect.

Hétérogénéité du signal

L'hétérogénéité des signaux mesurés sur une puce est souvent attribuée à des variations uniquement liées aux fluorochromes. Cette différence entre les signaux est décelable sur le nuage dont les points, représentant des spots, ont pour abscisses les moyennes entre les logarithmes des signaux verts et rouges (dites valeurs A) et pour ordonnées les différences entre logarithmes des signaux verts et rouges (dites valeurs M). Ce graphique, dit graphe MA, est obtenu par la fonction plotMA :

La fonction normalizeWithinArrays corrige cette hétérogénéité en ajustant un modèle des variations des valeurs M en fonction des valeurs A par une méthode de régression locale (dite méthode loess) :

Le résultat de la commande précédente, l'objet MA, ressemble en tous points à l'objet RG : seuls les tiroirs nommés R et G et contenant les signaux rouges et verts sont remplacés par d'autres, nommées M et A et contenant les valeurs du même nom. Le paramètre bc.method indique la méthode utilisée pour corriger le bruit de fond. La méthode par défaut est subtract qui soustrait les valeurs de bruit de fonds aux intensités. Il est parfois préférable, notamment pour l'analyse différentielle, d'utiliser la méthode normexp qui ajuste les intensités de manière adaptive au bruit de fond et mène à des valeurs d'intensité strictement positives. On peut en plus passer le paramètre offset qui sera ajouté à chaque intensité avant la transformation et qui permet de tirer vers 0 les log-ratios pour les faibles intensités. La fonction RG.MA permet d'obtenir les signaux verts et rouges normalisés. L'effet de la normalisation intra-puces est affiché par la commande suivante :

Le graphique obtenu montre bien un autre effet de la normalisation décrite ci-dessus : le centrage par puces des valeurs M. La commande suivante permet de visualiser la dispersion des valeurs de M dans les différentes puces.

La fonction plotDensities affiche les densités empiriques lissées pour les intensités rouges et vertes des puces. Sans normalisation, on observe une grande variabilité entre les canaux et entre les puces.

Après normalisation loess des valeurs M pour chaque puce, les distributions des intensités rouges et vertes devraient être essentiellement les mêmes pour chaque puce. Cependant, il peut être observé une grande variabilité entre les puces. La normalisaiton loess n'affecte pas les valeurs A. La fonction normalizeBetweenArrays permet de normaliser les distributions d'intensité entre les puces :

Essayez différentes méthodes de correction du bruit et de fond, et de normalisation (intra-puce), jusqu'à obtenir des distributions d'intensité satisfaisantes.

Création de la matrice gènesxpuces

La normalisation des données a conduit à la création d'une matrice gènes x puces de valeurs d'intensité :

Gènes différentiellement exprimés

Pour décider qu'un gène est différentiellement exprimé entre deux conditions, on va tester si les moyennes obtenues pour les différents réplicats sont égales au moyen du test de Student (t-test).

Evaluation du dispositif expérimental

Il est intéressant de voir quelle différence inter-groupe le dispositif expérimental permet de mettre en évidence. Pour cela, il est nécessaire de connaître l'ordre de grandeur des écart-types intra-groupes :

L'objet sdbyclass contient deux parties, nommées respectivement sdbyclass$R et sdbyclass$G. On obtient une rapide description de ces éléments de la manière suivante :

A partir de cette description, déterminez un seuil au dessous duquel on aura une grande majorité des écarts-types intra-groupe inférieure à ce seuil. Sur la base de ce seuil, la fonction power.t.test peut être utilisée pour calculer la différence de moyennes entre les deux groupes que l'on peut mettre en évidence avec un risque d'erreur de 5% :

Déterminez quel ordre de changement de moyenne il est possible de détecter avec un niveau de confiance de 95%.

Test de comparaison de deux moyennes

Le choix de la stratégie de test de comparaison de deux moyennes dépend de la validité de l'hypothèse d'égalité des variances. On va donc procéder à un test de comparaison de variances des expression géniques avant d'effectuer le test de comparaison de moyenne adéquat. On obtient la probabilité critique que les variances soit égales de la manière suivante :

Nous allons maintenant définir une fonction pour l'appliquer à chaque colonne de notre matrice gène x puce.

La fonction apply permet d'appliquer notre fonction compareMean à l'ensemble des lignes ou des colonnes d'un tableau. Le paramètre MARGIN permet de spécifier si on l'applique aux lignes (MARGIN=1) ou aux colonnes (MARGIN=2). On l'utilise donc de la manière suivante :

Une façon de visualiser les résultats de ces tests est d'afficher un volcano plot. Pour cela, il est nécessaire de calculer la moyenne des logarithmes de variation d'expression avant d'afficher le graphe :

On s'intéresse maintenant aux gènes deux fois plus ou deux fois moins exprimés ayant une p-valeur significative :

Analyse différentielle avec limma

limma permet aussi d'identifier les gènes différentiellement exprimés mais au moyen de régression linéaire et de modèle Bayésiens (entre autres). Vous allez comparer les résultats obtenus avec limma avec vos résultats précédents. Pour analyser vos puces avec limma, créez d'abord un fichier targets pour l'ensemble des puces puis utilisez les commandes suivantes :

Caractérisation d'un ensemble de gènes

Nous allons à présent analyser les gènes précédemment identifiés comme différentiellement exprimés. Pour cela, nous allons rechercher les annotations surreprésentées au sein de cet ensemble de gènes. Le principe est de tester si la fréquence d'une annotation est largement supérieure à la fréquence attendue (en générale la fréquence observée pour le génome) et d'effectuer ce test pour chaque annotation disponible.

A partir du mapping mots-clés SwissProt - gènes annotés, écrire un script Perl qui recherche les mots-clés surreprésentés au sein d'un ensemble de gènes. Dans un premier temps, vous utiliserez la loi hypergéométrique sans correction pour test multiple. Puis, vous ajouterez la correction de Benjamini et Hochberg (FDR). On vous donne les fonction pour le calcul de la p-valeur avec la loi hypergéomtétrique et la correction de ces p-valeurs avec la FDR suivantes :

sub cnp { my $n = shift; my $p = shift; if ($p>$n) { return 0; } elsif ($p==0) { return 1; } elsif ($p==1) { return $n; } elsif ($p==$n) { return 1; } else { if ($p>$n/2) { $p=$n-$p; } my $res=$n; for (my $i=2;$i<=$p;$i++) { $res = ($res * ($n-$i+1) )/$i; } return $res; } } sub hypergeometric { my $s = shift; # population size my $q = shift; # sample1 size my $t = shift; # sample2 size my $c = shift; # intersection size of sample1 and sample2 my $res = 0; my $min=$q; if ($t < $q) { $min=$t; } my $cnp_s_q = cnp($s,$q); for (my $i=$c ; $i <= $min ; $i++) { $res += (cnp($t,$i) * cnp($s-$t, $q-$i) / $cnp_s_q); } return $res; }

sub min3 { my $m = shift; my $a = shift; my $b = shift; if ($a<$m) { $m = $a; } if ($b<$m) { $m = $b; } return $m; } sub FDR { my $alpha = shift; my $h = shift; my @hits = @$h; my @shits = sort {$a->pvalue <=> $b->pvalue} @hits; my $n = scalar @shits; if ($n > 0) { $shits[$n-1]->pvalueadj($shits[$n-1]->pvalue); for (my $i=$n-1;$i>=1;$i--) { $shits[$i-1]->pvalueadj( min3($shits[$i]->pvalueadj, $n/$i*$shits[$i-1]->pvalue, 1) ); } } }

Analysez les ensembles de gènes différentiellement exprimés à 15, 30 et 60 minutes identifiés précédemment à l'aide du site YeastSpec/funSpec dédié aux annotations sur la levure.

Analyse de variance

On va maintenant s'intéresser aux gènes différentiellement exprimés dans au moins une des conditions. Téléchargez les données déjà normalisées et chargez-les dans R. Ces données ont été publiées avec un article dans lequel les auteurs ont montré que la flagelline et EF-Tu induisent un ensemble commun de voies de signalisation et de réponses immunitaires chez Arabidopsis thaliana. Pour cela, ils ont traité une souche sauvage d'arabette avec l'éliciteur EF-Tu et mesuré l'expression des gènes à 0, 30 et 60 minutes après ajout d'EF-Tu. Comme contrôle, ils ont utilisé le mutant fls2-17 qui n'a pas de récepteur à la flagelline et pour lequel le traitement avec EF-Tu ne devrait pas induire de réponse immunitaire.

Créez ensuite deux variables, wt et xt, auxquelles vous affecterez leurs niveaux d'expression respectifs. Ensuite, il ne reste plus qu'à effectuer une comparaison de moyennes. Comme nous disposons de trois conditions expérimentales (0, 30 et 60 minutes), nous allons effectuer une analyse de variances (anova). Définissez une variable classes afin de représenter les différents facteurs de wt et xt ainsi que la fonction suivante que vous passerez à la fonction apply pour l'appliquer à toutes les colonnes de wt et xt :

Combien de gènes différentiellement exprimés obtenez-vous ?

Correction des p-valeurs dans le cadre de tests mutliples

La probabilité de rejeter l'hypothèse nulle alors qu'elle est vraie augmente avec le nombre de tests effectués, c'est-à-dire que nous allons décider qu'un certain nombre de gènes sont différentiellement exprimés alors qu'ils ne le sont pas. Pour corriger les p-valeurs obtenues précédemment, nous allons utiliser la FDR (False Discovery Rate proposée par Benjamini et Hochberg en 95) qui est implémentée dans le package multtest qui fait partie de bioconductor. Charger la librairie avec la commande suivante :

Appliquez ensuite la FDR de la façon suivante :

Où testanova est la liste des p-valeurs obtenues par l'anova effectuée précédemment. Pour obtenir le nombre de gènes significatifs :

Une fois que vous avez la liste des gènes différentiellement exprimés chez le sauvage ainsi que chez le mutant, vous allez caractériser comme précédemment ces ensembles de gènes. Pour cela, il vous faudra convertir les identifiants d'oligonucléotides de la puce (Arabidopsis genome ATH1 121501) en identifiant de gènes. Vous pourrez utiliser l'outil de conversion en ligne proposé par le projet TAIR. Une fois les identifiants de gènes obtenus, utilisez l'outil FatiGO de la suite Babelomics pour analyser vos ensembles de gènes.

Partie 2 : Analyse d'un compendium de données de transcriptome

Dans cette partie, vous allez utiliser des données de transcriptome de levure obtenues pour environ 300 mutants (approche compendium). L'idée est d'effectuer un clustering hiérarchique des profils d'expression des gènes dans toutes ces conditions expérimentales afin d'identifier des ensembles de gènes co-exprimés et donc potentiellement co-régulés et qui pourraient donc participer à un même processus biologique. Une fois les gènes regroupés en cluster, le but est de caractériser ces ensembles de gènes par des méthodes de surreprésentation statistique pour identifier dans quel(s) processus biologique(s) ils sont impliqués. L'analyse plus fine d'un cluster de gènes pourra permettre par la suite d'attribuer un processus biologique potentiel pour un gène de fonction inconnue dont le profil d'expression est fortement corrélé avec ceux d'un groupe de gènes bien connus.

Les données que vous allez utiliser ont été rassemblées et analysées et ont fait l'objet d'une publication dont l'identifiant PubMed est 10929718. Retrouvez et lisez l'abstract sur PubMed. Téléchargez les données sur les ratios d'expression pour les 300 conditions et décompressez-les dans votre répertoire de travail. Importez-les ensuite dans MeV et effectuez un clustering hiérarchique en utilisant le coefficient de corrélation de Pearson. Remarque : vous pouvez ajuster la taille d'affichage des ratios dans le menu Display -> Set element size.

Avant de réellement analyser ces données, il est nécessaire d'appliquer un certain filtrage. On ne va retenir que les gènes au moins 2x plus ou 2x moins exprimés dans au moins 5 conditions ainsi que les conditions pour lesquelles on a au moins 5 gènes 2x plus ou 2x moins exprimés. Les données sont au format log10 des ratios d'expression. Dans un premier temps, réalisez un script Perl pour convertir ces données au format log2 des ratios d'expression. Ensuite, écrivez un second script qui filtrera ces données pour ne retenir que les gènes 2x plus surexprimés ou sousexprimés dans au moins 5 conditions, et les conditions qui contiennent au moins 5 gènes 2x plus surexprimés ou sousexprimés. Ensuite, importez ces données dans MeV et refaites le clustering hiérarchique.

A présent, sélectionnez, enregistrez et analysez les clusters qui vous semblent pertinents (au moins 5). Pour l'analyse, utilisez l'outil YeastSpec/funSpec en utilisant les sources de données suivantes : MIPS Functional Classification, MIPS Phenotypes, MIPS Subcellular Localization, MIPS Protein Complexes, et GO Biological Process. Interprétez les résultats obtenus.