Wieloetapowe dowody genetyczne, transkrypcyjne i jednokomórkowe wskazują na MAP1LC3A jako priorytetowy gen spośród genów związanych z ferroptozą w glejakach wielopostaciowych

Podsumowanie

Genetycznie zakotwiczone, wieloetapowe podejście integrujące randomizację mendlowską, transkryptomikę guzów oraz analizy pojedynczych komórek pozwoliło wytypować MAP1LC3A jako gen związany z ferroptozą, który wiąże się z podatnością na glejakoludzkie i stanowi kandydata do przyszłej walidacji eksperymentalnej.

Streszczenie

Glejaki wielopostaciowe (GBM) pozostają wysoce agresywnym nowotworem złośliwym, a wkład genów związanych z ferroptozą w podatność na tę chorobę nie jest jeszcze w pełni poznany. W celu priorytetyzacji genów związanych z ferroptozą i powiązanych z GBM zastosowano zakotwiczone genetycznie, wieloetapowe ramy analityczne. Spośród 483 genów wyselekcjonowanych z bazy FerrDb V2, 315 posiadało kandydujące locus cech ilościowych ekspresji cis (cis-eQTLs) w bazie eQTLGen, 250 zachowało co najmniej trzy niezależne instrumenty po grupowaniu nierównowagi sprzężeń (linkage disequilibrium clumping), a 226 dostarczyło wiarygodnych szacunków randomizacji mendlowskiej metodą ważenia odwrotnością wariancji (IVW), wykorzystując badanie asocjacyjne całego genomu (GWAS) dla GBM obejmujące 6 183 przypadki i 18 169 osób z grupy kontrolnej. Trzydzieści cztery geny spełniły eksploracyjne kryteria odkrycia P < 0,05 oraz stopień fałszywych odkryć Benjamini–Hochberga (BH-FDR) < 0,20, wykazując zgodność kierunkową w szacunkach bajesowskiej ważonej randomizacji mendlowskiej (BWMR). Randomizacja mendlowska na etapie replikacji z wykorzystaniem cis-eQTLs krwi pełnej z bazy GTEx V10 potwierdziła cztery geny: ATG7, RPTOR, MAP1LC3A oraz CHMP6. Ocena w trzech niezależnych kohortach transkryptomicznych (guz–kontrola) wykazała, że MAP1LC3A był konsekwentnie niedoregulowany w tkance nowotworowej i wykazał istotny połączony szacunek efektów losowych (log₂ fold change, −1,273; 95% przedział ufności, −1,625 do −0,920; stopień fałszywych odkryć = 0,016), podczas gdy pozostałe trzy geny nie wykazały porównywalnego statystycznego wsparcia międzykohortowego. Następnie przeprowadzono analizę wirtualnego nokautu na poziomie pojedynczych komórek w zbalansowanym podzbiorze 2 400 komórek złośliwych wybranych z 4 916 kwalifikujących się komórek z 20 guzków GBM u dorosłych o dzikim typie IDH. W pięciu niezależnie zainicjowanych uruchomieniach zidentyfikowano odpowiednio 3, 15, 4 i 7 silnych genów downstream dla ATG7, RPTOR, MAP1LC3A i CHMP6. Powstałe zbiory konsensusowe obejmowały 17 unikalnych genów, z których RND3 był wspólny dla wszystkich czterech celów. Analiza Gene Ontology wskazała na wzbogacenie procesów adhezji komórkowej i procesów powierzchni komórkowej, natomiast żadne ścieżki KEGG ani Reactome nie pozostały istotne po korekcie na wielokrotne testowanie. Łącznie wyniki te priorytetyzują MAP1LC3A do przyszłych badań eksperymentalnych, jednocześnie odróżniając asocjację genetyczną, zgodność ekspresji w guzach i perturbacje obliczeniowe od definitywnych dowodów przyczynowości lub mechanizmu.

Wprowadzenie

Glioblastoma (GBM) pozostaje paradygmatycznym przykładem nowotworu opornego na leczenie. Mimo coraz bardziej precyzyjnej klasyfikacji molekularnej i opieki multidyscyplinarnej, trwałe poprawy w wynikach leczenia pacjentów były ograniczone¹. W przypadku pacjentów w dobrym stanie ogólnym aktualne postępowanie obejmuje maksymalną bezpieczną resekcję, po której następuje radioterapia z jednoczesnym i uzupełniającym podawaniem temozolomidu – schemat ten został ustanowiony w przełomowym badaniu randomizowanym i utrzymany we współczesnych wytycznych klinicznych^1,2. Niemniej jednak, rozproszona infiltracja oraz rozległa heterogenność komórkowa i molekularna ograniczają trwałą kontrolę nad chorobą, a u większości pacjentów ostatecznie dochodzi do progresji lub nawrotu, na które nie istnieje uniwersalnie skuteczna standardowa metoda leczenia^1,3. Ta uporczywa luka między postępami w charakterystyce choroby a wynikami klinicznymi podkreśla potrzebę zidentyfikowania biologicznie istotnych zależności molekularnych, które mogą stać się podstawą nowych strategii terapeutycznych w leczeniu GBM.

Ferroptoza jest zależną od żelaza formą regulowanej śmierci komórkowej, charakteryzującą się niekontrolowaną peroksydacją fosfolipidów oraz niewydolnością komórkowych systemów obrony antyoksydacyjnej, co odróżnia ją mechanistycznie od apoptozy i innych kanonicznych programów śmierci komórkowej^4,5. Proces ten ma szczególne znaczenie w przypadku GBM, gdzie zmiany genetyczne i plastyczność metaboliczna przekształcają homeostazę żelaza, równowagę redoks oraz metabolizm lipidów. Zintegrowane profilowanie genomiczne i lipidomiczne wykazało, że delecja CDKN2A powoduje redystrybucję utlenialnych wielonienasyconych kwasów tłuszczowych, tworząc tym samym zależną od genotypu podatność na ferroptozę w modelach GBM⁶. Podobnie, analizy porównawcze guzów pierwotnych i nawrotowych pozwoliły zidentyfikować związane z nawrotem zmiany w GPX4, ACSL4 oraz innych regulatorach ferroptozy⁷. Wykazano również, że eksperymentalna modulacja szlaków obrony przed ferroptozą wpływa na odpowiedź komórek GBM i modeli ksenoprzeszczepów na temozolomid⁸. Zbiorczo wyniki te wskazują na ferroptozę jako biologicznie uzasadnioną słabość terapeutyczną w GBM. Odzwierciedlają one jednak przede wszystkim powiązania ze stanem guza lub kontekstowe obserwacje eksperymentalne i nie rozstrzygają, czy konstytutywna zmienność w ekspresji genów związanych z ferroptozą przyczynia się do dziedzicznej podatności na GBM.

Większość badań nad ferroptozą w glejakach u ludzi opierała się na analizach różnicowej ekspresji, modelowaniu przeżywalności oraz subtypowniu molekularnym z wykorzystaniem zbiorów danych TCGA, CGGA i GEO^9,10. Chociaż badania te wykazały istotność prognostyczną programów transkrypcyjnych związanych z ferroptozą, ich obserwacyjny charakter nie pozwala określić, czy zmieniona ekspresja genów przyczynia się do podatności na GBM, czy jest raczej konsekwencją rozwoju nowotworu. Randomizacja mendlowska w skali całego transkryptomu pozwoliła następnie na zidentyfikowanie regulowanych genetycznie, zależnych od tkanki genów związanych z ryzykiem glejaka, natomiast nowsze badania oparte na loci cech ilościowych ekspresji (eQTL) i białek (pQTL) zaczęły priorytetyzować potencjalne cele terapeutyczne dla GBM^11,12. Niemniej jednak wcześniejsze analizy zazwyczaj przyjmowały podejście obejmujące cały transkryptom lub zorientowane na cele leków, zamiast oceniać predefiniowany, kompleksowy zestaw genów związanych z ferroptozą. Co więcej, niewiele badań zintegrowało duże badanie asocjacyjne całego genomu (GWAS) dla GBM z randomizacją mendlowską na etapie replikacji przy użyciu niezależnego zasobu eQTL, a następnie przeprowadziło ocenę w wielu kohortach transkryptomicznych nowotwór-kontrola. Rozróżnienie to jest istotne, ponieważ genetyczna regulacja ekspresji genów różni się znacznie w zależności od tkanki, a asocjacje eQTL pochodzące z krwi nie muszą odzwierciedlać efektów regulacyjnych w obrębie guzów mózgu^11,13. W związku z tym konieczne jest opracowanie zintegrowanego modelu łączącego asocjację genetyczną, randomizację mendlowską na etapie replikacji, transkryptomikę nowotworów w wielu kohortach oraz funkcjonalną predykcję na poziomie pojedynczych komórek pochodzących od pacjentów, aby zidentyfikować geny związane z ferroptozą, których udział w GBM jest potwierdzony przez zbieżne linie dowodów.

W niniejszym badaniu zbadano, czy genetycznie regulowana ekspresja genów związanych z ferroptozą wiąże się z podatnością na GBM. Randomizacja mendlowska na etapie odkrywania i etapie replikacji została połączona z analizami ekspresji genów w trzech niezależnych kohortach transkrypcyjnych (guz–kontrola). Geny potwierdzone w obu etapach randomizacji mendlowskiej zostały następnie ocenione w danych transkrypcyjnych pojedynczych komórek pochodzących od pacjentów z wykorzystaniem wirtualnej perturbacji w celu charakterystyki przewidywanych odpowiedzi transkrypcyjnych w komórkach złośliwych. Zamiast priorytetyzować geny wyłącznie na podstawie sygnatur ekspresji w guzach, ten wieloetapowy schemat wykorzystał najpierw dziedziczoną zmienność genetyczną, a następnie ocenił wzorce ekspresji istotne dla choroby wraz z profilami perturbacji obliczeniowych w rozdzielczości pojedynczych komórek. Wynikowe dowody zbieżne posłużyły do wyłonienia priorytetowych genów związanych z ferroptozą do przyszłych badań eksperymentalnych w GBM.

Protokół

Niniejsze badanie zostało zwolnione z oceny etycznej przez Komisję Etyki Medycznej Pierwszego Szpitala Ludowego w Zhaoqing (numer referencyjny B2026-08-03). W badaniu wykorzystano retrospektywnie zgromadzone, zanonimizowane dane genetyczne i transkrypcyjne na poziomie zbiorczym, w tym dane z kontrolowanym dostępem zatwierdzone przez Komitet Dostępu do Danych z EGAD00010001657 oraz zestawy danych uzyskane z GEO, eQTLGen i GTEx zgodnie z ich odpowiednimi warunkami dostępu i użytkowania. Nie rekrutowano nowych uczestników, nie pobierano biospecymentów ani nie uzyskano dostępu do możliwych do zidentyfikowania danych na poziomie uczestnika. Zgoda etyczna i świadoma zgoda na badania pierwotne zostały uzyskane przez odpowiednich twórców danych, a dane z kontrolowanym dostępem były wykorzystywane zgodnie z obowiązującą Umową o Dostępie do Danych.

Projekt badania

W niniejszym badaniu zastosowano wieloetapowy schemat w celu priorytetyzacji genów związanych z ferroptozą, które korelują z podatnością na glejak wielopostaciowy (GBM), oraz oceny ich efektów transkrypcyjnych istotnych dla przebiegu choroby (Rycina 1). W pierwszej kolejności geny związane z ferroptozą, wyselekcjonowane z bazy FerrDb V2, oceniono za pomocą randomizacji mendlowskiej (MR) dla dwóch próbek, wykorzystując dane o cis-ekspresyjnych locus cech ilościowych (cis-eQTL) oraz szerokie badanie asocjacyjne całego genomu (GWAS) dla GBM. Podstawową metodą przesiewową była MR z ważeniem odwrotnością wariancji (IVW), bayesowska randomizacja mendlowska z ważeniem (BWMR) zapewniła uzupełniającą ocenę odporności wyników, a niezależny zbiór danych eQTL posłużył do MR na etapie replikacji. Po drugie, geny potwierdzone w analizach genetycznych oceniono w trzech niezależnych kohortach transkrypomicznych (guz–kontrola), a następnie przeprowadzono metaanalizę międzykohortową. Po trzecie, dane transkrypomiczne z pojedynczych komórek pochodzących od pacjentów wykorzystano do przeprowadzenia wirtualnej perturbacji genów w komórkach złośliwych oraz identyfikacji powtarzalnych odpowiedzi transkrypcyjnych downstream. Odpowiedzi te scharakteryzowano następnie poprzez analizy wzbogacenia funkcjonalnego oraz analizy sieci wspólnych. Ogólnie rzecz biorąc, analizy genetyczne zaprojektowano w celu identyfikacji genów związanych z podatnością na GBM, natomiast analizy transkrypomiczne i analizy pojedynczych komórek służyły ocenie zgodności biologicznej i generowaniu hipotez do późniejszej walidacji eksperymentalnej.

Źródła danych

Geny związane z ferroptozą pozyskano z bazy FerrDb V2, co po harmonizacji symboli genów i usunięciu duplikatów dało 483 unikalne ludzkie geny kodujące białka¹⁴. Statystyki podsumowujące cis-ilościowych cech ekspresyjnych (cis-eQTL) krwi pełnej z konsorcjum eQTLGen wykorzystano jako zbiór danych ekspozycji w etapie odkrywania w ramach randomizacji mendlowskiej (MR), natomiast dane cis-eQTL krwi pełnej z wersji GTEx V10 posłużyły jako niezależny zbiór danych ekspozycji w etapie replikacji MR^11,15. Powiązania z wynikiem GBM uzyskano ze statystyk podsumowujących badań asocjacyjnych całego genomu (GWAS) z kontrolowanym dostępem, dostępnych w European Genome-phenome Archive, obejmujących 6183 przypadki i 18 169 kontroli pochodzenia europejskiego¹⁶. Ekspresję na poziomie tkankowym genów o priorytecie genetycznym oceniono w trzech niezależnych kohortach Gene Expression Omnibus (GEO): GSE196533, obejmującej 61 próbek glejaka stopnia 4 opisanych w metadanych jako GBM oraz dziewięć nieżełudkowych próbek mózgu; GSE4290, obejmującej próbki GBM oraz nieguzowe próbki mózgu pochodzące od pacjentów z epilepsją; oraz GSE116520, zawierającej sparowane próbki rdzenia guza i obszarów okołoguzowych wraz z nieżełudkowymi kontrolami^17,18,19. Dane Smart-seq2 pochodzące od pacjentów z GSE131928 wykorzystano do analizy wirtualnej perturbacji komórek złośliwych²⁰. Charakterystyka zbiorów danych oraz ich odpowiednie role analityczne zostały podsumowane w Tabeli 1. Wszystkie analizy przeprowadzono z wykorzystaniem wcześniej zgromadzonych, zanonimizowanych zbiorów danych, dla których w oryginalnych badaniach uzyskano zgodę komisji etycznej oraz świadomą zgodę uczestników.

Wybór instrumentów genetycznych i harmonizacja danych

Instrumenty kandydujące ograniczono do cis-ekspresyjnych loci cech ilościowych (cis-eQTLs) powiązanych z ekspresją genów na poziomie istotności całego genomu (P < 5 × 10⁻⁸). Warianty pogrupowano (clumping) z wykorzystaniem panelu referencyjnego dla populacji europejskiej z projektu 1000 Genomes Project, stosując próg nierównowagi sprzężeń (LD) r² < 0.001 w oknie 10 000 kb. Geny, które po grupowaniu posiadały mniej niż trzy niezależne instrumenty, zostały wykluczone z głównej analizy randomizacji mendlowskiej (MR). Ze względu na eksploracyjny cel przesiewowy, z góry określono minimum trzech instrumentów, aby zachować geny ze słabym, ale silnym wsparciem cis-eQTL, umożliwiając jednocześnie estymację IVW przy użyciu wielu instrumentów. Próg ten został uzupełniony rygorystycznym grupowaniem LD oraz filtrowaniem statystyk F; estymacje oparte na zaledwie trzech lub czterech instrumentach interpretowano z ostrożnością, a analizy wrażliwości przeprowadzano tylko wtedy, gdy było to metodologicznie zasadne. Siłę instrumentu oceniano dla każdego wariantu za pomocą statystyki F (F = β²/SE²), gdzie β i SE reprezentują odpowiednio szacunek efektu cis-eQTL i jego błąd standardowy. Warianty z F < 10 zostały wykluczone, aby zminimalizować błąd wynikający ze słabych instrumentów^21,22. Zbiory danych ekspozycji i wyniku zharmonizowano poprzez wyrównanie alleli efektowych i kierunków efektu. Wykluczono warianty powielone, warianty niedostępne w zbiorze danych GBM GWAS oraz warianty z niekompatybilnym kodowaniem alleli. Ponieważ częstotliwości alleli efektowych nie były dostępne dla GBM GWAS, warianty palindromiczne o niejednoznacznej orientacji nici zostały usunięte zamiast być wnioskowane. Z tego samego powodu nie przeprowadzono formalnego testu kierunkowości Steigera.

Analizy randomizacji mendlowskiej

Związek między genetycznie przewidywaną ekspresją genów a podatnością na GBM oceniono za pomocą dwupróbkowej randomizacji mendlowskiej (MR). W etapie odkrywania uwzględniono wyłącznie geny posiadające co najmniej trzy niezależne instrumenty cis-eQTL (cis-expression quantitative trait locus), a głównym podejściem analitycznym była metoda ważonej odwrotnością wariancji (IVW). Szacunki efektów przedstawiono jako ilorazy szans (OR) z 95% przedziałami ufności (CI) na jednostkę wzrostu genetycznie przewidywanej ekspresji genów. Aby uwzględnić wielokrotne testowanie w obrębie ocenianych genów, wartości P dla metody IVW skorygowano przy użyciu procedury Benjamini-Hochberga²³. Jako kandydaci eksploracyjni zachowano geny z wartością P < 0,05 oraz wskaźnikiem fałszywych odkryć Benjamini-Hochberga (BH-FDR) < 0,20. Wybrano ten stosunkowo liberalny próg FDR, aby uniknąć przedwczesnego wykluczenia potencjalnie istotnych genów na etapie odkrywania; w związku z tym status kandydata interpretowano w powiązaniu z późniejszymi analizami, a nie jako dowód potwierdzający. Następnie do kandydatów z etapu odkrywania zastosowano bayesowską ważoną randomizację mendlowską (BWMR), wykorzystując te same zharmonizowane instrumenty²⁴. Zgodność między wynikami IVW i BWMR oceniano na podstawie istotności statystycznej oraz kierunku efektu. W przypadkach, gdy pozwalała na to liczba dostępnych instrumentów, przeprowadzono test Q Cochrana, test przecięcia MR-Egger, analizę MR-PRESSO oraz analizy leave-one-out w celu oceny heterogeniczności, plejotropii horyzontalnej, wpływowych wariantów odstających oraz wpływu poszczególnych polimorfizmów pojedynczego nukleotydu (SNP)²².

MR na etapie replikacji przeprowadzono z wykorzystaniem danych cis-eQTL z krwi pełnej z GTEx V10. W przypadku dostępności tylko jednego instrumentu zastosowano metodę stosunku Walda, natomiast w przypadku genów z dwoma lub więcej instrumentami wykorzystano metodę IVW. Replikację zdefiniowano jako P < 0,05, BH-FDR < 0,20 oraz kierunek efektu zgodny z odpowiadającą mu szacunkową wartością z etapu odkrycia. Ponieważ kilka szacunków replikacyjnych GTEx oparto na zaledwie jednym lub dwóch instrumentach, zinterpretowano je jako dowody wspierające replikację, a nie jako niezależne dowody przyczynowości.

Międzykohortowa ocena transkryptomiczna

Cztery geny wyłonione w analizach randomizacji mendlowskiej (MR) w fazie odkrywania i replikacji zostały ocenione w trzech niezależnych kohortach transkryptomicznych. W przypadku GSE196533 surowe dane liczbowe z sekwencjonowania RNA analizowano przy użyciu DESeq2²⁵. Wykluczono geny z liczbą zliczeń poniżej 10 we wszystkich próbach z wyjątkiem dwóch, natomiast cztery geny docelowe zachowano niezależnie od filtrowania ekspresji. Ekspresję różnicową oceniano pomiędzy 61 próbkami glejaka stopnia 4, oznaczonymi jako GBM w złożonych metadanych, a dziewięcioma niepotworowymi próbkami mózgu.

W przypadku zestawu GSE4290 wykluczono cztery próbki bez jednoznacznej diagnozy histopatologicznej, co pozostawiło 77 próbek GBM oraz 23 nieguzowe próbki mózgu. Przetworzone intensywności mikromacierzy poddano transformacji log₂, normalizacji kwantylowej oraz analizie przy użyciu odpornych liniowych modeli empirycznego Bayesa zaimplementowanych w limma²⁶. W sytuacjach, gdy wiele sond mapowało do tego samego genu, wybierano sondę o najwyższej średniej ekspresji we wszystkich uwzględnionych próbkach, niezależnie od istotności ekspresji różnicowej.

Zbiór GSE116520 obejmował sparowane próbki z rdzenia guza oraz obszarów peritumoranych pochodzące od 17 pacjentów, a także osiem kontroli nienowotworowych. Zdeponowane dane ekspresji, poddane transformacji logarytmicznej i normalizacji kwantylowej, zostały przeanalizowane przy użyciu pakietu limma. Korelacje wewnątrzpacjenckie między próbkami z rdzenia guza a próbkami peritumoranymi uwzględniono poprzez blokowanie na poziomie pacjenta oraz funkcję duplicateCorrelation. Porównanie rdzenia guza z grupą kontrolną stanowiło uprzednio określone zestawienie dla metaanalizy międzykohortowej, natomiast porównania peritumorane oraz trend uporządkowany kontrola–peritumoralny–rdzeń guza oceniano oddzielnie.

Specyficzne dla badania zmiany log₂ fold change oraz błędy standardowe zestawiono, wykorzystując model efektów losowych z ograniczoną maksymalną wiarygodnością i wnioskowaniem Hartunga–Knappa, zaimplementowany w pakiecie metafor. Heterogeniczność międzybadawczą oceniono za pomocą statystyki Q Cochrana oraz I². Połączone wartości P dla czterech genów docelowych skorygowano za pomocą procedury Benjamini–Hochberg. Silne wsparcie transkrypcyjne zdefiniowano jako stopień fałszywych odkryć (FDR) w metaanalizie < 0,05, istotność FDR na poziomie kohorty w co najmniej dwóch zbiorach danych oraz zgodny kierunek efektu we wszystkich trzech kohortach.

Analiza wirtualnego nokautu na poziomie pojedynczych komórek

Dane Smart-seq2 pochodzące od pacjentów z zestawu GSE131928 zostały wykorzystane do oceny czterech genów zreplikowanych w MR w komórkach złośliwych. Dorosłe komórki złośliwe zidentyfikowano zgodnie z adnotacjami zawartymi w oryginalnym badaniu, a z każdego kwalifikującego się pacjenta pobrano losowo taką samą liczbę komórek, aby zminimalizować dysproporcje w reprezentacji pacjentów. Wirtualny knockout przeprowadzono przy użyciu scTenifoldKnk i powtórzono w pięciu niezależnych seriach. Geny z wartością P < 0,05 po korekcie Benjamini–Hochberg w pojedynczej serii uznano za istotne. Geny downstream, które pojawiły się w co najmniej trzech z pięciu serii, zdefiniowano jako główny zestaw konsensusowy, natomiast do analizy czułości zastosowano bardziej rygorystyczne kryterium czterech z pięciu serii. Wyniki te zinterpretowano jako obliczeniowe przewidywania zaburzeń transkrypcyjnych, a nie jako dowód bezpośredniej regulacji molekularnej.

Analiza wzbogacenia funkcjonalnego i wspólnych sieci

Analizę wzbogacenia funkcjonalnego przeprowadzono z wykorzystaniem genów downstream specyficznych dla celu, które zostały powtarzalnie zidentyfikowane w co najmniej trzech z pięciu uruchomień wirtualnego nokautu. Wzbogacenie ścieżek Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) oraz Reactome oceniono za pomocą jednostronnych testów hipergeometrycznych, przyjmując 1 004 geny ujęte w wnioskowaniu sieci jako tło zestawu genów. Wartości P korygowano oddzielnie dla każdej bazy adnotacji zgodnie z procedurą Benjamini–Hochberga; za statystycznie istotne uznano skorygowaną wartość P < 0,05.

W celu przedstawienia relacji między czterema celami knockoutu a ich konsensusowymi genami downstream skonstruowano sieć bipartycyjną. Geny powiązane z wieloma celami zidentyfikowano na podstawie ich wspólnego stopnia, a stopień pokrycia zestawów genów specyficznych dla danych celów określono za pomocą liczby części wspólnych oraz indeksów Jaccarda. Krawędzie sieci reprezentują powiązania między powtarzalnymi perturbacjami obliczeniowymi i nie powinny być interpretowane jako dowód bezpośrednich oddziaływań molekularnych.

Analiza statystyczna i powtarzalność

O ile nie określono inaczej, testy statystyczne były dwustronne, a porównania wielokrotne kontrolowano przy użyciu procedury Benjamini–Hochberg. Kryteria istotności specyficzne dla poszczególnych analiz opisano w odpowiednich podrozdziałach. Wszystkie analizy przeprowadzono przy użyciu języków R lub Python. W procedurach randomizowanych zastosowano wcześniej określone ziarna (seeds), a kod analityczny, wersje oprogramowania oraz szczegółowe ustawienia parametrów zarchiwizowano w celu zapewnienia powtarzalności. Wszystkie zbiory danych zostały zebrane wcześniej i zanonimizowane; zgody etyczne oraz świadome zgody uczestników uzyskano w ramach oryginalnych badań.

Wyniki

Wybór genów związanych z ferroptozą oraz instrumentów genetycznych

Z bazy danych FerrDb V2 pozyskano łącznie 483 genów związanych z ferroptozą (Rysunek 1). Spośród nich 315 genów dopasowano do zbioru danych eQTLGen, a każdy z nich posiadał co najmniej jeden kandydujący cis-ilościowy cecha locus ekspresji (cis-eQTL). Po przeprowadzeniu klasteryzacji w oparciu o nierównowagę sprzężeń (linkage disequilibrium clumping), 250 genów zachowało co najmniej trzy niezależne instrumenty kandydujące. Po wyszukaniu wariantów wynikowych, harmonizacji alleli oraz kontroli jakości, 26 genów dostarczyło poprawnych szacunków metodą ważonej odwrotnością wariancji (IVW) i zostało włączonych do analizy randomizacji mendlowskiej (MR) na etapie odkrywania (Plik uzupełniający 1). Wszystkie 3 578 instrumentów zachowanych w analizie na etapie odkrywania posiadały statystykę F >10 (minimum 29,72; mediana 70,76), co wskazuje na brak dowodów na obciążenie wynikające ze słabych instrumentów. Wśród 34 kandydujących genów na etapie odkrywania mediana statystyki F wyniosła 67,22, przy wartości minimalnej 29,72.

Etap odkrycia MR identyfikuje geny związane z ferroptozą powiązane z podatnością na GBM

Spośród 26 genów, dla których uzyskano poprawne szacunki IVW, 34 spełniło z góry określone kryteria etapu odkrywania: IVW P < 0,05 oraz stopień fałszywych odkryć według metody Benjamini–Hochberga (BH-FDR) < 0,20, co obejmowało 19 korelacji odwrotnych i 15 dodatnich z podatnością na GBM (Rycyna 2A). Najsilniejsze dowody statystyczne zaobserwowano dla RPTOR (OR = 0,809, 95% CI 0,737–0,87; P = 7,02 × 10⁻⁶; BH-FDR = 0,012) oraz PLA2G6 (OR = 1,568, 95% CI 1,281–1,920; P = 1,08 × 10⁻⁵; BH-FDR = 0,0012). Szacunki bayesowskiej ważonej randomizacji mendlowskiej (BWMR) były nominalnie istotne i zgodne kierunkowo z szacunkami IVW dla wszystkich 34 genów kandydatów (Rycyna 2A). Testy punktu przecięcia MR-Egger nie dostarczyły dowodów na kierunkową plejotropię horyzontalną. Test Q Cochrana wykazał heterogeniczność jedynie dla MAP1LC3A (P = 0,043), natomiast globalne testy MR-PRESSO nie zidentyfikowały istotnych zakłóceń spowodowanych wartościami odstającymi wśród 3 ocenianych genów. Analizy MR-PRESSO nie można było przeprowadzić dla SLC7A1, ponieważ dostępne były tylko trzy instrumenty (Plik uzupełniający 1). Analizy leave-one-out dla poszczególnych genów, wykresy porównania metod oraz wykresy lejka dla czterech genów zreplikowanych w dalszej kolejności przedstawiono na Rycynie uzupełniającej 1. 34 kandydaci z etapu odkrywania zostali następnie ocenieni przy użyciu niezależnego zbioru danych eQTL. Spośród nich 26 posiadało wystarczającą liczbę instrumentów do MR na etapie replikacji, a cztery spełniły z góry określone kryteria replikacji.

Niezależna replikacja MR potwierdza czterech kandydatów z etapu odkrywania

Z 34 kandydatów na etapie odkrycia 26 posiadało co najmniej jeden kwalifikujący się instrument cis-eQTL w całości krwi w wersji GTEx V10 i zostało włączonych do analizy MR na etapie replikacji. Trzynaście genów reprezentowanych było przez pojedynczy instrument i zostało przeanalizowanych przy użyciu wskaźnika Walda, natomiast pozostałe 13 genów posiadało dwa lub więcej instrumentów i zostało przeanalizowanych metodą IVW. Cztery geny spełniły wcześniej określone kryteria replikacji: P < 0,05, BH-FDR < 0,20 oraz kierunek efektu zgodny z szacunkiem z etapu odkrycia (Ryc. 2B; Plik uzupełniający 1).

Wyższa genetycznie przewidywana ekspresja ATG7 (OR = 0.523, 95% CI 0.30–0.831; P = 0.061; BH-FDR = 0.0976), RPTOR (OR = 0.718, 95% CI 0.563–0.915; P = 0.075; BH-FDR = 0.0976) oraz MAP1LC3A (OR = 0.830, 95% CI 0.717–0.959; P = 0.017; BH-FDR = 0.1012) wiązała się ze zmniejszoną podatnością na GBM. W przeciwieństwie do tego, wyższa genetycznie przewidywana ekspresja CHMP6 wiązała się ze zwiększoną podatnością (OR = 1.378, 95% CI 1.032–1.838; P = 0.0295; BH-FDR = 0.1916). Kierunki efektów dla wszystkich czterech genów były zgodne z tymi zaobserwowanymi w analizie na etapie odkrywania. Nie wykryto istotnej heterogeniczności wśród genów, dla których można było obliczyć statystykę Q Cochrana (Supplementary File 1). Ponieważ większość szacunków replikacyjnych oparto tylko na jednym lub dwóch instrumentach, formalne testy plejotropii horyzontalnej i zniekształceń przez wartości odstające były możliwe do zastosowania jedynie dla ograniczonego podzbioru genów (Supplementary File 1). Odpowiednie wykresy diagnostyczne dla MAP1LC3A, RPTOR i CHMP6 przedstawiono na Supplementary Figure 2. Gen ATG7 nie kwalifikował się do wieloinstrumentalnych analiz diagnostycznych, ponieważ jego szacunek replikacyjny został wyprowadzony z jednowymiarowego stosunku Walda.

Międzykohortowa ocena transkryptomiczna wskazuje na priorytetowe znaczenie MAP1LC3A

Cztery geny potwierdzone zarówno w etapie odkrywania, jak i w etapie replikacji analiz MR zostały ocenione w trzech niezależnych kohortach transkrypcyjnych reprezentujących różne platformy ekspresji (Rycina 3; Tabela 2; Supplementary Figure 3; Supplementary File 1). Ekspresja MAP1LC3A była konsekwentnie obniżona w tkance nowotworowej we wszystkich trzech kohortach: GSE19653 (log₂FC = −1.53, FDR dla całego transkryptomu = 3.21 × 10⁻⁸), GSE4290 (log₂FC = −1.243, FDR = 3.5 × 10⁻¹²) oraz GSE16520 w rdzeniu guza w porównaniu z nieopiołową kontrolą (log₂FC = −1.204, FDR = 9.78 × 10⁻⁸). W kohorcie GSE116520 ekspresja MAP1LC3A była również niższa w tkance okołoguzowej niż w nieopiołowych kontrolach (log₂FC = −1.056, FDR = 2.85 × 10⁻⁶), wykazując znaczący trend spadkowy od kontroli, przez tkankę okołoguzową, aż do rdzenia guza (współczynnik trendu = −0.531, FDR = 8.59 × 10⁻⁶).

Metaanaliza z efektami losowymi potwierdziła istotnie niższą ekspresję MAP1LC3A w tkance nowotworowej (połączone log₂FC = −1.273, 95% CI −1.625 do −0.920; Hartung–Knapp P = 0.041; BH-FDR = 0.016), bez dowodów na heterogeniczność między badaniami (I² = 0%; Supplementary File 1). Ekspresja RPTOR była konsekwentnie niższa we wszystkich trzech kohortach i osiągnęła istotność w skali całego transkryptomu w GSE4290, choć jego połączona ocena nie była istotna statystycznie (log₂FC = −0.258, 95% CI −0.65 do 0.139; BH-FDR = 0.196; I² = 42.3%). Ekspresja CHMP6 była konsekwentnie wyższa w tkance nowotworowej i osiągnęła istotność w GSE4290, podczas gdy połączona ocena pozostała nieistotna (log₂FC = 0.150, 95% CI −0.130 do 0.431; BH-FDR = 0.196; I² = 52.0%). ATG7 wykazał niewielkie, niespójne kierunkowo różnice między kohortami oraz brak istotnej połączonej asocjacji (log₂FC = 0.036, 95% CI −0.073 do 0.146; BH-FDR = 0.291; I² = 0%). Zatem spośród czterech genów powielonych metodą MR, MAP1LC3A wykazał najsilniejsze i najbardziej spójne dowody na różnicową ekspresję związaną z nowotworem.

Wirtualny knockout pojedynczych komórek ujawnia powtarzalne, specyficzne dla celu zaburzenia transkrypcyjne

Cztery geny zreplikowane metodą MR oceniono w 4 916 kwalifikujących się komórkach złośliwych pochodzących z 20 guzów GBM typu IDH-wild-type u dorosłych w zbiorze GSE131928. Geny ATG7, RPTOR, MAP1LC3A i CHMP6 wykryto odpowiednio w 42,78%, 45,89%, 46,89% i 31,90% kwalifikujących się komórek złośliwych, co uzasadniało ich włączenie do analizy wirtualnego nokautu (Rysunek dodatkowy 4; Plik dodatkowy 1). Aby zminimalizować brak równowagi w reprezentacji pacjentów, z każdego guza losowo pobrano 120 komórek, co pozwoliło uzyskać zbalansowany względem pacjentów zbiór danych obejmujący 2 40 komórek złośliwych. Każdy gen docelowy oceniano w pięciu niezależnie zainicjowanych uruchomieniach, co zaowocowało przeprowadzeniem 20 analiz wirtualnego nokautu.

Stosując z góry określone kryterium BH-FDR < 0,05 w co najmniej trzech z pięciu powtórzeń, wirtualny nokaut zidentyfikował trzy stabilne geny downstream dla ATG7, 15 dla RPTOR, cztery dla MAP1LC3A oraz siedem dla CHMP6 (Rysunek 4A; Rysunek uzupełniający 5). Zbiór konsensusowy dla RPTOR obejmował RND3, NKAIN4, CHI3L1, CDKN1A, BCAN, PDGFRA, OLIG1, LHFPL3, ENO2, HILPDA, LGALS3, ANXA1, CNTN1, NAMPT oraz SCRG1. Zbiór konsensusowy dla MAP1LC3A zawierał RND3, CD24, BCAN i S100B, natomiast zbiory konsensusowe dla ATG7 i CHMP6 zawierały odpowiednio trzy i siedem genów. Zastosowanie bardziej rygorystycznego kryterium istotności w co najmniej czterech z pięciu powtórzeń zredukowało zbiory konsensusowe do dwóch genów powiązanych z ATG7, dziewięciu genów powiązanych z RPTOR, jednego genu powiązanego z MAP1LC3A oraz czterech genów powiązanych z CHMP6. Łącznie analizy te pozwoliły zidentyfikować powtarzalne, specyficzne dla celu zaburzenia transkrypcyjne w ramach wnioskowanej sieci regulacyjnej komórek nowotworowych.

Analizy wzbogacenia funkcjonalnego i sieci wspólnych identyfikują zbieżne odpowiedzi związane z adhezją

Cztery specyficzne dla celów zestawy konsensusowe obejmowały 17 unikalnych genów downstream. Analiza sieciowa zidentyfikowała RND3 jako jedyny gen wspólny dla wszystkich czterech wirtualnych nokautów, podczas gdy BCAN, CD24 i NKAIN4 były wspólne dla trzech z czterech wirtualnych nokautów. CHI3L1, LHFPL3 i PDGFRA były wspólne dla dwóch celów, natomiast pozostałe dziesięć genów było specyficznych dla danego celu (Rysunek 4B,C). Największa bezwzględna część wspólna w parach wystąpiła między RPTOR a CHMP6, które miały sześć wspólnych genów downstream. Na podstawie podobieństwa Jaccarda największa proporcjonalna część wspólna została zaobserwowana między ATG7 a CHMP6 (indeks Jaccarda = 0.429), a następnie między RPTOR–CHMP6 i MAP1LC3A–CHMP6 (w obu przypadkach 0.375).

Analiza Gene Ontology dla połączonego konsensusowego zestawu 17 genów zidentyfikowała dziesięć istotnie wzbogaconych terminów po korekcie Benjamini–Hochberg (Rysunek 4D; Rysunek uzupełniający 6). Wzbogacone terminy procesów biologicznych obejmowały adhezję komórkową (BH-FDR = 0.028), pozytywną regulację proliferacji populacji komórek (BH-FDR = 0.028), odpowiedź zapalną (BH-FDR = 0.0139), pozytywną regulację kaskady ERK1/ERK2 (BH-FDR = 0.0165) oraz adhezję międzykomórkową (BH-FDR = 0.0196). Istotne terminy komponentów komórkowych obejmowały powierzchnię komórki, region zewnątrzkomórkowy, błonę plazmatyczną oraz macierz zewnątrzkomórkową, natomiast wiązanie węglowodanów było jedynym istotnie wzbogaconym terminem funkcji molekularnej. Adhezja komórkowa pozostała istotnie wzbogacona, gdy analiza została ograniczona do genów wspólnych dla co najmniej dwóch celów oraz gdy zastosowano bardziej rygorystyczne kryterium konsensusu czterech z pięciu serii. Żaden szlak KEGG ani Reactome nie pozostał istotny po korekcie BH.

Wzbogacenie specyficzne dla celu było najbardziej rozległe w przypadku RPTOR, którego 15-genowy zestaw konsensusowy został wzbogacony o cztery terminy procesów biologicznych i cztery terminy komponentów komórkowych (Rycina uzupełniająca 7). Zestaw konsensusowy MAP1LC3A był wzbogacony o adhezję komórkową (BH-FDR = 8,74 × 10⁻⁴) oraz rozwój ośrodkowego układu nerwowego (BH-FDR = 0,0364), natomiast zestaw konsensusowy CHMP6 był wzbogacony o adhezję komórkową (BH-FDR = 0,075). Żaden termin Gene Ontology nie osiągnął wartości BH-FDR < 0,05 dla trzygenowego zestawu konsensusowego ATG7.

DOSTĘPNOŚĆ DANYCH:

Publiczne zestawy danych transkrypomicznych są dostępne w bazie GEO pod numerami dostępu GSE196533, GSE4290, GSE16520 i GSE131928. Statystyki podsumowujące wyniki GBM zostały zdeponowane z kontrolowanym dostępem w Europejskim Archiwum Genomów i Fenomów (EGA), zestaw danych EGAD0101657 (https://ega-archive.org/datasets/EGAD0101657). Dostęp jest zarządzany przez odpowiedni Komitet Dostępu do Danych i wymaga zatwierdzonego wniosku oraz Umowy o Dostępie do Danych. Zgodnie z obowiązującą umową, autorzy nie są uprawnieni do redystrybucji plików ani ich deponowania w publicznym repozytorium. Podsumowania danych eQTL są dostępne z Konsorcjum eQTLGen oraz GTEx zgodnie z ich odpowiednimi warunkami dostępu i użytkowania. Skrypty analityczne wspierające niniejsze badanie zostały udostępnione jako Supplementary Coding File 1 oraz Supplementary Coding File 2.

Rycina 1: Projekt badania i schemat integracji dowodów dla priorytetyzacji genów związanych z ferroptozą w glejakach wielopostaciowych w oparciu o zakotwiczenie genetyczne. Geny związane z ferroptozą, wyselekcjonowane z bazy FerrDb V2, przypisano do zbioru eQTLGen, przesiano w celu znalezienia niezależnych instrumentów cis-lokusa cechy ilościowej ekspresji (cis-eQTL), a następnie oceniono za pomocą randomizacji mendlowskiej (MR) na etapie odkrywania. Z 483 wyselekcjonowanych genów 315 posiadało co najmniej jeden kandydujący cis-eQTL, 250 zachowało co najmniej trzy niezależne instrumenty po grupowaniu w oparciu o nierównowagę sprzężeń (LD), a 26 dostarczyło poprawnych szacunków metodą ważenia odwrotnością wariancji (IVW) po weryfikacji wariantów efektu końcowego i harmonizacji alleli. Trzydzieści cztery geny spełniło kryteria etapu odkrywania, po czym zastosowano bayesowską ważoną randomizację mendlowską (BWMR) w celu oceny odporności wyników. Dwadzieścia sześć genów można było następnie ocenić w MR na etapie replikacji, wykorzystując cis-eQTL krwi pełnej z bazy GTEx V10. Cztery geny (ATG7, RPTOR, MAP1LC3A i CHMP6) spełniły kryteria replikacji i zostały poddane dalszej ocenie w trzech niezależnych kohortach transkrypomicznych oraz za pomocą wirtualnej perturbacji w pochodzących od pacjentów komórkach nowotworowych. Integracja tych komplementarnych analiz pozwoliła na wytypowanie genu MAP1LC3A do dalszych badań. BWMR = bayesowska ważona randomizacja mendlowska; eQTL = lokus cechy ilościowej ekspresji; IVW = ważenie odwrotnością wariancji; LD = nierównowaga sprzężeń; MR = randomizacja mendlowska. Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Rycina 2: Solidność etapu odkrywania oraz niezależna replikacja genetycznie przewidywanych efektów genów związanych z ferroptozą w ryzyku wystąpienia glejaka wielopostaciowego. (A) Sparowane wykresy leśne porównujące oszacowania uzyskane metodą ważenia odwrotnością wariancji (IVW) oraz bayesowską ważoną randomizację mendlowską (BWMR) dla 34 genów spełniających kryteria etapu odkrywania metody IVW P < 0,05 oraz współczynnik fałszywych odkryć Benjamini–Hochberga (BH-FDR) < Kwadraty poprzedzające nazwy genów oznaczają geny, które zostały następnie potwierdzone w niezależnej analizie replikacyjnej. Trójkąt identyfikuje LPIN1, dla którego szacunki IVW i BWMR wykazały rozbieżne kierunki efektu. (B) Wykres leśny (forest plot) 26 genów ocenianych w zbiorze danych etapu replikacji. Przedstawiono szacunki IVW dla genów z co najmniej dwoma instrumentami, natomiast szacunki stosunku Walda (Wald-ratio) dla genów z pojedynczym instrumentem. Symbole wypełnione na pomarańczowo identyfikują geny ATG7, RPTOR, MAP1LC3A oraz CHMP6, które spełniły kryteria replikacji (P < 0,05 i BH-FDR < 0,20). Punkty i linie poziome reprezentują odpowiednio ilorazy szans (OR) oraz 95-procentowe przedziały ufności (CI); pionowa przerywana linia oznacza OR = 1. Ilorazy szans przedstawiono w skali logarytmicznej. GBM = glejak wielopostaciowy. Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Rysunek 3Transkryptomiczna ocena czterech genów zreplikowanych w badaniu MR w różnych kohortach. Ekspresja ATG7, RPTOR, MAP1LC3A i CHMP6 w GSE19653 (61 próbek glejaka stopnia 4, opisanych w złożonych metadanych jako GBM, oraz dziewięć niepowrotnych próbek mózgu); (A) GSE4290 (77 próbek GBM i 23 nieguzowych próbek mózgu); (B) oraz GSE116520 (17 próbek z rdzenia guza, 17 dopasowanych do pacjenta próbek okołoguzowych oraz osiem nieopuchowcowych próbek kontrolnych); (C) Pudełka wskazują medianę oraz rozstęp międzykwartylny (IQR), wąsy sięgają 1,5 × IQR, a kropki reprezentują poszczególne próbki. (D) Specyficzne dla badania zmiany krotności (log₂ fold changes) oraz metaanaliza z efektami losowymi porównująca tkankę guza lub rdzenia guza z nienowotworową tkanką mózgu. Punkty i linie poziome oznaczają szacunki specyficzne dla danego badania oraz 95% przedziały ufności, natomiast romby reprezentują połączone szacunki metodą ograniczonej maksymalnej wiarygodności (REML) z wnioskowaniem w modelu Hartunga–Knappa. Wartości dodatnie wskazują na wyższą ekspresję w tkance guza. FDR = stopień fałszywych odkryć (false discovery rate); GBM = glejak wielopostaciowy (glioblastoma); MR = randomizacja mendlowska (Mendelian randomization). Kliknij tutaj, aby wyświetlić powiększoną wersję tej ilustracji.

Rycina 4Konsensus międzysezyjny i konwergencja funkcjonalna wirtualnych knockoutów jednokomórkowych w złośliwych komórkach glejaka wielopostaciowego. (A) Liczba zidentyfikowanych stabilnych genów downstream dla każdego celu, określona przy użyciu z góry zdefiniowanego kryterium istotności w co najmniej trzech z pięciu powtórzeń oraz bardziej rygorystycznego kryterium czułości w czterech z pięciu powtórzeń. (B) Parzyste nakładanie się stabilnych genów downstream; komórki wskazują liczbę wspólnych elementów oraz współczynniki podobieństwa Jaccarda. (C) Sieć dwudzielna łącząca cztery cele wirtualnego nokautu (romby) z odpornymi genami downstream (koła). Kolory krawędzi oznaczają zaburzony cel, natomiast rozmiar kół i intensywność ich koloru wskazują liczbę celów współdzielących każdą odpowiedź downstream. Krawędzie reprezentują powiązania między powtarzalnymi zaburzeniami obliczeniowymi, a nie bezpośrednie oddziaływania molekularne. (D) Istotne wzbogacenie ontologii genów (Gene Ontology) dla połączonego konsensusowego zestawu 17 genów. Długość słupka reprezentuje −log₁₀(BH-FDR), linia przerywana wskazuje próg istotności (BH-FDR = 0,05), a kolory oznaczają proces biologiczny (BP), komponent komórkowy (CC) i funkcję molekularną (MF). Analiza wzbogacenia funkcjonalnego została przeprowadzona z wykorzystaniem tła sieci regulacyjnej zrównoważonej dla pacjentów obejmującej 104 geny. Żaden szlak KEGG ani Reactome nie pozostał istotny po korekcie Benjamini–Hochberg. Skróty: BH-FDR = współczynnik fałszywych odkryć Benjamini–Hochberga; BP = proces biologiczny; CC = komponent komórkowy; GO = Gene Ontology; KEGG = Kyoto Encyclopedia of Genes and Genomes; MF = funkcja molekularna. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

Tabela 1: Przegląd źródeł danych i ich ról analitycznych w badaniu. Liczba próbek reprezentuje obserwacje uwzględnione w niniejszych analizach. BH-FDR = stopa fałszywych odkryć Benjamini–Hochberg; cis-eQTL = cis-ilościowe cechy loci ekspresji; EGA = European Genome-phenome Archive; GBM = glejak wielopostaciowy; GTEx = Genotype-Tissue Expression; GWAS = badanie asocjacyjne całego genomu; IV = zmienna instrumentalna; MR = randomizacja mendlowska; RNA-seq = sekwencjonowanie RNA. Kliknij tutaj, aby pobrać tę tabelę.

Tabela 2: Międzykohortowe dowody transkryptomiczne dla czterech genów zreplikowanych metodą MR. Wartości reprezentują log₂ krotności zmian (fold change) dla GBM lub tkanki rdzenia guza w stosunku do nienowotworowej tkanki mózgu. Połączone szacunki uzyskano przy użyciu modeli efektów losowych z ograniczoną maksymalną wiarygodnością (REML) z wnioskowaniem Hartunga–Knappa. CI = przedział ufności; FDR = stopa fałszywych odkryć. Kliknij tutaj, aby pobrać tę tabelę.

Rysunek uzupełniający 1: Analizy wrażliwości randomizacji mendlowskiej na etapie odkrywania dla czterech powtórzonych genów. Dla każdego z genów MAP1LC3A, ATG7, RPTOR i CHMP6 przedstawiono odpowiednio: (A) analizę leave-one-out; (B) wykres rozrzutu porównujący metody; oraz (C) wykres lejkowy.Kliknij tutaj, aby pobrać ten plik.

Rysunek uzupełniający 2: Wykresy diagnostyczne randomizacji mendlowskiej na etapie replikacji dla trzech zreplikowanych genów z wieloma instrumentami. Dla genów MAP1LC3A, RPTOR i CHMP6 przedstawiono odpowiednio: (A) wykres rozrzutu porównujący metody oraz (B) wykres lejkowy. Wartość dla ATG7 została oszacowana przy użyciu wskaźnika Walda z pojedynczym instrumentem, w związku z czym gen ten nie kwalifikował się do wykresów diagnostycznych dla wielu instrumentów.Kliknij tutaj, aby pobrać ten plik.

Rysunek uzupełniający 3: Analiza głównych składowych trzech niezależnych kohort transkryptomicznych. (A) Kohorta sekwencjonowania RNA GSE19653. (B) Kohorta Affymetrix GPL570 GSE4290. (C) Kohorta Illumina GPL1058 GSE16520. Analizę głównych składowych przeprowadzono z wykorzystaniem 50 genów lub sond o największej wariancji wewnątrz kohorty. Każdy punkt reprezentuje próbkę biologiczną; kolory wskazują grupy tkanek; a etykiety osi informują o wariancji wyjaśnionej przez każdą główną składową.Kliknij tutaj, aby pobrać ten plik.

Rysunek uzupełniający 4: Wykrywalność czterech genów zreplikowanych w MR w dorosłych złośliwych komórkach GBM. (A) Ogólne wskaźniki wykrywalności ATG7, RPTOR, MAP1LC3A oraz CHMP6 wśród 4 916 złośliwych komórek z 20 guzów GBM typu dzikiego IDH u dorosłych w zbiorze GSE131928/SCP393. (B) Wskaźniki wykrywalności tych samych genów na poziomie pacjenta. Kolor wskazuje procent złośliwych komórek z TPM > 0.Kliknij tutaj, aby pobrać ten plik.

Rysunek uzupełniający 5: Powtarzalność sygnałów downstream w wirtualnych knockoutach dla różnych ziarn losowości (cross-seed). (A) Liczba genów downstream istotnych statystycznie (BH-FDR) zidentyfikowanych w pięciu niezależnych uruchomieniach dla każdego celu. Punkty reprezentują ziarna losowości, a poziome paski wskazują mediany. (B) Geny downstream istotne w co najmniej trzech z pięciu uruchomień. Oś X przedstawia liczbę istotnych uruchomień, kolory identyfikują zaburzony cel, a rozmiar punktu reprezentuje medianę statystyki Z scTenifoldKnk. Sam gen docelowy został wykluczony.Kliknij tutaj, aby pobrać ten plik.

Rysunek uzupełniający 6: Analizy wrażliwości wzbogacenia dla konsensusu połączonego, współdzielonego i o rygorystycznym progu. Wzbogacenie funkcjonalne (A) połączonego konsensusu zdefiniowanego przez istotność w co najmniej trzech z pięciu powtórzeń; (B) genów współdzielonych przez co najmniej dwa cele przy kryterium trzech z pięciu; (C) połączonego rygorystycznego konsensusu zdefiniowanego przez istotność w co najmniej czterech z pięciu powtórzeń; oraz (D) genów współdzielonych przez co najmniej dwa cele przy kryterium czterech z pięciu. Oś x przedstawia −log₁₀(nominal P), wielkość punktu odzwierciedla liczbę pokrycia, a kolory oznaczają bazę danych adnotacji. Wypełnione punkty osiągnęły BH-FDR < 0.05, natomiast puste punkty oznaczają terminy eksploracyjne z nominal P < 0.05. Przez cały czas stosowano tło sieci regulacyjnej obejmujące 1 04 geny.Prosimy kliknąć tutaj, aby pobrać ten plik.

Rysunek uzupełniający 7: Celowo swoista analiza wzbogaceń funkcjonalnych dla silnych odpowiedzi po wirtualnym nokaucie. Wzbogacenie genów downstream po wirtualnym nokaucie (A) ATG7; (B) RPTOR; (C) MAP1LC3A oraz (D) CHMP6. Oś x przedstawia −log₁₀(nominalne P), rozmiar punktu odzwierciedla liczbę wspólnych elementów, a kolory oznaczają GO: BP, GO: CC, GO: MF, KEGG lub Reactome. Punkty wypełnione osiągnęły BH-FDR < 0.05, natomiast punkty puste oznaczają terminy eksploracyjne z nominalnym P < 0.05. Tło dla analizy wzbogaceń stanowiło 1 04 geny tworzące zbalansowaną względem pacjentów sieć regulacyjną.Kliknij tutaj, aby pobrać ten plik.

Plik uzupełniający 1: Tabele uzupełniające wspierające wieloetapową priorytetyzację genów związanych z ferroptozą w kontekście podatności na glejakowioma. Ten plik uzupełniający zawiera wszystkie tabele uzupełniające wspierające analizy genetyczne, transkrypcyjne oraz analizy w pojedynczych komórkach. Obejmuje on screening i selekcję genów związanych z ferroptozą oraz instrumentów genetycznych; pełne wyniki randomizacji mendlowskiej z etapu odkrywania i replikacji wraz z analizami wrażliwości, w tym ocenami heterogeniczności, plejotropii horyzontalnej i analizami MR-PRESSO; charakterystykę kohort, analizy różnicowej ekspresji oraz metaanalizę międzykohortową genów priorytetyzowanych genetycznie; a także analizy wirtualnego nokautu w pojedynczych komórkach, oceny reprodukowalności, analizy wzbogacenia funkcjonalnego oraz wyniki sieci wspólnych.Kliknij tutaj, aby pobrać ten plik.

Uzupełniający plik z kodem 1: Skrypty w językach R i Python wykorzystane do przeprowadzenia randomizacji mendlowskiej, analiz transkryptomicznych, wirtualnego nokautu w pojedynczych komórkach, wzbogacenia funkcjonalnego oraz analiz sieci opisanych w niniejszym badaniu.Aby pobrać ten plik, kliknij tutaj.

Uzupełniający plik z kodem 2: Pomocnicze skrypty analityczne, rutyny wykresowania i narzędzia przepływu pracy wykorzystane do wygenerowania wyników badania, rycin oraz wyników uzupełniających.Aby pobrać ten plik, kliknij tutaj.

Dyskusja

W niniejszym badaniu zintegrowano analizy asocjacji genetycznych, randomizację mendlowską (MR) na etapie replikacji, transkryptomikę guzów oraz modelowanie w pojedynczych komórkach pochodzących od pacjentów, aby wytypować geny związane z ferroptozą, które są powiązane z podatnością na glejak wielopostaciowy (GBM). Przesiew 483 wyselekcjonowanych genów w oparciu o badanie asocjacji całego genomu (GWAS) dla GBM, obejmujące 6183 przypadki i 18 169 kontroli, pozwolił na zidentyfikowanie 34 kandydatów na etapie odkrywczym, z których czterech — ATG7, RPTOR, MAP1LC3A oraz CHMP6 — został potwierdzonych w analizie replikacyjnej z wykorzystaniem niezależnego zasobu loci cech ilościowych ekspresji (eQTL). Dowody stawały się coraz bardziej selektywne poza analizami genetycznymi. Gen MAP1LC3A był konsekwentnie wyciszany w trzech niezależnych kohortach guzów i pozostał istotny w metaanalizie międzykohortowej. Ten wzorzec ekspresji w guzie uzupełniał ochronne powiązanie zaobserwowane w analizach MR, mimo że oba podejścia odnoszą się do odmiennych aspektów biologii choroby. RPTOR i CHMP6 wykazały spójne kierunkowo, lecz mniej przekonujące dowody transkryptomiczne, natomiast dla ATG7 zabrakło powtarzalnego wsparcia na poziomie tkankowym. Wirtualny nokaut ujawnił ponadto specyficzne dla celu, choć częściowo nakładające się odpowiedzi transkrypcyjne w komórkach złośliwych. Łącznie te sekwencyjne warstwy analityczne doprecyzowały wstępne wyniki MR poprzez odróżnienie kandydatów o różnym stopniu poparcia istotnego dla przebiegu choroby, przy czym MAP1LC3A okazał się najsilniejszym ogólnym kandydatem.

Większość istniejących dowodów dotyczących ludzi, łączących ferroptozę z glejakiem, pochodzi z badań nad ekspresją w guzach. Analizy zbiorów danych TCGA, CGGA oraz innych publicznych kohort wielokrotnie identyfikowały sygnatury związane z ferroptozą, które korelują z przeżywalnością, stopniem złośliwości guza, charakterystyką molekularną i cechami immunologicznymi^9,27. Choć badania te potwierdziły znaczenie kliniczne programów transkrypcyjnych związanych z ferroptozą, profile ekspresji uzyskane z rozwiniętych guzów nie pozwalają odróżnić genów związanych z podatnością na chorobę od zmian transkrypcyjnych powstających podczas progresji nowotworu lub odzwierciedlających różnice w składzie komórkowym. Analizy genetyczne zapewniają uzupełniającą perspektywę. Robinson i współpracownicy zintegrowali dane GWAS dla glejaków z zestawami danych eQTL dla mózgu i krwi pełnej, wykorzystując MR oraz kolokalicację, co pozwoliło na wytypowanie domniemanych genów podatności o efektach zależnych od tkanki i wykazało ograniczoną zgodność między szacunkami pochodzącymi z krwi i mózgu¹¹. Niedawno zintegrowane badanie MR oparte na eQTL i pQTL połączyło dowody genetyczne z analizami różnicowej ekspresji i kolokalicacji, aby wytypować GPX7 i CXCL10 do dalszej oceny w GBM12. W przeciwieństwie do nich, niniejsze badanie rozpoczęło się od zdefiniowanego zestawu genów związanych z ferroptozą i oceniało genetycznie wytypowanych kandydatów poprzez MR na etapie replikacji, transkryptomikę guza oraz obliczeniową perturbację w rozdzielczości pojedynczych komórek. Stopniowe zawężanie wyników z 34 powiązań na etapie odkrywania do czterech zreplikowanych genów, a ostatecznie do MAP1LC3A jako jedynego genu wykazującego statystycznie istotną różnicową ekspresję w różnych kohortach, ilustruje wartość dyskryminacyjną integracji wielu komplementarnych podejść analitycznych. Co istotne, analizy transkryptomiczne nie miały na celu walidacji instrumentów genetycznych pochodzących z krwi, lecz ustalenie, czy geny wytypowane genetycznie wykazują również powtarzalne wzorce ekspresji istotne dla choroby.

MAP1LC3A budzi szczególne zainteresowanie, ponieważ niniejsze wyniki rozszerzają jego dotychczasową charakterystykę jako markera związanego z nowotworem oraz markera prognostycznego. Wcześniejsze bioinformatyczne badanie wielkohortowe włączyło MAP1LC3A do sygnatury sześciu genów związanych z przeżywalnością i nawrotami GBM, a także zgłosiło zmienioną metylację MAP1LC3A, choć jego wpływ na podatność na chorobę pozostał nieustalony²⁸. W niniejszej pracy wyższa genetycznie przewidywana ekspresja MAP1LC3A była konsekwentnie związana z niższą podatnością na GBM w obu etapach MR. Co więcej, MAP1LC3A był powtarzalnie obniżony w trzech niezależnych kohortach nowotworowych pomimo różnic w platformach ekspresyjnych, składzie próbek i metodologiach analitycznych, a zbiorcza ocena metaanalizy nie wykazała wykrywalnej heterogeniczności między badaniami. Chociaż wyniki te nie dowodzą, że zredukowana ekspresja MAP1LC3A inicjuje GBM, dostarczają one silniejszych dowodów łączących ten gen z podatnością na chorobę niż same analizy różnicowej ekspresji w nowotworach. MAP1LC3A koduje izoformy LC3A w obrębie ssaczej rodziny białek ATG8. Bai i współpracownicy wykazali, że wariant 1 LC3A ulega konjugacji z fosfatydyloetanolaminą w celu wytworzenia LC3A-II i lokalizuje się w autofagosomach podczas indukowanej autofagii²⁹. Wykazano również, że autofagiczna wymiana ferrytyny wpływa na wrażliwość na ferroptozę w komórkach GBM, w tym w warunkach deprywacji cystyny oraz w modelach zależnych od ALDH1A3^30,31. Jednakże badania te analizowały głównie całkowite LC3-II lub LC3B, a nie konkretnie MAP1LC3A. W przeprowadzonych analizach pojedynczych komórek wirtualna perturbacja MAP1LC3A wywołała powtarzalne odpowiedzi downstream, w których wzbogacone były procesy związane z adhezją komórkową. Razem spostrzeżenia te wskazują MAP1LC3A jako ukierunkowanego kandydata do badania, w jaki sposób regulacja związana z autofagią, podatność na ferroptozę i zachowanie komórek złośliwych przenikają się w GBM.

Pozostałe geny zreplikowane metodą MR otrzymały różny stopień poparcia w kolejnych analizach. Wyższa genetycznie przewidywana ekspresja RPTOR wiązała się z niższą podatnością na GBM w obu etapach MR, a jej poziom był konsekwentnie niższy we wszystkich trzech kohortach nowotworowych, choć połączona ocena nie osiągnęła istotności statystycznej. Wirtualny nokaut RPTOR wygenerował największy zestaw powtarzalnych zmian transkrypcyjnych w dół, z wzbogaceniem obejmującym sygnalizację ERK, odpowiedzi zapalne, proliferację komórek oraz adhezję komórkową. Choć wyniki te są zgodne z uznaną rolą RPTOR jako rusztowania mTORC1, skala odpowiedzi transkrypcyjnej nie powinna być interpretowana jako dowód na silniejszy efekt przyczynowy³². CHMP6 wykazał analogicznie zgodne powiązania MR w obu etapach, gdzie wyższa genetycznie przewidywana ekspresja wiązała się ze zwiększoną podatnością na GBM. Mimo że ekspresja CHMP6 była konsekwentnie podwyższona we wszystkich trzech kohortach nowotworowych, połączony przedział ufności obejmował wartość zerową, a heterogenność między badaniami była umiarkowana. Dowody eksperymentalne wykazujące, że zależna od CHMP6 naprawa błon ESCRT-III hamuje ferroptotyczną śmierć komórek, dostarczają prawdopodobnego kontekstu mechanistycznego, choć wyniki te uzyskano poza modelami GBM³³. Z kolei ATG7 wykazał zreplikowane ochronne powiązanie genetyczne, ale brak powtarzalnego wzorca ekspresji w guzach. Wirtualny nokaut zidentyfikował tylko trzy solidne geny w dół, a żadna kategoria funkcjonalna nie pozostała istotna po korekcie na wielokrotne testowanie. Poprzednie badania eksperymentalne wskazały na zaangażowanie zależnej od ATG7 autofagii w adaptację GBM i odpowiedź na leczenie^34,35, lecz obserwacje te nie wyjaśniają stosunkowo słabego poparcia międzyplatformowego zaobserwowanego w niniejszej pracy. W związku z tym RPTOR, CHMP6 i ATG7 pozostają prawdopodobnymi kandydatami drugorzędnymi, podczas gdy MAP1LC3A wykazał najsilniejszą zbieżność w analizach genetycznych, transkryptomicznych oraz w badaniach perturbacji obliczeniowych.

Analizy wirtualnych perturbacji nie wykazały pojedynczej ścieżki downstream wspólnej dla wszystkich czterech priorytetowych genów. Zamiast tego, powtarzalne odpowiedzi transkrypcyjne wykazały jedynie częściową zgodność, przy czym RND3 był jedynym genem downstream wspólnym dla wszystkich czterech sieci specyficznych dla celów. Najwyraźniejsza zbieżność funkcjonalna dotyczyła adhezji komórkowej oraz procesów zewnątrzkomórkowych lub powierzchniowych komórki, a wzbogacenie w adhezję komórkową pozostało istotne przy zastosowaniu bardziej rygorystycznego kryterium cross-seed. Żadna ścieżka KEGG ani Reactome nie pozostała istotna po korekcie na wielokrotne testowanie. Obserwacja ta jest istotna, ponieważ mimo że geny kandydackie wybrano z kuratorowanego zbioru genów związanych z ferroptozą, ich przewidywane efekty downstream w złośliwych komórkach GBM nie były zdominowane przez kanoniczne ścieżki ferroptozy. Raczej ich wkład w podatność GBM może obejmować szersze procesy komórkowe, w ramach których działa mechanizm związany z ferroptozą. Niniejsze analizy nie ustanawiają wspólnego mechanizmu molekularnego ani nie identyfikują RND3 jako mediatora przyczynowego. Zamiast tego wskazują one ograniczony zestaw programów komórek złośliwych, w szczególności tych związanych z adhezją komórkową, które wymagają dalszych badań eksperymentalnych.

Niniejsze badanie należy interpretować jako etapowe ramy priorytetyzacji, a nie jako ostateczne przypisanie genów przyczynowych. Żadna pojedyncza warstwa analityczna nie została uznana za rozstrzygającą; zamiast tego asocjacje na etapie odkrywania były sukcesywnie oceniane przy użyciu BWMR, niezależnego zasobu eQTL, trzech kohort transkrypcyjnych oraz modelowania regulacyjnego komórek złośliwych pochodzących od pacjentów. Należy uznać kilka ograniczeń. Po pierwsze, próg odkrywania BH-FDR < 0,20 był przeznaczony do przesiewania kandydatów, a nie do wnioskowania potwierdzającego, i tylko 26 z 34 kandydatów na etapie odkrywania można było poddać replikacji. Po drugie, kilka genów było reprezentowanych przez stosunkowo niewielką liczbę instrumentów genetycznych i nie przeprowadzono formalnych analiz mocy na poziomie genów; w konsekwencji słabe lub nieistotne asocjacje należy interpretować z ostrożnością. Próg kwalifikowalności trzech instrumentów zwiększył pokrycie genów, ale ograniczył zakres i stabilność analiz wrażliwości dla genów reprezentowanych przez zaledwie trzy lub cztery warianty. Chociaż wszystkie zachowane instrumenty na etapie odkrywania przekroczyły konwencjonalny próg F > 10, a kandydaci zostali dodatkowo ocenieni za pomocą BWMR i niezależnej replikacji, zabezpieczenia te nie kompensują w pełni rzadkości instrumentów; takie szacunki powinny zatem pozostać charakterem eksploracyjnym. Po trzecie, oba zasoby eQTL pochodziły z krwi pełnej i mogą nie oddawać dokładnie efektów regulacyjnych specyficznych dla mózgu lub guza. Po czwarte, dostępne statystyki podsumowujące GWAS dla GBM nie zawierały informacji wymaganych do testowania kierunkowości Steigera oraz formalnych analiz kolokalizacji. W konsekwencji nie jest pewne, czy sygnały asocjacji eQTL i GBM w każdym locus pochodzą z tego samego wariantu przyczynowego, czy z odrębnych wariantów w nierównowadze sprzężeń. Chociaż metoda BWMR została zaprojektowana tak, aby uwzględniać rozległą plejotropię horyzontalną i instrumenty odstające, zgodność między IVW a BWMR nie może wykluczyć plejotropii resztkowej ani zastąpić formalnych analiz kolokalizacji. Ponadto kohorty transkrypcyjne oceniały ustalone guzy, a nie podatność na chorobę, a jedna z kohort składała się z próbek glejaków stopnia 4, a nie wyłącznie z GBM typu dzikiego IDH. Wreszcie, analizy jednokomórkowe ograniczono do komórek złośliwych z jednego zestawu danych i modelowano zaburzenia regulacyjne obliczeniowo, a nie eksperymentalnie; nie oceniano komórek niemalignicznych w mikrośrodowisku guza ani nie odtworzono bezpośrednio zaburzeń genów in vitro lub in vivo. W związku z tym warianty przyczynowe, mechanizmy specyficzne dla typu komórek i konsekwencje biologiczne pozostają do ustalenia.

Wśród czterech powielonych genów, MAP1LC3A wykazał najbardziej spójne wsparcie w analizach genetycznych, transkrypcyjnych oraz analizach pojedynczych komórek. RPTOR, CHMP6 i ATG7 zachowały dowody z dwuetapowych analiz MR, lecz wykazały mniej spójne wsparcie w późniejszych analizach transkrypcyjnych i analizach perturbacyjnych. W związku z tym MAP1LC3A należy uznać za priorytetowego kandydata do dalszych badań, a nie za udowodniony gen przyczynowy lub cel terapeutyczny. Przyszłe badania powinny w pierwszej kolejności określić, czy sygnały asocjacji eQTL i GBM kolokalizują, wykorzystując pełne zestawy danych na poziomie locus wraz z zasobami regulacyjnymi specyficznymi dla mózgu lub guza. Następnie dwukierunkowe badania perturbacyjne w modelach GBM pochodnych od pacjentów mogłyby zbadać wrażliwość na ferroptozę, peroksydację lipidów, przeżywalność komórek oraz programy transkrypcyjne związane z adhezją, zidentyfikowane w analizach obliczeniowych. Takie eksperymenty będą niezbędne do odróżnienia efektów wpływających na uwarunkowaną genetycznie podatność na chorobę od tych, które wpływają na zachowanie utrwalonych komórek nowotworowych, oraz do bezpośredniego przetestowania zbieżnych asocjacji zidentyfikowanych w niniejszym badaniu.

Materiały

Lista materiałów użytych w tym artykule

Nazwa	Firma	Numer katalogowy	Komentarze
BWMR	pakiet R	BWMR	Bayesowska ważona randomizacja mendlowska
DESeq2	Bioconductor	Version 1.46.0	Analiza różnicowej ekspresji w sekwencjonowaniu RNA
FerrDb V2	FerrDb	Version 2	Źródło 483 wyselekcjonowanych genów związanych z ferroptozą
Glioblastoma bulk microarray	NCBI Gene Expression Omnibus	GSE4290	Kohorta oceny transkryptomicznej
Glioblastoma GWAS summary statistics	European Genome-phenome Archive	EGAD00010001657	Dane o wynikach z kontrolowanym dostępem; 6 183 przypadki i 18 169 kontrole
Glioblastoma Smart-seq2 single-cell RNA sequencing	NCBI Gene Expression Omnibus	GSE131928	Analiza wirtualnego nokautu komórek złośliwych
Grade 4 glioma bulk RNA sequencing	NCBI Gene Expression Omnibus	GSE196533	Kohorta oceny transkryptomicznej
GTEx whole-blood cis-eQTL summary statistics	Genotype-Tissue Expression project	GTEx V10	Dane ekspozycji w etapie replikacji
limma	Bioconductor	Version 3.62.2	Analiza różnicowej ekspresji w mikromierzach
metafor	pakiet R	Version 4.8-0	Metaanaliza z efektami losowymi
R	R Foundation for Statistical Computing	Version 4.4.2	Środowisko obliczeń statystycznych
scTenifoldKnk	pakiet R	Version 1.0.3	Analiza wirtualnego nokautu w pojedynczych komórkach
Tumour-core and peritumoural bulk microarray	NCBI Gene Expression Omnibus	GSE116520	Kohorta oceny transkryptomicznej
TwoSampleMR	pakiet R	Version 0.6.29	Randomizacja mendlowska dla dwóch próbek
Whole-blood cis-eQTL summary statistics	eQTLGen Consortium	eQTLGen	Dane ekspozycji w etapie odkrywania