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

Zestaw GSE99671 zidentyfikował 105 różnicowo wyrażonych genów związanych ze starzeniem komórkowym według bazy CellAge
GSE99671 obejmował 36 próbek z 18 par tkanek. Po filtrowaniu niskich wartości liczników zachowano 16 683 geny. Przy skorygowanej wartości P < 0,05 stwierdzono różnicową ekspresję 2 248 genów. Przy bardziej rygorystycznym 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 (Ryc. 1A,B). Przecięcie 2 248 genów różnicowo wyrażonych z 866 genami z bazy CellAge pozwoliło wyłonić 105 różnicowo wyrażonych genów związanych ze starzeniem komórkowym (Ryc. 1C). Ekspresja PPARG była obniżona 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, przy wartości P dla testu Wilcoxona dla par = 0,0294 (Ryc. 1D).

Wielomodelowy screening prognostyczny wskazał PPARG jako kluczowy gen kandydujący
W badaniu TARGET-OS uwzględniono 85 pacjentów, u których odnotowano 27 zgonów. Zintegrowany screening 105 genów związanych z procesem starzenia wykazujących różnicową ekspresję, z wykorzystaniem jednowymiarowej regresji Coxa, analizy Kaplana-Meiera, analizy ROC przeżywalności, regresji LASSO, powtarzalnej regresji LASSO oraz modelowania za pomocą losowego lasu przeżywalności (random survival forest), doprowadził do wyłonienia 18 kandydujących genów centralnych (hub genes) przed korektą kliniczną (Ryc. 2A). W jednowymiarowej analizie Coxa PPARG był powiązany z przeżywalnością całkowitą (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 Kaplana-Meiera porównująca grupy o wysokiej i niskiej ekspresji dała wynik P = 0.00784 (Ryc. 2B). Zależne od czasu wartości AUC dla 1, 3 i 5 lat wyniosły odpowiednio 0.603, 0.760 i 0.776 (Ryc. 2C). PPARG wykazał częstotliwość wyboru w powtarzalnym LASSO na poziomie 0.920 oraz istotność w losowym lesie przeżywalności wynoszącą 0.0398 (Ryc. 2D–F). Początkowa punktacja hub score dla PPARG wyniosła 6/6, ponieważ spełnił on wszystkie sześć określonych wcześniej kryteriów screeningu.

Korekta kliniczna potwierdziła wartość prognostyczną PPARG
Po uwzględnieniu klinicznych współzmiennych PPARG pozostało istotnie powiązane 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 wykazał C-index równy 0.707 oraz AIC równy 89.921 (Tabela uzupełniająca 1). Dodanie PPARG zwiększyło C-index do 0.829, zmniejszył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 wykluczeniu operacji radykalnej utrzymała ochronny związek PPARG (HR = 0.249, P = 0.00185; Tabela uzupełniająca 3, Rysunek 2J). PPARG uzyskało najwyższy zintegrowany wynik kliniczny równy 13 (początkowy wynik hub 6 plus siedem punktów integracji klinicznej) i wykazało końcowe wartości AUC w 3. i 5. roku wynoszące odpowiednio 0.770 i 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 w zestawach: 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 z TARGET-OS PPARG nie korelował istotnie z ogólnym wynikiem starzenia CellAge (Spearman ρ = 0.022, P = 0.837; Rysunek 3B), ale korelował z kilkoma pojedynczymi genami CellAge (Rysunek 3C). Analiza mikrośrodowiska immunologicznego wykazała pozytywne korelacje między PPARG a makrofagami (ρ = 0.485, FDR = 2.7 × 10^-5), limfocytami T CD8 (ρ = 0.410, FDR = 5.88 × 10^-4), sygnaturą przypominającą osteoklasty (ρ = 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 wykazywały wyższe sygnatury komórek przypominających osteoklasty, makrofagów, limfocytów T CD8, komórek dendrytycznych, monocytów, neutrofili, komórek NK i 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 zawierał 68 336 komórek i 32 297 genów. Ekspresja PPARG różniła się znacząco pomiędzy typami komórek (Rysunek 4A). Najwyższą średnią ekspresję zaobserwowano w komórkach śródbłonka (średnia ekspresja = 0,540; odsetek dodatnich = 44,33%), perycytach (średnia ekspresja = 0,439; odsetek dodatnich = 41,61%), makrofagach/monocytach (średnia ekspresja = 0,363; odsetek dodatnich = 31,62%) oraz komórkach zrębu związanych z guzem (średnia ekspresja = 0,361; odsetek dodatnich = 44,10%; Rysunek 4B–E). Podgrupa złośliwych komórek kostniakomięsaka wykazywała ekspresję PPARG (średnia ekspresja = 0,163; odsetek dodatnich = 16,70%), jednak ekspresja w złośliwych komórkach kostniakomięsaka nie była znacząco wyższa niż w innych komórkach (FDR = 0,151). Wyniki te sugerują, że ekspresja PPARG w kostniakomięsaku odzwierciedlała przede wszystkim stany mikrośrodowiskowe naczyniowe, mieloidalne i zrębowe, a nie była ograniczona do komórek złośliwych (Rysunek 4F).

Transkryptomika przestrzenna powiązała PPARG z przestrzennymi stanami związanymi ze starzeniem oraz niszami naczyniowymi
Po kontroli jakości zachowano 4 572 przestrzenne punkty (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). Osobno, spośród 866 genów CellAge, 845 dopasowano w macierzy ekspresji przestrzennej (97,58%). PPARG wykazywało ogniskową ekspresję przestrzenną (Rysunek 5C). Przestrzenny wynik starzenia 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 wynikiem śródbłonkowym (ρ = 0,0433, FDR = 0,00592) oraz wynikiem perycytoidnym (ρ = 0,0367, FDR = 0,0181), podczas gdy PPARG korelowało ujemnie z wynikiem złośliwego kostniakomięsaka (ρ = -0,0592, FDR = 0,000219) oraz wynikiem zrazowym guza (ρ = -0,0531, FDR = 0,000774; Tabela uzupełniająca 5, Rysunek 5E). Analiza transferu etykiet (label transfer) wykazała analogicznie dodatnie korelacje z wynikiem predykcji śródbłonka (ρ = 0,0507, FDR = 0,00120) i wynikiem predykcji perycytów (ρ = 0,0394, FDR = 0,0123), wraz z ujemną korelacją z wynikiem predykcji komórek złośliwego kostniakomięsaka (ρ = -0,0699, FDR = 1,1 × 10^-5; Rysunek 5F–H).

Zewnętrzna ekspresja i walidacja eksperymentalna potwierdziły obniżenie poziomu PPARG
Zbiór GSE36001 obejmował 19 próbek kostniakomięsaka oraz sześć prawidłowych kontroli. Poziom PPARG był znacząco obniżony w kostniakomięsaku (logFC = -1.429, P = 0.00730, skorygowana wartość P = 0.0435; Rysunek 6A). W walidacji na liniach komórkowych, 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; Rysunek 6B). Ekspresja białka PPARG była również znacząco zmniejszona w komórkach 143B, co potwierdzono metodą Western blotting (P < 0.01; Rysunek 6C,D). Wyniki uzyskane z zewnętrznej kohorty oraz na poziomie mRNA i białka konsekwentnie 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 dostarczył jedynie zewnętrznej walidacji ekspresji, a nie niezależnej walidacji prognostycznej.

DOSTĘPNOŚĆ DANYCH:
Wszystkie zbiory danych wykorzystane w niniejszym badaniu są publicznie dostępne. Dane GSE99671 i GSE36001 pobrano z bazy danych 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). Dane transkryptomiczne i kliniczne TARGET-OS pobrano z UCSC Xena (https://xena.ucsc.edu/). Geny związane ze starzeniem komórkowym pobrano 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 oraz transkryptomiki przestrzennej 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 użyte do reprodukcji przedstawionych wyników zostały zestawione i przekazane jako Supplementary File 1.

Rycina 1: Identyfikacja genów różnicowo wyrażonych oraz kandydackich genów związanych ze starzeniem pochodzących z bazy CellAge w kostniakomiaku. (A) Wykres wulkaniczny przedstawiający geny różnicowo wyrażone pomiędzy tkankami kostniakomiaka a dopasowanymi kontrolnymi tkankami nietumorowymi w zbiorze danych GSE99671. Istotnie upregulowane i downregulowane geny zostały zaznaczone zgodnie z przyjętymi 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 oraz genów związanych ze starzeniem z bazy CellAge. (D) Porównanie ekspresji PPARG w parach tkanek kostniakomiaka i dopasowanych kontrolnych tkanek nietumorowych w GSE99671. Kliknij tutaj, aby zobaczyć powiększoną wersję tej ryciny.

Rysunek 2: Analizy przeżycia oparte na uczeniu maszynowym i skorygowane klinicznie identyfikują PPARG jako kluczowy gen hub związany z prognostyką starzenia w kostniakomiaku. (A) Wykres leśny przedstawiający wyniki jednowymiarowej regresji Coxa dla kandydujących genów związanych ze starzeniem w kohorcie TARGET-OS. (B) Krzywa przeżycia Kaplana–Meiera porównująca przeżycie całkowite pomiędzy pacjentami 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 Coxa dla selekcji kandydujących genów prognostycznych. (E) Analiza stabilności powtarzanej metody LASSO pokazująca częstotliwości wyboru lambda.min w 300 pięciokrotnych powtórzeniach. (F) Analiza losowego lasu przeżycia (Random survival forest) przedstawiająca wyniki istotności zmiennych z 1000 drzew. (G) Wykres leśny przedstawiający wyniki skorygowanej klinicznie regresji Coxa dla kandydujących genów hub. (H) Zmiany AIC po dodaniu poszczególnych genów hub do modelu klinicznego. (I) Poprawa wskaźnika C-index po dodaniu poszczególnych genów hub do modelu klinicznego. (J) Końcowy ranking zintegrowanego wyniku klinicznego (zakres 0–13) 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. Kliknij tutaj, aby wyświetlić większą wersję tego rysunku.

Rysunek 3: Analiza wzbogacenia funkcjonalnego i 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 ogólnym wynikiem starzenia CellAge w danych bulk 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 za pomocą współczynnika korelacji rang Spearmana (ρ), a skorygowane wartości P obliczano metodą Benjamini-Hochberga. Skróty: GSEA = analiza wzbogacenia zbiorów genów; NF-κB = jądrowy czynnik kappa-B; JAK-STAT = kinaza Janusa-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 mięsaka kostniaka. (A) Wizualizacja UMAP głównych typów komórek w ludzkim zbiorze danych transkryptomicznych pojedynczych komórek mięsaka kostniaka po uproszczonej adnotacji ręcznej. (B) FeaturePlot przedstawiający globalny rozkład ekspresji PPARG w pojedynczych komórkach. (C) DotPlot przedstawiający ekspresję PPARG w głównych typach komórek. (D) Wykres skrzypcowy (violin plot) przedstawiający poziomy ekspresji PPARG w różnych typach komórek. (E) Wykres słupkowy przedstawiający proporcję komórek PPARG-dodatnich w każdym głównym typie komórek. (F) Wizualizacja UMAP przedstawiająca ekspresję PPARG w złośliwych komórkach mięsaka kostniaka. Skrót: UMAP = uniform manifold approximation and projection. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

Rycina 5: Przestrzenna lokalizacja transkryptomiczna PPARG oraz przestrzennych cech związanych ze starzeniem w osteosarcoma. (A) Rozkład przestrzenny klastrów zdefiniowanych na podstawie transkryptomu w przekroju z przestrzennej transkryptomiki osteosarcoma SP_BS3. (B) Rozkład przestrzenny wykrytych genów na spot, przedstawiony przez nFeature_Spatial. (C) Przestrzenny wzorzec ekspresji PPARG w spotach SP_BS3. (D) Rozkład przestrzenny wyniku starzenia (senescence score) wyznaczonego na podstawie CellAge. (E) Analiza korelacji między ekspresją PPARG a przestrzennym wynikiem starzenia wyznaczonym na podstawie CellAge lub wynikami nisz ekologicznych komórek. (F) Mapa predykcji z przenoszenia etykiet (label-transfer) pokazująca dominujący typ komórek wyznaczony z analizy pojedynczych komórek dla każdego spotu przestrzennego. (G) Analiza korelacji między ekspresją PPARG, wynikiem starzenia wyznaczonym na podstawie CellAge a wynikami predykcji typów komórek z przenoszenia etykiet. (H) Rozkład przestrzenny spotów z wysokim poziomem PPARG (PPARG-high) i niskim poziomem PPARG (PPARG-low). Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Rycina 6: Zewnętrzna ekspresja i walidacja eksperymentalna obniżenia poziomu PPARG w osteosarkomie. (A) Wykres pudełkowy przedstawiający poziomy ekspresji PPARG w próbkach osteosarkomy (n = 19) i normalnych próbkach kontrolnych (n = 6) w zbiorze danych GSE36001. (B) Analiza qRT-PCR ekspresji mRNA PPARG w ludzkich komórkach osteosarkomy 143B i ludzkich komórkach kontrolnych osteoblastów. (C) Reprezentatywny Western blot wykazujący ekspresję białka PPARG i GAPDH w ludzkich komórkach kontrolnych osteoblastów oraz komórkach osteosarkomy 143B. GAPDH zastosowano jako kontrolę załadunku. (D) Ilościowa analiza densytometryczna prążków Western blot wykazująca względne poziomy białka PPARG znormalizowane do GAPDH. W punktach (B) i (D) dane przedstawiono jako średnia ± SD z trzech niezależnych eksperymentów. P < 0,01 i 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 po odwrotnej transkrypcji; 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. Indeks C, AIC oraz podsumowanie modelu Coxa zbudowanego z wykorzystaniem wyłącznie zmiennych klinicznych, w tym płci, wieku, stanu choroby w momencie diagnozy, pierwotnej lokalizacji guza, konkretnego regionu guza oraz statusu operacji radykalnej. Skrót: AIC = kryterium informacyjne Akaike. Kliknij tutaj, aby pobrać ten plik.

Tabela uzupełniająca 2: Porównanie modeli Coxa opartych wyłącznie na danych klinicznych oraz modeli łączących dane kliniczne i genetyczne. Wyniki porównania modeli po dodaniu poszczególnych kandydujących genów węzłowych (hub genes) do modelu klinicznego, obejmujące wskaźnik C-index, AIC, statystyki testu ilorazu wiarygodności oraz mierniki poprawy modelu. Skrót: AIC = kryterium informacyjne Akaike. Prosimy kliknąć tutaj, aby pobrać ten plik.

Tabela uzupełniająca 3: Analiza wrażliwości po usunięciu zmiennej dotyczącej operacji radykalnej. Wyniki wrażliwości regresji Coxa oceniające, czy powiązania prognostyczne kandydujących genów hubowych, w szczególności PPARG, pozostały stabilne po wykluczeniu zmiennej 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ównania grup między ekspresją PPARG a sygnaturami ssGSEA związanymi z układem odpornościowym, zrębem, naczyniami, stanem zapalnym i SASP, w tym 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 z senescencją; ssGSEA = analiza wzbogacenia zestawów genów dla pojedynczej próbki. Kliknij tutaj, aby pobrać ten plik.

Tabela uzupełniająca 5: Analiza korelacji transkryptomicznych przestrzennie dla PPARG w SP_BS3. Wyniki korelacji między ekspresją PPARG a przestrzenną oceną starzenia pochodzącą z CellAge, ocenami niszy ekologicznej komórek oraz ocenami przewidywania typu komórek pochodzącymi z transferu etykiet w przekroju transkryptomicznym przestrzennie 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.