Der konzeptionelle Rahmen von GREIN ist in Abb. 1 dargestellt. Individuelle RNA-seq-Datensätze werden von der GREP2-Pipeline verarbeitet und lokal als R Expression Sets gespeichert. Benutzer können über die grafische Benutzeroberfläche (GUI) von GREIN auf vorverarbeitete Datensätze zugreifen und diese analysieren oder noch nicht verarbeitete Datensätze zur Verarbeitung einreichen. GUI-gesteuerte Arbeitsabläufe erleichtern die Untersuchung und Visualisierung von Daten, die statistische Analyse, die Erstellung von Transkriptionssignaturen und die systembiologische Interpretation von differentiell exprimierten Genen (DE). Sowohl GREIN als auch die Back-End-Pipeline (GREP2) sind in R geschrieben und als Docker-Container bzw. R-Paket veröffentlicht. Die grafischen Benutzeroberflächen für GREIN sind in Shiny16 implementiert, einem Web-Framework zur Erstellung dynamischer Webanwendungen in R. Die Webinstanz auf https://shiny.ilincs.org/grein wird über einen robusten Docker-Schwarm von Shiny-Servern mit Lastausgleich bereitgestellt. Die gesamte GREIN-Infrastruktur, einschließlich der Verarbeitungspipeline, wird über Docker-Container bereitgestellt.
Schematischer Arbeitsablauf von GREP2, Web-Interface und Ausgaben von GREIN. GEO-Datensätze werden systematisch mit der GREP2-Pipeline verarbeitet und in der Back-End-Datensatzbibliothek gespeichert. GUI-gesteuerte GREIN-Workflows erleichtern die umfassende Analyse und Visualisierung der verarbeiteten Datensätze.
Benutzerfreundliche GUI-gesteuerte Workflows in GREIN erleichtern typische Wiederverwendungsszenarien für RNA-seq-Daten, wie z. B. die Untersuchung von Qualitätskontrollmaßnahmen und die Visualisierung von Expressionsmustern im gesamten Datensatz, Stichprobengrößen- und Power-Analysen zur Information über die Versuchsplanung künftiger Studien, statistische differentielle Genexpression, Genlistenanreicherung und Netzwerkanalysen. Neben dem standardmäßigen Zwei-Gruppen-Vergleich unterstützt das Modul zur differenziellen Genexpressionsanalyse auch die Anpassung eines verallgemeinerten linearen Modells, das Kovariaten oder Chargeneffekte berücksichtigt. Zu den implementierten interaktiven Visualisierungs- und Explorationswerkzeugen gehören Clusteranalyse, interaktive Heatmaps, Hauptkomponentenanalyse (PCA), t-distributed stochastic neighbor embedding (t-SNE) usw. (Ergänzende Tabelle S1). Der Benutzer kann auch nach ontologischen Annotationen von menschlichen RNA-seq-Proben und Datensätzen suchen, die vom MetaSRA-Projekt13 bereitgestellt werden. Jede verarbeitete menschliche RNA-seq-Probe ist mit der MetaSRA-Zuordnung biomedizinischer Ontologien gekennzeichnet, darunter Disease Ontology, Cell Ontology, Experimental Factor Ontology, Cellosaurus und Uberon. Die biologische Interpretation der differentiellen Genexpressionen wird durch direkte Links zu anderen Online-Tools unterstützt, mit denen typische Post-hoc-Analysen durchgeführt werden können, wie z. B. die Analyse der Genliste und der Anreicherung von Signalwegen sowie die Netzwerkanalyse differentiell exprimierter Gene (DE-Gene). Die Verbindung zu diesen analytischen Webdiensten wird hergestellt, indem die differenzielle Genexpressionssignatur (d. h. die Liste der durchschnittlichen Veränderungen in der Genexpression und die zugehörigen p-Werte für alle analysierten/aufwärts/abwärts regulierten Gene) an iLINCS17 (Integrative LINCS) übermittelt wird. iLINCS bietet auch die Signaturen-Konnektivitätsanalyse für die kürzlich veröffentlichten Connectivity Map L1000-Signaturen18. Detaillierte Schritt-für-Schritt-Anleitungen zu den GREIN-Analyse-Workflows finden Sie im ergänzenden Material und im Abschnitt „Hilfe“ in GREIN.
Schlüsselfunktionen
Suchen oder zur Verarbeitung einreichen
Der Benutzer kann entweder in der Box „Suche nach GEO-Serien (GSE)-Zugang“ nach einem bereits verarbeiteten GEO-Datensatz suchen oder einen Datensatz zur Verarbeitung einreichen, wenn der Datensatz noch nicht verarbeitet ist (ergänzende Abb. S2). Zum jetzigen Zeitpunkt ist die überwiegende Mehrheit der GEO-RNA-seq-Datensätze für Menschen, Mäuse und Ratten bereits vorverarbeitet und die Einreichung von GEO-Datensätzen zur Verarbeitung durch den Benutzer ist nur noch gelegentlich erforderlich. Der Benutzer kann den Verarbeitungsstatus des angeforderten Datensatzes auf der Registerkarte „Verarbeitungskonsole“ überprüfen (ergänzende Abb. S3). Zu den weiteren Suchoptionen gehören die Suche nach Schlüsselwörtern in den Metadaten der Datensätze und die Suche nach Proben in biomedizinischen Ontologien über die ontologischen MetaSRA-Annotationen.
Explore dataset
GREIN ermöglicht den Zugriff auf rohe und normalisierte (counts per million und transcript per million) Gen- und Transkriptdaten. GREIN verfügt über mehrere interaktive und anpassbare Tools zur Visualisierung von Expressionsmustern, wie z. B. interaktive Heatmaps von geclusterten Genen und Proben, Dichteplots für alle oder eine Teilmenge von Proben, Variabilitätsanalysen zwischen und innerhalb von Gruppen durch 2D- und 3D-Dimensionalitätsreduktionsanalysen und Visualisierungen wie PCA und t-SNE (Abb. 2). Der Benutzer kann auch das Expressionsprofil jedes Gens separat visualisieren (Supplementary Fig. S6).
Explorative Analyseplots in GREIN. (A) Die Korrelations-Heatmap zeigt eine hohe Korrelation innerhalb der Zelllinien und eine geringe Korrelation zwischen den Zelllinien. Generell deuten hohe Korrelationen innerhalb jeder Zelllinie auf eine hohe Qualität der Transkriptionsprofile hin. (B) Hierarchisches Clustering auf der Grundlage der Pearson-Korrelation der 500 variabelsten Gene auf der Grundlage der mittleren absoluten Abweichung als Maß für die Variabilität. Die Daten sind normalisiert und auf den Mittelwert zentriert. (C) Dreidimensionales Hauptkomponentenanalyse-Diagramm der Zelllinien. (D) Zweidimensionales t-SNE-Diagramm der Behandlungsbedingung und der Zelllinie zeigt eine klare Trennung der Zelllinien und dann der RNA-Fraktionen, was auf zwei dominante Quellen der Variabilität zwischen den RNA-seq-Profilen hinweist.
Qualitätskontrolle
Die Qualität der RNA-seq-Daten in öffentlichen Repositories ist weiterhin ein großes Problem. In einer kürzlich durchgeführten Studie von Deelen et al.19 musste mehr als die Hälfte von 65.000 verarbeiteten öffentlichen RNA-seq-Proben aufgrund von Qualitätsproblemen entfernt werden. Anstatt Proben zu entfernen, bietet GREIN einen umfassenden Bericht zur Qualitätskontrolle (QC) der Rohsequenzdaten und der Sequenzzuordnung für jede Probe (ergänzende Abb. S7) und ermöglicht es dem Benutzer, eine Entscheidung darüber zu treffen, welche Proben von nachgelagerten Analysen ausgeschlossen werden sollten.
Statistische Power-Analyse
Das Power-Analyse-Modul in GREIN erleichtert die Berechnung und Visualisierung der statistischen Power für den Nachweis unterschiedlich exprimierter Gene in zukünftigen Studien mit ähnlichen biologischen Proben. Die Schätzung einer angemessenen Stichprobengröße für künftige Studien mit ähnlichen biologischen Proben ist oft der wichtigste Motivationsfaktor bei der erneuten Analyse von RNA-seq-Daten. Die Power-Analyse erleichtert auch die Post-hoc-Analyse der Falsch-Negativ-Raten im aktuellen Datensatz. Fehlende statistische Leistung und Unterschiede in der statistischen Leistung zwischen Genen können zu falsch negativen Ergebnissen und damit zu falschen Schlussfolgerungen führen20. Das Segment „Leistungskurve“ liefert Leistungsschätzungen für eine unterschiedliche Anzahl von Proben auf der Grundlage eines einzelnen Gens (Abb. 3A). Der Benutzer kann die Standardwerte der Parameter ändern. Das Diagramm „Detektierbarkeit von Genen“ visualisiert die Leistungsschätzung für jedes der Gene auf der Grundlage der ausgewählten Gruppen und der genweisen Streuung (Abb. 3B). Die mittlere Abdeckung der Gene wird gegen ihre biologische Variabilität aufgetragen und in zwei Gruppen auf der Grundlage ihres Erkennbarkeitsstatus (Power ≥0,8 und Power < 0,8) dargestellt.
Power-Analyse zur Bewertung transkriptioneller Veränderungen in der nicht bösartigen MCF10A-Zelllinie. (A) Auf einem einzelnen Gen basierende Power-Schätzungen für eine unterschiedliche Anzahl von Proben in jeder Gruppe mit einer minimalen Fold-Change von 2 und statistischer Signifikanz α = 0,01. (B) Genweise Erkennbarkeit auf der log2CPM-BCOV-Ebene mit FDR ≤0,1 und zwei Proben in jeder Gruppe.
Differenzielle Genexpression
Die Erstellung und Interpretation differenzieller Genexpressionssignaturen ist ein typisches Analyseszenario in RNA-seq-Experimenten. Mit GREIN kann der Benutzer eine Signatur erstellen, indem er die Genexpression zwischen zwei Gruppen von Proben mit oder ohne Anpassungen für experimentelle Kovariaten oder Chargeneffekte vergleicht. GREIN kann mit komplexen Versuchsplänen umgehen, indem es die Flexibilität bietet, Gruppen und Untergruppen neu zu ordnen oder bestimmte Proben auszuwählen. Die differenzielle Expressionssignatur kann über interaktive Grafiken visualisiert werden, darunter eine Heatmap der am stärksten differenziell deregulierten Gene (ergänzende Abb. S15), geordnet nach der Falschentdeckungsrate (FDR), ein Diagramm der logarithmischen Faltenveränderung im Vergleich zur logarithmischen durchschnittlichen Expression (MA) (ergänzende Abb. S16) und ein Diagramm der Gennachweisbarkeit (ergänzende Abb. S17). Differentielle Expressionssignaturen, mit oder ohne Berücksichtigung potenziell falsch negativer Ergebnisse, können direkt in iLINCS exportiert werden, um Anreicherungs- und Konnektivitätsanalysen durchzuführen.
Anwendungsfall: Analyse der transkriptionellen und translationalen Regulation von Hypoxie in nicht-malignen Brustepithel- und dreifach-negativen Brustkrebszelllinien
Wir demonstrieren die Verwendung von GREIN anhand einer Neuanalyse von kürzlich veröffentlichten GEO RNA-seq-Daten (GSE104193). Sesé et al.21 untersuchten die Transkriptions- und Translationsregulation des hormonrefraktären dreifach-negativen Brustkrebs-Subtyps (TNBC) unter einer Kombination aus Hypoxie und mTOR-Inhibitor-Behandlung (mechanistic target of rapamycin). Die Autoren analysierten insbesondere die Expressionsprofile von TNBC-Zellen (MDA-MB-231) und nicht bösartigen Brustepithelzellen (MCF10A), die normoxischen (21 % O2) und hypoxischen (0,5 % O2) Bedingungen ausgesetzt und/oder mit dem mTORC1- und -2-Inhibitor PP242 behandelt wurden. Jede der Proben wurde für die gesamte (T) und polysomal gebundene (P) mRNA sequenziert. Der Datensatz enthält 32 Proben, die zwei biologische Wiederholungen für jede Kombination von Zelllinie, Sauerstoffgehalt, Behandlungsstatus und mRNA-Fraktion darstellen.
Die explorative Analyse des verarbeiteten Datensatzes in GREIN (Abb. 2) zeigt, dass die stärkste Variationsquelle zwischen den Proben die Unterschiede zwischen den beiden Zelllinien sind. Dies wird durch die Korrelationsanalyse der vollständigen Expressionsprofile (Abb. 2A), das hierarchische Clustering der 500 wichtigsten hochvariablen Gene auf der Grundlage der mittleren absoluten Abweichung (Abb. 2B), die 3D-PCA-Darstellung der Proben (Abb. 2C) und die 2D-t-SNE-Darstellung (Abb. 2D) bestätigt. Darüber hinaus weisen hohe Korrelationen zwischen Expressionsprofilen für dieselbe Zelllinie (Abb. 2A) auf ein gutes Signal-Rausch-Verhältnis bei den Genexpressionsmessungen hin. Die zusätzliche Substruktur der Daten, die durch den 2D t-SNE Plot angezeigt wird, wurde untersucht, indem die Proben nach verschiedenen Attributen gezeichnet wurden (ergänzende Abb. S1). Diese Analyse ergab, dass die Unterschiede innerhalb der einzelnen Zelllinien durch verschiedene mRNA-Fraktionen und dann durch Unterschiede zwischen den Versuchsbedingungen verursacht werden.
Als Nächstes verwendeten wir GREIN, um eine statistische Leistungsanalyse auf der Grundlage des Musters der in diesem Datensatz beobachteten biologischen Variabilität durchzuführen. Wir betrachteten die Transkriptionsprofile jeder Zelllinie, die der Hypoxie ausgesetzt war und mit oder ohne PP242 behandelt wurde, was zu vier Vergleichen führt. Unter der Annahme eines mindestens zweifachen Expressionsunterschieds zwischen den Gruppen bei einer statistischen Signifikanz von α = 0,01 und mit nur zwei Wiederholungen in jeder Gruppe liegt die statistische Aussagekraft eines Gens, das als differenziell exprimiert erkannt wird, bei allen Vergleichen unter 0,55 (Tabelle 2). Unsere Analyse zeigt, dass man vier Wiederholungen pro Gruppe benötigt, um eine 80%ige Power zu erreichen, die eine zweifache Veränderung der Expression nachweist (Tabelle 2 und Abb. 3A). In einem typischen RNA-seq-Experiment ist eine Sequenzierungstiefe von 20-30 Millionen ausreichend, um die Genexpression für fast alle Gene zu quantifizieren4,22 , was auch in diesem Datensatz deutlich wird. Wir bewerteten auch die statistische Leistung jedes Gens, das als differenziell exprimiert erkannt werden kann, anhand des Plots „Erkennbarkeit von Genen“. Die durchschnittlichen logarithmischen Werte der Zählungen pro Million (CPM) der Gene wurden gegen den genweisen biologischen Variationskoeffizienten (BCOV) aufgetragen, und die Leistung wurde für die entsprechenden Gene berechnet (Abb. 3B). Zur Schätzung der statistischen Signifikanz wurde eine kontrollierte Falschentdeckungsrate von 0,05 und ein erwarteter Prozentsatz wahrer Positiver von 10 % verwendet. Wir definieren ein Gen als nachweisbar als differenziell exprimiert unter hypoxischen Bedingungen, wenn seine Power 0,8 oder mehr beträgt. Wie erwartet, besteht eine inverse Beziehung zwischen BCOV und Power (Abb. 3B). Außerdem steigt die Aussagekraft zur Erkennung der differentiellen Expression eines Gens mit einem höheren log CPM oder einer höheren Effektgröße.
Eines der Ziele der Studie war die Analyse der transkriptionellen Veränderungen unter hypoxischen und normoxischen Bedingungen mit und ohne PP242-Behandlung in den Zelllinien MCF10A und MDA-MB-231. Wir erstellten Transkriptionssignaturen von Hypoxie und Hypoxie + PP242 in der gesamten mRNA durch differenzielle Expressionsanalyse zwischen Hypoxie- bzw. Hypoxie + PP242-Proben im Vergleich zu den Kontrollproben, wobei wir den Chargeneffekt berücksichtigten, indem wir „Replikat“ als Kovariate behandelten, und zwar für jede Zelllinie separat. Sowohl unter Hypoxie als auch unter Hypoxie + PP242 (Abb. 4A) fanden wir in der MCF10A-Zelllinie im Vergleich zu MDA-MB-231 eine höhere Anzahl unterschiedlich exprimierter Gene (DE), was darauf hindeutet, dass die Tumorzelllinie möglicherweise besser mit Hypoxie umgehen kann. Diese Analyse zeigte auch, dass die meisten nicht differenziell exprimierten Gene ebenfalls nicht nachweisbar sind, was darauf hindeutet, dass es sich um falsch negative Ergebnisse handeln könnte. Dies steht im Einklang mit der Power-Analyse, die zeigt, dass 4 Proben pro Gruppe erforderlich wären, um differenziell exprimierte Gene mit durchschnittlichem BCOV konsistent zu identifizieren. Um geringer exprimierte Gene zu identifizieren, wäre ein noch höherer Probenumfang erforderlich.
Differenzielle Expression und Nachweisbarkeit der Gene. (A) Die Anzahl der Gene (log10-Skala), die nicht differenziell exprimiert und nicht nachweisbar (NDE&NDT), differenziell exprimiert (DE) und nicht differenziell exprimiert, aber nachweisbar (NDE&DT) in den Vergleichen mit Normoxie für die gesamte mRNA-Fraktion. Wir bezeichnen ein Gen als nachweisbar (DT), wenn seine Leistung ≥0,8 ist, und als differenzierbar, wenn die FDR < 0,05 ist. (B) Das Gendetektionsdiagramm für den ersten Vergleich (MCF10A und Hypoxie), das die oben genannte Liste von Genen zusammen mit ihren jeweiligen Fold Changes (FC) visualisiert.
Um die differentiell exprimierten Gene im Hinblick auf die betroffenen biologischen Pfade zu interpretieren, haben wir die differentiellen Genexpressionssignaturen von Hypoxie über iLINCS an Online-Anreicherungstools (DAVID23, ToppGene24, Enrichr25 und Reactome26) übermittelt. Die eingereichten Signaturen enthielten eine kombinierte Liste von DE- und NDE&DT-Genen, die wahrscheinlich echte positive und echte negative Gene repräsentieren. Die Gene wurden auf der Grundlage eines Grenzwertes von 0,7 und 0,01 für die statistische Aussagekraft bzw. FDR ausgewählt. Abbildung 5 zeigt die mit ToppGene erzielten Anreicherungsergebnisse für die MCF10-Hypoxie-Signatur. Signifikant angereichert (FDR < 0.05) Zu den Top-10-Gene-Ontology (GO)-Kategorien aus ToppGene und dem funktionalen Annotationstool DAVID gehören Reaktion auf Hypoxie, Reaktion auf verminderte Sauerstoffkonzentration, Angiogenese, Regulierung der Zellproliferation, Oxidations-Reduktionsprozesse und Reaktion auf abiotische Reize, die in beiden Zelllinien vorkommen (ergänzende Tabelle S2 und ergänzende Tabelle S3). Die meisten dieser Kategorien stimmen mit der ursprünglichen Studie überein. Darüber hinaus identifizierte die ToppGene-Suite das Hypoxie-induzierte Faktor (HIF-1-alpha)-Transkriptionsfaktor-Netzwerk, das in beiden Zelllinien aktiviert wurde (ergänzende Tabelle S4 und ergänzende Tabelle S5).
Snapshot einiger signifikanter Pfad- und Gen-Ontologie (GO)-Kategorien aus ToppGene über iLINCS. Diese Kategorien wurden beim Vergleich zwischen Hypoxie und Normoxie in der MCF10A-Zelllinie anhand einer kombinierten Liste von DE- und NDE&DT-Genen gefunden. Die rote vertikale Linie ist der gewählte Grenzwert von 0,05.
Schließlich nutzten wir die GREIN-Verbindung mit iLINCS, um die hochgeladene Signatur mit den LINCS27-Konsensgen-Knockdown-Signaturen (CGS) zu „verbinden „18. Wir fanden 3.727 LINCS-Konsens-Gen-Knockdown-Signaturen, die signifikant (pValue < 0,05) mit unserer hochgeladenen Signatur verbunden waren. Die Zielgene der 100 am stärksten verknüpften Signaturen wurden für eine weitere Anreicherungsanalyse ausgewählt. Wir fanden die zelluläre Reaktion auf Hypoxie und die Regulierung des Hypoxie-induzierbaren Faktors (HIF) durch Sauerstoff in der Liste der 10 wichtigsten aktivierten Signalwege in beiden Zelllinien (ergänzende Tabelle S6 und ergänzende Tabelle S7). Während diese Analyse ähnliche angereicherte funktionelle Kategorien wie die ursprüngliche Anreicherungsanalyse ergibt, ergänzt sie die ursprüngliche Analyse, indem sie mehrere Zielgene aufzeigt, die nicht unterschiedlich exprimiert werden, obwohl sie ausreichend hoch exprimiert sind, um nach unseren Leistungsanalysen nachweisbar zu sein. Die Verknüpfung dieser beiden Ergebnisse deutet darauf hin, dass diese Gene potenzielle übergeordnete Regulatoren der Reaktion auf Hypoxie sind.