Zintegrowane analizy transkryptomiczne w ujęciu bulk, pojedynczych komórek oraz przestrzenne wskazują PPARG jako kandydujący gen prognostyczny związany ze starzeniem w kostniakomiaku

Podsumowanie

Wykorzystując analizy transkryptomiczne typu bulk, pojedynczych komórek oraz analizy przestrzenne, wraz z modelowaniem przeżycia opartym na uczeniu maszynowym i walidacją eksperymentalną, niniejsze badanie wskazuje PPARG jako kandydata na gen prognostyczny związany z procesem starzenia w osteosarcoma oraz wiąże jego obniżoną ekspresję ze złym rokowaniem przeżycia w bazie TARGET-OS oraz z cechami naczyniowymi i mikrootoczeniowymi.

Streszczenie

Kostniak złośliwy pozostaje wyzwaniem w przypadkach przerzutowych, nawrotowych lub opornych na terapię. Celem niniejszego badania była identyfikacja genów prognostycznych związanych z procesem starzenia się komórek oraz charakterystyka ich kontekstów przestrzennych. W zbiorze GSE99671 przeprowadzono analizę różnicowej ekspresji par za pomocą DESeq2, a następnie wyznaczono część wspólną z zestawem genów starzenia CellAge. Dane transkryptomiczne z kohorty TARGET-OS pobrano z bazy UCSC Xena. Geny kandydackie oceniono za pomocą jednowymiarowej regresji Coxa, analizy Kaplana-Meiera, analizy charakterystyki operacyjnej odbiornika (ROC) zależnej od czasu, regresji Coxa LASSO, powtarzalnej analizy LASSO oraz modelowania mgaławki przeżycia (random survival forest), przy czym współzmienne kliniczne włączono do skorygowanych modeli Coxa. Do charakterystyki wykorzystano analizę wzbogacenia funkcjonalnego, analizę mikrośrodowiska immunologicznego, transkryptomikę jednokomórkową, transkryptomikę przestrzenną SP_BS3, walidację ekspresji w GSE36001 oraz walidację za pomocą qRT-PCR i Western blot w komórkach kostniaka 143B oraz komórkach osteoblastów. W GSE99671 zidentyfikowano 2 248 genów o różnicowej ekspresji przy skorygowanym P < 0,05, a ich część wspólna z 866 genami CellAge dała 105 różnicowo wyrażonych genów związanych ze starzeniem. W TARGET-OS niższa ekspresja PPARG była powiązana z wyższym ryzykiem śmiertelności (jednowymiarowe HR = 0,603, 95% CI = 0,454–0,802, P = 0,000494; skorygowane HR = 0,224, 95% CI = 0,085–0,589, P = 0,00241). Dodanie PPARG do modelu klinicznego zwiększyło wskaźnik C-index z 0,707 do 0,829. PPARG był obniżony zarówno w GSE99671, jak i w GSE36001, a qRT-PCR i Western blotting potwierdziły niższą ekspresję mRNA i białka PPARG w komórkach kostniaka 143B niż w komórkach osteoblastów. Analiza jednokomórkowa zlokalizowała PPARG w komórkach śródbłonka, perycytach, makrofagach/monocytach oraz w komórkach zrębu związanych z nowotworem. Analiza transkryptomiki przestrzennej wykazała słabe, ale istotne dodatnie korelacje między ekspresją PPARG a wskaźnikami starzenia CellAge, wskaźnikami związanymi ze śródbłonkiem oraz wskaźnikami związanymi z perycytami. Wyniki te wskazują na PPARG jako kandydata na biomarker prognostyczny pochodzący z CellAge, związany z niekorzystnym przeżyciem i cechami naczyniowego mikrośrodowiska w kostniaku złośliwym, oraz wspierają badanie jego potencjalnego znaczenia dla stratyfikacji ryzyka i mikrośrodowiska nowotworowego związanego ze starzeniem komórkowym. Zbiór GSE36001 posłużył wyłącznie do zewnętrznej walidacji ekspresji; nie przeprowadzono niezależnej walidacji przeżycia.

Wprowadzenie

Kostniak osteogenny jest najczęstszym pierwotnym złośliwym nowotworem kości u dzieci, młodzieży i młodych dorosłych¹. Chociaż wielolekowa chemioterapia w połączeniu z chirurgią poprawiła wyniki leczenia w przypadku choroby zlokalizowanej^2,3, pacjenci z przerzutowym, nawracającym lub opornym na terapię kostniakiem osteogennym nadal wykazują niską przeżywalność długoterminową^3,4. Pojawiające się dowody wskazują, że indukowane stresem oksydacyjnym przebudowanie epigenetyczne może ułatwiać adaptację przerzutową i progresję nowotworu, co podkreśla złożoną plastyczność molekularną leżącą u podstaw agresywnych fenotypów raka⁵. Niezawodne biomarkery, które byłyby jednocześnie interpretowalne klinicznie i informatywne biologicznie, pozostają nieliczne. Dlatego kluczowe pozostaje zidentyfikowanie cech molekularnych, które oddają heterogeniczność kostniaka osteogennego i ryzyko prognostyczne na wielu poziomach danych.

Starzenie się komórek (senescencja) to stabilny program zatrzymania cyklu komórkowego indukowany przez dysfunkcję telomerów, uszkodzenia DNA, stres oksydacyjny, aktywację onkogenów oraz presję terapeutyczną⁶. Senescencja może ograniczać aberracyjną proliferację; jednakże komórki senescentne mogą również przekształcać mikrośrodowisko guza poprzez programy zapalne, wydzielanie chemokin, czynników wzrostu oraz przebudowę macierzy zewnątrzkomórkowej^7,8. W przypadku osteosarkoma geny związane z senescencją mogą odzwierciedlać zarówno wewnętrzne stany stresowe komórek nowotworowych, jak i niemaligniczne komponenty mikrośrodowiska, jednak ich znaczenie prognostyczne i organizacja przestrzenna nie zostały systematycznie ocenione.

PPARG koduje receptor gamma aktywowany przez proliferatory peroksysomów, który jest ligandowo aktywowanym receptorem jądrowym zaangażowanym w metabolizm lipidów, regulację zapalną, różnicowanie komórkowe oraz modulację odpornościową⁹. Rola PPARG w nowotworach zależy od kontekstu¹⁰. W niektórych przypadkach PPARG wiąże się z procesami różnicowania i stanami przeciwzapalnymi, podczas gdy w innych może wspierać adaptacyjne programy nowotworowe lub zrębowe. Niemniej jednak, wzorzec jego ekspresji, wartość prognostyczna oraz lokalizacja komórkowa i przestrzenna w kostniakomiaku nie zostały jeszcze w pełni scharakteryzowane.

W niniejszym badaniu zidentyfikowano geny o różnym poziomie ekspresji w zbiorze GSE99671, a następnie wyznaczono ich część wspólną z zestawem genów starzenia CellAge, co pozwoliło na wyłonienie 105 genów związanych ze starzeniem wykazujących różnicową ekspresję. Wykorzystując dane o przeżywalności z bazy TARGET-OS, wiele podejść opartych na uczeniu maszynowym i modelowaniu przeżywalności oraz korektę kliniczną, wskazano PPARG jako gen kluczowy. Dalsza charakterystyka PPARG została przeprowadzona przy użyciu analiz funkcjonalnych i mikrośrodowiska immunologicznego w badaniach bulk, transkryptomiki jednokomórkowej, transkryptomiki przestrzennej, zewnętrznej walidacji ekspresji z GSE36001 oraz walidacji metodą qRT-PCR i Western blot w linii komórkowej kostniakomięsaka 143B oraz w ludzkich osteoblastach.

Protokół

W niniejszym badaniu do analizy i walidacji wykorzystano publiczne zbiory danych oraz linie komórkowe; nie zaangażowano w nie uczestników będących ludźmi ani klinicznych próbek tkanek, w związku z czym zgoda komisji etycznej nie była wymagana.

Analiza różnicowej ekspresji w zbiorze GSE99671 i przecięcie z bazą CellAge
Surowe dane licznikowe (raw count data) oraz informacje o grupowaniu próbek dla GSE99671 pobrano z GEO^11,12. Przeanalizowano osiemnaście par próbek z kostniakomięsaka i dopasowanych próbek prawidłowych. Analizę różnicowej ekspresji par przeprowadscono przy użyciu pakietu DESeq2 z formułą projektu ~ pair_id + condition, gdzie pair_id uwzględniał efekty parowania osobników, a condition porównywał tkankę nowotworową z tkanką prawidłową¹³. Zachowano geny, których liczba odczytów wynosiła co najmniej 10 w co najmniej trzech próbkach. Różnicową ekspresję zdefiniowano jako skorygowane P < 0,05, a do wizualizacji zastosowano bardziej rygorystyczny próg: skorygowane P < 0,05 oraz |log2FC| ≥ 1. Geny różnicowo eksponowane przecięto z listą 866 genów starzenia z bazy CellAge po konwersji symboli genów na wielkie litery¹⁴. Analizę różnicowej ekspresji przeprowadzono w środowisku R (wersja 4.3.2) z użyciem DESeq2 (wersja 1.40.2), a skorygowane wartości P obliczono metodą Benjamini-Hochberga.

Kohorta TARGET-OS i modelowanie prognostyczne
Transkrypcyjne i kliniczne dane z projektu TARGET-OS pozyskano z bazy UCSC Xena¹⁵. Wyekstrahowano wartości ekspresji dla 105 genów o różnej ekspresji związanych z procesem starzenia komórkowego. Uwzględniono osiemdziesięciu pięciu pacjentów z pełnymi danymi dotyczącymi czasu przeżycia, statusu przeżycia oraz ekspresji genów kandydackich, w tym 27 zdarzeń zgonu. Standaryzowane wartości ekspresji genów przeanalizowano za pomocą jednowymiarowej regresji Coxa¹⁶, analizy przeżycia Kaplana-Meiera oraz analizy charakterystyki operacyjnej odbiornika (ROC) zależnej od czasu¹⁷. Do oceny stabilności wyboru i istotności zmiennych wykorzystano regresję Coxa metodą LASSO¹⁸, powtórną analizę stabilności LASSO oraz modelowanie za pomocą losowego lasu przeżycia (random survival forest)¹⁹. Zintegrowany wynik hub (hub score) oraz zintegrowany ranking kliniczny obliczono zgodnie z poniższymi jednoznacznymi kryteriami binarnymi. Analizy przeżycia przeprowadzono w środowisku R przy użyciu pakietów survival (wersja 3.5-7), timeROC (wersja 0.4), glmnet (wersja 4.1-8) oraz randomForestSRC (wersja 3.2.2). W przypadku jednowymiarowej analizy przesiewowej Coxa dla 105 genów kandydackich zastosowano korektę współczynnika odkryć fałszywych (FDR) według metody Benjamini-Hochberga; geny z FDR < 0.05 uznano za istotne statystycznie.

Preprocessing modelu, grupowanie PPARG i zależna od czasu analiza ROC
W grupie 85 pacjentów, wśród których odnotowano 27 zgonów, wykluczono geny o zerowej wariancji, brakujące wartości ekspresji genów kandydujących uzupełniono za pomocą imputacji medianą, a każdy gen kandydujący poddano standaryzacji z-score. Do analizy Kaplana-Meiera ekspresję zdychotomizowano na podstawie mediany kohorty: wartości ściśle powyżej mediany przypisano do grupy wysokiej ekspresji, a wartości równe lub poniżej mediany do grupy niskiej ekspresji (PPARG: n = 42 wysoka i n = 43 niska). Testy log-rank były dwustronne. W analizach ROC zależnych od czasu wykorzystano pakiet timeROC z parametrem event cause = 1, brzegowym ważeniem odwrotnością prawdopodobieństwa cenzurowania, czasami ewaluacji wynoszącymi 365, 1095 i 1825 dni oraz iid = FALSE. Aby zapewnić, że wyższe wartości markerów konsekwentnie wskazują na wyższe ryzyko, w przypadku genów z dodatnimi współczynnikami Coxa stosowano wystandaryzowane wartości ekspresji, natomiast dla genów ochronnych z ujemnymi współczynnikami wartości te mnożono przez -1.

LASSO i powtórzone LASSO
Model Cox LASSO dopasowano za pomocą pakietu glmnet, przyjmując family = "cox", alpha = 1, wcześniejszą standaryzację z-score (zatem standardize = FALSE), pięciokrotną walidację krzyżową, type.measure = "deviance" oraz losowy seed 123. Jako główne rozwiązanie współczynników wykorzystano lambda.min. Analiza stabilności polegała na 300-krotnym powtórzeniu tej samej pięciokrotnej walidacji krzyżowej; powtórzenie b wykorzystywało seed 1000 + b (b = 1,...,300). Dla każdego genu częstotliwość selekcji określono jako proporcję powtórzeń z niezerowym współczynnikiem przy lambda.min; odnotowano również niezerową selekcję przy lambda.1se.

Losowy las przeżycia
Las przeżycia dopasowano do wszystkich 105 znormalizowanych genów kandydujących przy użyciu randomForestSRC (wersja 3.2.2), stosując ziarno 123, ntree = 1,000, importance = TRUE oraz na.action = "na.impute". Dla danych przeżycia zachowano ustawienia domyślne pakietu: podział log-rank, mtry = 11 (sufiks pierwiastka kwadratowego z 105 predyktorów), minimalny rozmiar węzła końcowego = 15, nsplit = 10 losowych punktów podziału, próbkowanie bez zastawiania z frakcją próbkowania 0.632 oraz istotność zmiennych opartą na anti-split.

Zintegrowany wynik hubu
Każdy ze 105 kandydujących genów o zróżnicowanej ekspresji związanych z CellAge otrzymał jeden punkt za każde z sześciu kryteriów binarnych: (1) przynależność do części wspólnej zróżnicowanej ekspresji/CellAge (wszyscy kandydaci otrzymali ten punkt, ponieważ wymaganą wartość przed punktacją było skorygowane P < 0,05); (2) nominalne jednowymiarowe P Coxa < 0,05; (3) P testu log-rank Kaplan-Meier < 0,05; (4) średnie zależne od czasu AUC dla 3 i 5 lat ≥ 0,65; (5) częstotliwość wyboru lambda.min w powtarzalnym LASSO na poziomie 70. percentyla lub powyżej wśród kandydatów oraz > 0; oraz (6) istotność w random-survival-forest na poziomie 70. percentyla lub powyżej wśród kandydatów oraz > 0. Wszystkie kryteria miały taką samą wagę jednostkową, co dało wynik hubu od 0 do 6; geny z wynikiem ≥ 4 zakwalifikowały się do korekty klinicznej (18 genów). Raportowano również wartości FDR dla jednowymiarowego modelu Coxa, a FDR < 0,05 służyło do oznaczenia istotności przy wielokrotnym testowaniu, jednak wstępnie określony wskaźnik punktowy opierał się na nominalnym P < 0,05.

Ranking zintegrowany klinicznie
Po połączeniu danych o ekspresji i rekordów klinicznych, skorygowane analizy objęły 40 pacjentów z pełnymi zapisami współzmiennych oraz 13 zgonów. Końcowy wynik stanowił suma początkowego wyniku hub plus jeden punkt za każde z siedmiu kryteriów: skorygowane P Cox < 0,05, skorygowane P Cox < 0,10, czułe P Cox < 0,05 po usunięciu operacji radykalnej, czułe P Cox < 0,10, średnia AUC dla 3 i 5 lat ≥ 0,65, P testu ilorazu wiarygodności < 0,10 dla modelu klinicznego z genami w porównaniu z modelem tylko klinicznym oraz delta AIC < 0. Ponieważ progi 0,05 i 0,10 były zagnieżdżone, wartość P < 0,05 zapewniała dwa punkty, co nadawało większą wagę wyraźnie istotnym dowodom z analizy Coxa (skorygowanej i czułej). Całkowity zakres wyników wynosił od 0 do 13; remisy rozstrzygano na podstawie mniejszej wartości skorygowanego P Cox, a następnie większej średniej AUC dla 3/5 lat. PPARG otrzymał wszystkie sześć punktów początkowych oraz wszystkie siedem punktów z integracji klinicznej (13/13), zajmując pierwsze miejsce. Poprawa wskaźnika C-index została opisana opisowo i nie przypisano jej punktacji.

Korekta kliniczna
Kandydackie geny hub zostały zintegrowane ze zmiennymi klinicznymi z bazy TARGET-OS, w tym z płcią, wiekiem, stanem choroby w momencie diagnozy, miejscem wystąpienia guza pierwotnego, konkretnym regionem guza oraz operacją radykalną. Do analiz z korektą kliniczną włączono czterdziestu pacjentów z pełną dokumentacją ekspresji i danymi klinicznymi, w tym 13 zgonów. Modele Coxa zawierające wyłącznie zmienne kliniczne porównano z modelami uwzględniającymi zmienne kliniczne oraz ekspresję genów. Do oceny poprawy modelu wykorzystano wskaźnik C-index, kryterium informacyjne Akaike (AIC) oraz wartości P z testu ilorazu wiarygodności. Analizę wrażliwości przeprowadzono po usunięciu zmiennej dotyczącej operacji.

Analiza wzbogacenia funkcjonalnego i mikrośrodowiska odpornościowego
Próbki TARGET-OS zostały podzielone na warstwy zgodnie z ekspresją PPARG. Różnice w ekspresji między grupami PPARG-high a PPARG-low wykorzystano do wygenerowania rankingowych list genów do analizy wzbogacenia zbiorów genów (GSEA)²⁰. Przedstawione ścieżki to: Nemeth Inflammatory Response LPS Up, Burton Adipogenesis 5, Burton Adipogenesis 6, Krieg KDM3A Targets Not Hypoxia, Reactome: Transcriptional Regulation By TP53, Fulcher Inflammatory Response Lectin Vs LPS Dn, Hollmann Apoptosis Via CD40 Dn, Zhou Inflammatory Response Live Dn, WP: Fatty Acids And Lipoproteins Transport In Hepatocytes, Sweet Lung Cancer KRAS Up, KEGG Medicus Pathogen HIV Tat To TLR2/4 NF-kB Signaling Pathway oraz Reactome: Fatty Acids. Obliczono całkowity wynik starzenia CellAge, a analizę korelacji Spearmana zastosowano do oceny powiązań między PPARG a genami starzenia lub sygnaturami mikrośrodowiska odpornościowego²¹. Różnice w wynikach mikrośrodowiska między grupami PPARG-high a PPARG-low oceniano za pomocą testów nieparametrycznych z korekcją na wielokrotne testowanie.

Analiza transkryptomiczna pojedynczych komórek
Opublikowany zestaw danych transkryptomicznych pojedynczych komórek z ludzkiego kostniakomięsaka przeanalizowano przy użyciu wstępnie przetworzonego obiektu, dla którego przeprowadzono już kontrolę jakości, redukcję wymiarowości, klasteryzację oraz adnotację ręczną^22,23,24. Uwzględniono łącznie 68 336 komórek i 32 297 genów. Na potrzeby interpretacji tekstu głównego adnotacje uproszczono do 13 głównych typów komórek: limfocytów B, CAF, komórek w cyklu komórkowym, komórek śródbłonka, komórek erytroidalnych, makrofagów/monocytów, złośliwych komórek kostniakomięsaka, komórek miogennych, neutrofili, komórek podobnych do osteoklastów, perycytów, limfocytów T/NK oraz stromalnych komórek związanych z guzem. Do wizualizacji lokalizacji PPARG wykorzystano wykresy redukcji wymiarowości, ekspresji cech, wykresy kropkowe (dot plots) oraz skrzypcowe (violin plots). Komórki PPARG-dodatnie zdefiniowano jako komórki o ekspresji większej od zera. Różnice między typami komórek oceniano za pomocą testów Kruskala-Wallisa i Wilcoxona z korekcją Benjamini-Hochberga. Analizy pojedynczych komórek przeprowadzono w środowisku R z wykorzystaniem pakietu Seurat (wersja 5.0.1).

Analiza transkryptomiki przestrzennej
Próbkę transkryptomiki przestrzennej SP_BS3 wykorzystano do stworzenia obiektu ekspresji przestrzennej^25,26. Progi kontroli jakości ustalono na poziomie nFeature_Spatial ≥ 200 oraz percent.mt ≤ 30, co pozwoliło na pozostawienie 4 572 spotów do analizy. Dane zostały znormalizowane, wyselekcjonowano 3 000 genów o wysokiej zmienności, a następnie przeprowadzono skalowanie danych, analizę głównych składowych, konstrukcję grafu sąsiedztwa, klastrowanie spotów przestrzennych oraz redukcję wymiarowości. Wynik starzenia przestrzennego (spatial CellAge senescence score) obliczono po usunięciu PPARG ze zbioru genów, aby uniknąć korelacji kołowej. Skonstruowano i oceniono sygnatury śródbłonka, perycytów, makrofagów/monocytów, stromy związanej z guzem, złośliwego osteosarcoma oraz komórek podobnych do osteoklastów. Do oceny powiązań między ekspresją PPARG a wynikami przestrzennymi wykorzystano analizę korelacji Spearmana. Przeniesienie etykiet (label transfer) wykonano, przyjmując zbiór danych z pojedynczych komórek jako referencyjny, a zbiór danych przestrzennych jako zapytanie, w celu wywnioskowania przewidywanych wyników typów komórek dla każdego spotu^23,24. Analizy transkryptomiki przestrzennej przeprowadzono w środowisku R przy użyciu pakietu Seurat (wersja 5.0.1), a do wszystkich wartości P korelacji przestrzennych zastosowano korekcję FDR według metody Benjamini-Hochberga.

Zewnętrzna walidacja ekspresji w GSE36001
Zbiór danych GEO GSE36001 został wykorzystany wyłącznie jako niezależna kohorta walidacyjna ekspresji; ze względu na brak dostępnych danych dotyczących przeżywalności, nie został on użyty do walidacji prognostycznej^11,27. Zbiór danych obejmował 19 próbek kostniakomięsaka oraz sześć normalnych próbek kontrolnych. Do konwersji identyfikatorów sond na symbole genów wykorzystano adnotacje platformy GPL6102. W przypadku przypisania wielu sond do tego samego genu, zachowano sondę o najwyższej średniej ekspresji. Różnicową ekspresję między grupą nowotworową a normalną oceniono za pomocą pakietu limma²⁸. Analizy przeprowadzono w środowisku R z użyciem pakietu limma (wersja 3.56.2), a skorygowane wartości P obliczono metodą Benjamini-Hochberga.

Walidacja metodą qRT-PCR i Western blot
Walidację eksperymentalną przeprowadzono z wykorzystaniem ludzkiej linii komórkowej osteosarkoma 143B oraz ludzkich komórek osteoblastów. Komórki osteosarkoma hodowano w pożywce Dulbecco's modified Eagle's uzupełnionej o 10% płodowej surowicy bydlęcej i 1% penicyliny-streptomycyny w temperaturze 37 °C w wilgotnej atmosferze zawierającej 5% CO₂, a pasażowanie wykonywano przy użyciu 0,25% trypsyny-EDTA przy 80%–90% konfluencji. Ludzkie komórki osteoblastów utrzymywano w zalecanych warunkach hodowli. Potwierdzono, że wszystkie linie komórkowe są wolne od zanieczyszczenia mycoplasma. Do analiz qRT-PCR całkowite RNA wyekstrahowano przy użyciu odczynnika do ekstrakcji RNA opartego na fenolu i guanidynium, a stężenie i czystość RNA oceniono spektrofotometrycznie. Jeden mikrogram całkowitego RNA poddano odwrotnej transkrypcji przy użyciu odczynnika do odwrotnej transkrypcji zgodnie z zalecanym protokołem. qRT-PCR przeprowadzono z wykorzystaniem chemii opartej na fluorescencyjnym barwniku wiążącym DNA w następujących warunkach cykliczności: wstępna denaturacja w 95 °C przez 30 s, a następnie 40 cykli w 95 °C przez 5 s i 60 °C przez 30 s, wraz z analizą krzywych topnienia w celu potwierdzenia specyficzności amplifikacji. Każda reakcję przeprowadzono w trzech powtórzeniach technicznych i wykonano trzy niezależne eksperymenty biologiczne. Jako kontrolę wewnętrzną zastosowano GAPDH, a względną ekspresję PPARG obliczono metodą 2^-ΔΔCt29. Starter forward dla PPARG miał sekwencję 5'-CGAAGACATTCCATTCACAAGAACAG-3', a starter reverse 5'-AGATGCAGGCTCCACTTTGATTG-3'.

Wykonano analizę Western blotting w celu zbadania ekspresji białka PPARG. Komórki poddano lizie w buforze do radioimmunoprecypitacji uzupełnionym o inhibitory proteaz, a stężenia białek określono za pomocą testu z kwasem bicinchoninowym. Naniesiono równe ilości białka (30 µg (na ścieżkę) rozdzielono metodą elektroforezy w żelu poliakryloamidowym z dodatkiem 10% dodecylosiarczanu sodu, a następnie przeniesiono na membrany z difluorowodoru poliwinylydenu. Po blokowaniu 5% mlekiem beztłuszczowym przez 1 h w temperaturze pokojowej, membrany inkubowano przez noc w 4 °C z przeciwciałami pierwszorzędowymi przeciwko PPARG (1:1 000) i GAPDH (1:5 000), a następnie inkubowano z przeciwciałem wtórnym sprzężonym z peroksydazą chrzanową (1:5 000) przez 1 h w temperaturze pokojowej. Prążki białkowe wizualizowano za pomocą detekcji chemiluminescencyjnej; przeprowadzono trzy niezależne eksperymenty. Intensywność prążków oznaczono ilościowo przy użyciu oprogramowania do analizy obrazu.³⁰Różnice międzygrupowe analizowano przy użyciu obustronnego nieparzystego testu t Studenta t testy. Dane przedstawiono jako średnią ± odchylenie standardowe (SD) z trzech niezależnych eksperymentów oraz P < Przyjęto, że wartość 0,05 jest statystycznie istotna. Analizy statystyczne danych eksperymentalnych przeprowadzono przy użyciu oprogramowania do analizy statystycznej (wersja 9.0).

Analiza statystyczna
O ile nie zaznaczono inaczej, wszystkie analizy bioinformatyczne przeprowadzono w programie R (wersja 4.3.2). Za statystycznie istotne uznano dwustronne wartości P < 0,05. Korelacje oceniano przy użyciu współczynnika korelacji rang Spearmana (ρ). W stosownych przypadkach korektę dla wielokrotnych testów wykonano metodą FDR Benjamini-Hochberga. Dane eksperymentalne przedstawiono jako średnia ± SD i porównano za pomocą dwustronnych nieparzystych testów t Studenta. Eksperymentalne analizy statystyczne przeprowadzono przy użyciu oprogramowania do analizy statystycznej (wersja 9.0).

Wyniki

Zbiór GSE99671 zidentyfikował 105 genów związanych ze starzeniem komórkowym (CellAge) o różnym poziomie ekspresji
Zbiór GSE99671 obejmował 36 próbek z 18 par tkanek. Po odfiltrowaniu genów o niskiej liczbie odczytów zachowano 16 683 geny. Przy skorygowanym P < 0,05 stwierdzono różną ekspresję dla 2 248 genów. Przy zastosowaniu bardziej rygorystycznego progu skorygowanego P < 0,05 i |log2FC| ≥ 1, istotnych było 594 geny, w tym 102 geny o zwiększonej ekspresji i 492 geny o zmniejszonej ekspresji w guzach (Rycina 1A,B). Przecięcie zbioru 2 248 genów o różnej ekspresji z 866 genami CellAge pozwoliło wyodrębnić 105 genów związanych ze starzeniem komórkowym o różnej ekspresji (Rycina 1C). Gen PPARG wykazywał obniżoną ekspresję w GSE99671, przy log2FC = -0,644, P = 0,00451 i skorygowanej wartości P = 0,0309. W 13 z 18 par ekspresja PPARG była wyższa w tkance prawidłowej niż w tkance nowotworowej, z wartością P testu Wilcoxona dla par = 0,0294 (Rycina 1D).

Wielomodelowy screening prognostyczny wskazał PPARG jako główny gen kandydacki
W badaniu TARGET-OS uwzględniono 85 pacjentów, w tym 27 zdarzeń związanych ze zgonem. Zintegrowany screening 105 genów różnicowo wyrażonych, związanych ze starzeniem, z wykorzystaniem jednowymiarowej regresji Coxa, analizy Kaplan-Mierera, analizy ROC przeżycia, LASSO, powtarzalnego LASSO oraz modelowania losowego lasu przeżycia (random survival forest) pozwolił na wyłonienie 18 genów hubów kandydackich przed korektą kliniczną (Rysunek 2A). Gen PPARG był powiązany z przeżyciem całkowitym w jednowymiarowej analizie Coxa (HR = 0,603, 95% CI = 0,454–0,802, P = 0,000494, FDR = 0,0447), co wskazuje, że wyższa ekspresja PPARG wiązała się z niższym ryzykiem śmiertelności. Analiza Kaplan-Mierera porównująca grupy o wysokiej i niskiej ekspresji dała wynik P = 0,00784 (Rysunek 2B). Zależne od czasu wartości AUC w 1., 3. i 5. roku wyniosły odpowiednio 0,603, 0,760 i 0,776 (Rysunek 2C). PPARG wykazał częstotliwość selekcji w powtarzalnym LASSO na poziomie 0,920 oraz istotność w losowym lesie przeżycia wynoszącą 0,0398 (Rysunek 2D–F). Początkowy wynik hub dla PPARG wyniósł 6/6, ponieważ spełnił on wszystkie sześć z góry określonych kryteriów screeningu.

Korekta kliniczna potwierdziła wartość prognostyczną PPARG
Po uwzględnieniu klinicznych zmiennych towarzyszących, PPARG pozostał istotnie powiązany z przeżyciem całkowitym (skorygowane HR = 0,224, 95% CI = 0,085–0,589, P = 0,00241; Rysunek 2G). Model oparty wyłącznie na danych klinicznych uzyskał C-index wynoszący 0,707 oraz AIC wynoszący 89,921 (Tabela uzupełniająca 1). Dodanie PPARG zwiększyło C-index do 0,829, obniżyło AIC do 78,466 i istotnie poprawiło dopasowanie modelu zgodnie z testem ilorazu wiarygodności (P = 0,000244; Tabela uzupełniająca 2, Rysunek 2H,I). Analiza wrażliwości po usunięciu operacji radykalnej zachowała ochronny charakter powiązania z PPARG (HR = 0,249, P = 0,00185; Tabela uzupełniająca 3, Rysunek 2J). PPARG osiągnął najwyższy zintegrowany wynik kliniczny równy 13 (początkowy wynik hub 6 plus siedem punktów za integrację kliniczną) i wykazał końcowe wartości AUC dla 3 i 5 lat wynoszące odpowiednio 0,770 oraz 0,813 (Rysunek 2K).

Cechy funkcjonalne i immunologiczne mikrośrodowiska związane z PPARG
Analiza GSEA porównująca grupy z wysoką i niską ekspresją PPARG wykazała wzbogacenie zestawów: Nemeth Inflammatory Response LPS Up, Burton Adipogenesis 5, Burton Adipogenesis 6, Krieg KDM3A Targets Not Hypoxia, Reactome: Transcriptional Regulation By TP53, Fulcher Inflammatory Response Lectin Vs LPS Dn, Hollmann Apoptosis Via CD40 Dn, Zhou Inflammatory Response Live Dn, WP: Fatty Acids And Lipoproteins Transport In Hepatocytes, Sweet Lung Cancer KRAS Up, KEGG Medicus Pathogen HIV Tat To TLR2/4 NF-kB Signaling Pathway oraz Reactome: Fatty Acids (Rysunek 3A). W danych bulk TARGET-OS PPARG nie korelował istotnie z ogólnym wskaźnikiem starzenia CellAge (Spearman ρ = 0.022, P = 0.837; Rysunek 3B), ale wykazał korelację z kilkoma pojedynczymi genami CellAge (Rysunek 3C). Analiza mikrośrodowiska immunologicznego wykazała dodatnie korelacje między PPARG a makrofagami (ρ = 0.485, FDR = 2.7 × 10^-5), limfocytami T CD8 (ρ = 0.410, FDR = 5.88 × 10^-4), sygnaturą osteoklastopodobną (ρ = 0.383, FDR = 0.00120), neutrofilami (ρ = 0.376, FDR = 0.00120) oraz komórkami dendrytycznymi (ρ = 0.370, FDR = 0.00123). Guzy z wysoką ekspresją PPARG wykazały wyższe sygnatury komórek osteoklastopodobnych, makrofagów, limfocytów T CD8, komórek dendrytycznych, monocytów, neutrofili, komórek NK oraz komórek śródbłonka po korekcie FDR (Tabela uzupełniająca 4, Rysunek 3D,E).

Transkryptomika pojedynczych komórek zlokalizowała PPARG w przedziałach naczyniowych i mikrośrodowiskowych
Zbiór danych z pojedynczych komórek obejmował 68 336 komórek i 32 297 genów. Ekspresja PPARG różniła się znacząco między typami komórek (Rysunek 4A). Najwyższa średnia ekspresja została zaobserwowana w komórkach śródbłonkowych (średnia ekspresja = 0,540; odsetek pozytywnych = 44,33%), perycytach (średnia ekspresja = 0,439; odsetek pozytywnych = 41,61%), makrofagach/monocytach (średnia ekspresja = 0,363; odsetek pozytywnych = 31,62%) oraz komórkach zrazowych związanych z guzem (średnia ekspresja = 0,361; odsetek pozytywnych = 44,10%; Rysunek 4B–E). Podgrupa złośliwych komórek osteosarcoma wykazywała ekspresję PPARG (średnia ekspresja = 0,163; odsetek pozytywnych = 16,70%), jednak ekspresja w tych komórkach nie była znacząco wyższa niż w innych komórkach (FDR = 0,151). Wyniki te sugerują, że ekspresja PPARG w osteosarcoma odzwierciedlała przede wszystkim stany naczyniowe, mieloidalne i mikrośrodowiska zrazowego, a nie była ograniczona do komórek złośliwych (Rysunek 4F).

Transkryptomika przestrzenna powiązała PPARG z przestrzennymi stanami związanymi z senescencją oraz niszami naczyniowymi
Po kontroli jakości zachowano 4 572 przestrzennych punktów (spots) SP_BS3, które pogrupowano w siedem klastrów przestrzennych (Rysunek 5A). Rysunek 5B przedstawia rozkład przestrzenny nFeature_Spatial (liczba wykrytych genów na punkt). Oddzielnie, z 866 genów CellAge, 845 dopasowano w macierzy ekspresji przestrzennej (97,58%). PPARG wykazał ogniskową ekspresję przestrzenną (Rysunek 5C). Przestrzenny wskaźnik senescencji CellAge, obliczony po usunięciu PPARG, wykazał słabą, ale statystycznie istotną dodatnią korelację z ekspresją PPARG (ρ = 0,0692, P = 3,0 × 10^-6, FDR = 1,9 × 10^-5; Rysunek 5D). Ocena nisz przestrzennych wykazała dodatnie korelacje między PPARG a wskaźnikiem śródbłonka (ρ = 0,0433, FDR = 0,00592) i wskaźnikiem perycytów (ρ = 0,0367, FDR = 0,0181), podczas gdy PPARG korelował ujemnie ze wskaźnikiem złośliwego kostniakomiąka (ρ = -0,0592, FDR = 0,000219) oraz wskaźnikiem zrębu guza (ρ = -0,0531, FDR = 0,000774; Tabela uzupełniająca 5, Rysunek 5E). Analiza transferu etykiet (label transfer) wykazała analogicznie dodatnie korelacje z predykcyjnym wskaźnikiem śródbłonka (ρ = 0,0507, FDR = 0,00120) i predykcyjnym wskaźnikiem perycytów (ρ = 0,0394, FDR = 0,0123), wraz z ujemną korelacją z predykcyjnym wskaźnikiem komórek złośliwego kostniakomiąka (ρ = -0,0699, FDR = 1,1 × 10^-5; Rysunek 5F–H).

Ekspresja zewnętrzna i walidacja eksperymentalna potwierdziły obniżenie poziomu PPARG
Zbiór GSE36001 obejmował 19 próbek kostniakomięsaka oraz sześć prawidłowych próbek kontrolnych. PPARG był znacząco obniżony w kostniakomięsaku (logFC = -1,429, P = 0,00730, skorygowana wartość P = 0,0435; Rycina 6A). W walidacji komórkowej ekspresja mRNA PPARG była znacząco niższa w komórkach kostniakomięsaka 143B niż w ludzkich komórkach osteoblastów, co wykazano za pomocą qRT-PCR (P < 0,001; Rycina 6B). Ekspresja białka PPARG była również znacząco zredukowana w komórkach 143B, co potwierdzono metodą Western blotting (P < 0,01; Rycina 6C,D). Wyniki uzyskane z zewnętrznej kohorty oraz na poziomie mRNA i białka spójnie potwierdziły obniżoną ekspresję PPARG w kostniakomięsaku. Zbiór GSE36001 nie zawierał danych dotyczących przeżywalności, w związku z czym zapewnił jedynie zewnętrzną walidację ekspresji, a nie niezależną walidację prognostyczną.

DOSTĘPNOŚĆ DANYCH:
Wszystkie zbiory danych wykorzystane w niniejszym badaniu są publicznie dostępne. Dane GSE99671 oraz GSE36001 pobrano z bazy Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi/acc=GSE99671; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi/acc=GSE36001). Transkryptomiczne i kliniczne dane TARGET-OS pobrano z UCSC Xena (https://xena.ucsc.edu/). Geny związane z procesem starzenia uzyskano z CellAge: The Database of Cell Senescence Genes, będącej częścią Human Ageing Genomic Resources (https://genomics.senescence.info/cells/). Zbiory danych transkryptomicznych pojedynczych komórek i danych przestrzennych dla ludzkiego kostniakomięsaka pobrano z opublikowanego atlasu oraz powiązanego z nim repozytorium GitHub (https://github.com/zhengxj1/A-Single-Cell-and-Spatially-Resolved-Atlas-of-Human-Osteosarcomas). Zbiory genów do analizy wzbogacenia pobrano z MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/). Przetworzone dane wygenerowane w tym badaniu oraz skrypty analityczne wykorzystane do odwzorowania raportowanych wyników zostały zebrane i przesłane jako Plik uzupełniający 1.

Rysunek 1: Identyfikacja genów różnicowo wyrażonych oraz potencjalnych genów związanych ze starzeniem się komórek wyodrębnionych z bazy CellAge w kostniakomiaku. (A) Wykres wulkaniczny przedstawiający geny różnicowo wyrażone pomiędzy tkankami kostniakomiaka a dopasowanymi kontrolnymi tkankami nienowotworowymi w zbiorze danych GSE99671. Istotnie zwiększone i zmniejszone poziomy ekspresji genów zostały wyróżnione zgodnie z predefiniowanymi kryteriami odcięcia. (B) Mapa ciepła przedstawiająca wzorce ekspresji reprezentatywnych genów różnicowo wyrażonych w próbkach kostniakomiaka i dopasowanych próbkach kontrolnych w GSE99671. (C) Diagram Venna przedstawiający część wspólną genów różnicowo wyrażonych w GSE99671 i genów związanych ze starzeniem się z bazy CellAge. (D) Porównanie ekspresji PPARG w parach tkanek kostniakomiaka i dopasowanych kontrolnych tkanek nienowotworowych w GSE99671. Kliknij tutaj, aby wyświetlić większą wersję tego rysunku.

Rycina 2: Analizy przeżycia z wykorzystaniem uczenia maszynowego i korekty klinicznej identyfikują PPARG jako kluczowy prognostyczny gen hub związany z senescencją w osteosarcoma. (A) Wykres leśny przedstawiający wyniki jednowymiarowej regresji Coxa dla kandydujących genów związanych z senescencją w kohorcie TARGET-OS. (B) Krzywa przeżycia Kaplana-Meiera porównująca przeżycie całkowite pacjentów z wysoką i niską ekspresją PPARG. (C) Zależne od czasu krzywe ROC oceniające wydajność predykcyjną PPARG w odniesieniu do przeżycia całkowitego. (D) Krzywa walidacji krzyżowej regresji LASSO Cox dla selekcji kandydujących genów prognostycznych. (E) Powtórna analiza stabilności LASSO pokazująca częstotliwość wyboru lambda.min w 300 powtórzeniach pięciokrotnych. (F) Analiza losowego lasu przeżycia (Random Survival Forest) przedstawiająca wyniki istotności zmiennych z 1 000 drzew. (G) Wykres leśny przedstawiający wyniki klinicznie skorygowanej regresji Coxa dla kandydujących genów hub. (H) Zmiany AIC po dodaniu poszczególnych genów hub do modelu klinicznego. (I) Poprawa indeksu C po dodaniu poszczególnych genów hub do modelu klinicznego. (J) Końcowy ranking zintegrowanego wyniku klinicznego (zakres 0–13) dla kandydującego genu hub. (K) Zależne od czasu krzywe ROC dla PPARG w podzbiorze analizy klinicznej obejmującym 40 pacjentów, pokazujące AUC dla 3 i 5 lat; AUC dla 1 roku nie była możliwa do oszacowania. Proszę kliknąć tutaj, aby zobaczyć powiększoną wersję tej ryciny.

Rysunek 3: Analiza wzbogacenia funkcjonalnego oraz mikrośrodowiska immunologicznego powiązanego z PPARG. (A) Wykres bąbelkowy GSEA porównujący grupy PPARG-high i PPARG-low, przedstawiający Nemeth Inflammatory Response LPS Up, Burton Adipogenesis 5, Burton Adipogenesis 6, Krieg KDM3A Targets Not Hypoxia, Reactome: Transcriptional Regulation By TP53, Fulcher Inflammatory Response Lectin Vs LPS Dn, Hollmann Apoptosis Via CD40 Dn, Zhou Inflammatory Response Live Dn, WP: Fatty Acids And Lipoproteins Transport In Hepatocytes, Sweet Lung Cancer KRAS Up, KEGG Medicus Pathogen HIV Tat To TLR2/4 NF-kB Signaling Pathway oraz Reactome: Fatty Acids. (B) Korelacja między PPARG a całkowitym wynikiem starzenia CellAge w zbiorczych danych TARGET-OS. (C) Korelacje między PPARG a poszczególnymi genami CellAge. (D) Korelacje między PPARG a sygnaturami mikrośrodowiska immunologicznego. (E) Różnice w wynikach mikrośrodowiska między grupami PPARG-high i PPARG-low. Korelacje oceniano przy użyciu współczynnika korelacji rang Spearmana (ρ), a skorygowane wartości P obliczono metodą Benjamini-Hochberga. Skróty: GSEA = analiza wzbogacenia zbiorów genów; NF-κB = jądrowy czynnik kappa-B; JAK-STAT = kinaza Janus i przekaźnik sygnału i aktywator transkrypcji; IL-12 = interleukina-12; FDR = wskaźnik fałszywych odkryć. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

Rysunek 4: Lokalizacja PPARG w przedziałach komórkowych w danych transkryptomicznych pojedynczych komórek kostniakomięsaka. (A) Wizualizacja UMAP głównych typów komórek w ludzkim zbiorze danych transkryptomicznych pojedynczych komórek kostniakomięsaka po uproszczonej adnotacji ręcznej. (B) FeaturePlot pokazujący globalny rozkład ekspresji PPARG w pojedynczych komórkach. (C) DotPlot pokazujący ekspresję PPARG w głównych typach komórek. (D) Wykres skrzypcowy pokazujący poziomy ekspresji PPARG w różnych typach komórek. (E) Wykres słupkowy pokazujący proporcję komórek PPARG-dodatnich w każdym głównym typie komórek. (F) Wizualizacja UMAP pokazująca ekspresję PPARG w złośliwych komórkach kostniakomięsaka. Skrót: UMAP = uniform manifold approximation and projection. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

Rysunek 5: Przestrzenna lokalizacja transkrypcyjna PPARG i cech przestrzennych związanych ze starzeniem w osteosarcoma. (A) Rozkład przestrzenny klastrów zdefiniowanych na podstawie transkryptomu w przekroju przestrzennym transkryptomiki osteosarcoma SP_BS3. (B) Rozkład przestrzenny wykrytych genów na punkt, przedstawiony przez nFeature_Spatial. (C) Przestrzenny wzorzec ekspresji PPARG w punktach SP_BS3. (D) Rozkład przestrzenny wskaźnika starzenia wyznaczonego na podstawie CellAge. (E) Analiza korelacji między ekspresją PPARG a przestrzennym wskaźnikiem starzenia wyznaczonym na podstawie CellAge lub wskaźnikami ekologicznej niszy komórkowej. (F) Mapa predykcji z transferu etykiet pokazująca dominujący typ komórek pochodzący z analizy single-cell dla każdego punktu przestrzennego. (G) Analiza korelacji między ekspresją PPARG, wskaźnikiem starzenia wyznaczonym na podstawie CellAge a wynikami predykcji typu komórek z transferu etykiet. (H) Rozkład przestrzenny punktów z wysokim i niskim poziomem PPARG. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

Rycina 6: Zewnętrzna analiza ekspresji i walidacja eksperymentalna obniżenia poziomu PPARG w osteosarcoma. (A) Wykres pudełkowy przedstawiający poziomy ekspresji PPARG w próbkach z osteosarcoma (n = 19) oraz w prawidłowych próbkach kontrolnych (n = 6) w zbiorze danych GSE36001. (B) Analiza qRT-PCR ekspresji mRNA PPARG w ludzkich komórkach osteosarcoma 143B oraz w ludzkich komórkach kontrolnych osteoblastów. (C) Reprezentatywny western blot wykazujący ekspresję białek PPARG i GAPDH w ludzkich komórkach kontrolnych osteoblastów oraz komórkach osteosarcoma 143B. GAPDH posłużył jako kontrola ładowania. (D) Kwantyfikacja densytometryczna prążków western blot przedstawiająca względne poziomy białka PPARG znormalizowane do GAPDH. Na wykresach (B) i (D) dane przedstawiono jako średnia ± SD z trzech niezależnych eksperymentów. P < 0,01 oraz P < 0,001 w stosunku do grupy kontrolnej ludzkich osteoblastów, określone za pomocą dwustronnego nieparzystego testu t Studenta. Skróty: qRT-PCR = ilościowa reakcja polimerazy w łańcuchu z odwrotną transkrypcją; GAPDH = gliceraldehyd-3-fosforan dehydrogenaza; SD = odchylenie standardowe. Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Tabela uzupełniająca 1: Wydajność modelu Coxa opartego wyłącznie na danych klinicznych w kohorcie TARGET-OS. Współczynnik C-index, AIC oraz podsumowanie modelu Coxa zbudowanego z wykorzystaniem wyłącznie zmiennych klinicznych, w tym płci, wieku, stanu choroby w momencie diagnozy, lokalizacji guza pierwotnego, specyficznego regionu guza oraz stanu operacji radykalnej. Skrót: AIC = kryterium informacyjne Akaike. Aby pobrać ten plik, kliknij tutaj.

Tabela uzupełniająca 2: Porównanie modeli Coxa opartych wyłącznie na danych klinicznych oraz modeli łączących dane kliniczne z genetycznymi. Wyniki porównania modeli po dodaniu poszczególnych kandydowanych genów hub do modelu klinicznego, obejmujące C-index, AIC, statystyki testu ilorazu wiarygodności oraz wskaźniki poprawy modelu. Skrót: AIC = kryterium informacyjne Akaike. Aby pobrać ten plik, kliknij tutaj.

Tabela uzupełniająca 3: Analiza wrażliwości po usunięciu zmiennej dotyczącej operacji radykalnej. Wyniki analizy wrażliwości regresji Coxa oceniające, czy powiązania prognostyczne kandydatów na geny hubowe, w szczególności PPARG, pozostały stabilne po wykluczeniu zmiennej dotyczącej operacji radykalnej z dostosowanego modelu klinicznego. Kliknij tutaj, aby pobrać ten plik.

Tabela uzupełniająca 4: Sygnatury mikrośrodowiska immunologicznego i zrębowego powiązane z PPARG w TARGET-OS. Wyniki korelacji i porównań grupowych pomiędzy ekspresją PPARG a sygnaturami ssGSEA związanymi z układem odpornościowym, zrębem, naczyniami, stanem zapalnym oraz SASP, obejmujące współczynniki korelacji Spearmana, wartości P, skorygowane wartości P oraz porównania grupy PPARG-high względem PPARG-low. Skróty: SASP = fenotyp wydzielniczy związany ze starzeniem; ssGSEA = analiza wzbogacenia zestawów genów dla pojedynczej próbki. Kliknij tutaj, aby pobrać ten plik.

Tabela uzupełniająca 5: Analiza korelacji przestrzennej transkryptomiki PPARG w SP_BS3. Wyniki korelacji między ekspresją PPARG a wynikami starzenia opracowanymi za pomocą CellAge, wynikami nisz ekologicznych komórek oraz wynikami przewidywania typów komórek na podstawie transferu etykiet w przekroju transkryptomicznym przestrennym osteosarcoma SP_BS3. Kliknij tutaj, aby pobrać ten plik.

Dyskusja

W niniejszym badaniu, poprzez integrację analizy różnicowej ekspresji, przecięcia zbiorów genów CellAge, modelowania przeżywalności TARGET-OS, korekty klinicznej, lokalizacji multiomicznej oraz walidacji eksperymentalnej, wskazano PPARG jako potencjalny gen prognostyczny związany ze starzeniem komórkowym w kostniakomiaku. Przyjęte ramy analityczne były zgodne z obecnym naciskiem na badania nad kostniakomiakiem oparte na informacjach molekularnych oraz kuratorowanymi zasobami genów starzenia w celu interpretacji biologii związanej ze starzeniem komórkowym^14,31. Ekspresja PPARG była obniżona w kostniakomiaku w porównaniu z tkanką prawidłową, a niższa ekspresja PPARG wiązała się z gorszą przeżywalnością całkowitą w kohorcie TARGET-OS. Wyniki te sugerują, że PPARG nie tylko ulega zmianom transkrypcyjnym w kostniakomiaku, ale może również nieść klinicznie istotne informacje prognostyczne. Jednakże, ponieważ wybór kandydata i ocena modelu przeprowadzono w tej samej kohorcie TARGET-OS (85 pacjentów, 27 zgonów), zaobserwowana poprawa wydajności modelu była podatna na nadmierny optymizm i przeuczenie; w związku z tym wartość prognostyczną PPARG należy traktować jako generującą hipotezę do czasu jej walidacji w niezależnej kohorcie przeżywalności w kostniakomiaku. Niemniej jednak biologię PPARG należy interpretować w sposób zależny od kontekstu, ponieważ eksperymentalne badania nad kostniakomiakiem zgłaszały zarówno przeciwnowotworowe efekty modulacji receptora jądrowego PPAR-gamma, jak i związane z PPARG programy osteoklastyczne, które mogą wspierać progresję choroby^32,33,34.

Ważnym niuansem było to, że PPARG nie powinno być interpretowane jako prosty substytut całkowitego wyniku CellAge. W danych bulk z bazy TARGET-OS PPARG nie korelował istotnie z globalnym wynikiem starzenia CellAge, podczas gdy w transkryptomice przestrzennej PPARG wykazało słabą, ale istotną korelację z wynikiem CellAge obliczonym po wyłączeniu PPARG. Różnica ta mogła wynikać z efektów składu komórkowego w danych bulk, wielofunkcyjnej natury zestawów genów starzenia oraz lokalnego wzbogacenia nisz mikrośrodowiskowych w plamkach przestrzennych. Badania konsensusowe i transkryptomiczne podkreśliły, że starzenie komórkowe jest heterogeniczne, dynamiczne i zależne od typu komórki, czynnika stresowego oraz kontekstu tkankowego, podczas gdy programy SASP mogą wywierać przeciwstawne efekty podczas progresji nowotworu^25,35,36,37. PPARG zdefiniowano zatem konserwatywnie jako gen prognostyczny związany ze starzeniem wywiedziony z CellAge, a nie jako udowodniony czynnik napędzający starzenie. W związku z tym oznaczenie PPARG jako związane ze starzeniem odzwierciedlało przynależność do zestawu genów CellAge, a nie wykazane zaangażowanie mechanistyczne w starzenie, i PPARG nie powinno być stosowane jako ilościowy substytut ogólnej aktywności starzenia.

W szerszym ujęciu wyniki przestrzenne te zilustrowały rozróżnienie między istotnością statystyczną a istotnością biologiczną. Przy 4 572 punktach przestrzennych nawet bardzo słabe korelacje mogły przekraczać konwencjonalne progi istotności; na przykład korelacja między PPARG a przestrzennym wskaźnikiem starzenia CellAge (ρ = 0,0692) wyjaśniała jedynie około 0,48% wariancji, a mimo to osiągnęła wartość P równą 3,0 × 10^-6, ponieważ duża liczba punktów zapewniła znaczną moc statystyczną. Powiązania na poziomie punktów należy zatem traktować jako sygnały wykrywalne statystycznie, lecz biologicznie umiarkowane i służące do generowania hipotez, a interpretację biologiczną powinny kierować wielkości efektów, a nie same wartości P. Istotności statystycznej w dużych zbiorach danych na poziomie punktów nie należy utożsamiać z silnymi efektami biologicznymi.

Analizy pojedynczych komórek oraz analizy przestrzenne dostarczyły kontekstu biologicznego dla PPARG. PPARG był wzbogacony w komórkach śródbłonka, perycytach, makrofagach/monocytach oraz w komórkach zrębowych związanych z nowotworem, a przestrzennie wiązał się z wynikami dla śródbłonka i perycytów. Wyniki te sugerują, że informacje prognostyczne niesione przez PPARG mogą być powiązane z komponentami naczyniowymi i mikrośrodowiskowymi w kostniakomiaku. Interpretacja ta była zgodna z niedawnymi atlasami pojedynczych komórek i analizami przestrzennymi kostniakomiaka oraz z szerszymi dowodami na to, że nisze immunologiczne, naczyniowe i zrębowe kształtują heterogenność guza, odpowiedź na leczenie i progresję^25,38,39. Ponieważ złośliwe komórki kostniakomiaka wykazywały ekspresję PPARG tylko w podzbiorze komórek i nie wykazywały znacząco wyższej ekspresji niż inne komórki, interpretacja ograniczona wyłącznie do komórek nowotworowych byłaby niepełna. Z drugiej strony, obniżenie poziomu ekspresji PPARG w analizach bulk zaobserwowane w GSE99671 i GSE36001 mogło częściowo odzwierciedlać różnice w składzie komórek zrębowych, naczyniowych, szpiku kostnego, adipogennych lub immunologicznych pomiędzy tkankami nowotworowymi a prawidłowymi, a nie rzeczywisty spadek ekspresji w złośliwych komórkach kostniakomiaka. Ponieważ w analizach bulk nie skorygowano wprost czystości guza ani liczebności poszczególnych typów komórek, możliwości tej nie można było wykluczyć, co uzasadniało przeprowadzenie ukierunkowanych badań.

Analizy funkcjonalne i mikrośrodowiska immunologicznego były zgodne z tą interpretacją. PPARG korelował z sygnaturami makrofagów, limfocytów T CD8, komórek dendrytycznych, neutrofili, komórek NK i komórek śródbłonka, a analiza GSEA wyróżniła szlaki związane ze stanem zapalnym, NF-κB, sygnalizacją JAK-STAT/IL-12, regulacją TP53 oraz punktami kontrolnymi uszkodzeń DNA. Wspólnie wyniki te sugerowały, że PPARG może wyznaczać złożony stan mikrośrodowiskowy obejmujący stres związany ze starzeniem, infiltrację immunologiczną oraz komponenty naczyniowe i zrębowe. Interpretacja ta była biologicznie prawdopodobna, ponieważ PPAR-gamma pełni uznane role w hamowaniu aktywacji zapalnej makrofagów/monocytów, w tym wywierając wpływ na programy transkrypcyjne związane z AP-1, STAT i NF-κB, a mikrośrodowisko immunologiczne osteosarkoma zawiera elementy mieloidalne, limfoidalne i naczyniowe o funkcjach zarówno pro- i antynowotworowych^38,40,41. Te obserwacje z poziomu pojedynczych komórek oraz analizy przestrzenne miały charakter opisowy i służyły generowaniu hipotez; same w sobie nie ustanowiły mechanizmów niszy naczyniowej, programów starzenia ani szlaków prognostycznych.

Należy wskazać na kilka ograniczeń. Po pierwsze, główna analiza prognostyczna opierała się na retrospektywnej publicznej kohorcie TARGET-OS, z ograniczoną wielkością próby i liczbą zdarzeń; w związku z tym wartość prognostyczna PPARG powinna zostać zwalidowana w niezależnych kohortach zgodnie z przyjętymi zasadami raportowania i walidacji markerów nowotworowych⁴². W szczególności przesiew kandydatów przeprowadzono u 85 pacjentów z grupy TARGET-OS, wśród których odnotowano 27 zgonów, natomiast ocena modelu skorygowanego klinicznie wykorzystała nakładający się podzbiór 40 pacjentów z 13 zgonami. Ponieważ obie analizy pochodziły z tej samej kohorty źródłowej, raportowana poprawa wskaźnika C-index (z 0,707 do 0,829) oraz redukcja AIC są prawdopodobnie zbyt optymistyczne. Brak było niezależnej kohorty przeżywalności do zewnętrznej walidacji wartości prognostycznej PPARG; zbiór GSE36001 został wykorzystany wyłącznie do walidacji ekspresji w tkance nowotworowej w porównaniu z prawidłową. Ponadto zintegrowany wynik hub score był eksploracyjną wewnętrzną heurystyką rankingową, a nie zwalidowanym narzędziem prognostycznym. Po drugie, analiza transkryptomiki przestrzennej opierała się na pojedynczej próbce SP_BS3; choć korelacje były istotne statystycznie, wielkości efektów były małe i wymagają walidacji na dodatkowych próbkach przestrzennych. W szczególności zwiększenie liczby próbek transkryptomiki przestrzennej od niezależnych pacjentów będzie niezbędne do uzyskania bardziej wiarygodnych szacunków tych słabych powiązań, a uzasadnione jest przeprowadzenie badań opartych na większych kohortach przestrzennych. Transkryptomika przestrzenna dostarcza cennego kontekstu molekularnego in situ, jednak interpretacja pozostaje pod wpływem rozdzielczości platformy, strategii pobierania próbek, jakości tkanki oraz wyborów dotyczących integracji obliczeniowej⁴³. Co więcej, analiza przestrzenna opierała się na wstępnie przetworzonym publicznym obiekcie jednokomórkowym z uproszczonymi adnotacjami oraz na obliczeniowym punktowaniu i transferze etykiet; biorąc pod uwagę minimalne wielkości efektów, dane te wspierały opisowe stwierdzenia dotyczące lokalizacji, a nie twierdzenia mechanistyczne dotyczące nisz naczyniowych lub programów starzenia. Po trzecie, walidacja eksperymentalna opierała się na linii komórkowej osteosarcoma 143B oraz ludzkich osteoblastach; niezbędne są dodatkowe linie komórkowe osteosarcoma oraz próbki kliniczne, a walidacja w jednej linii komórkowej nie może potwierdzić specyficzności typu komórkowego, klinicznej istotności prognostycznej ani biologii starzenia. Po czwarte, badanie to wykazało korelację, a nie przyczynowość. Funkcjonalna perturbacja PPARG w komórkach osteosarcoma i modelach mikrośrodowiskowych będzie potrzebna, aby określić, czy PPARG bezpośrednio reguluje programy związane ze starzeniem, nisze naczyniowe lub progresję nowotworu. Po piąte, powiązanie ze starzeniem zdefiniowano poprzez nakładanie się ze zbiorem genów CellAge, a PPARG nie korelował z ogólnym wynikiem starzenia CellAge w TARGET-OS. Ponadto porównania tkanki nowotworowej z prawidłową w analizie bulk nie zostały skorygowane o czystość nowotworu ani skład typów komórek, więc obserwowana downregulation może częściowo odzwierciedlać różnice w składzie mikrośrodowiska, a nie zmiany wewnątrzkomórkowe w komórkach złośliwych.

PPARG jest genem pochodzącym z bazy CellAge, którego niższa ekspresja wiąże się z niekorzystnym przeżyciem całkowitym w kohorcie TARGET-OS. Jego obniżona ekspresja w porównaniach masowych tkanek kostniakomięsaka, wraz z wzbogaceniem w przedziałach naczyniowych, mieloidalnych i zrębowych, sugeruje, że całkowite poziomy PPARG mogą częściowo odzwierciedlać skład komórkowy mikrośrodowiska, a nie ekspresję wewnątrzpochodną komórek złośliwych. Wyniki analizy pojedynczych komórek oraz analizy przestrzennej mają charakter opisowy, a powiązanie prognostyczne z przeżyciem nie zostało niezależnie potwierdzone i wymaga walidacji w zewnętrznej kohorcie z uwzględnieniem wyników przeżycia. Wyniki te wskazują na PPARG jako potencjalny biomarker do oceny prognostycznej w kostniakomięsaku oraz do badań nad mikrośrodowiskiem związanym z senescencją, w oczekiwaniu na walidację zewnętrzną i badania funkcjonalne. W transkryptomice przestrzennej korelacje między PPARG a wynikami nisz naczyniowych lub związanych z senescencją charakteryzowały się słabą wielkością efektu, choć były istotne statystycznie ze względu na dużą liczbę punktów przestrzennych, i powinny być zatem interpretowane z ostrożnością.

Materiały

Lista materiałów użytych w tym artykule

Nazwa	Firma	Numer katalogowy	Komentarze
Przeciwciało pierwszorzędowe anty-GAPDH	Proteintech Group, Wuhan, China	10494-1-AP	Przeciwciało pierwszorzędowe stosowane do wykrywania GAPDH jako kontroli załadunku w technice western blotting; rozcieńczenie 1:5,000.
Przeciwciało pierwszorzędowe anty-PPARG	Proteintech Group, Wuhan, China	16643-1-AP	Przeciwciało pierwszorzędowe stosowane do wykrywania białka PPARG w technice western blotting; rozcieńczenie 1:1,000.
Zestaw do oznaczania białka BCA	Beijing Solarbio Science & Technology Co., Ltd., Beijing, China	PC0020	Test kolorymetryczny stosowany do określania całkowitego stężenia białka przed elektroforezą.
DESeq2	Bioconductor	Version 1.40.2	Pakiet R stosowany do analizy różnicowej ekspresji genów w danych transkryptomicznych opartych na liczbie odczytów.
Podłoże DMEM (Dulbecco's modified Eagle's medium)	Beijing Solarbio Science & Technology Co., Ltd., Beijing, China	11995	Podstawowe podłoże hodowlane stosowane do utrzymania komórek osteosarkoma 143B.
Odczynnik do detekcji ECL	Beijing Solarbio Science & Technology Co., Ltd., Beijing, China	PE0010	Substrat chemiluminescencyjny stosowany do wykrywania przeciwciał sprzężonych z HRP w technice western blotting.
Płodowa surowica bydlęca (FBS)	Zhejiang Tianhang Biotechnology Co., Ltd. (Sijiqing), Huzhou, Zhejiang, China	11011-8611	Surowica dodawana do podłoża hodowlanego w celu wspierania wzrostu i żywotności komórek.
glmnet	CRAN	Version 4.1-8	Pakiet R stosowany do analiz regresji z regularyzacją, w tym modelowania LASSO i elastic-net.
GraphPad Prism	GraphPad Software, San Diego, CA, USA	Version 9.0	Oprogramowanie stosowane do analizy statystycznej, generowania wykresów i wizualizacji danych eksperymentalnych.
Sprzężone z HRP przeciwciało wtórne	Proteintech Group, Wuhan, China	SA00001-2	Przeciwciało wtórne sprzężone z peroksydazą chrzanową stosowane do detekcji w western blot; rozcieńczenie 1:5,000.
Ludzkie komórki osteoblastów	Cell Applications, Inc., San Diego, CA, USA	406-05A	Pierwotne ludzkie komórki osteoblastów stosowane jako niezłośliwy typ komórek porównawczych/kontrolnych.
Ludzka linia komórkowa osteosarkoma 143B	American Type Culture Collection (ATCC), Manassas, VA, USA	CRL-8303	Ludzka linia komórkowa osteosarkoma stosowana do eksperymentów walidacyjnych in vitro oraz analiz molekularnych.
ImageJ	National Institutes of Health (NIH), USA	Version 1.53	Oprogramowanie do analizy obrazu stosowane do ilościowej analizy zdjęć eksperymentalnych.
limma	Bioconductor	Version 3.56.2	Pakiet R stosowany do analizy różnicowej ekspresji i testowania statystycznego opartego na modelach liniowych.
Penicylina-streptomycyna	Beijing Solarbio Science & Technology Co., Ltd., Beijing, China	P1400	Suplement antybiotykowy stosowany w podłożu hodowlanym w celu ograniczenia kontaminacji bakteryjnej.
Startery PPARG i GAPDH	Sangon Biotech (Shanghai) Co., Ltd., Shanghai, China	Custom synthesis; sequences provided in the Methods	Dedykowane startery oligonukleotydowe stosowane do analizy qPCR ekspresji PPARG i GAPDH.
Membrany PVDF, 0.45 µm	Beijing Solarbio Science & Technology Co., Ltd., Beijing, China	YA1701	Membrany stosowane do transferu białek podczas western blotting.
Oprogramowanie statystyczne R	R Foundation for Statistical Computing, Vienna, Austria	Version 4.3.2	Środowisko obliczeń statystycznych stosowane do analiz bioinformatycznych, konstrukcji modeli i wizualizacji.
randomForestSRC	CRAN	Version 3.2.2	Pakiet R stosowany do modelowania za pomocą random survival forest i analizy istotności cech.
Zestaw do odwrotnej transkrypcji	Beyotime Biotech Inc., Shanghai, China	D7168M	Stosowany do syntezy komplementarnego DNA (cDNA) z wyizolowanego RNA przed ilościowym PCR.
Bufor do lizy RIPA	Beijing Solarbio Science & Technology Co., Ltd., Beijing, China	R0010	Bufor do ekstrakcji białek stosowany do lizy komórek w analizie western blot.
Seurat	Satija Laboratory	Version 5.0.1	Pakiet R stosowany do przetwarzania, integracji, klastrowania i wizualizacji danych z sekwencjonowania RNA pojedynczych komórek.
survival	CRAN	Version 3.5-7	Pakiet R stosowany do analizy przeżywalności, w tym modelowania proporcjonalnych hazardów Coxa.
SYBR Green qPCR master mix	Beijing Solarbio Science & Technology Co., Ltd., Beijing, China	SR1110	Fluorescencyjny master mix stosowany do ilościowej amplifikacji PCR w czasie rzeczywistym.
timeROC	CRAN	Version 0.4	Pakiet R stosowany do generowania zależnych od czasu krzywych charakterystyki operacyjnej odbiornika (ROC) i obliczania wydajności predykcyjnej w czasie.
Odczynnik TRIzol	Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA	15596026CN	Odczynnik stosowany do ekstrakcji całkowitego RNA z hodowli komórkowych.
Trypsyna-EDTA, 0.25%	Beijing Solarbio Science & Technology Co., Ltd., Beijing, China	T1300	Odczynnik do dysocjacji komórek stosowany do pasażowania i zbierania komórek adherentnych.