Le schéma conceptuel de GREIN est montré dans la Fig. 1. Les ensembles de données RNA-seq individuels sont traités par le pipeline GREP2 et stockés localement en tant que R Expression Sets. L’utilisateur peut accéder aux ensembles de données prétraitées et les analyser via l’interface utilisateur graphique (GUI) de GREIN ou soumettre pour traitement des ensembles de données qui n’ont pas encore été traités. Les flux de travail pilotés par l’interface graphique facilitent l’examen et la visualisation des données, l’analyse statistique, la construction de signatures transcriptionnelles et l’interprétation en biologie systémique des gènes différentiellement exprimés (DE). GREIN et le pipeline dorsal (GREP2) sont tous deux écrits en R et publiés respectivement sous la forme d’un conteneur Docker et d’un paquetage R. Les interfaces utilisateur graphiques de GREIN sont implémentées dans Shiny16, un cadre web permettant de créer des applications web dynamiques dans R. L’instance web à https://shiny.ilincs.org/grein est déployée via un essaim Docker robuste de serveurs Shiny à charge équilibrée. L’infrastructure complète de GREIN, y compris le pipeline de traitement, est déployée via des conteneurs Docker.
Flux schématique de GREP2, interface web et sorties de GREIN. Les jeux de données GEO sont systématiquement traités à l’aide du pipeline GREP2 et stockés dans la bibliothèque de jeux de données back-end. Les flux de travail GREIN pilotés par interface graphique facilitent l’analyse complète et la visualisation des ensembles de données traités.
Les flux de travail conviviaux pilotés par interface graphique dans GREIN facilitent les scénarios de réutilisation typiques des données RNA-seq, tels que l’examen des mesures de contrôle de la qualité et la visualisation des modèles d’expression dans l’ensemble de l’ensemble de données, l’analyse de la taille de l’échantillon et de la puissance dans le but d’informer la conception expérimentale des études futures, l’expression génétique différentielle statistique, l’enrichissement de la liste des gènes et l’analyse de réseau. Outre la comparaison standard entre deux groupes, le module d’analyse de l’expression différentielle des gènes prend également en charge l’ajustement d’un modèle linéaire généralisé qui tient compte des covariables ou des effets de lot. Les outils de visualisation et d’exploration interactifs mis en œuvre comprennent l’analyse en grappes, les cartes thermiques interactives, l’analyse en composantes principales (ACP), l’intégration stochastique des voisins distribuée en t (t-SNE), etc. (tableau supplémentaire S1). L’utilisateur peut également rechercher les annotations ontologiques des échantillons d’ARN-seq humain et des ensembles de données fournis par le projet MetaSRA13. Chaque échantillon d’ARN-seq humain traité est étiqueté avec la cartographie MetaSRA des ontologies biomédicales, y compris l’ontologie des maladies, l’ontologie des cellules, l’ontologie des facteurs expérimentaux, Cellosaurus et Uberon. L’interprétation biologique des expressions génétiques différentielles est facilitée par des liens directs vers d’autres outils en ligne permettant d’effectuer des analyses post-hoc typiques telles que l’analyse d’enrichissement des listes de gènes et des voies et l’analyse de réseau des gènes différentiellement exprimés (DE). La connexion à ces services Web analytiques est mise en œuvre en soumettant la signature d’expression génique différentielle (c’est-à-dire la liste des changements moyens dans l’expression génique et les valeurs p associées pour tous les gènes analysés régulés à la hausse ou à la baisse) à iLINCS17 (Integrative LINCS). iLINCS fournit également l’analyse de connectivité des signatures pour les signatures Connectivity Map L1000 récemment publiées18. Des instructions détaillées, étape par étape, sur les flux de travail d’analyse GREIN sont fournies dans le matériel supplémentaire et dans la section » Aide » de GREIN.
Fonctionnalités clés
Recherche ou soumission pour traitement
L’utilisateur peut soit rechercher un ensemble de données GEO déjà traité dans la case » Recherche de séries GEO (GSE) accession « , soit soumettre un ensemble de données pour traitement si l’ensemble de données n’est pas déjà traité (figure supplémentaire S2). À l’heure actuelle, la grande majorité des ensembles de données RNA-seq GEO d’humains, de souris et de rats ont été prétraités et la soumission par l’utilisateur d’ensembles de données GEO pour traitement ne sera requise qu’occasionnellement. L’utilisateur peut vérifier l’état du traitement de l’ensemble de données demandé dans l’onglet « Processing console » (figure supplémentaire S3). D’autres options de recherche comprennent la recherche par mots-clés à travers les métadonnées des ensembles de données et la recherche d’échantillons à travers les ontologies biomédicales via les annotations ontologiques MetaSRA.
Explorer les ensembles de données
GREIN permet d’accéder aux données brutes et normalisées (comptes par million et transcription par million) au niveau des gènes et des transcriptions. GREIN est livré avec plusieurs outils interactifs et personnalisables pour visualiser les profils d’expression tels que des cartes thermiques interactives de gènes et d’échantillons groupés, des graphiques de densité pour tous les échantillons ou un sous-ensemble d’échantillons, une analyse de la variabilité entre et au sein d’un groupe par le biais d’analyses de réduction de la dimensionnalité en 2D et 3D et de visualisations telles que PCA et t-SNE (Fig. 2). L’utilisateur peut également visualiser le profil d’expression de chaque gène séparément (figure supplémentaire S6).
Plots d’analyse exploratoire dans GREIN. (A) La carte thermique de corrélation montre une corrélation plus élevée au sein des lignées cellulaires et une faible corrélation entre les lignées cellulaires. En général, les corrélations élevées au sein de chaque lignée cellulaire indiquent une qualité élevée des profils transcriptionnels. (B) Regroupement hiérarchique basé sur la corrélation de Pearson des 500 gènes les plus variables en fonction de l’écart absolu médian comme mesure de variabilité. Les données sont normalisées et centrées sur la moyenne. (C) Graphique tridimensionnel de l’analyse en composantes principales des lignées cellulaires. (D) Tracé bidimensionnel t-SNE de la condition de traitement et de la lignée cellulaire montre une séparation claire des lignées cellulaires, puis des fractions d’ARN indiquant deux sources dominantes de la variabilité entre les profils RNA-seq.
Contrôle de qualité
La qualité des données RNA-seq dans les dépôts publics continue d’être un problème majeur. Dans une étude récente de Deelen et al.19, plus de la moitié des 65 000 échantillons RNA-seq publics traités ont dû être retirés en raison de problèmes de CQ. Plutôt que de retirer des échantillons, GREIN fournit un rapport de contrôle de qualité (CQ) complet des données de séquence brutes et de la cartographie de séquence pour chaque échantillon (figure supplémentaire S7), et permet à l’utilisateur de prendre une décision sur les échantillons à exclure des analyses en aval.
Analyse de puissance statistique
Le module d’analyse de puissance de GREIN facilite le calcul et la visualisation de la puissance statistique de détection des gènes exprimés de manière différentielle dans les études futures utilisant des échantillons biologiques similaires. L’estimation de la taille d’échantillon appropriée pour les études futures avec des échantillons biologiques similaires est souvent le facteur de motivation clé dans la réanalyse des données RNA-seq. L’analyse de puissance facilite également l’analyse post-hoc des taux de faux négatifs dans l’ensemble de données actuel. Le manque de puissance statistique et les différences de puissance statistique entre les gènes peuvent produire des résultats faux négatifs conduisant à des conclusions erronées20. Le segment « Courbe de puissance » fournit des estimations de puissance pour différents nombres d’échantillons basés sur un seul gène (Fig. 3A). L’utilisateur peut modifier les valeurs par défaut des paramètres. Le graphique « Détectabilité des gènes » visualise l’estimation de la puissance de chacun des gènes en fonction des groupes sélectionnés et de la dispersion par gène (Fig. 3B). La couverture moyenne des gènes est tracée par rapport à leur variabilité biologique et est affichée en deux ensembles basés sur leur statut de détectabilité (puissance ≥0,8 et puissance < 0,8).
Analyse de puissance pour évaluer les changements transcriptionnels dans la lignée cellulaire MCF10A non maligne. (A) Estimations de puissance basées sur un seul gène pour différents nombres d’échantillons dans chaque groupe avec un changement de pli minimal de 2 et une signification statistique α = 0,01. (B) Détectabilité par gène sur le plan log2CPM-BCOV avec FDR ≤0,1 et deux échantillons dans chaque groupe.
Expression génique différentielle
Créer et interpréter une signature d’expression génique différentielle est un scénario d’analyse typique dans les expériences RNA-seq. Avec GREIN, l’utilisateur peut créer une signature en comparant l’expression génique entre deux groupes d’échantillons avec ou sans ajustements pour les covariables expérimentales ou les effets de lot. GREIN peut gérer des plans expérimentaux complexes en offrant la flexibilité de réarranger les groupes et sous-groupes ou de sélectionner des échantillons spécifiques. La signature d’expression différentielle peut être visualisée via des graphiques interactifs qui incluent une carte thermique des gènes les plus déréglés de manière différentielle (figure supplémentaire S15) classés par taux de fausse découverte (FDR), un graphique de changement de pli logarithmique par rapport à l’expression moyenne logarithmique (MA) (figure supplémentaire S16) et un graphique de détectabilité des gènes (figure supplémentaire S17). La signature d’expression différentielle, avec ou sans prise en compte des résultats potentiellement faux négatifs, peut être directement exportée vers iLINCS pour une analyse d’enrichissement et de connectivité.
Cas d’utilisation : Analyse de la régulation transcriptionnelle et translationnelle de l’hypoxie dans des lignées cellulaires de cancer du sein épithélial non malin et triple négatif
Nous démontrons l’utilisation de GREIN en réanalysant une donnée GEO RNA-seq récemment publiée (GSE104193). Sesé et al.21 ont examiné la régulation transcriptionnelle et translationnelle du sous-type de cancer du sein triple négatif (TNBC) hormono-réfractaire sous une combinaison d’hypoxie et de traitement par inhibiteur mTOR (cible mécaniste de la rapamycine). Les auteurs ont notamment analysé les profils d’expression de cellules TNBC (MDA-MB-231) et de cellules épithéliales mammaires non malignes (MCF10A) exposées à des conditions normoxiques (21% O2) et hypoxiques (0,5% O2) et/ou traitées avec un inhibiteur de mTORC1 et -2, le PP242. Chacun des échantillons a été séquencé pour les ARNm totaux (T) et liés aux polysomes (P). L’ensemble de données contient 32 échantillons, représentant deux réplicats biologiques pour chaque combinaison de lignée cellulaire, de niveau d’oxygène, de statut de traitement et de fraction d’ARNm.
L’analyse exploratoire de l’ensemble de données traitées dans GREIN (Fig. 2) montre que la plus forte source de variation entre les échantillons provient des différences entre les deux lignées cellulaires. Ceci est confirmé par l’analyse de corrélation des profils d’expression complets (Fig. 2A), le regroupement hiérarchique des 500 gènes les plus variables sur la base de l’écart absolu médian (Fig. 2B), le graphique PCA 3D des échantillons (Fig. 2C) et le graphique t-SNE 2D (Fig. 2D). En outre, les corrélations élevées entre les profils d’expression pour la même lignée cellulaire (Fig. 2A) indiquent un bon rapport signal/bruit dans les mesures d’expression génique. La sous-structure supplémentaire des données indiquée par le graphique 2D t-SNE a été examinée en peignant les échantillons en fonction de différents attributs (figure supplémentaire S1). Cette analyse a révélé que les séparations au sein de chaque lignée cellulaire sont induites par différentes fractions d’ARNm, puis par des différences entre les conditions expérimentales.
Puis, nous avons utilisé GREIN pour effectuer une analyse de puissance statistique basée sur le modèle de variabilité biologique observé dans cet ensemble de données. Nous avons considéré les profils transcriptionnels de chaque lignée cellulaire exposée à l’hypoxie et traitée avec ou sans PP242, ce qui conduit à quatre comparaisons. En supposant une différence d’expression d’au moins deux fois entre les groupes, à la signification statistique de α = 0,01, et avec seulement deux répétitions dans chaque groupe, la puissance statistique d’un gène à être détecté comme différentiellement exprimé est inférieure à 0,55 dans toutes les comparaisons (tableau 2). Notre analyse indique qu’il faudrait quatre répliques par groupe pour atteindre une puissance de 80 % en détectant un changement de deux fois dans l’expression (tableau 2 et figure 3A). Dans une expérience typique de RNA-seq, une profondeur de séquençage de 20 à 30 millions est suffisante pour quantifier l’expression des gènes pour presque tous les gènes4,22, ce qui est également évident dans cet ensemble de données. Nous avons également évalué le pouvoir statistique de chaque gène d’être détecté comme différentiellement exprimé à partir du graphique « Détectabilité des gènes ». Le logarithme moyen des valeurs de comptage par million (CPM) des gènes a été représenté par rapport au coefficient de variation biologique (BCOV) par gène et la puissance a été calculée pour les gènes correspondants (Fig. 3B). Un taux de fausses découvertes contrôlé de 0,05 et un pourcentage attendu de vrais positifs de 10 % ont été utilisés pour estimer la signification statistique. Nous définissons un gène comme étant détectable comme étant différentiellement exprimé en condition hypoxique si sa puissance est de 0,8 ou plus. Comme prévu, il existe une relation inverse entre le BCOV et la puissance (Fig. 3B). De plus, la puissance pour détecter l’expression différentielle d’un gène augmente avec un log CPM ou une taille d’effet plus élevée.
L’un des objectifs de l’étude était d’analyser les changements transcriptionnels dans des conditions hypoxiques et normoxiques avec et sans traitement par PP242 dans les lignées cellulaires MCF10A et MDA-MB-231. Nous avons créé des signatures transcriptionnelles de l’hypoxie et de l’hypoxie + PP242 dans l’ARNm total par une analyse d’expression différentielle entre les échantillons d’hypoxie et d’hypoxie + PP242 respectivement par rapport aux échantillons de contrôle tout en ajustant l’effet de lot en traitant la » réplique » comme une covariable, pour chaque lignée cellulaire séparément. Nous avons trouvé un nombre plus élevé de gènes exprimés de manière différentielle (DE) dans les lignées cellulaires MCF10A par rapport à MDA-MB-231 à la fois en hypoxie et en hypoxie + PP242 (Fig. 4A), ce qui indique que la lignée cellulaire tumorale est peut-être mieux équipée pour faire face à l’hypoxie. Cette analyse a également montré que la plupart des gènes exprimés de manière non différentielle ne sont pas non plus détectables, ce qui indique qu’ils peuvent représenter des résultats faussement négatifs. Ceci est conforme à l’analyse de puissance montrant que 4 échantillons par groupe seraient nécessaires pour identifier de manière cohérente les gènes exprimés de manière différentielle avec une BCOV moyenne. Pour identifier les gènes moins exprimés, une taille d’échantillon encore plus élevée serait nécessaire.
Expression différentielle et détectabilité des gènes. (A) Le nombre de gènes (échelle log10) non différentiellement exprimés et non détectables (NDE&NDT), différentiellement exprimés (DE), et non différentiellement exprimés mais détectables (NDE&DT) dans les comparaisons avec la normoxie pour la fraction totale d’ARNm. Nous appelons un gène détectable (DT) si sa puissance ≥0,8 et différenciable si FDR < 0,05. (B) Le tracé de détectabilité des gènes pour la première comparaison (MCF10A et hypoxie) qui visualise la liste des gènes susmentionnés ainsi que leurs fold changes (FC) respectifs.
Pour interpréter les gènes différentiellement exprimés en termes de voies biologiques affectées, nous avons soumis les signatures d’expression génique différentielle de l’hypoxie à des outils d’enrichissement en ligne (DAVID23, ToppGene24, Enrichr25 et Reactome26) via iLINCS. Les signatures soumises comprenaient une liste combinée de gènes DE et NDE&DT représentant les vrais positifs et les vrais négatifs probables. Les gènes ont été sélectionnés sur la base d’un seuil de 0,7 et 0,01 pour la puissance statistique et le FDR respectivement. La figure 5 illustre les résultats d’enrichissement obtenus à partir de ToppGene pour la signature d’hypoxie de MCF10. Les gènes significativement enrichis (FDR < 0.05) Les 10 principales catégories d’ontologie génique (GO) de ToppGene et de l’outil d’annotation fonctionnelle DAVID comprennent la réponse à l’hypoxie, la réponse à une diminution des niveaux d’oxygène, l’angiogenèse, la régulation de la prolifération cellulaire, le processus d’oxydoréduction et la réponse à un stimulus abiotique qui sont communs aux deux lignées cellulaires (tableau supplémentaire S2 et tableau supplémentaire S3). La plupart de ces catégories sont cohérentes avec l’étude originale. En outre, la suite ToppGene a identifié le réseau de facteurs de transcription du facteur induit par l’hypoxie (HIF-1-alpha) qui était activé dans les deux lignées cellulaires (tableau supplémentaire S4 et tableau supplémentaire S5).
Snapshot de certaines des catégories significatives de voies et d’ontologie génique (GO) de ToppGene via iLINCS. Ces catégories sont trouvées dans la comparaison entre l’hypoxie et la normoxie dans la lignée cellulaire MCF10A en utilisant une liste combinée de gènes DE et NDE&DT. La ligne verticale rouge représente le seuil sélectionné de 0,05.
Enfin, nous avons utilisé la connexion GREIN avec iLINCS pour » connecter » la signature téléchargée avec les signatures LINCS27 consensus (CGS) knockdown de gènes18. Nous avons trouvé 3 727 signatures LINCS consensus knockdown gène qui ont été significativement (pValue < 0,05) connecté avec notre signature téléchargée. Les gènes cibles des 100 premières signatures connectées ont été sélectionnés pour une analyse d’enrichissement supplémentaire. Nous avons trouvé la réponse cellulaire à l’hypoxie et la régulation du facteur inductible de l’hypoxie (HIF) par l’oxygène dans la liste des 10 principales voies activées dans les deux lignées cellulaires (tableau supplémentaire S6 et tableau supplémentaire S7). Bien que cette analyse produise des catégories fonctionnelles enrichies similaires à celles de l’analyse d’enrichissement initiale, elle complète cette dernière en impliquant plusieurs gènes cibles qui ne sont pas exprimés de manière différentielle bien qu’ils soient suffisamment exprimés pour être détectables selon nos analyses de puissance. L’association de ces deux résultats implique ces gènes comme régulateurs potentiels de plus haut niveau de la réponse à l’hypoxie.