GREIN: An Interactive Web Platform for Re-analyzing GEO RNA-seq Data

Zarys koncepcyjny GREIN jest przedstawiony na Rys. 1. Poszczególne zestawy danych RNA-seq są przetwarzane przez potok GREP2 i przechowywane lokalnie jako R Expression Sets. Użytkownik może uzyskać dostęp i analizować wstępnie przetworzone zbiory danych poprzez graficzny interfejs użytkownika GREIN (GUI) lub przesłać do przetworzenia zbiory danych, które nie zostały jeszcze przetworzone. Przepływy pracy oparte na GUI ułatwiają analizę i wizualizację danych, analizę statystyczną, konstrukcję sygnatur transkrypcyjnych oraz interpretację genów ulegających różnicy ekspresji (DE) w ramach biologii systemów. Zarówno GREIN jak i back-end pipeline (GREP2) zostały napisane w języku R i wydane odpowiednio jako kontener Docker i pakiet R. Graficzne interfejsy użytkownika dla GREIN są zaimplementowane w Shiny16, frameworku sieciowym do budowania dynamicznych aplikacji internetowych w R. Instancja sieciowa w https://shiny.ilincs.org/grein jest rozmieszczona poprzez solidny rój Docker z równoważonymi obciążeniowo serwerami Shiny. Kompletna infrastruktura GREIN, w tym potok przetwarzania, jest wdrażana za pomocą kontenerów Docker.

Przyjazne dla użytkownika przepływy pracy oparte na GUI w GREIN ułatwiają typowe scenariusze ponownego wykorzystania danych RNA-seq, takie jak badanie środków kontroli jakości i wizualizacja wzorców ekspresji w całym zbiorze danych, analiza wielkości próby i mocy w celu informowania o projektach eksperymentalnych przyszłych badań, statystyczna ekspresja różnicowa genów, wzbogacanie listy genów i analiza sieciowa. Poza standardowym porównaniem dwóch grup, moduł analizy różnicowej ekspresji genów wspiera również dopasowanie uogólnionego modelu liniowego, który uwzględnia zmienne lub efekty wsadowe. Zaimplementowane narzędzia interaktywnej wizualizacji i eksploracji obejmują analizę skupień, interaktywne mapy cieplne, analizę składowych głównych (PCA), t-dystrybuowane stochastyczne osadzanie sąsiadów (t-SNE), itp. (Supplementary Table S1). Użytkownik może również wyszukiwać adnotacje ontologiczne dla próbek ludzkiego RNA-seq i zbiorów danych dostarczonych przez projekt MetaSRA13. Każda przetworzona próbka ludzkiego RNA-seq jest opatrzona etykietą MetaSRA mapującą ontologie biomedyczne, w tym Disease Ontology, Cell Ontology, Experimental Factor Ontology, Cellosaurus i Uberon. Biologiczna interpretacja różnicowej ekspresji genów jest wspomagana przez bezpośrednie linki do innych narzędzi online do wykonywania typowych analiz post-hoc, takich jak analiza wzbogacenia listy genów i ścieżek oraz analiza sieciowa genów różnie wyrażonych (DE). Połączenie z tymi analitycznymi usługami internetowymi jest realizowane poprzez przesłanie sygnatury różnicowej ekspresji genów (tj. listy średnich zmian w ekspresji genów i związanych z nimi wartości p dla wszystkich analizowanych genów regulowanych w górę/w dół) do iLINCS17 (Integrative LINCS). iLINCS zapewnia również analizę łączności sygnatur dla ostatnio wydanych sygnatur Connectivity Map L100018. Szczegółowe instrukcje krok po kroku dotyczące przepływów pracy analizy GREIN znajdują się w Materiałach uzupełniających i sekcji 'Help’ w programie GREIN.

Kluczowe funkcjonalności

Wyszukaj lub prześlij do przetwarzania

Użytkownik może albo wyszukać już przetworzony zestaw danych GEO w polu 'Search for GEO series (GSE) accession’, albo przesłać zestaw danych do przetwarzania, jeśli zestaw danych nie jest jeszcze przetworzony (Supplementary Fig. S2). W tym momencie zdecydowana większość zbiorów danych GEO human, mouse i rat RNA-seq została wstępnie przetworzona, a przesłanie przez użytkownika zbiorów danych GEO do przetworzenia będzie wymagane tylko sporadycznie. Użytkownik może sprawdzić status przetwarzania żądanego zbioru danych w zakładce „Processing console” (Supplementary Fig. S3). Inne opcje wyszukiwania obejmują wyszukiwanie słów kluczowych w metadanych zbiorów danych oraz wyszukiwanie próbek w ontologiach biomedycznych poprzez adnotacje ontologiczne MetaSRA.

Explore dataset

GREIN umożliwia dostęp zarówno do surowych jak i znormalizowanych (counts per million i transcript per million) danych na poziomie genów i transkryptów. GREIN wyposażony jest w kilka interaktywnych i konfigurowalnych narzędzi do wizualizacji wzorców ekspresji, takich jak interaktywne mapy cieplne zgrupowanych genów i próbek, wykresy gęstości dla wszystkich lub podzbiorów próbek, analiza zmienności pomiędzy i wewnątrz grupy poprzez dwu- i trójwymiarowe analizy redukcji wymiarowości i wizualizacje takie jak PCA i t-SNE (Rys. 2). Użytkownik może również wizualizować profil ekspresji każdego genu osobno (Supplementary Fig. S6).

Kontrola jakości

Jakość danych RNA-seq w publicznych repozytoriach nadal stanowi poważny problem. W ostatnim badaniu przeprowadzonym przez Deelen et al.19, ponad połowa z 65 000 przetworzonych publicznych próbek RNA-seq musiała zostać usunięta z powodu problemów z QC. Zamiast usuwać próbki, GREIN dostarcza kompleksowy raport kontroli jakości (QC) surowych danych sekwencji i mapowania sekwencji dla każdej próbki (Supplementary Fig. S7), i pozwala użytkownikowi podjąć decyzję, które próbki powinny być wykluczone z dalszych analiz.

Statystyczna analiza mocy

Moduł analizy mocy w GREIN ułatwia obliczanie i wizualizację statystycznej mocy wykrywania różnie wyrażonych genów w przyszłych badaniach wykorzystujących podobne próbki biologiczne. Oszacowanie odpowiedniej wielkości próby dla przyszłych badań z podobnymi próbkami biologicznymi jest często kluczowym czynnikiem motywującym do ponownej analizy danych RNA-seq. Analiza mocy ułatwia również analizę post-hoc wskaźników fałszywie ujemnych w bieżącym zbiorze danych. Brak mocy statystycznej i różnice w mocy statystycznej między genami mogą dawać wyniki fałszywie negatywne, prowadzące do błędnych wniosków20. Segment „Krzywa mocy” dostarcza oszacowań mocy dla różnej liczby próbek opartych na pojedynczym genie (Rys. 3A). Użytkownik może modyfikować domyślne wartości parametrów. Wykres „Wykrywalność genów” wizualizuje szacunkową moc dla każdego z genów w oparciu o wybrane grupy i rozproszenie genów (Rys. 3B). Średnie pokrycie genów jest wykreślone względem ich zmienności biologicznej i jest wyświetlane w dwóch zestawach w oparciu o ich status wykrywalności (moc ≥0,8 i moc < 0,8).

Różnicowa ekspresja genów

Tworzenie i interpretacja sygnatur różnicowej ekspresji genów jest typowym scenariuszem analizy w eksperymentach RNA-seq. Dzięki GREIN, użytkownik może stworzyć sygnaturę porównując ekspresję genów pomiędzy dwoma grupami próbek z lub bez korekty dla zmiennych eksperymentalnych lub efektów partii. GREIN może obsługiwać złożone projekty eksperymentów, zapewniając elastyczność w rearanżacji grup i podgrup lub wybierając konkretne próbki. Sygnatura ekspresji różnicowej może być wizualizowana za pomocą interaktywnej grafiki, która zawiera mapę cieplną genów podlegających deregulacji (Supplementary Fig. S15) uszeregowanych według współczynnika fałszywego odkrycia (FDR), wykres log fold change vs. log średniej ekspresji (MA) (Supplementary Fig. S16) oraz wykres wykrywalności genów (Supplementary Fig. S17). Sygnatury ekspresji różnicowej, z lub bez uwzględnienia potencjalnie fałszywie ujemnych wyników, mogą być bezpośrednio wyeksportowane do iLINCS do analizy wzbogacania i łączności.

Przypadek użycia: Analiza transkrypcyjnej i translacyjnej regulacji hipoksji w niezłośliwych liniach komórkowych nabłonka piersi i potrójnie ujemnego raka piersi

Demonstrujemy użycie GREIN poprzez ponowną analizę niedawno opublikowanych danych GEO RNA-seq (GSE104193). Sesé i wsp.21 badali transkrypcyjną i translacyjną regulację podtypu TNBC (triple-negative breast cancer) opornego na leczenie hormonalne w połączeniu z hipoksją i inhibitorem mTOR (mechanistic target of rapamycin). W szczególności autorzy analizowali profile ekspresji komórek TNBC (MDA-MB-231) i komórek nabłonka niezłośliwego raka piersi (MCF10A) poddanych działaniu warunków normoksycznych (21% O2) i hipoksycznych (0,5% O2) i/lub traktowanych inhibitorem mTORC1 i -2 PP242. Każda z próbek została poddana sekwencjonowaniu dla całkowitego (T) i związanego z polisomami (P) mRNA. Zestaw danych zawiera 32 próbki, reprezentujące dwie repliki biologiczne dla każdej kombinacji linii komórkowej, poziomu tlenu, statusu leczenia i frakcji mRNA.

Eksploracyjna analiza przetworzonego zestawu danych w GREIN (Rys. 2) pokazuje, że najsilniejsze źródło zmienności między próbkami pochodzi z różnic między dwiema liniami komórkowymi. Potwierdza to analiza korelacji pełnych profili ekspresji (Ryc. 2A), hierarchiczne grupowanie 500 najbardziej zmiennych genów na podstawie mediany odchylenia bezwzględnego (Ryc. 2B), wykres 3D PCA próbek (Ryc. 2C) oraz wykres 2D t-SNE (Ryc. 2D). Ponadto, wysokie korelacje pomiędzy profilami ekspresji dla tej samej linii komórkowej (Rys. 2A) wskazują na dobry stosunek sygnału do szumu w pomiarach ekspresji genów. Dodatkowa podstruktura danych wskazana przez wykres 2D t-SNE została zbadana poprzez malowanie próbek według różnych atrybutów (Supplementary Fig. S1). Analiza ta ujawniła, że separacje w obrębie każdej linii komórkowej są indukowane przez różne frakcje mRNA, a następnie różnice między warunkami eksperymentalnymi.

Następnie użyliśmy GREIN do przeprowadzenia statystycznej analizy mocy w oparciu o wzór biologicznej zmienności obserwowanej w tym zbiorze danych. Wzięliśmy pod uwagę profile transkrypcyjne każdej linii komórkowej narażonej na hipoksję i traktowanej z lub bez PP242, co prowadzi do czterech porównań. Zakładając co najmniej dwukrotną różnicę w ekspresji między grupami, przy istotności statystycznej α = 0,01, i przy tylko dwóch powtórzeniach w każdej grupie, moc statystyczna genu do wykrycia jako różnie wyrażonego jest poniżej 0,55 we wszystkich porównaniach (Tabela 2). Nasza analiza wskazuje, że potrzebne byłyby cztery replikacje na grupę, aby osiągnąć 80% mocy wykrywania dwukrotnej zmiany w ekspresji (Tabela 2 i Rys. 3A). W typowym eksperymencie RNA-seq, głębokość sekwencjonowania 20-30 milionów jest wystarczająca do ilościowego określenia ekspresji genów dla prawie wszystkich genów4,22 , co jest również widoczne w tym zbiorze danych. Oceniliśmy również moc statystyczną każdego genu do wykrycia jako różnie wyrażonego z wykresu „Wykrywalność genów”. Średnie wartości log log zliczeń na milion (CPM) genów zostały wykreślone w stosunku do biologicznego współczynnika zmienności (BCOV) dla każdego genu, a moc została obliczona dla odpowiednich genów (Fig. 3B). Do oszacowania istotności statystycznej użyto kontrolowanego współczynnika fałszywego odkrycia 0,05 i oczekiwanego odsetka wyników prawdziwie pozytywnych 10%. Zdefiniowaliśmy gen jako wykrywalny jako różnie wyrażony w warunkach hipoksji, jeśli jego moc wynosi 0,8 lub więcej. Zgodnie z oczekiwaniami, istnieje odwrotna zależność pomiędzy BCOV i mocą (Rys. 3B). Również moc do wykrycia różnicowej ekspresji genu wzrasta wraz z wyższym log CPM lub wielkością efektu.

Jednym z celów badania była analiza zmian transkrypcyjnych w warunkach hipoksji i normoksji z i bez leczenia PP242 w obu liniach komórkowych MCF10A i MDA-MB-231. Stworzyliśmy sygnatury transkrypcyjne hipoksji i hipoksji + PP242 w całkowitym mRNA poprzez analizę ekspresji różnicowej pomiędzy próbkami hipoksji i hipoksji + PP242 odpowiednio w stosunku do próbek kontrolnych, korygując efekt serii poprzez traktowanie „replikacji” jako kowariantu, dla każdej linii komórkowej osobno. Stwierdziliśmy większą liczbę genów ulegających różnej ekspresji (DE) w liniach komórkowych MCF10A w porównaniu do MDA-MB-231 zarówno w hipoksji jak i hipoksji + PP242 (Rys. 4A) wskazując, że być może nowotworowa linia komórkowa jest lepiej przygotowana do radzenia sobie z hipoksją. Analiza ta wykazała również, że większość genów nie wykazujących ekspresji różnicowej również nie jest wykrywalna, co wskazuje, że mogą one reprezentować wyniki fałszywie negatywne. Jest to zgodne z analizą mocy pokazującą, że 4 próbki na grupę byłyby potrzebne do konsekwentnej identyfikacji różnie wyrażonych genów ze średnią BCOV. Aby zidentyfikować geny o niższej ekspresji, wymagana byłaby jeszcze większa wielkość próbki.

Aby zinterpretować różnie wyrażone geny w kategoriach dotkniętych ścieżek biologicznych, przesłaliśmy sygnatury różnicowej ekspresji genów hipoksji do narzędzi wzbogacania online (DAVID23, ToppGene24, Enrichr25 i Reactome26) za pośrednictwem iLINCS. Przesłane sygnatury zawierały połączoną listę genów DE i NDE&DT reprezentujących prawdopodobne wyniki prawdziwie pozytywne i prawdziwie negatywne. Geny zostały wybrane na podstawie wartości odcięcia 0,7 i 0,01 odpowiednio dla mocy statystycznej i FDR. Rycina 5 ilustruje wyniki wzbogacenia uzyskane z ToppGene dla sygnatury hipoksji MCF10. Istotnie wzbogacone (FDR < 0.05) 10 najlepszych kategorii ontologii genów (GO) z ToppGene i narzędzia funkcjonalnej anotacji DAVID obejmuje odpowiedź na hipoksję, odpowiedź na obniżony poziom tlenu, angiogenezę, regulację proliferacji komórek, proces utleniania-redukcji oraz odpowiedź na bodziec abiotyczny, które są wspólne dla obu linii komórkowych (Tabela Uzupełniająca S2 i Tabela Uzupełniająca S3). Większość z tych kategorii jest zgodna z badaniami oryginalnymi. Dodatkowo, pakiet ToppGene zidentyfikował sieć czynnika transkrypcyjnego indukowanego hipoksją (HIF-1-alfa), który był aktywowany w obu liniach komórkowych (Tabela uzupełniająca S4 i Tabela uzupełniająca S5).

Wreszcie, wykorzystaliśmy połączenie GREIN z iLINCS, aby „połączyć” przesłaną sygnaturę z LINCS27 konsensusowymi (CGS) sygnaturami knockdown genów18. Znaleźliśmy 3,727 LINCS konsensusowych sygnatur knockdown genów, które były znacząco (pValue < 0.05) połączone z naszą załadowaną sygnaturą. Geny docelowe 100 najlepszych połączonych sygnatur zostały wybrane do dalszej analizy wzbogacania. Znaleźliśmy odpowiedź komórkową na hipoksję i regulację czynnika indukującego hipoksję (HIF) przez tlen na liście 10 najlepszych aktywowanych ścieżek w obu liniach komórkowych (Tabela uzupełniająca S6 i Tabela uzupełniająca S7). Podczas gdy ta analiza daje podobne wzbogacone kategorie funkcjonalne jak początkowa analiza wzbogacenia, uzupełnia ona oryginalną analizę poprzez wskazanie kilku genów docelowych, które nie ulegają różnej ekspresji, chociaż są wystarczająco wysoko wyrażone, aby być wykrywalne zgodnie z naszymi analizami mocy. Powiązanie tych dwóch wyników razem implikuje te geny jako potencjalne regulatory wyższego poziomu odpowiedzi na hipoksję.