Une étude comparée de listes de gènes différentiels générées par séquençage haut-débit [454 Roche] et micro-réseaux d’AD

1. Statistiques pour la Biologie Intégrative [groupe SIB] Une étude comparée de listes de gènes différentiels générées par séquençage haut-débit [454Roche] et micro-réseaux d’ADN [Agilent] G. Guernec1, J. Montfort1, A. Lecam1, R. Reinhart2, P. Prunet1, Y. Guiguen1,P.Y Rescan1, D. Power3 1 INRA UR1037 SCRIBE, Campus de Beaulieu, 35000 - Rennes, FRANCE 2 Max Planck Institute for Molecular Genetics, htpt group, Ihnestr. 63, 14195 - Berlin, ALLEMAGNE 3Universidade do Algarve, Campus de Gambelas, 8005-139, Faro, PORTUGAL V2: Rennes – Groupe SIB – le 19 avril 2011 V1: Versailles – Rencontres Bio-informaticiens et Statisticiens de l’INRA- le 25 mars 2011

2. Octobre 2009 : Acquisition du premier jeu de données de comptage généré par technologie de séquençage haut débit … à l’INRA Scribe … Travaux inclus dans le projet Aquagenome 2008 • Questionnement légitime des biologistes : • «Qu’en est-il de la fiabilité des résultats obtenus à partir de ces nouveaux supports émergents? » • … Report à la bibliographie : • Marioni & Al, 2008 [1] / Comparaison [Illumina] vs [Affymetrix] • Qu’en est-il avec la technologie 454 [Roche] ? • … Vise à répondre à une question type de transcriptomique : • « Recherche des gènes actifs dans le muscle à un stade de développement donné chez la truite arc en ciel » • Analyse différentielle sur les transcrits d’ARN • - 2 stades de développement : Larvaire [J] vs Adulte [B] Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

3. Détails du protocole expérimental • Bonne gestion et prise en charge des données issues de microréseaux d’ADN par le plateau transcriptomique du Scribe. • Mise en place d’une étude comparative entre une technologie microréseaux [Agilent] et une technologie de séquençage haut-débit [454 Roche Titanium] Connaissances a priori • Analyses préalables en microarrays et biblio. : 2 stades très marqués chez le muscle • Planification expérimentale en microarrays : 5 répétitions par condition [Agilentmonocouleur] • … fixe le nombre de répétitions (nombre de chemins) pour le RNA seq (454 -peu d’informations) Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

4. Les grandes étapes d’une analyse de séquençage haut- débit … • …GTGCCCGATCTCGATCAGGTAATCGTATATTAAAGTGCCCGATCTCGATCAGGT… • Départ : Constitution d’une banque d’ADNc pour une espèce donnée en 4 étapes • 1) Fragmenter l’ADNc d’un génome en plusieurs milliers (milliards) de petits morceaux • 2) Lire la séquence de chaque petit fragment (read) (Détection par fluorescence : Sanger) • PB : Les séquenceurs automatiques ne donnent pas tjs des lectures exactes des séquences • Le taux d’erreur dépend: • De la pureté et de l’homogénéïté des échantillons d’ADN initiaux • Séquence spécifique des paires de bases dans l’échantillon d’ADN • Solution : Multiples lectures indépendantes (10*) pour chaque pb du génôme Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

5. Les grandes étapes d’une analyse de séquençage haut- débit … 3) Faire coïncider par informatique, les fragments où les séquences identiques : Alignement (Mapping)des séquences (Tophat1.2…) Zoom sur l’alignement de 2 séquences : Cas de figure possible AGTGCCCGATTATA [0 mismatch] AGTACCCGATTATA [1 mismatch] AGTACCCGACTATA [2 mismatchs] AGTGCCCG-TTATA TATATTACAGTGCCCGATTATATTAAGTG … 0, 1 ou 2 mismatches tolérés 4) Assemblage des lectures (cufflinks…) ayant des séquences chevauchantes: Théorie du contigage - Compare chaque séquence avec chacune des autres séquences - Génère une liste des régions ayant certains critères de similarités de séquences Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

6. Les grandes étapes d’une analyse de séquençage haut- débit … Zoom sur la théorie du contigage Couverture 543210 Contig ou consensus READS Remarque :«L’obtention d’une séquence consensus par assemblage des contigs est ici une préoccupation de second ordre » Exemple de couverture d’un contig … Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

7. Les grandes étapes d’une analyse de séquençage haut- débit … Les principaux supports Définitions : 1) « Un run (réalisation d’un processus complet par la machine) produit un grand nombre de lectures (reads) correspondant à des séquences d’ADN ou d’ARN de l’espèce étudiée. La capacité de la machine se mesure en nombre total de bases séquencées » 2) Library size : nombre total de reads alignés pour un échantillon donné Problèmes associés à l’assemblage : Faible couverture Erreurs de séquençage Séquences répétées Michael Metker. Nature ReviewsGenetics, 2010 [3] Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

8. Détails du protocole expérimental Muscle pool – StadeLarvaire (J) Muscle – StadeAdulte (B) Technologie 454 Roche Technologie micro-réseaux Agilent • Séquençage en simultané de millions de petits fragments d’ADN • (330 pb de moyenne) • Désigné par Salem & Al, 2008 [4*44K] • Oligonucléotides de 60 mers obtenus à partir de processus de synthèse • 45220 sondes incluant 37394 oligos uniques (+6409 oligos dupliqués aléatoirement) • Technique de F. Sanger utilisée en parallèle • Chaque échantillon séquencé sur 1 des 10 pistes • 55793 contigs disponibles Plate-forme GENOTOUL, Toulouse Plateau transcriptomique Scribe, Rennes Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

9. Le jeu de données de comptage généré [454] • Aucun contrôle possible des variabilités engendrées par les étapes d’alignement et d’assemblage Comparaison directe des données d’expression (agilent) et de comptage (454) générés Stade J (*5) Stade B (*5) Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

10. Descriptif des données d’expression de micro-réseaux Départ : 37454 sondes différentes (oligos + contrôles) Boxplot des log(intensités) normalisées • Etape de filtrages classiques (Agilent) en 3 parties : • - uniformité du signal • - intensité du spot > sign. intensité BG • - intensité < seuil de saturation • + Règle de décision / sonde • Suppression des contrôles Nouveau : 26496 sondes conservées (72%) • Normalisation par la médiane des 10 échantillons Gamme dynamique plus importante au stade larvaire Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

11. Descriptif données de comptage [séquençage 454] Départ : 55793 contigs Nombre total de lectures : 1.069.535 reads (compte-tenue de la qté d’ARN de départ) Remarques : « Près de 35% des contigs comptent 0 lecture à chacune des 10 pistes étudiées et ne sont donc pas conservés pour la suite de l’étude » 24592 contigs (>67%) comptent 0 lecture à un des 2 stades [contigs spécifiques] Problème lors de l’estimation des Fold Changes : FC = Jest/ Best  0 ? … ou inversement …. Nouveau : 36411 contigs Evaluation de la variabilité globale et intra - conditions • Elaboration d’un score /20 d’estimation de la qualité des répétitions biologiques (ie inter-runs) • 20 tests exact de Fisher (loi hypergéométrique) avec une erreur de 1% (+ ajustement Bonferroni) • - 5 runs « stade larvaire » = 10 tests « runs 2 à 2 » • - 5 runs « stade adulte » = 10 tests « runs 2 à 2 » Pour un gène … Plus le score est faible plus la répétabilité inter-runs est bonne Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

12. Evaluation de la variabilité globale et intra – conditions Départ : 36411 contigs Nombre total de lectures : 1.069.535 reads (compte-tenue de la qté d’ARN de départ) Remarques : « 92% des contigs testés ont une répétabilité globale bonne voire excellente » On observe une répétabilité inter-runs plus homogène au stade adulte qu’au stade larvaire Mise en place de l’étude comparative … … A partir de 2 fichiers de même dimension Recherche de la correspondance oligos<-> contigs à partir d’un fichier d’annotations swissprot Suppression des oligosagilent pour lesquels il n’y a pas de correspondance Conservation des oligos pour laquelle la correspondance est unique Final : Etude comparative à partir de 2 fichiers de 7010 gènes / identifiant unique Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

13. Normalisation des données de comptage (N=7010) • Etude comparative proposée par [Marioni & Al. ,2008] : PAS de normalisation des données de comptage • Depuis : « Importance de la normalisation sur le résultat final » [Bullard et al, 2010] [2] • Choix de la technique de normalisation ? ….. Encore flou ……. Démarche proposée : Ajout d’un 1 à l’ensemble des reads (cellules du fichier) : Pour s’affranchir du pb des 0s (cf estimation des Fold Changes) Ajustement par le total des reads sur le run : Xij(k) : Nbr de reads gène i / échantillon (run) j / condition k Nj(k) : Production totale reads sur run j (library size) / condition k 3) Choix de la constante : Cste = max(Nj(k)) …. Et arrondis …. Pourquoi ? Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

14. Effets de la normalisation des données de comptage (N=7010) Remarques : Variabilité inter – runs reste plus importante au stade J avec N = 7010 Amélioration significative de la répétabilité après normalisation (tests de Mac Némar 5%) … Semble contradictoire avec les résultats de microarrays : Répétabilité entre les échantillons globalement meilleure au stade J Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

15. Descriptifs et comparaisons (N=7010) 4) Gamme dynamique plus importante au stade J (idem microarrays) 5) La corrélation intra-groupe reste globalement très forte en microarrays *: unequalssubgroupsbecause of ex-aecquos 6) Elle est d’autant plus forte que les gènes observés comptent beaucoup de lectures Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

16. Analyse différentielle – Démarche proposée Les analyses se font gène par gène puis gestion des tests multiples par correction FDR … Pour les micro-réseaux tests limma (G. K. Smyth – BioConductor) … Pour le séquençage 454 • Régression de Poisson [contrainte forte : E(Y)=VAR(Y)=λ] : « Le nombre de lectures observées pour une piste et un contig donné est représenté par la variable Y de paramètre λ = (λJ; λB) telle que : Avec : et i=1, .., 5 • Estimation des Fold Changes (max. de vraisemblance) : • Le problème de la surdispersion: …Biais d’optimisme possible sur les p-values Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

17. Analyse différentielle – Démarche proposée … Pour le séquençage 454 • Détection de la surdispersion : test de chi2 de Pearson • 256 modèles / 7010 en surdispersion • Recours à la régression négative binômiale : Inclusion d’un terme de bruit … • Test du rapport de vraisemblance pour déterminer la significativité de la variable « stade » • « Sous H0: La statistique de test suit une loi de chi2 à 1 ddl » Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

18. N= 7010 - 1% FDR Analyse différentielle – Résultats Comment expliquer les divergences des listes de gènes obtenues à partir de ces 2 supports ? Remarques : 1) Proportion importante de gènes différentiels en microarrays (plus de 50%) 2) … spécifiques aux données de microarrays. Pourquoi ? Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

19. Analyse différentielle – Une explication des divergences … • En supposant que la répétabilité inter –runs est parfaite : Sommation des reads intra-conditions Exemple : gène X -> 0 0 1 1 1 (reads / run : stade J) 4 5 4 5 5 (stade B) J B Tableau de contingence 2*2 J B Gène X 3 23 Gène X 0 0 1 1 1 4 5 4 5 5 ∑ ∑ TOTAL 546812 528436 TOTAL • Cadre classique : test exact de Fisher 1% (loi hypergéométrique) • H0 : Prop. Estimée (stadeJ) = Prop. Population (stade J) Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

20. Analyse différentielle – Une explication des divergences … • Intervalles de confiance (99%) des proportions relatives à toutes les combinaisons possibles pour chacun des 7010 gènes • Exemple ci-dessous : 2 gènes comptant 95 (à gauche) et 4 lectures (à droite) au total 99% CI Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

21. Visualisation ordonnée des graphes « On ne peut conclure qu’un gène est différentiel dès lors qu’il est associé à moins de 5 lectures au total » -> Manque de puissance Analyse différentielle – Reprise des diagrammes • On constate une baisse très nette du nombre de gènes différentiels « spécifiques » microarrays : Près de 3000 gènes différentiels [Agilent] n’ont pas bénéficié d’une couverture suffisante pour sortir différentiels en RNAseq … • Quelle interprétation pour les 1460 gènes différentiels restants ? Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

22. Analyse différentielle – Une interprétation des divergences … Nouveau zoom sur la répétabilité inter – runs : Remarques : Les gènes différentiels en microréseaux spécifiquement ont une variabilité inter-runs significativement plus importante que dans les autres groupes Evocation d’un problème de l’adéquation du modèle aux données … Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

23. Analyse différentielle – Comparaison des Fold Changes Remarques : 1) On retrouve la caractéristique principale de la comparaison Affymetrix <-> Agilent « The set of differentially expressed genes that show the strongest correlation (…) seems to be those that are mapped by many reads » [ J.C Marioni (Genome Research -2008) ] 2) Proportion conséquente de fold changes atypiques (changement de signes) Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

24. Analyse différentielle – Interprétation des autres groupes spécifiques Proportion des Fold Changes atypiques en fonction su seuil Proportion des gènes spécifiquement différentiels en séquençage en fonction du seuil Remarques : La proportion des fold changes « atypiques » diminue avec le seuil de test Les p-values [microarray] des gènes spécifiques « RNA seq» restent en bonne partie voisine du seuil choisi Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

25. Analyse différentielle – Influence de la taille des contigs ? • Répartition de l’ensemble des contigs de départ en fonction de leurs tailles (nombre de bases de nucléotides) • Répartition des effectifs des groupes de gènes différentiels (%) par taille de transcrits Remarque(s) : « On observe une légère tendance mais pas d’effets notables de la taille des transcrits sur l’expression des contigs » Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

26. Analyse différentielle – Influence des répétitions biologiques en 454 ? • Remarque 1: « Dérive actuelle » - Ne séquencer qu’un échantillon d’ARN par condition puis tests exacts de Fisher et corrections FDR • Quel est le pouvoir de généralisation de résultats ainsi obtenus ? • Remarque 2: A- t’on besoin d’autant de répétitions biologiques en 454 qu’en Agilent pour obtenir des listes de gènes différentielles quasi - similaires ? Principe Approcher la liste de gènes différentiels obtenue par micro-réseau en faisant varier le nombre de réplicats biologiques (runs) du 454 Venn ’s group II [N= 2681] Remarque(s) : « Amélioration notable de l’intersection et des FC » Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

27. Analyse différentielle – Comparaison des normalisations proposées pour le 454 Panorama des normalisations existantes (04/2011) • Total read Count normalization (TC) • Hypothèse : Comptages proportionnels au niveau d’expression et à la profondeur de séquençage • Plus la profondeur de la librairie est importante, plus le comptage associé au transcrit est grand (indépendamment de la taille du transcrit et du niveau d’expression) • Full Quantile normalization (Bolstad & Al, 2002) (FQ) [4] • Equivalent « normalisation par quantile » en microréseaux • Implémentation limma 3.6.9 [5] • Hypothèse : Distribution identique des comptages d’une librairie à l’autre • DESEQ normalization (Anders & Huber, 2010) (DESEQ) [6] • Implémentation package DESEQ • Xij~NB(μij,σ2ij) • μij= qij(k)*sj • …avec qij(k) : Concentration « vraie » gène i /échantillon j • Adjusted Total read Count normalization (MAX) • Hypothèse : Augmentation proportionnel du nombre de comptages par gène pour suggérer une couverture suffisante par gène a • Scalingnormalization (Robinson & Al, 2010) (TMM) [7] • Principe : • fj= si / sj • Mrij : log ratio gène i entre échantillon j et r • wij: Correction sur la variance • Li : Taille du gène i • Single Quantile Normalization (Bullard & al, 2010) (MED, UQ) [2] • Hypothèse : Nombre total de reads dépend surtout de qqs contigs fortement exprimés Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

28. Analyse différentielle – Comparaison des normalisations proposées pour le 454 [N= 2681] Résultats : Tableau (ci-dessus): Pourcentage de gènes différentiels communs en faisant varier le type de normalisation pour les données de 454 2) Graphe (ci-contre) : Similarités des listes de gènes différentiels obtenues à partir des données de 454 en faisant varier le type de normalisation Clustering hiérarchique Distance de corrélation et méthode de Ward Conclusion : Diagrammes de Venn très similaires d’une normalisation à l’autre Particularité de la normalisation FQ Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

29. Conclusions & Perspectives [Technologie 454] Etude essentiellement basée sur la comparaison de diagrammes de Venn MAIS mesurait-on vraiment la même chose au départ ? 2) L’intérêt de la technologie 454 est liée à la garantie pour l’utilisateur d’une couverture suffisante a priori sur l’ensemble des contigs étudiés : « Solutions bioinformatiques et/ ou planifications expérimentales adaptées passent par une coordination soutenue des compétences bio-informatiques et statistiques » 3) Nécessité de « garder le contact » avec les biologistes. Importance du travail de vulgarisation scientifique par la mise à disposition d’outils intermédiaires adaptés : - Score d’estimation de la variabilité des répétitions biologiques et/ou techniques - Outils d’estimation de seuil de lectures « limite » 4) Les modèles d’analyse différentielle utilisés en séquençage 454 reposent encore aujourd’hui sur des hypothèses fortes (loi hypergéométrique), et nécessitent d’être améliorés … Après normalisation des données de comptage, certains gènes différentiels en microarrays ne le sont plus en séquençage du fait d’une mauvaise prise en compte de la variabilité qui perdure … L’adéquation du choix de l’association « normalisation + modèle » se doit d’être rigoureusement déterminée Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

30. Références [liste non exhaustive] • Metzker ML. Sequencing technologies –the next generation. Nature Reviews Genetics, 2010 Jan; 11(1):31-46 [3] • • Bullard JH, Purdom E, Hansen KD, Dudoit S. Evaluation of statisticalmethods for normalization and differential expression in mRNA-seqexperiments. BMC Bioinformatics2010, 11:94. [2] • • Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mappingandquantifyingmammaliantranscriptomes by RNA-seq. Nature Methods, 2008 Jul; 5(7):621-628. • • Audic, S. and Claverie, J-M. The significance of digital gene expression profiles. Genome research 1997; 7:986-995 • • Cloonan, N et al. Stem celltranscriptomeprofiling via massive-scalemRNAsequencing. Nature Methods2008; 5(7):613-619. • • Robinson MD, Smyth, GK. Moderatedstatistical tests for assessingdifferences in tag abundance. Bioinformatics23(21):2881-2887. [7] • • Robinson MD, McCarthy DJ, Smyth, GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics2009. • • Marioni JC, Mason CE et al. RNA-seq: An assessmentoftechnicalreproducibility and comparisonwithgene expression arrays. GenomeResearch2008, 18:1509-1517. [1] • • Robinson MD, Oshlack A. A scalingnormalizationmethod for differential expression analysis of RNA-seqdata. GenomeBiology2010, 11:R25 • • Anders, S and Huber, W. Differential expression analysis for sequence count data. Nature Precedings 2010, march [6] • Torres T, Metta M, Ottenwälder B, and Schoötterer C. 2008 : Gene expression profiling by massively parallel sequencing, Genome Research, 1: 172-177 Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011

31. Références [suite] • Bolstad B, Irizarry R, Astrand M, and Speed T. A comparison of normalization methods for high density nucleotide array data based on variance and bias. Technical report, UC Berkeley, 2002 [4] • Gordon K. Smyth. Limma : linearmodels for microarray data. In R. Gentleman, V. Carey, S. Dudoit, R. Irizarry, and W. Huber, editors, Bioinformatics and Computational Biology Solutions using R and Bioconductor, pages 397, 420. Springer, New York, 2005 [5] Statistiques pour la Biologie Intégrative [Groupe SIB] – Rennes, le 19 avril 2011